Validation of a stereological method for estimating particle size and density from 2D projections with high accuracy

Stereological methods for estimating the 3D particle size and density from 2D projections are essential to many research fields. These methods are, however, prone to errors arising from undetected particle profiles due to sectioning and limited resolution, known as ‘lost caps’. A potential solution by Keiding et al. (1972) accounts for lost caps by quantifying the smallest detectable profiles in terms of their limiting section angle (ϕ). However, this simple solution has not been widely adopted nor validated. Rather, model-independent design-based stereological methods, which do not explicitly account for lost caps, have come to the fore. Here, we provide the first experimental validation of the Keiding model by quantifying ϕ of synaptic vesicles using electron-tomography 3D reconstructions. This analysis reveals a Gaussian distribution for ϕ rather than a single value. Nevertheless, curve fits of the Keiding model to the 2D diameter distribution accurately estimate the mean ϕ and 3D diameter distribution. While systematic testing using Monte Carlo simulations revealed an upper limit to determining ϕ, our analysis shows that mean ϕ can be used to estimate the 3D particle density from the 2D density under a wide range of conditions, and this method is potentially more accurate than minimum-size-based lost-cap corrections and disector methods. We applied the Keiding model to estimate the size and density of somata, nuclei and vesicles in rodent cerebella, where high packing density can be problematic for design-based methods.


36
Estimating the size and density of particles from their orthogonal projection, such as a 2D 37 image, is a common stereological endeavour in the fields of biosciences, petrography, 38 materials science and astronomy 1 . This approach is particularly valuable in the field of 39 biosciences where the size and density of biological structures, such as cells and organelles, 40 are often compared before and after drug perturbations, between normal and disease 41 conditions, species, ages or critical periods of development 2 . Moreover, measures of particle 42 size and density, or the equivalent measure of volume fraction (VF; see Definition of Key 43 Terms), form the basis of our understanding of a wide range of biological phenomena. For 44 example, the density of synaptic vesicles near the active zone has been related to measures 45 of synaptic plasticity 3-5 , the density of cerebellar granule cells (GCs) and mossy fiber 46 terminals (MFTs) has been used to estimate the amount of information transferred across 47 the input layer of the cerebellum 6 and the amount of energy expenditure at the cellular and 48 subcellular level 7 . Stereological measures are also commonly used to assay disease states, 49 such as that of the lung, liver, kidney and brain [8][9][10][11] . Hence, stereological methods for 50 estimating particle size and density have wide application and are of great practical utility. 51 Recent advances in high-resolution volumetric imaging have significantly improved the 52 morphological information available about cells and tissue structure 12-15 making them ideal 53 for 3D analysis. However, these technologies are expensive and full reconstructions are 54 both labour and computationally intensive. The use of stereological methods for analysing 55 a relatively small sample of 2D projections is therefore still the most time efficient and 56 practical solution for most laboratories. 57 Stereological methods for estimating size and density have developed along two distinct 58 approaches: a model-based approach that makes basic assumptions about the geometry of 59 the particle of interest, e.g. Wicksell's transformation 16 and the Abercrombie correction 17 that 60 'caps'. Besides distorting the diameter distribution, caps also introduce a distortion of the 91 apparent density, an effect known as overprojection, the split-cell error or the Holmes 92 effect 17,40-43 . 93 A fundamental limitation of the Wicksell 16 and Bach 38 models, however, is that they assume 94 all caps are resolvable. While this might be true for the largest caps with diameters on the 95 order of F(d), the smallest caps are usually unresolvable, falling below the limits of resolution 96 and contrast or blending in with their surrounding environment 44 It was not until the 1970s that a key innovation for accounting for lost caps was developed 105 by Keiding et al. 49 whereby lost caps are defined with respect to a lower limit of the angle 106 subtended by the observable caps (cap-angle limit, ϕ; Figure 1). In this conceptual model, ϕ 107 is independent of particle diameter, which is important since a distribution of limiting cap 108 sizes can arise when particle size varies while specimen contrast, rather than microscope 109 resolution, limits cap detection. Incorporating ϕ into the Wicksell-Bach model, Keiding et al. 110 derived the following relationship between F(d) and G(d) for spherical particles: 111 where d is the 2D particle diameter, y is the variable of integration and ϕ can vary from 0°, 113 where no caps are lost, to 90°, where all caps are lost. Here, ζ is the mean axial length 114 spanning from below to above the section that contains the center points of those particles 115 observed within the projection (Figure 1), defined as follows: 116 ζ = T + μD⋅cosϕ (2) 117 As expected, Equation 1 reduces to Bach's analytical solution when ϕ = 0° and to Wicksell's 118 analytical solution when ϕ = 0° and T = 0. 119 Another innovation of Keiding et al. 49 was to estimate both F(d) and ϕ from G(d) using a 120 maximum likelihood estimation (MLE) algorithm, rather than using an unfolding algorithm, 121 thereby providing a better quantification of ϕ. Knowing ϕ is particularly useful since it can be 122 used via Equation 2 to estimate the 3D particle density (λ3D) from the measured 2D particle 123 density (λ2D) as follows 50 : 124 λ3D = λ2D / ζ (3) 125 (Appendix A). This 'correction' method for estimating λ3D is potentially more accurate than 126 using the classic Abercrombie correction 17 , which assumes no caps are lost (i.e. ϕ = 0°), or 127 the Floderus 42 and Konigsmark 47 corrections that use the minimum cap penetration depth 128 (hmin) or equivalent minimum cap diameter (dmin), both of which are likely to be outlier 129 measures. Unfortunately, despite potentially providing the most accurate description of the 130 lost-cap distribution via ϕ, the Keiding model has not been widely adopted nor validated. 131 G(d), especially the cap tail, depends on the relative size of the particles compared to the 212 section thickness and the imaging method used to acquire the projections. 213  To quantify the relationship between G(d), F(d) and lost caps, we computed the distribution 235 of lost caps (L(d)) for ϕ = 5-90° and compared L(d) to F(d) (Figure 3F). Results showed L(d) 236 had a negatively skewed Gaussian shape with tail descending to 0 u.d. The upper limit of 237 L(d) showed a dispersion, rather than a hard cutoff, since a variation in particle size in 238 combination with a fixed-ϕ limit resulted in limiting caps of different size. For ϕ < 55°, there 239 was a clear separation between L(d) and F(d), i.e. nearly all lost caps had a diameter smaller 240 than those that define F(d). For ϕ > 55°, on the other hand, L(d) overlapped F(d), in which 241 case there was a large number of lost caps with diameters equivalent to those that define 242 F(d). At the most extreme condition where all caps are lost (i.e. ϕ = 90°) L(d) spanned from 243 0 to the largest diameter of F(d). Hence, at these larger ϕ, L(d) and F(d) cannot be delineated 244 by diameter size. 245 2.3 Estimating the 3D particle size and cap-angle limit using the Keiding 246 model 247 Next, we tested the Keiding model's capacity to estimate F(d) and ϕ from G(d). To do this, 248 we used a 3D Monte Carlo simulation package to compute virtual projections of spherical 249 particles for planar and thick sections, creating lost caps according to the Keiding model (i.e. 250 a cap was removed from the projection when its θ < ϕ; Figure S2; Materials and Methods). 251

Exploration of the effects of section thickness and lost caps on G(d)
To make a direct comparison to our numerical solutions of Equation 1, we matched F(d) of 252 the simulated particles to that used in the numerical solutions. Moreover, we set the number 253 of measured 2D diameters per projection (~500) to match the average sample size of our 254 experimental G(d) of GC somata and nuclei, the two datasets with the largest number of 255 measured diameters (Table 2). From the 2D diameters of the simulated particles, we 256 On the other hand, when ϕ > 55°, estimates of μD and σD were less accurate, with positive 262 and negative biases, respectively, and estimates of ϕ often had a large negative bias (Figure 263 4B and C). In this case, the LSE routine had difficulty estimating true ϕ since G(d) ≈ F(d) as 264 ϕ approached 90°. The similarity between G(d) and F(d) at true ϕ > 55° was greatest for 265 thick sections in which case fitting a simple Gaussian function (Equation 5) to G(d), which is 266 equivalent to curve fitting the Keiding model with ϕ = 90°, where all caps are lost, resulted 267 in similar estimates of μD and σD, as did simply using the 2D measures μd and σd as 268 estimates. Moreover, the estimates of μD and σD of the three approaches all showed 269 relatively small biases (< 2 and 5%, respectively, for T = 1 u.d.). An overall comparison 270 revealed that thick sections were moderately better for estimating μD and σD, and planar 271 sections were moderately better for estimating ϕ. Repeating the error analysis for G(d) 272 computed from ~2000 diameters gave qualitatively similar results, except μD, σD and ϕ had 273 smaller biases and confidence intervals, where the confidence intervals followed a 1/√n 274 relation ( Figure S5). 275 While the above results show the Keiding model works well in estimating μD, σD and ϕ, 276 estimates were most accurate when ϕ < 55° ( Figure 4C). However, this upper limit of ϕ 277 (ϕcutoff) is dependent on the spread of F(d) and the number of measured diameters, both of 278 which were set in our Monte Carlo simulations to match our experimental data (CVD = 0.09, 279 ~500 diameters). For a wider F(d) and/or smaller number of diameters, accurate parameter 280 estimation will be limited to a smaller range of true ϕ, and vice versa ( Figures 4C and S5 test the accuracy of estimated ϕ using estimated μD and σD, since the estimation errors of 286 μD and σD are relatively small even when true ϕ > ϕcutoff ( Figures 4C and S5). Results showed 287 that indeed this was possible with relatively small adjustments to the sinϕcutoff relations 288 (Equation 9). Hence, if one requires an accurate estimate of ϕ, for example when computing 289 particle density as discussed further below, and one only has estimates of μD and σD, then 290 Equation 9 can be used as an accuracy test, i.e. by requiring estimated ϕ < estimated ϕcutoff 291 (Figure S6). Although the fit error of ϕ can also be used as a measure of accuracy, ~10% of 292 our curve fits to simulated G(d) for T = 0 u.d. and true ϕ > 55° showed small fit errors (< 1°) 293 but large estimation errors (|Δϕ| > 5°). Hence, a combination of the fit error and Equation 9  294 can be used as an accuracy test for ϕ. 295

296
A potential source of error of the Bach 38 and Keiding 49 models for thick sections is that the 297 models assume particles near the bottom of the section can be identified and outlined in a 298 projection as well as those toward the top of the section. While this may be true for thin 299 sections, it is unlikely to be true for thick sections with overlapping particle projections 2,57 . 300 To investigate what effects overlapping projections might have on the ability to accurately 301 estimate F(d) and ϕ from G(d), we used simulations to compute the sum of 2D projection 302 overlaps between a given particle and particles higher in a thick section, where particles had 303 a random distribution and a modest to high VF (0.15-0.45). Results of the analysis showed 304 that, as expected, particles toward the bottom of the section experienced more overlaps in 305 a projection, and this effect was greater for thicker sections and higher particle density 306 (Figure 5A). To investigate what effect projection overlaps might have on the shape of G(d), 307 we set an upper limit to the amount of projection overlaps an observable particle can have, 308 i.e. particles with overlaps above the set limit (ψ) were considered hidden from view, i.e. 309 lost, and excluded from the projection. Results of these 'semi-transparent' particle 310 simulations showed that, for thick sections and high particle density, a large number of 311 particles near the bottom of the section were excluded from the projection, thereby reducing 312 the effective section thickness (T) for a given projection ( Figure 5B; ψ = 0.25). Nevertheless, 313 excluding bottom-dwelling particles from the projection had minimal effect on the shape of 314 G(d), and therefore had little effect on estimates of F(d), as long as the top-dwelling particles 315 were simulated as transparent, i.e. small degrees of projection overlap did not interfere with 316 counting the projections or computing their size ( Figure 5C). This is in contrast to when 317 particles were simulated as opaque such that overlapping projections reduced the number 318 of observed projections and increased their size; in this case, opaque particles created a 319 positive skew in G(d) ( Figure 5D). However, the distortions in G(d) produced only modest 320 changes in estimates of F(d). Hence, these results indicate overlapping projections of semi-321 transparent and opaque particles in thick sections create only small biases in estimates of 322 F(d). However, the overlapping projections do have the potential to reduce the effective 323 section thickness (Figure 5B) and therefore affect the estimates of the 3D density (see 324 below (i.e. stochasticity) in particle location and size, they lacked variation normally associated with 330 experimental data such as that due to finite spatial resolution, limited contrast, irregular 331 particle shape and blending with surrounding material (e.g. intracellular/extracellular 332 proteins). Hence, to address this shortcoming, we tested the accuracy of the Keiding model 333 using a volumetric ET z-stack of clusters of MFT vesicles, since this allowed us to directly 334 measure F(d) and ϕ, as well as G(d). First, F(d) was measured by outlining individual 335 vesicles across multiple planes of the z-stack ( Figure 6A; ET11). From the outlines of a given 336 vesicle, equivalent radii were computed as a function of z-depth and fit to a circle (Equation 337 10), resulting in estimates for the vesicle's 3D diameter (D) and z-axis center point (z0) 338 (Figure S7). Aligning the vesicle profiles at their centers revealed the profiles overlaid each 339 other and fit well to a semi-circle ( Figure 6B). From the measured 3D diameters, we 340 computed μD ± σD = 42.7 ± 3.4 nm and F(d) ( Figure 6D; n = 234). To determine the shape 341 of F(d), we fit the distribution to Gaussian, chi and gamma functions (Equations 5, 6, 7) and 342 found closely overlapping fits ( Figure S8). Moreover, we found parameter f > 70 for the chi 343 and gamma fits indicating Gaussian-like distributions. Repeating the same analysis for F(d) 344 derived from another ET z-stack (ET10) and another study of MFT vesicles 5 produced 345 similar results. Hence, F(d) of the MFT vesicles was well described by a Gaussian function, 346 which is consistent with the findings of a 3D analysis of vesicles in hippocampal CA1 347 excitatory synapses 58 . 348 Next, we measured ϕ by computing the axial extent of each vesicle profile, which only 349 partially extended to the north and south poles ( Figures 6B and S7A). The missing profiles 350 at the pole regions, i.e. lost caps, are likely due to limited resolution and poor contrast of 351 orthogonally oriented membranes 45 . Plotting the minimum measured diameter (δmin) for the 352 vesicle poles interior to the z-stack versus their 3D diameter revealed a linear relation that 353 was well described by the Keiding model, i.e. δmin = D⋅sinϕ ( Figure 6C1). In contrast, the 354 minimum of all measured diameters (dmin), which has been used as a correction for lost caps 355 (Appendix A), was a poor fit to the δmin-D relation. Converting all δmin to ϕ values revealed a 356 distribution that spanned 27-65° ( Figure 6C2) with μϕ ± σϕ = 41.4 ± 6.8° (CVϕ = 0.2; n = 397) 357 and was well described by a Gaussian distribution ( Figure 6D, inset). Hence, this analysis 358 revealed a source of variation not accounted for in the Keiding model, which assumes all 359 particles have the same ϕ. The ramifications of a variation in ϕ are explored in the next 360 section. 361 Next, we computed G(d) from all measured 2D diameters ( Figure 6D). Comparison of G(d) 362 to F(d) revealed a negative skew in G(d) as expected for planar sections (estimated T ≈ 0.1 363 u.d.; Table 1). Finally, using the measured G(d), F(d) and ϕ, we tested the Keiding model C and S12). The ϕ-accuracy test computed via Equation 9 showed estimated ϕ < estimated 431 ϕcutoff (~44-50°) for all but two fits where estimated ϕ was ~3-4° above its estimated ϕcutoff. 432 For the two fits that failed the ϕ-accuracy test, we recomputed the fits while fixing their ϕ to 433 that estimated from the other fits of the same preparation (this made little difference in 434 estimates of μD and σD). Averaging results across the 3 rats gave μD = 5.78 ± 0.16 μm, σD = 435 0.41 ± 0.03 μm and ϕ = 37 ± 2° (±SEM; Table 2 Having an accurate estimate of the number of particles per unit volume is essential to many 480 fields of science. Historically, estimates of 3D particle density (λ3D) were obtained by 481 computing the number of particles per unit area (λ2D; Figure S1) observed in a projection 482 divided by the section thickness (T). However, particle caps on the top and bottom of the 483 section inflate the particle count with respect to T, i.e. they create an overprojection 17,40-43 . 484 To correct for overprojection, one can use ϕ to estimate λ3D via Equation 3 (Appendix A). 485 As a first test of the validity of Equation 3, we used the Monte Carlo simulations in Figure  486 4C to investigate the relation between the measured λ2D and true λ3D by computing ζ = λ2D / 487 λ3D. Comparison of this measured 'true' ζ to the expected ζ computed via Equation 2, using 488 true T, μD and ϕ, showed a close agreement ( Figure 10A1). Hence, our simulations of the 489 Keiding model were well described by Equation 3. Next, we tested the ability of the Keiding 490 model to accurately estimate λ3D via Equation 3 by computing Δλ3D for the same simulations 491 using parameters μD and ϕ estimated via Keiding-model curve fits to G(d) ( Figure 4C). 492 Results showed accurate estimates of λ3D, except when true ϕ > 45°, in which case the 493 estimation errors of μD and ϕ translated into estimation errors of λ3D ( Figure 10A2). Here, 494 the finding that ϕcutoff for estimated λ3D is less than that for estimated μD, σD and ϕ (~55°) is 495 consistent with an increase in variability from using estimated μD and ϕ to compute ζ, and 496 therefore consistent with our ϕ-accuracy test, where estimated ϕcutoff ≈ 43° (Equation 9). 497 Hence, these results demonstrate the ability of Equation 3 to accurately estimate λ3D using 498 Keiding-model estimates of μD and ϕ as long as estimated ϕ < estimated ϕcutoff. Comparing 499 results for planar and thick sections showed qualitatively similar results, except errors were 500 smaller for thick sections due to smaller errors in μD (but see next paragraph for caveats of 501 using thick sections). Repeating the error analysis for G(d) computed from ~2000 diameters 502 gave qualitatively similar results compared to those for ~500 diameters, but λ2D and λ3D had 503 smaller biases and confidence intervals, where the confidence intervals followed a 1/√n 504 relation for ϕ < 45° ( Figure S17) similar to that of estimated μD and ϕ ( Figure S5). 505 For thick sections, a major caveat of estimating λ3D from λ2D via Equation 3 is that the 506 calculation assumes one has a good estimate of ζ. As shown in Figure 5B, however, 507 overlapping projections of semi-transparent and opaque particles may preclude counting 508 particles at the bottom of the section, thereby reducing ζ. To demonstrate this, we computed 509 Δζ for simulations similar to those in Figure 5 and found large positive biases for both semi-510 transparent and opaque particles, i.e. estimated ζ was larger than true ζ ( Figure 10A1). As 511 expected, the overestimation of ζ translated into large negative biases in estimates of λ3D 512 ( Figure 10A2). Given this caveat, therefore, the best approach to estimate λ3D from λ2D for 513 particles with a high density is to use planar sections, in which case there will be little to no 514 interference in counting particles from overlapping particle projections. 515 As a second independent method for computing the 3D density, we used the particle area 516 fraction (AF) to VF relation of Weibel and Paumgartner 60 (VF = Kv⋅AF) to compute the VF of 517 our simulations, where AF is the sum of all projection areas divided by the total projection 518 area (Areaxy) and Kv is a proportional scale factor that is a function of T and μD, but modified 519 to be a function of ϕ rather than hmin (Appendix B). Results showed excellent agreement 520 between the estimated and true VF of our simulations ( Figure 10A3). However, for the 521 simulations of overlapping projections of semi-transparent and opaque particles in thick 522 sections, there was a significant underestimation of the VF, as expected. ), which was also nearly the same to λ3D computed via a 3D measurement 536 (Δλ3D = -3%). In contrast, when we used dmin of the z-stack analysis (19 nm) to estimate ζ 537 (using measured μD) rather than ϕ (where dϕ = 28 nm), we computed λ3D = 9688 μm -3 . 538 Hence, in this example, using dmin to estimate λ3D caused significant error (Δλ3D = -15%). 539 Next, we investigated the range of expected error for the z-stack density analysis, as well 540 as that from assuming a fixed ϕ ( Figure S10), by computing λ2D and λ3D for Monte Carlo 541 simulations that mimicked the 3D analysis. Results gave ranges of Δλ3D that were consistent 542 with the above measured Δλ3D ( Figure S18) and showed the assumption of a fixed ϕ 543 introduces only small errors for estimating λ3D from λ2D via Equation 3 for ϕ < 50° (Figure 544 S19). Again, using dmin to estimate λ3D caused significant error (Δλ3D = -19.2 ± 1.9%; n = 545 100 simulation repetitions). Finally, using the VF = Kv⋅AF relation (Equations B.1, B.2c) we 546 computed a VF via the 2D analysis (0.46) that matched that computed via the 3D analysis 547 (0.47). Repeating the same 2D versus 3D density analysis for our other ET z-stack (ET10) 548 produced similar results ( Figure 11B; Table 3). However, unlike ET11, the density analysis 549 of ET10 had to be confined within a subregion of the z-stack due to a nonhomogeneous 550 distribution of vesicles in the axial axis. Overall, our 3D analysis of MFT vesicles shows that 551 Equation 3 accurately estimated λ3D from λ2D with only small error when true ϕ < ϕcutoff, 552 confirming our previous results from simulations (Figure 10). 553 2.12 Estimating the 3D particle density via the disector method (a 554 comparison) 555 The disector method 18,52 is a popular method for estimating λ3D of particles with arbitrary 556 geometry using two adjacent sections from a z-stack, referred to as the 'reference' and 557 'lookup' section. This method counts the number of particles that appear within the reference 558 section but do not appear in the lookup section, i.e. the method counts a particle only if its 559 leading edge appears within the reference section. To compute λ3D, one determines λ2D from 560 the leading-edge count and divides λ2D by the distance between the reference and lookup 561 sections (i.e. the section thickness for adjacent pairs). Although theoretically, lost caps 562 should not introduce error into this counting method 52 , underestimation errors due to lost 563 caps have been reported 61 . Given the popularity of the disector method, as well as the 564 potential source of error due to lost caps, we decided to use the disector method to estimate 565 λ3D of our simulated projections where ϕ = 10-80° and compare these results to our analysis 566 in Figure 10, where λ3D was estimated via Equation 3, i.e. the ϕ-correction method. Using T 567 = 0.3 u.d. as recommended 18 , we found no bias due to lost caps for all ϕ tested (Figure 12, 568 black squares; ~500 particles per projection). This result is consistent with historical thinking 569 that lost caps do not create a bias in counting: while the effect of lost caps on any given 570 reference section is to reduce the particle count, the effect of lost caps on the lookup section 571 is to increase the particle count for that reference section, the two effects thereby 572 cancelling 52 . In essence, lost caps simply shift the reference section to which particles are 573 counted and therefore have no effect on the final particle count. However, these simulations 574 did not take into account a potentially significant bias pointed out by Hedreen 27 : fewer caps 575 are likely to be identified in the reference section compared to the lookup section since the 576 former is done blind and the latter is not (the researcher searches for caps in the lookup 577 section only after identifying a particle in the reference section). In fact, our analysis of ET 578 z-stacks of MFT vesicles revealed a 17-22° bias in ϕ between a blind and nonblind vesicle 579 detection (Figures 7). To investigate the effects of such a bias, we added a bias to our 580 disector simulations by increasing ϕ of the reference section with respect to the lookup 581 section (ϕbias = 5-20°) thereby reducing cap identification in the reference section (i.e. 582 increasing the number of lost caps) with respect to the lookup section. Interestingly, results 583 of these simulations showed large underestimation errors even for small ϕbias (Figure 12; -584 5% < Δλ3D < -20% for ϕbias = 10°). Hence, the blind-versus-nonblind bias in particle detection 585 has the potential to create a significant underestimation of λ3D. Repeating the disector bias 586 analysis using our experimental ET z-stack data gave similar results, showing a large error 587 (Δλ3D = -41 ± 5%) for the average measured ϕbias = 20°. 588 Besides the potential error due to a blind-versus-nonblind bias in particle detection, our 589 simulations show that the disector method is ~2 to 3-fold less accurate at estimating λ3D (i.e. 590 Δλ3D has a larger ±σ) compared to the ϕ-correction method ( Figure S17; ±σ panel) due to a 591 lower particle count per section. Similar results were found using the disector method to 592 compute λ3D of our ET z-stack of MFT vesicles, where Δλ3D = ±6% for the disector method 593 (±σ; ϕbias = 0°) compared to Δλ3D = ±2% for λ3D computed via Equation 3 ( Figure S18). 594 Hence, the inherent smaller count per region of interest (ROI) of the disector method leads 595 to a larger uncertainty in estimates of λ3D. to the estimation of λ3D from λ2D for our sample of GC somata. Using the same confocal 600 images used to compute the 9 G(d) of the GC somata of 3 rats, we computed λ2D = 17.7-601 22.6 × 10 3 mm -2 ( Figure S1A; 8792-29,845 μm 2 ROI area per image, 196-528 somata per 602 ROI, 3 images per rat). Using estimated μD and ϕ for each G(d), we then computed ζ = 5.5-603 7.6 μm via Equation 2 (estimated T = 1.4-2.4 μm, scaled for z-shrinkage; Table 1). Finally, 604 we computed λ3D = λ2D/ζ = 2.7-4.1 × 10 6 mm -3 with mean λ3D = 3.2 ± 0.2 × 10 6 mm -3 and VF 605 = 0.32 ± 0.01 (±SEM; Table 2; Equation 4). To verify these results using the relation VF = 606 Kv⋅AF, we computed Kv = 0.58-0.77, AF = 0.39-0.50 and VF = Kv⋅AF = 0.27-0.37, with 607 mean VF = 0.32 ± 0.02 and λ3D = 3.1 ± 0.1 × 10 6 mm -3 (±SEM) which is not significantly 608 different to λ3D computed via Equation 3 (p = 0.4, paired t-test). 609 Next, we estimated λ2D of GC nuclei from the TEM images used to compute G(d) (λ2D = 610 17.1-34.2 × 10 3 mm -2 ; Figure S1B; 2046-5747 μm 2 ROI area per image, 35-158 nuclei per 611 ROI, 6-7 images per mouse). Using estimated μD and ϕ for each G(d), we then computed ζ 612 = 4.2-4.8 μm via Equation 2 (T = 0.06 μm) and λ3D = λ2D/ζ = 3.5-8.1 × 10 6 mm -3 , with mean 613 λ3D = 5.9 ± 0.4 × 10 6 mm -3 and VF = 0.34 ± 0.01 (±SEM; Table 2 Unlike our estimates for F(d) of GCs in rats and mice, which show no difference, comparison 624 of the above estimates for λ3D in rats and mice show a 2-fold difference: estimated λ3D of the 625 rat dataset (3.2 × 10 6 mm -3 ) is significantly smaller than estimated λ3D of the mouse dataset 626 (5.9 × 10 6 mm -3 ; p = 0.004). While this difference could reflect a true difference between 627 species, a previous comparative study of the cerebellum, which used the same methodology 628 across species, reported similar densities within the GC layer of rats and mice 62 . Hence, it 629 is more likely that the 2-fold difference in estimated λ3D is due to one or more differences in 630 methodology. Although there are perhaps too many differences in methodology to make a 631 decisive conclusion, notwithstanding the confounding problem of the 'reference trap' 63 , it is 632 still instructive to consider them here. The first difference in methodology is section 633 thickness: thick versus planar. Because the rat dataset was derived from thick sections, 634 overlapping opaque projections could have created an underestimate of λ3D compared to 635 the mouse dataset ( Figure 10). The second difference in methodology is imaging 636 technology: confocal versus TEM. Because of the lower contrast and resolution of the 637 confocal images, it was considerably harder to identify and delineate GC somata profiles 638 compared to the GC nuclei profiles. The poor delineation of GCs would create an 639 undercount. The third difference in methodology is tissue preparation: chemical versus cryo 640 fixation. A change in the extracellular volume due to chemical fixation of the rat tissue 641 preparation could have created a lower λ3D. However, this explanation is unlikely since 642 chemical fixation tends to reduce the extracellular volume 59 which would result in a higher 643 rather than lower λ3D. The final difference in methodology is the sampling space: ROIs were 644 larger in the rat dataset compared to the mouse dataset. Because GCs are not distributed 645 uniformly, but rather form high-density clumps interspersed by MFTs and blood vessels, 646 there could be a bias towards a larger λ3D in the mouse dataset due to smaller ROIs. 647 However, our analysis of the two datasets using different sized ROIs (not shown) indicates 648 the difference in ROI size cannot account for the 2-fold difference in estimated λ3D. Hence, 649 we suspect our estimate of λ3D in the rat dataset is underestimated, and this is most likely 650 due to overlapping projections in thick sections and poor delineation of GC profiles in the 651 confocal images. Because of the tight packing of GC somata, estimates of λ3D via Equation 652 3 are best achieved by using planar sections, a higher contrast preparation and superior 653 microscope resolution. We therefore believe our most accurate estimate of the GC λ3D is 654 that of our mouse dataset. Although our estimated λ3D of GC somata in rats is likely to be 655 underestimated, it is still 1.7-fold larger than our previous estimate from the same confocal 656 images 6 . Because the latter estimate of λ3D was computed via the disector method, we 657 suspect it is underestimated due to the blind-versus-nonblind bias in particle detection 658 discussed in the previous section. Hence, this result suggests the local MFT-to-GC 659 expansion ratio is larger than previously estimated 6 . 660 ). Finally, we computed λ3D = 3100-7200 667 μm -3 via Equation 3. Not surprisingly, this range of λ3D is considerably smaller than that for 668 our 3D ET analysis of MFT vesicles (λ3D = 9,000-11,000 μm -3 ). Again, similar to our analysis 669 of GCs in confocal images, we suspect λ3D is underestimated due to overlapping projections 670 in thick sections (Figure 10). To estimate the effective section depth of those vesicles 671 counted, we assumed λ3D of the TEM dataset was the same as that of our ET dataset and 672 computed ζ = λ3D/λ2D = 30-43 nm. This range of ζ is ~2.0-fold smaller than that estimated 673 via Equation 2, indicating vesicles at the bottom of the tissue sections were not counted 674 ( Figure 5B), creating an underestimation of λ3D. Hence, this analysis lends support to our 675 conclusion that the optimal method for estimating λ3D from λ2D for particles with a high 676 density is to use planar sections, in which case there will be no interference of counting from 677 overlapping projections. Since the thinnest possible tissue section created via an 678 ultramicrotome is currently on the order of a vesicle, the best option for computing vesicle 679 density is via ET reconstructions (Figure 11). 680

681
Stereological methods for estimating the 3D size distribution (F(d)) and density (λ3D) of a 682 collection of particles from their 2D projection are essential tools in many fields of science. 683 These methods, however, inevitably contain sources of error, one being the unresolved or 684 nonexistent profiles known as lost caps. Surprisingly, the simple solution for lost caps 685 developed by Keiding et al. 49 , which defines lost caps of spherical particles with respect to 686 a single (i.e. fixed) cap-angle limit (ϕ), has not been widely adopted and has never been 687 validated. Here, we provide the first experimental validation of the Keiding model by 688 quantifying ϕ of unresolved vesicle caps within 3D reconstructions. While this analysis 689 reveals a Gaussian distribution for ϕ rather than a single value, curve fits of the Keiding 690 model to the 2D diameter distribution (G(d)) nevertheless accurately estimate the mean ϕ, 691 as well as F(d). Parameter space evaluation with Monte Carlo simulations revealed that the 692 estimates are most accurate when ϕ falls below a specific value (ϕcutoff). Hence, our 693 experimental and theoretical analyses reveal that, if one only wishes to estimate F(d) from 694 a 2D projection, then the Keiding model can be called to task, whether one is using planar 695 or thick sections. On the other hand, if one wishes to estimate both F(d) and λ3D from a 2D 696 projection, then one will need an accurate estimate of ϕ ( Figure 13). As we discuss below, 697 obtaining an accurate estimate of ϕ for some preparations may require optimising 698 experimental conditions. 699

700
There are five basic assumptions of the Keiding model 49 that one must keep in mind when 701 applying it to the estimation of particle size and density via Equations 1-3. 702 The first assumption is that the particles of interest are approximately spherical, i.e. they are 703 convex with the average shape of a sphere in rotation, which includes elliptical 16 . The 704 assumption of a spherical shape is usually valid for vesicles, vacuoles and nuclei 1,17,44,49,64 , 705 but also for large structures such as follicles 16 and glomeruli 65 . 706 The second assumption is that F(d) is well described by a probability density function (PDF), 707 e.g. a Gaussian, chi or gamma distribution. The assumption that F(d) was a simple Gaussian 708 distribution worked well for our analysis of GC somata and nuclei and MFT vesicles, and for 709 the liver cell nuclei of Keiding et al. 49 (Figure S3). Similarly, the assumption that F(d) was a 710 chi distribution worked well for Wicksell's corpuscles ( Figure S4). If uncertain of the shape 711 of F(d), one can compare fits to distributions obtained from thick sections (T > D) where it is 712 (Figures 2C2 and S3). 713 The third assumption is that the particles have a random distribution; or if the particles are 714 clustered, then the clusters have a random distribution. It is not necessary that the spatial 715 distribution is perfectly random (i.e. Poisson) but rather there is no order to the distribution. 716 For example, a lattice structure (e.g. hexagonal packing) would be problematic since the 717 particles in a 2D projection would have a discrete size distribution. 718 The fourth assumption is that ϕ is the same (fixed) for all particles. Because our 3D analysis 719 of MFT vesicles revealed not a single value for ϕ, but a range of ϕ well described by a 720 Gaussian distribution (CVϕ = 0.2), this assumption could be problematic. However, we found 721 curve fits of the Keiding model to G(d) accurately estimated the mean of this distribution 722 ( Figures 6D and S9D). Moreover, by incorporating a Gaussian-ϕ model into our Monte Carlo 723 simulations, we were able to measure the bias introduced by the fixed-ϕ assumption over a 724 range of ϕ distributions (i.e. μϕ) and found the bias was relatively small for μϕ < ϕcutoff (Figures 725 S11 and S19). Because vesicles are at the lower limit of resolution, we suspect the spread 726 of our measured ϕ distribution, i.e. CVϕ = 0.2, may be close to a worst-case scenario. 727 The fifth assumption pertains to thick sections: that they are perfectly transparent so that 728 one observes 2D projections of particles at the bottom of the section as well as those of 729 particles at the top of the section. For opaque particles with high density this assumption is 730 clearly problematic 2,57 . For semi-transparent particles with high density this assumption is 731 less problematic since overlapping projections do not necessarily preclude drawing their 732 outline ( Figure 2C1) or counting them ( Figure S1C). Yet, one can imagine a section of 733 enough thickness that overlapping projections preclude outlining or counting the particles at 734 the bottom of the section. Using Monte Carlo simulations, we tested the effect of such a 735 scenario by removing particles from simulated projections if they experienced a total sum of 736 overlaps greater than a set limit. Interestingly, results of these simulations showed only a 737 small effect on the shape of G(d) and therefore only small biases for estimated F(d) (Figure 738  imaging or ET). Second, the lateral resolution of the microscope (ρxy) and image (Rxy) should 761 be high with respect to the size of the particle. Third, the efficiency of particle staining and 762 contrast between the particles and their surrounding environment should be high. Fourth, 763 surfaces of the tissue sections should be avoided when creating images, e.g. using guard 764 zones in the axial axis, to avoid lost caps of the nonexistent type, i.e. caps that either fall off 765 the surfaces of tissue sections or never existed, i.e. the microtome failed to transect the 766 particles during sectioning 27,46-48 . Fifth, if the images to analyse exist within a z-stack, cap 767 identification should be performed in a nonblind manner, i.e. by tracking particles through 768 adjacent planes of the z-stack. 769 Given these considerations, our most suitable preparation for computing G(d) was that of 770 the GC nuclei in high-resolution TEM images, where there were few lost caps and a small 771 ϕ = 20°. For this dataset, the tissue sections were planar (T ≈ 0.01 u.d.) and the lateral 772 resolution of the microscope was high with respect to the nuclei (ρxy = 1 × 10 -4 u.d.; Table  773 1). Moreover, the GC nuclei were easy to identify and delineate due to their dark spotted 774 appearance. In contrast, our analysis of GC somata in confocal images gave a larger 775 number of lost caps and larger ϕ (37°). For this dataset, the optical sections were not planar 776 (T ≈ 0.3 u.d), the lateral resolution of the microscope was comparatively low (ρxy = 0.06 u.d.) 777 and, due to the opaque immunolabeling and dense packing of the somata, the task of 778 delineating between adjacent somata was more difficult compared to that of the nuclei in 779 TEM images. and H is the distance between optical sections. Hence, to obtain an accurate estimate of 809 λ3D, one must consider obtaining an accurate measure of T or H 2,25,26,45,66 . However, 810 distortions of particle density along the axial axis of tissue sections must also be 811 considered 2,28,30,48,67,68 . These include uniform shrinkage and differential deformations along 812 the axial axis, and lost caps at the surfaces of the sections (i.e. nonexistent caps). 813 The challenges of estimating T and avoiding axial distortions of particle density also pertain 814 to serial 3D reconstructions. Hence, the reason 3D reconstructions should not automatically 815 be assumed to be the 'gold standard'. For the 3D reconstructions used in this study, we 816 avoided lost caps of the nonexistent type by avoiding the section surfaces. Moreover, the 817 spherical shape of the vesicles allowed us to estimate the axial tissue shrinkage ( Figure S7) 818 which for EM can be considerable depending on the amount of electron beam exposure 69 . 819 Finally, plots of vesicle density as a function of z-depth allowed us to verify that λ3D was 820 computed within a homogeneous vesicle distribution ( Figure 11A2 and B). Hence, the effort 821 to estimate λ3D from a 3D reconstruction in this study was nontrivial. 822 With these axial-axes difficulties in mind, it becomes clear that the 2D methods have a 823 distinct advantage over the 3D methods in that they can effectively remove T from the 824 estimation of λ3D by using planar sections. Under this condition, any bias due to an 825 inaccurate estimate of T or uniform tissue shrinkage along the axial axis would be small. 826 Hence, the 2D ϕ-correction method for estimating λ3D, when applied to data derived from 827 planar sections, is potentially the most accurate method of estimating λ3D. However, planar 828 sections come with their own challenges. First, planar sections can be technically 829 challenging and costly to achieve, perhaps requiring TEM or ET. Second, due to the use of 830 high-resolution microscopy, planar sections are likely to have a significantly smaller field of 831 view, potentially creating bias if particles have a nonhomogeneous distribution. A smaller 832 field of view, however, is not strictly prohibitive since it can be counteracted by analysing 833 more images that have been acquired using a design-based random sampling strategy 70,71 . 834 Our recommendation for the use of planar sections with 2D methods is counter to previous 835 recommendations of using thick sections 25,26,47 . The reasoning put forth for using thick 836 sections is that any bias introduced by lost caps (which affects the magnitude of μD in 837 Equation 3) will be relatively small in comparison to T. However, this approach will only be 838 valid if particles have a low density, one can reliably count particles at the bottom of the 839 sections, one has a good measure of T and one can correct for any distortions of particle 840 density along the axial axis, as discussed above. Moreover, our analysis indicates that when 841 sections are thick, most caps are likely to be lost, in which case ϕ will be close to 90° and 842 indeterminable ( Figure 9 and Table 4). In the case of an indeterminable ϕ, one can use a 843 range of ϕ (e.g. ϕcutoff-90°) to estimate λ3D. This would be an improvement to using the 844 Abercrombie correction 17 that assumes no caps are lost (ϕ = 0°) or dmin correction 47 (or 845 equivalent hmin correction 42 ) since dmin is not a good measure of the lost-cap distribution 846 when ϕ > 20° (Appendix A). 847

848
Although ϕ was never reported for Wicksell's spleen corpuscles 16 , our curve-fit analysis of 849 Wicksell's G(d), which produced an estimated F(d) that matched that of Wicksell, resulted 850 in ϕ = 25° with a small fit error of 3° ( Figure S4). In the Wicksell study, the corpuscles were 851 large (μD = 323 μm) in comparison to the tissue section (18 μm), in which case T ≈ 0.06 u.d. 852 Hence, the small estimated ϕ for Wicksell's dataset is consistent with that of our GC nuclei 853 in planar sections where ϕ = 20°. 854 Another study of GC nuclei 64 computed a lost-caps correction factor f = 0.978. If one 855 expresses f as a cap-angle limit, then ϕ = cos -1 (f) = 12°. This ϕ is smaller than that computed 856 for the GC nuclei in this study (20°), which is unexpected given Harvey and Napper used 857 thin sections (T ≈ 0.3) and light microscopy compared to planar sections and TEM. It seems 858 likely that the smaller ϕ of Harvey and Napper is due to their use of dmin to estimate f, in 859 which case f is likely to be underestimated since dmin < dϕ with high probability (Appendix A; 860  problem. To improve our estimate for MFT vesicles, we performed a nonblind cap detection 868 in planar sections of ET z-stacks to obtain ϕ < ϕcutoff (Figures 6 and S9).  (Table 4) with a large estimation error (their Table 3 reference and lookup sections were equally probable, since any bias due to lost caps on the 916 reference section cancelled that on the lookup section. However, nearly all disector analyses 917 are performed sequentially, with a blind particle detection on the reference section followed 918 by nonblind particle detection on the lookup section, in which case there is an increased 919 probability of identifying caps on the lookup section (a scenario proposed by Hedreen 27 ). 920 When such an asymmetrical bias was added to our simulations, there was an 921 underestimation error of λ3D that was large even for small degrees of bias. Using the blind-922 versus-nonblind bias in particle detection measured from our analysis of MFT vesicles (ϕbias 923 = 20°), for example, we found a large underestimation error (Δλ3D = -40%; Figure 12). 924 Interestingly, an underestimation counting error of -15% for adjacent reference and lookup 925 sections has been previously reported and attributed to lost caps 61 . To remove the blind-926 versus-nonblind bias in the disector method, Hedreen 27 suggests using a third section 927 immediately below the reference section to guide the identification of caps in the reference 928 section, i.e. a nonblind-nonblind particle detection. 929 Besides the potential underestimation error due to the blind-versus-nonblind bias in 930 identifying caps, the low particle count of the disector method makes it inherently less 931 accurate in estimating λ3D compared to the ϕ-correction method. Our simulations indicate 932 that, given the same number of particles per projection, the disector method is ~2 to 3-fold 933 less accurate than the ϕ-correction method due to the smaller counts per section. To get the 934 same level of accuracy one would need to increase the cross-sectional area of the section 935 (or ROI) by ~4-fold. Hence, the ϕ-correction method for estimating λ3D is potentially more 936 efficient and accurate than the disector method. 937 Finally, for particles with a high density, especially those that are touching one another (e.g. 938 cerebellar GCs and MFT vesicles), the disector method is not recommended 18 . In this case, 939 the ϕ-correction method can be used in conjunction with planar sections. That said, it is 940 important to keep in mind that the ϕ-correction method requires an accurate estimate of ϕ, 941 which is not always possible for a given particle and imaging technique, and is designed for 942 particles with a spherical geometry. it assumes all particles are the same size (Appendix A). The use of dmin (or hmin) to correct 956 for lost caps is also problematic since this measure is susceptible to being an outlier (e.g. a 957 single false-positive measurement). Moreover, those using the Abercrombie and Floderus 958 corrections typically use approximate measures of μD, potentially adding an additional bias. 959 As discussed above, the design-based disector method has been deemed unbiased at its 960 conception 52 ; however, our analysis of this method has confirmed that a blind-versus-961 nonblind bias in particle detection, as proposed by Hedreen 27 , leads to an underestimation 962 of λ3D. Moreover, biases due to inaccurate measures of section thickness (T) and z-axis 963 tissue distortions can lead to significant biases in estimated λ3D for both design-based and 964 model-based methods, as discussed above. Hence, it is not surprising that comparisons of 965 estimated λ3D between model-based and design-based methods show large discrepancies. 966 Our analysis supports the conclusion that neither model-based or design-based methods 967 should automatically be assumed unbiased 25,26 and both methods need to be 968 verified/calibrated, preferably via 3D reconstructions 24,28,29,51 . Yet, 3D reconstructions should 969 not automatically be assumed to be the 'gold standard' due to potential z-axis distortions of 970 tissue sections, as discussed above. 971 Here, we validated a model-based ϕ-correction method for estimating λ3D that has been long 972 overlooked since the 1980s 50 , most likely since design-based methods have become the 973 defacto tools in modern stereology. Results of the validation show the ϕ-correction method 974 can estimate λ3D with high accuracy. A high accuracy is achieved via a superior model of 975 the lost-cap distribution, i.e. the Keiding model, that represents the mean cap-angle limit of 976 a population of spherical particles. Moreover, use of an LSE routine (or MLE routine) allows 977 one to use all measured 2D diameters (i.e. G(d)) to estimate ϕ, which is a significant 978 improvement to using dmin, a single measure that is likely to be an outlier, and also gives an 979 accurate estimate of μD. Comparison of estimated λ3D computed via the ϕ-correction method 980 to that computed via the Abercrombie and Floderus correction methods show biases 981 (underestimations) in the latter corrections by as much as ~20% for our MFT vesicle dataset 982 (Table 6). 983 To estimate particle size, we used outlines to compute the cross-sectional area of a particle's 984 projection in 2D images, which is equivalent to a high-resolution design-based point-grid 985 method, since pixels define a grid 70 . Because particles are never perfect spheres, cross-986 sectional area is a better measure of 2D size than the commonly used diameter line-987 segment measures dlong and dshort (Table 5) and 5% CO2, 325 mOsm) using a Leica microsystems vibratome (VT1200S) as previously 1010 described 3 . Sections were allowed to recover in high-sucrose ACSF at 35°C for 30-45 min, 1011 then in normal ACSF (125 mM NaCl, 25 mM NaHCO3, 25 mM D-glucose, 2.5 mM KCl, 1.25 1012 mM NaH2PO4, 2 mM CaCl2, and 1 mM MgCl2, equilibrated with 5% CO2 and 95% O2) at 1013 room temperature (~23°C). 1014 To prepare the cerebellar sections for high-pressure freezing, sections were heated to 37°C 1015 for 5-10 minutes in ACSF, then mounted into a sample 'sandwich' on a table maintained at 1016 37°C. The sample sandwich was assembled by placing a 6 mm sapphire disk on the middle 1017 plate of a transparent cartridge system, followed by a spacer ring, the section, a drop of 1018 ACSF containing 15% of polyvinylpyrrolidone for cryoprotection and adhesion, another 1019 sapphire disk and finally a spacer ring. The sample sandwich was frozen via a Leica EM 1020 ICE high-pressure freezing machine. 1021 Freeze substitution of the frozen samples was performed using an AFS1 or AFS4 Leica 1022 system equipped with an agitation module 76 . While in liquid nitrogen, frozen samples were 1023 transferred from storage vials to freeze-substitution vials containing 0.1% tannic acid and 1024 acetone, previously frozen in liquid nitrogen. Vials were transferred to the AFS1/AFS4 1025 system and shaken for 22-24 hours at -90°C. Inside the AFS1/AFS4 system, samples were 1026 washed for 10 minutes in pre-chilled acetone at -90°C for 3-4 repetitions. Next, a contrasting 1027 cocktail with 2% osmium and 0.2% uranyl acetate in acetone was chilled to -90°C and added 1028 to each vial. The temperature of the vials was kept at -90°C for 7-10 hours, raised to -60°C 1029 within 2 hours (15°C/hour), kept at -60°C for 3.5 hours, raised to -30°C within 4 hours 1030 (7.5°C/hour), kept at -30°C for 3.5 hours, raised to 0°C within 3 hours (10°C/hour), kept at 1031 0°C for ~10 min, then transferred to ice where samples were washed with acetone (3 × 10 1032 min). Samples were transferred from the vials to glass dishes containing acetone at room 1033 temperature and inspected for intactness and proper infiltration. Samples were washed with 1034 propylene oxide (2 × 10 min) and infiltrated with Durcupan resin at 2:1, 1: (±σ). During the fit routine, parameter T was fixed at its estimated value. The initial guess 1091 for ϕ was set to θmin, where θmin = sin -1 (dmin/μD) and dmin is the smallest non-zero diameter 1092 bin of G(d). Initial guesses for μD and σD were set to μ and σ of G(d), where μ of G(d) was 1093 computed as the sum of d⋅G(d)⋅dx over all bins and σ 2 was computed as the sum of (d -1094 μ) 2 ⋅G(d)⋅dx. For a small number of fits to the simulated G(d), usually for conditions of true ϕ 1095 > ϕcutoff, initial guesses had to be adjusted to get a successful fit. No parameters were 1096 constrained during the fits (e.g. 0 ≤ ϕ ≤ 90°) since testing of the LSE routine using a variety 1097 of datasets showed such constraints were never active or violated during the test fits. To 1098 validate the LSE routine, the MLE fits of Keiding et al. 49 were replicated, showing nearly 1099 identical results ( Figure S3; Table 4). Likewise, an LSE fit to Wicksell's G(d) of spleen 1100 corpuscles 16 resulted in an estimated F(d) that was nearly the same to that of Wicksell's 1101 unfolding solution ( Figure S4). 1102 The distribution of lost caps, L(d), was computed from G(d) via Equation 1 for T = 0 u.d. as 1103 follows: L(d, ϕ) = G(d, ϕ = 0°) -G(d, ϕ), where G(d, ϕ = 0°) and G(d, ϕ) were computed over 1104 the range d = 0-3 u.d. and G(d, ϕ) was normalised so that its last data point at d = 3 u.d. 1105 equaled that of G(d, ϕ = 0°). 1106 To compute λ2D from a 2D image using Fiji, a rectangular ROI was defined within a 1107 distribution of the particles of interest and two adjacent borders were designated as inclusive 1108 and the other two as exclusive ( Figure S1). Particles were counted if they touched the 1109 inclusive borders or were completely contained within the ROI, and not counted if they 1110 touched the exclusive border 78 . λ2D was computed as the particle count (N2D) divided by the 1111 ROI area. Using λ2D, μD and ϕ, λ3D was estimated via Equation 3, in which case it was 1112 important that λ2D was computed from the same image (or z-stack) from which μD and ϕ 1113 were estimated, since there was variation in μD and ϕ between sections. To compute the 1114 particle VF for a given λ3D and F(d), the following expression was used: 1115 The xy-square dimensions of the cuboid were adjusted to accommodate the required 1149 number of particles per projection, and the z-dimension was adjusted to accommodate the 1150 required number of projections. The particle VF = 0.40 unless specified. 3D particle 1151 diameters (D) were randomly drawn from a Gaussian distribution for a given F(d). 1152 Projections in the xy-plane were computed by identifying those particles with their center 1153 point located within a given section (interior particles), and those with their center point 1154 located above or below the section (caps) at a distance dz, where dz < D/2. The 2D projected 1155 diameter (d) was computed as d = D for an interior particle and d = (D 2 -4dz 2 ) ½ for a cap, 1156 derived from the trigonometric relation: (½D) 2 = (½d) 2 + dz 2 . A particle's cap angle was 1157 computed as θ = sin -1 (d/D). Other measures computed for each particle were the particle's 1158 distance from the section surface and the sum of overlaps within the xy-projection (Ω; circle-1159 circle overlaps expressed as a fraction) between the particle and particles higher in the 1160 section. To simulate lost caps, particles were excluded from the projection if θ was less than 1161 a fixed lower limit (ϕ), as in the Keiding model 49 . For a few simulations, however, ϕ was not 1162 fixed but variable, in which case particles were assigned a ϕ randomly drawn from a 1163 Gaussian distribution (μϕ ± σϕ) for a given CVϕ. To simulate the inability to observe particles 1164 deep in the section due to overlapping projections, particles were excluded from the 1165 projection if their Ω was greater than a fixed upper limit (ψ). To simulate the merging of 1166 circular projections for opaque particles, the xy-distance between two projections was 1167 computed according to the α-parameter of Hilliard 57 : α = (½d1⋅d2)(4d12 2 -d1 2 -d2 2 ), where d1 1168 and d2 are the projection diameters and d12 is the distance between the projection center 1169 points; if -1 < α < 0, then the projections were merged into one, resulting in a decrease in 1170 projection count and increase in projection size. The procedure for merging projections 1171 began with the particle closest to the section surface, which was then merged with other 1172 particles if -1 < α < 0. The merging procedure continued with the next particle closest to the 1173 section surface, and so on. To compute λ2D, the number of particles in a projection was 1174 divided by the geometry Areaxy. Because periodic boundary conditions were used in the 1175 simulations ( Figure S2), this λ2D is equivalent to one computed using inclusive/exclusive 1176 rectangular borders for counting (Figures S1). 1177 To simulate the disector method of computing density 18,52 , particles with a Gaussian F(d) 1178 were randomly distributed within a cuboid geometry whose xy-square dimensions were 1179 adjusted to accommodate ~500 particles per projection, and the z-dimension was adjusted 1180 to accommodate 100 sections with T = 0.3 u.d. To compute the density of a given section, 1181 a projection was computed for that section (the reference projection) as well as an adjacent 1182 section of equal thickness (the lookup projection). Particles were counted if they appeared 1183 within the reference projection but not the lookup projection. λ2D was computed as particle 1184 count per Areaxy and λ3D = λ2D/T. Simulations included parameter ϕ, i.e. particles were 1185 excluded from a given projection if their θ < ϕ. To simulate a blind-versus-nonblind bias in 1186 vesicle detection, ϕ was defined separately for the reference section (ϕref) and lookup section 1187 (ϕlookup) such that ϕref ≥ ϕlookup, with their difference (bias) defined as ϕbias = ϕlookupϕref. 1188 To quantify ϕcutoff, the estimation error Δϕ (estimated ϕ -true ϕ) was computed from curve 1189 where CVD is computed using true μD and σD. To investigate whether this expression can 1200 be used to test the accuracy of estimated ϕ, ϕcutoff was computed using estimated μD and σD 1201 of each simulation (rather than true μD and σD) and this 'estimated' ϕcutoff was compared to 1202 the corresponding estimated ϕ. Results showed that the use of estimated μD and σD to 1203 compute CVD translated into negative offsets in ϕcutoff, ranging from -17° to -10° for n = 200 1204 to 2000 diameters, respectively. To account for these offsets, diameters in the dcutoff matrix 1205 were adjusted to remove the offsets and the matrix was refit to the bivariate polynomial, 1206 resulting in: 1207 estimated ϕcutoff ≈ sin -1 (0.987 -2.071⋅CVD + 0.124/√n -35.059⋅CVD/√n) (9) 1208 where CVD is computed using estimated μD and σD. In conjunction with the fit error of ϕ, this 1209 expression was used as an accuracy test of estimated ϕ, i.e. estimated ϕ was considered 1210 accurate if it was less than estimated ϕcutoff. For the analysis of ET z-stacks (ET10 and ET11; 1211 where α is the maximum scan angle 80 . For this study, α = 65° in which case exz = 1.4. Hence, 1244 ρz = ρx⋅exz = 6.5 nm. For ET11, where Ttissue = 133 nm, ρx = 3.5 nm and ρz = 4.9 nm. Next, 1245 ρz was estimated from experimental data by curve fitting Equation 1 to G(d) computed from 1246 the MFT vesicle analysis (Figures 6D and S9D) while fixing μD, σD and ϕ to their 'true' values 1247 measured from the 3D analysis and leaving T as the one free parameter. Results gave 1248 estimated T = 5.3 ± 0.8 nm for ET10 and 0.2 ± 0.4 nm for ET11. Hence, estimates of T (i.e. 1249 ρz) of the ET z-stacks spanned 0-6 nm. Given the small range of estimated T, and 1250 comparatively large dimensions of the MFT vesicles, the ET analysis was simplified by 1251 assuming T = 0 nm. To test what effect this assumption might have estimates of μD, σD and 1252 ϕ, the curve fits in Figures 6D and S9D were recomputed assuming T = 6 nm for ET10 and 1253 T = 5 nm for ET11 and found ΔμD, ΔσD and Δϕ were similar to those for assuming T = 0 nm: 1254 To estimate vesicle density within a cluster for ET10, it was necessary to confine the density 1257 analysis to a subregion of the original z-stack along the axial axis since the vesicle density 1258 as a function of z-depth was nonhomogeneous, being smaller at the top and bottom of the 1259 stack ( Figure 11B). Within this subregion, we estimated λ2D by counting the number of 1260 outlines that fell within a ROI (Areaxy = 0.039 μm 2 ) using inclusive/exclusive borders (Figure 1261 S1C) at the center of the vesicle cluster for 10 z-stack images spaced 10-15 nm apart along 1262 the axial axis, giving λ2D = 304.2 ± 15.5 μm -2 (±SEM). There was an average of 12 vesicles 1263 per ROI, which is 7-fold larger than the theoretical optimal number of particles for the disector 1264 the VF estimated via the 2D analysis is similar to that estimated via the 3D analysis. 1273 To estimate vesicle λ3D for ET11 via the 'physical' disector method, a reference section with 1274 T = 12.8 nm (0.3 u.d.) was randomly located within the center of the z-stack and a 1275 corresponding adjacent lookup section with the same T was defined. Vesicles that appeared 1276 in the reference section (i.e. vesicles that had one or more of their 2D outlines from the 1277 nonblind analysis in Figure 6 appear in the reference section) but not the adjacent lookup 1278 section were counted and used to compute λ3D = count/(Areaxy⋅T), where Areaxy = 0.144 1279 μm 2 . This analysis resulted in ~20 vesicles per section, or λ3D ≈ 11,000 μm -3 . To simulate a 1280 bias between a blind reference vesicle detection and nonblind lookup vesicle detection 1281 (ϕbias), vesicle outlines from the nonblind analysis were used for the lookup vesicle detection 1282 (mean ϕlookup = 42°) and a copy of the same outlines for the reference vesicle detection, but 1283 modified to have a larger ϕ (ϕref = ϕlookup + ϕbias) by deleting the necessary number of extreme 1284 outlines from the negative and positive pole regions to achieve the desired ϕref. Setting ϕbias 1285 = 17° ( Figure 7A) resulted in ~14 counts per section, or λ3D ≈ 7000 μm -3 , with estimation 1286 error Δλ3D = -32%. To compare these results to those of the Monte Carlo disector 1287 simulations, results of 7-9 reference sections as just described were combined to give a 1288 total of ~500 vesicles per reference section for a given ϕbias, and an average Δλ3D was 1289 computed from 100 such reference sections. where noted (unpaired two-tailed equal-variance, unless specified differently; significant p < 1294 0.05; F-test used to verify equal variance). Errors reported in the text and graphs 1295 (bars/shading) indicate the standard deviation (±σ), except in a few instances they indicate 1296 the standard error of the mean (±SEM) which is noted. 1297 The estimation error (Δ) of parameters μD, σD, λ2D or λ3D was computed as the percent 1298 difference between a parameter's estimated (ε) and true (t) value [Δ = 100(ε -t)/t], except 1299 for ϕ, which was computed as a difference (Δ = ε -t) since division by ϕ caused distortion at 1300 small ϕ. For Monte Carlo simulations with multiple repetitions, the mean and standard 1301 deviation of a parameter's estimation error (μΔ ± σΔ) is referred to as the bias and (68%) 1302 confidence interval, respectively. For simulations with true ϕ < ϕcutoff, the Δ-distributions were 1303 typically normal. However, for simulations with true ϕ ≥ ϕcutoff, the Δ-distributions were often 1304 skewed (i.e. absolute skew > 0.5); in this case, μΔ was computed as the median of the Δ-1305 distribution and +σΔ and -σΔ were computed separately above and below μΔ. 1306 particles are the same size. If there is a distribution of particle size, then each particle will 1334 have its own minimum observable diameter (δmin) and dmin will be the minimum of all δmin 1335 ( Figures 6C1 and S9C1). Depending on the spread of the δmin distribution, the difference 1336 between dmin and the mean δmin can be large, in which case the dmin correction will create a 1337 large underestimation of λ3D. Moreover, there is a high probability that dmin is an outlier, i.e. 1338 an unusually small value; this could happen when conditions for identifying small caps are 1339 better than average, or when there are false positive measurements. 1340 Our 3D analysis of MFT vesicles showed the Keiding model accurately describes the δmin-1341 D relation of the vesicles (Figure 6C1 and S9C1) and therefore gives an accurate estimate 1342 of λ3D via Equation 3 (Figure 11). This is in contrast to the dmin correction which 1343 underestimated λ3D for the same datasets by 13-20% and the Abercrombie correction which 1344 underestimated λ3D by 21-23% (Table 6; Figure S3). Note, estimated ϕ > 1726 ϕcutoff, where ϕcutoff was computed via Equation 9 using estimated μD and σD (n = 500). Hence, 1727 ϕ is considered indeterminable. 1728 1729 GC soma Confo +15% *** -13% *** -0% +1% * GC nucleus TEM +17% *** -21% *** -5% *** -2% *** MFT vesicle TEM +2% ** -12% *** -5% *** -5% *** MFT vesicle ET11 +0% -20% *** -11% *** -10% *** Confo: confocal. Percent difference (Δ) measured with respect to darea. Significance 1731 measured via paired t-tests in comparison to darea (*p < 0.05, **p < 0.01, ***p < 0.001). There 1732 were 50 particles per sample. 1733  Tables 2 and 3. θmin = sin -1 (dmin/μD), used in Equations 2 and 1736 3 as a substitute for ϕ. Right column: Abercrombie correction 17 , which assumes no lost caps 1737 (θmin = 0°). Δ = 100[X(θmin) -X(ϕ)]/X(ϕ) where X = ζ and λ3D. For vesicles, comparisons are 1738 with respect to 'true' ζ and λ3D computed via 3D reconstructions. See Appendix A. 1739 For thick sections and a high particle density, there are usually projection overlaps that make 1758 counting/outlining the projections more difficult; moreover, AF > VF, a condition known as 1759 overprojection. 1760