A rapid workflow for neuron counting in combined light sheet microscopy and magnetic resonance histology

Information on regional variation in cell numbers and densities in the CNS provides critical insight into structure, function, and the progression of CNS diseases. However, variability can be real or a consequence of methods that do not account for technical biases, including morphologic deformations, errors in the application of cell type labels and boundaries of regions, errors of counting rules and sampling sites. We address these issues in a mouse model by introducing a workflow that consists of the following steps: 1. Magnetic resonance histology (MRH) to establish the size, shape, and regional morphology of the mouse brain in situ. 2. Light-sheet microscopy (LSM) to selectively label neurons or other cells in the entire brain without sectioning artifacts. 3. Register LSM volumes to MRH volumes to correct for dissection errors and both global and regional deformations. 4. Implement stereological protocols for automated sampling and counting of cells in 3D LSM volumes. This workflow can analyze the cell densities of one brain region in less than 1 min and is highly replicable in cortical and subcortical gray matter regions and structures throughout the brain. This method demonstrates the advantage of not requiring an extensive amount of training data, achieving a F1 score of approximately 0.9 with just 20 training nuclei. We report deformation-corrected neuron (NeuN) counts and neuronal density in 13 representative regions in 5 C57BL/6J cases and 2 BXD strains. The data represent the variability among specimens for the same brain region and across regions within the specimen. Neuronal densities estimated with our workflow are within the range of values in previous classical stereological studies. We demonstrate the application of our workflow to a mouse model of aging. This workflow improves the accuracy of neuron counting and the assessment of neuronal density on a region-by-region basis, with broad applications for studies of how genetics, environment, and development across the lifespan impact cell numbers in the CNS.

1. Introduction

Accurate counts of neurons are essential metrics for understanding the structure and function of the brain (Bonthius et al., 2004; Herculano-Houzel and Lent, 2005). This is particularly important for pre-clinical studies of human disease, as it provides a necessary starting point for assessing neurodegenerative changes on a regional basis (Price et al., 2001; Rosen and Williams, 2001; Bandeira et al., 2009). In addition to absolute counts, measuring neuron density can provide important insights into the structural complexity of local circuits. The distribution of neurons in the mouse brain on a region-by-region basis is also crucial for understanding variation between individuals and strains. The density of neurons can reflect the distinctive cytoarchitecture of brain regions, including laminar organization, size and shape of constituent neurons, and the volume and composition of associated neuropil (Wree et al., 1982; Spocter et al., 2012; Kasthuri et al., 2015). Furthermore, assessing neuronal densities in targeted brain regions under varying conditions may reveal impacts that change the composition of neuropil with or without associated loss (or gain) of neuronal cell bodies (Amunts et al., 1996; Selemon and Goldman-Rakic, 1999). The heterogeneity of neuronal density within a given region can provide insights into the complexity, developmental history, and functional diversity of the region (Peters et al., 1991). Establishing such datasets enables pre-clinical studies of the impact of genetics, environment, and development across the lifespan on brain structure.

To investigate cells numbers and densities in CNS, researchers have developed various quantitative methods, including several stereological approaches (Williams and Rakic, 1988; West, 1999; von Bartheld, 2002; Schmitz and Hof, 2005; Deniz et al., 2018) and approaches that first homogenize brain tissue and dissociate brain cells (Herculano-Houzel and Lent, 2005; Collins et al., 2010; Young et al., 2012). Stereology is a quantitative method for estimating the three-dimensional characteristics of biological structures using two-dimensional histological sections or images and systematic random sampling within delineated brain regions. In traditional approaches to stereology, the specimen is preserved by means of chemical fixation, the brain is removed from skull, sectioned into thin slices (usually 5–50 μm thick), and stained with specific dyes to highlight different cell types. The sections are then viewed under a microscope, and measurements are taken to estimate the cell numbers within delineated brain regions. Both 3D counting and the optical fractionation are alternatives to traditional stereological methods introduced by Abercrombie (1946, reviewed in Rosen and Williams, 2001) that is design-based in its approach to sampling a volume of tissue, counting targeted objects (e.g., neurons), and generating statistics that assess the precision of the counts (West, 1999). These methods are considered unbiased, but still can suffer from bias due to differential z-axis shrinkage and variable penetration of stains. By providing a less biased estimate of structural parameters, these modern stereological methods have contributed significantly to our understanding of many biological systems. However, there are some limitations. Stereological methods require extensive preparation of tissue and the use of microtomes to produce thin tissue sections. This can introduce damaging artifacts, lost sections, and distortions, even in the hands of experts. The results are sensitive to a range of factors, including the choice of chemical fixation, the post-fixation treatment of the specimens, the specific histological protocols used, and the specific features of the design-based approach to sampling and counting. Thus, it has be difficult to determine the accuracy of the measurements obtained in tissue that was subject to such deformation and the range of operational variables employed. Finally, stereological analysis can be a laborious and time-consuming process (Bonthius et al., 2004). When multiple brain regions are analyzed, the time and effort required can increase significantly, and the manual aspect of the process can introduce additional sources of variability.

In contrast, the isotropic fractionator (IF) involves homogenizing a small tissue sample from a region of interest—or the entire brain —and processing it into a uniform suspension of cells. The cells are then stained with a dye and placed into a special counting chamber, where the number of cells in a known volume is counted under a microscope (Herculano-Houzel and Lent, 2005), or by the flow cytometry (Collins et al., 2010; Young et al., 2012). The total number of cells in the brain region can then be estimated by extrapolating the cell density from the counting chamber to the reported or measured volume of the entire brain region. The downsides to isotropic fractionator are similar to stereology. One limitation is that both IF and stereology use the brain which is taken out of the skull and processed, introducing swelling and/or shrinkage, which can affect the accuracy of volume estimation.

Light sheet microscopy (LSM) has emerged as a powerful tool for 3D visualization of targeted populations of cells in intact, cleared tissue samples using fluorescent immunocytochemistry. Recent advances in tissue clearing (Ueda et al., 2020) have enabled the capture of high-quality, countable images in the whole-brain using LSM (Hillman et al., 2019). These intact 3D volumes can then be used for accurate counting of neurons and assessing neuronal densities, providing a more comprehensive and representative view of the brain region of interest. Typically, the whole brain images are registered to a standard reference atlas, such as the Waxholm Space Atlas (WHS) (Johnson et al., 2010) or the Allen Brain Atlas (ABA) (Wang et al., 2020). The combination of the digital images and the standard reference atlas facilitates the segmentation of regions and correction of altered brain morphology, allowing for comparison of neuron numbers and density within and across studies.

Researchers have employed this approach of registering LSM to WHS or ABA atlases in a number of recent studies (Susaki et al., 2014; Menegas et al., 2015; Renier et al., 2016; Zhang et al., 2017; Krupa et al., 2021). Susaki et al. (2014) mapped LSM to WHS space but did not proceed with cytometric quantitative analysis. The remaining studies mapped the LSM images to ABA space for analysis. ClearMap, as described by Renier et al. (2016), utilizes a peak detection algorithm with a threshold determined by comparing manual and machine counting, which may not yield optimal results when applied to complex and heterogeneous neuron distributions. Zhang et al. (2017) employed L1 minimization (Lasso), a machine learning regression method, to detect neurons in the whole brain. However, this method assumes sparsity of objects and a linear relationship between input features and output variables, with all features equally important, which may not hold true in the case of heterogeneous neuron distributions. Krupa et al. (2021) used a 3D Unet for cell detection, but these deep learning methods require large inputs and considerable training time compared to other machine learning methods. Menegas et al. (2015) utilized Ilastics¹ to segment cells in eight brain regions, with separate training for each region.

Registration of LSM images obtained from mouse brains to ABA space for correcting deformation and quantitatively segmenting brain structures and regions has one fundamental deficiency. The ABA atlas was generated from multiple specimens that were imaged after the brains were removed from the skulls (Wang et al., 2020). The lack of skull support and the effects of processing will deform the brain. Furthermore, ABA assembled their atlas from 1,600 animals, which were registered into a volume that does not conform to the size and morphology of the in situ mouse brain—as we demonstrate in this work.

We present a novel approach that is highly reproducible and allows efficient counting across the whole brain. Our method utilizes magnetic resonance histology (MRH), an extension of MRI to microscopic resolution of fixed tissue specimens (Johnson et al., 1993). MRH images are acquired with the brain inside the skull, resulting in a representation that more closely approximates the size and morphology of the in vivo brain (Johnson et al., 2023), as compared to reconstructions from histological methods that require cranial dissection, serial sectioning of the brain, and chemical treatment of thin brain sections. As a result, MRH images can serve as the “gold standard” for correcting ex situ brain morphology. Our workflow offers a powerful combination of morphological correction based on MRH, machine learning-based neuron classification, and post-processing. It provides accurate and detailed neuron counts and calculations of neuronal density, which can be readily derived from different brain regions and specimens. Our method features a large, well-defined field of view in 3D, which allows for a comprehensive analysis of an entire subvolume under study. Moreover, our workflow applies the principles of optical fractionation to systematically sample and count neuron numbers in 3D volumes, which reduces the computational costs and represents a novel digital approach not previously explored.

2. Materials and methods

2.1. Specimens

Two groups of animals were used to test the methods (see Table 1). The first group included 5 C57BL/6 J mice (4 male and 1 female) that were sacrificed at 90 ± 3 days to test the consistency of the counts and provide measures that could be related to existing literature. A second experiment with BXD89 mice included two male specimens at 111 and 687 days to test the sensitivity of the method to changes in neuronal density arising from aging. The mice were obtained from independent litters. The BXD strains are a set of well-characterized recombinant inbred mouse strains, making them a valuable tool for systems genetics studies (Ashbrook et al., 2021).

TABLE 1

Table 1. Test specimens for neuron counting.

2.2. MRH, label map and LSM scanning

All animal procedures were carried out in accordance with guidelines approved by the Duke Institutional Animal Care and Use Committee. Specimens were perfusion-fixed using an active staining method as previously described (Johnson et al., 2019). Briefly, warm saline was perfused through a catheter in the left ventricle to exsanguinate, followed by ~5 min of perfusion (50 mL) of a mixture of 10% buffered formalin and 10% Prohance (Gadoteridol) to reduce the spin lattice relaxation time (T1) of the tissue for accelerated scanning. MRH scanning was performed on a 9.4 T vertical bore magnet with a Resonance Research, Inc., gradient coil yielding peak gradients up to 2,500 mT/m, controlled by an Agilent console running VnmrJ 4.0. The MRH scanning was performed using a Stesjkal Tanner spin echo sequence with b values of 3,000 s/mm² and 108 angular samples spaced uniformly on the unit sphere. Compressed sensing (Lustig et al., 2007) was used with a compression factor of 8X (Wang et al., 2018), resulting in a large (252 GB) 4D volume with isotropic resolution of 15 µm (Johnson et al., 2023).

The label map in this study is based on a modified version of the Common Coordinate Framework (CCFv3) (Johnson et al., 2023) from ABA (Wang et al., 2020). The CCFv3 atlas includes 461 carefully curated regions of interest (ROIs). Our workflow relies on an initial mapping of labels from our canonical MRH atlas of a 90-day male C57BL/6 J mouse to the MRH of the specimen under study using a pipeline built around Advanced Normalization Tools (ANTs) (Avants et al., 2011; Anderson et al., 2019). Many of the CCFv3 ROIs are less than 1 mm³, including numerous subdivisions of cortical areas (laminae) and subcortical structures (subnuclei) that are of limited value within the present scope of our analyses. Accordingly, we reduced the label set to 360 ROI (180 per hemisphere) gray matter and white matter structures. This modified atlas, which we refer to as the reduced CCFv3 (r1CCFv3), consolidates these smaller subregions in CCFv3. The r1CCFv3 provides a label set that registers to the MRH volumes reproducibly enabling precise neuron counting within a large number of structures. The details of this registration process are described in Johnson et al. (2023).

Following the MRH scans, the brains were carefully removed from the skulls. Paraformaldehyde-fixed samples were preserved with using SHIELD reagents (LifeCanvas Technologies) using the manufacturer’s instructions (Park et al., 2019). Samples were delipidated using Clear+ delipidation reagents. Following delipidation samples were labeled using eFLASH technology which integrates stochastic electrotransport (Kim et al., 2015) and SWITCH (Murray et al., 2015). The samples were then washed in PBS for 7-8 h before overnight fixation in 4% paraformaldehyde followed by incubation in secondary labeling buffer at 37°C with two refreshes over the course of 7–8 h before secondary labeling in the SmartLabel device. For each brain, 10 μg of rabbit anti-Iba1 (Cell Signaling Technologies 17198S*) primary antibody, 10 μg mouse anti-NeuN (Encor MCA-1B7), 6 μg rabbit anti-NeuN (Cell Signaling Technologies 24307S*), or 20 μg mouse anti-MBP (Encor MCA-7G7). Secondary antibodies were used at a 2:1 Secondary:Primary molar ratio. After immunolabeling, samples were incubated in 50% EasyIndex (RI = 1.52, LifeCanvas Technologies) overnight at 37°C followed by 1 d incubation in 100% EasyIndex for refractive index matching. After index matching the samples were imaged using a SmartSPIM axially-swept light sheet microscope using a 3.6x objective (0.2 NA) (LifeCanvas Technologies, Cambridge, MA). The processing streams for MRH and LSM are depicted in Figure 1A.

FIGURE 1

Figure 1. Overview of the workflow for assessing the density of neurons. (A) The mouse brain is imaged using two modalities: MRH imaging while the brain is in the skull, followed by LSM after the brain is removed from the skull and subjected to tissue clearing. (B) The LSM data are pre-processed by registering to MRH correcting the deformation in brain morphology. (C) The automated workflow locates the region with the label from r1CCFv3 and generates random subvolumes within that region to sample, applying the design-based principles of optical fractionation. Neurons in each subvolume are identified via a random forest algorithm followed by 3D watershed and volume filters and counted.

2.3. Automated neuron counting

2.3.1. Overview of the algorithm

To address the significant distortion introduced during skull removal and chemical processing in the preparation of LSM samples, the LSM images are first registered to the MRH volumes of the same specimen (Figure 1B; Tian et al., 2023). Section 3.1, accompanied by Figure 2, provides a detailed explanation of the impact of correcting the morphology and extraction of the cytometric statistics. We examine the effects of the preprocessing steps on the neuron density measurements, providing insights into the importance of these steps for accurate data analysis.

FIGURE 2

Figure 2. Effect of morphology correction on light-sheet datasets. (A) A single slice from a whole-brain LSM before correction. (B) The same slice as in (A) after correction using MRH. Magnified axial views of auditory areas and hippocampus before (C) and after (D) correction, and coronal views of the anterior temporal and posterior diencephalic region (showing the amygdaloid complex) before (E) and after (F) correction. The specimen is the same on both sides, and the colormap is used to illustrate the contrast between before and after correction.

The workflow is visualized in Figure 1C. After warping the LSM and label map to the MRH space, the user selects a region for analysis by inputting the corresponding label from the label map into Fiji. This generates a surface defining the region of the selected label, as well as a group of subvolumes within the surface. These subvolumes are distributed in a uniform but random manner, which improves the accuracy of neuron counting and facilitates analysis of the heterogeneity of neuron distribution throughout the selected region. Please note that subvolumes located on boundaries and containing broken tissues will be automatically excluded by screening the label value of the subvolume. For each subvolume (Figure 1C), we apply a lightly supervised algorithm called random forest neuron segmentation, which utilizes an ensemble of decision trees to generalize the projection between graphical features and labels. This algorithm is commonly used in image processing tasks and requires only a small amount of training data to perform well (Pavlov, 2000). Prior to segmentation, the random forest classifier is pre-trained using training data consisting of binary signatures, where neurons are labeled as 1 and non-neurons as 0, as well as feature vectors extracted from various visual characteristics of the image, such as texture, shape, and size. The training process optimizes the relationship between graphical features of an image and the assigned neuron labels. The labeling example can be seen in Supplementary Figure S1. In general, labeling approximately 20 neurons with both the cell body and background enables the algorithm to learn the essential features of neurons, resulting in a commendable F1 score (Supplementary Section S2) of around 0.9, as illustrated in Supplementary Figure S2. After the random forest segmentation is applied to each subvolume, the workflow will perform post-processing on the objects classified as neurons. This post-processing uses 3D watershed (showcased in Supplementary Figure S3) and volume filters to further isolate the neurons as discrete objects. The final count of the neurons in each subvolume is generated based on these refined objects. The density is calculated by averaging densities obtained from subvolumes, and the standard deviation of the density provides insight into the variability of the neuron distribution within the region.

2.3.2. Big data environment

Whole brain LSM images are large – typically ~300 GB. To enable efficient, interactive development and execution of the workflow we have assembled a hardware/software environment suited for these large data. The pipeline has been implemented on two high performance Dell servers: Dell E5-2400, NVIDIA Tesla V100 GPU and a Dell E52670, NVIDIA V100 GPU. Each server has 1.5 TB of memory.

Several software packages have been integrated into the big data environment. Imaris² is a commercial software package designed for interactive 3D/4D viewing of large microscopy data sets. It accommodates simultaneous visualization of multiple light sheet volumes using a hierarchical data format. Imaris has been developed to allow memory sharing with external packages. One of these is Fiji,³ a distribution of Image J that has been developed for applications such as that envisioned here through the extensive use of plugins. One crucial plugin for our work is BigDataViewer (Pietzsch et al., 2015), an open source solution for accommodating large volumes in Fiji and supporting plugins for post processing. This includes LabKit (Arzt et al., 2022) a user friendly Fiji plugin for microscopy segmentation using the random forest algorithm.

2.3.3. Compiled algorithms

To improve the accuracy and efficiency of neuron counting, we developed two pipelines for different purposes: one for visualized validation and the other for production use. The initial pipeline (available on GitHub⁴) provides visualization of classification and segmentation performance. This algorithm is built on the coding interface of the image rendering software, Imaris.⁵ The procedures, from applying classifiers to obtaining statistics, are executed in the GUI of Imaris. Although the workflow cannot be fully automated through the Imaris coding interface, it can be automated through recording applications such as OS Automator.⁶ The second algorithm (available on GitHub⁷) provides an automated workflow constructed using Python and macros. This method utilizes Fiji macro GUIs to automate segmentation, 3D watershed, and region counting plugins. The watershed uses morphological erosion to identify the center of each object, followed by calculating a distance map from the object center points to the object edges. The resulting topological map is then filled with imaginary water, and dams are built at locations where two watersheds meet to separate them. The region counting method assigns unique labels to each unconnected component to facilitate subsequent counting. The plugin used for constructing the random forest model is Labkit.⁸ The watershed method and the 3D connected region counting method in the first algorithm are available in Imaris. The 3D watershed and 3D connected region counting method in the second algorithm are the Fiji/binary/watershed function and from the open-source plugin MorphoLibj.⁹

2.3.3.1. Imaris pipeline for visualization

The user identifies the structure to be analyzed and loads the surface defining that structure. A Python script generates a collection of 3D subregions that are randomly placed within the volume of the structure defined by the surface. Each counting subregion is a 100 μm × 100 μm × 100 μm cube. To avoid oversampling or undersampling in different regions and to accurately capture the heterogeneity, we analyzed 15 subvolumes per brain region while keeping the number of subvolumes consistent across regions to facilitate statistical analysis. The user selects “Labkit” in the Fuji extension. The GUI displays the subregion’s volume, and the user either imports a pre-trained classifier or starts labeling neurons and background with binary labels. Supplementary Figure S1 shows representative images of the training process. Once a classifier has been trained, it can be saved for future use in other ROI or other specimens, with some limitations discussed below. The classification is then sent to Imaris. Imaris applies thresholding to the binary labels and generates individual surfaces based on the results. The watershed algorithm uses seed points to split touching surfaces. After approximating the size of neurons through visualization, the volume filter is applied, and the count is generated.

2.3.3.2. Fiji pipeline for high throughput

A Python script generates sub-regions inside the surface of the structure being analyzed and saves these subvolumes as individual TIFF files. The same Python script creates and executes a FIJI macro that reads the TIFF files in batches, applies the pre-trained random forest classifier, and performs 3D watershed. Components in the resulting image are labeled with different colors using the FIJI function “connected components labeling” and saved as a TIFF file. A second Python script reads these TIFF files in batches, applies volume filters, and generate counts of the components. The details of the volume filters are provided in Supplementary Section S6.

2.4. Data availability statement

Original datasets are available in a repository. The access will be granted upon request.

3. Results

3.1. Visualization of morphology correction

Figure 2 displays the comparison between images from an uncorrected and deformation-corrected whole brain light sheet volume. Much of the existing literature for stereology has employed optical light microscopy and confocal microscopy with higher resolution and smaller field of view than that in the whole brain LSM images. However, the deformation, which is dependent on the fixation and staining method, is likely to be similar to that presented in Figures 2A,C,E. The compensation of this deformation is usually done through the application of an isotropic correction factor, but this factor may be insufficient to fully address the nonuniform deformation arising from tissue processing. Figure 2B shows the LSM image with morphology correction by MRH, which significantly reduces the irregular distortion present in the raw LSM images. The magnified views of the posterior cerebrum (showing the auditory cortex and hippocampus (Figure 2D) and LGN (Figure 2F)) illustrate the benefits of the morphology correction method in improving the quality of imaging data and segmentation label maps, compared to Figures 2C,E. Supplementary Figure S4 illustrates the differences in neuron density between raw LSM and deformation corrected LSM for 6 brain regions. The smallest percent difference (calculated as the difference between raw and corrected density over the corrected density), observed in the primary visual cortex, is 7%, while the largest, observed in the subiculum, is 64%. These results demonstrate that deformation is not uniform throughout the brain.

Figures 3A,C,E displays the label-map CCFv3. Data for the ABA were obtained by averaging 2-photon microscopy images of 1,600 animals. The distortions from tissue handling are not consistent across all 1,600 animals. Wang et al. (2020) report variability in the volume of the regions in their atlas ranging from 6 to 80% so tissue handling introduces variation. Figures 3B,D,F displays the r1CCFv3 labels with morphology correction by MRH. The resulting corrected label map provides more accurate representation of brain regions (Figure 3).

FIGURE 3

Figure 3. Comparison of similar levels from the ABA (A,C,E) and a single specimen imaged by MRH. (B,D,F) The colormaps are different because the ABA atlas is shown with the full complement of 461 (CCFv3) labels and the MRH is shown with the reduced (r1CCFv3) labels in which some ROI have been combined.

3.2. Comparison of machine counting and human counting

To assess the accuracy and precision of our workflow, we chose five brain regions with cell densities we determined would be countable: Dorsal part of the lateral geniculate complex (LGd), Auditory area (AUD), retrosplenial area (RSP), orbital area (Orb), and subicular region (SUBR). We created 15–20 3D subvolumes for sampling within these regions with dimensions of 1,000 μm × 1,000 μm × 12 μm (556 × 556 × 3 voxels) for visualization in Figures 4A,B, and 100 μm × 100 μm × 100 μm (56 × 56 × 25 voxels) for routine implementation of the workflow. The random forest classifier, trained on a random sub volume in the auditory cortex of a specimen not involved in the experiment, was applied automatically to identify the objects to be included. We randomly selected 10–15 subvolumes within the same datasets, and an experienced researcher manually counted the neurons within them by labeling neurons in subvolumes across 3–5 specimens. We compared the results of both methods and found that both manual counting and the machine had overall performance that was comparable (Figures 4A,B). One consideration to note is that during the manual labeling process, the researcher was involved in selecting the subvolumes, which introduced potential bias. Although efforts were made to minimize bias and select subregions randomly, there may have been a tendency to choose subregions with visually distinguishable neurons, leading to an uneven distribution. That means the manual classification has a tendency to underestimatesthe heterogenieoty of the neuron distribution. However, in the subsequent machine classifier analysis, subvolumes were automatically generated without bias. As a result, the machine counts in this panel exhibit greater variability in comparison. A statistical comparison is shown in Figure 4D.

FIGURE 4

Figure 4. Comparison of neuron density between the machine workflow (A) and manual counting (B). Note that both the labeling will go through watershed illustrated in (C) to ensure the connected surfaces are separated. Standard deviation is shown in (D), representing statistical variations across five different 90-day C57BL/6J specimens.

Supplementary Table S2 presents a comparison between the proposed method and other cell counting software solutions.

3.3. Neuron counting in C57BL/6J mouse brain

We demonstrated the use of the pipeline in counting neurons in 13 different brain regions across 7 specimens. We choose representative regions to demonstrate the capability of our workflow in neocortex, hippocampus, amygdala, thalamus, and brainstem. These regions include (1) neocortex: orbital area (Orb), auditory area (AUD), retrosplenial area (RSP), primary visual area (VISp); (2) hippocampal region: entorhinal area (ENT), subicular region (SUBR), field CA1 (CA1), field CA3 (CA3); (3) amygdala: basolateral amygdalar nucleus (BLA); (4) thalamus: dorsal part of the lateral geniculate complex (LGd), and the entire thalamus (Th); (5) brainstem: facial motor nucleus (VII), principal sensory nucleus of the trigeminal (PSV). Initial measurements were performed on 5, 90-day old male C57BL/6J specimens. The results show that neuron density in some regions (LGd, AUD, RSP, Orb) is quite consistent across specimens, and much more variable while there is a high standard deviation in others (PSV, VISp, TH, BLA), as indicated by the relatively high standard deviation of the means (Figure 5). Supplementary Figure S5 and Supplementary Table S1 present the comparison between the neuron density obtained from our workflow and previous studies.

FIGURE 5

Figure 5. Demonstration of the variation of neuron density (A) and neuron number (B) across brain regions, as well as variations across specimens. Error bars are standard deviations of the means. The dots indicate neuron counts from individual specimens. Region: LGd: Dorsal lateral geniculate nucleus, AUD: auditory area, RSP: retrosplenial area, Orb: orbital area, SUBR: subiculum, VII: facial motor cortex, PSV: trigeminal, ENT: entorhinal area, CA1: field CA1, CA3: field CA3, VISp: primary cortex, TH: thalamus, BLA: Basolateral amygdala.

To assess the variability of our measurements, we calculated the coefficient of variation (CV) for the density and number of neurons in several brain regions across different specimens (Supplementary Table S3). The CV for neuron density ranged from 0.062 in CA1 to 0.333 in VII, with an average of 0.139. The CV for neuron number ranged from 0.06 in CA1 to 0.44 in PSV, with an average of 0.17. These results indicate that there is considerable variation in neuron density and number across individual animals, even within the same brain region. Notably, the regions with the highest CV for neuron density and number were VII (0.333) and PSV (0.44), respectively. This is likely to be induced by the small sizes of these regions. The higher coefficient of variation for small regions may be due to minor displacements between the placement of individual delineations from the labelmap and the actual neuroanatomical structures, as discerned in LSM, especially for brainstem regions where deformation was significant. Conversely, the regions with the lowest CV for neuron density and number were CA1 (0.062), RSP (0.073), AUD (0.076), Orb (0.092), CA3 (0.091), and SUBR (0.097), respectively. These findings underscore the importance of accounting for inter-individual variability when analyzing neuronal parameters and highlight the need for careful consideration of sample size and statistical power in studies of brain structure.

3.4. Illustration of aging effects on the neuron density and number

Figure 6 presents a comparison between one young (111 day) and one old (687 day) BXD89 specimen to demonstrate the application of our workflow to a different strain and its potential to be used to examine the impact of aging on neuron density. It is important to note that only two animals are being assessed in this Figure and the standard deviation thus reflects the heterogeneity of neuron distribution within each brain region, which was obtained by analyzing mean neuronal counts within the subvolumes. Our analysis revealed that in certain regions, such as the orbital area (Orb), the observed decrease in neuron density was primarily due to a reduction in the total number of neurons, while in other regions such as the entorhinal area (ENT) and basolateral amygdala (BLA), the observed decrease in density was largely attributable to regional enlargement of the brain during aging. These findings underscore regional variation and the complexity of the aging process and highlight the need for careful consideration of volume and cellular number when studying brain development across the lifespan.

FIGURE 6

Figure 6. Comparison of the neuron density (A) and number (B) in multiple brain regions of a young and an old BXD89 mouse, demonstrating the application of the workflow to a study of the impact of aging on neuronal populations in a different strain. The error bars in this figure represent the standard deviation of neuron distribution within each region, obtained from subvolumes.

4. Discussion

We present a workflow that corrects the deformation introduced by the chemical treatment of brain tissue for LSM and uses an automated algorithm for neuron segmentation and counting through machine learning and post-processing. To improve the accuracy of the results, we use the MRH of the same specimen as a reference for LSM image morphology, as MRH provides in-skull images of the brain that closely approximate the compartmental configuration and accurate size of the living brain. Our method offers a reliable means of acquiring neuron number and density and requires less human labor and subjective judgment than previous methods. The method exhibits high flexibility, as it requires a much lower amount of pre-training data and demonstrates efficient computational time and performance compared to deep learning methods, as depicted in Supplementary Figure S6. It is also much easier to replicate across multiple brain regions and specimens. By correcting raw images and conducting region segmentation, our method provides improved estimates over previous counting methods. Furthermore, by incorporating the principles of optical fractionation, our workflow provides an unbiased sampling method with equal probability to sample for each section with a relatively low sampling ratio (Supplementary Table S4), resulting in reduced computational costs and representing a novel digital approach that has not been explored before.

To validate the method, we used the workflow to analyze 13 regions in different divisions of 5 C57BL/6 J specimens. The results, as shown in Figure 5, demonstrate the variability in neuron density across regions within the same animal and across animals for the same regions. In the Supplementary materials, we conducted a systematic search of the literature providing neuron counts of the same regions for comparison with our method. As our method produces regional neuron density as its direct output, we assessed the accuracy of our results in terms of density. For most cortical regions, our findings were consistent with those of previous studies or fell within the range of previously reported results. However, in some regions, e.g., LGd, PSV, VISp, TH, and BLA, our results showed a higher neuron density than what was reported in earlier studies. This is likely due to our utilization of the morphology-corrected space for calculating neuron density. This space closely approximates the volume and morphology of the brain confined within the skull, which may have contributed to the observed differences in density measurements. In certain regions where our density measurements were markedly different from those reported in previous studies, such as the facial motor nucleus, we manually counted the sub-regions to validate the neuron distribution. Our findings indicated that the neuron distribution in this region was indeed sparse, with a density measurement about six times lower than that reported in the Blue Brain Map (BBM) (Ero et al., 2018). We hypothesize that this discrepancy may be due to differences in data quality or acquisition methods between studies. BBM estimated neuron counts using a transfer function and Monte Carlo method applied to the ABA Nissl atlas. The authors only used whole-brain values from literature, such as the total number of cells and neurons in the mouse brain, to constrain their estimates of cell densities in each brain region. Further investigation is required to elucidate the causes of the observed discrepancy. Nonetheless, our method, which directly measures regional neuron density, provides a reliable alternative to indirect methods such as those used by BBM, and can contribute to a more accurate appreciation of brain structure.

Our method has two limitations. In those regions in which the cell density is so high that the resolution of the LSM we acquired is not sufficient to differentiate countable cells (e.g., dentate gyrus, granule cell layer of cerebellar cortex), the random forest classifier and watershed algorithms will fail, and hence the performance hinges on the volume filter, which approximates the neuron count by partitioning the volumes of nested neurons based on the average neuron size (Supplementary Section S7). A second limitation is that performance is dependent on the quality of the raw data. In cases where the tissue preparation and scanning are suboptimal, the method will not function properly. Although the data provided by LifeCanvas, which we used for our study, are generally of high quality, there are occasionally regions with poor image quality (as shown in Supplementary Figure S7), which can limit the method’s reliability. Advanced imaging techniques such as expansion microscopy, as demonstrated by Wassie et al. (2019) or acquisition with higher power objectives will address resolution questions and are likely to produce reliable and accurate counts of neurons using our workflow.

The application of stereology to counting cells has a venerable history. As acquisition methods have improved the methods to count cells have been adapted to the newer data. The advent of tissue clearing and volume light sheet imaging yield opportunities for significant improvement in precision and accuracy. This is particularly important in our studies of the mouse brain as we seek apply this workflow to quantitative studies of the mophological and cytological impacts of neurodegeneration and aging. Correction for tissue changes, use of a more accurate label set, and automation described in this work will provide crucial tools in quantitative neuropathology.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement

The animal study was approved by Duke Institutional Animal Care and Use Committee. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

YT: Designed and coded the workflow, conducted the evaluation, and wrote the manuscript. GJ: Conceived the project, guided YT, and contributed to manuscript writing. RW: Offered valuable expertise and guidance in stereology. LW: Provided valuable expertise and guidance into neuron density, neuroscience, and strategy. GJ and RW: Provided funding and contributed to manuscript writing. All authors participated in manuscript revisions, reviewed, and approved the final submission.

Funding

This work was supported by National Institute of Aging R01AG070913 to GJ and RW.

Acknowledgments

We are grateful to James J. Cook for technical support and assistance in software development. We are grateful to Tatiana Johnson for assistance in editing the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2023.1223226/full#supplementary-material

Footnotes

1. ^https://www.ilastik.org

2. ^https://imaris.oxinst.com

3. ^https://fiji.sc/

4. ^https://github.com/YuqiTianCIVM/NeuronCounting/blob/main/Algor_forVisualization.py

5. ^https://imaris.oxinst.com

6. ^https://support.apple.com/guide/automator/intro-to-automator-aut6e8156d85/mac

7. ^https://github.com/YuqiTianCIVM/NeuronCounting/blob/main/MainAlgor.py

8. ^https://imagej.net/plugins/labkit/

9. ^http://imagej.net/MorphoLibJ

References

Abercrombie, M. (1946). Estimation of nuclear population from microtome sections. Anat. Rec. 80, 191–202. doi: 10.1002/ar.1090940210

CrossRef Full Text | Google Scholar

Amunts, K., Schlaug, G., Schleicher, A., Steinmetz, H., Dabringhaus, A., Roland, P. E., et al. (1996). Asymmetry in the human motor cortex and handedness. NeuroImage 4, 216–222. doi: 10.1006/nimg.1996.0073

CrossRef Full Text | Google Scholar

Anderson, R. J., Cook, J. J., Delpratt, N., Nouls, J. C., Gu, B., McNamara, J. O., et al. (2019). Small animal multivariate brain analysis (SAMBA) – a high throughput pipeline with a validation framework. Neuroinformatics 17, 451–472. doi: 10.1007/s12021-018-9410-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Arzt, M., Jr, D., Schmied, C., Pietzsch, T., Schmidt, D., Tomancak, P., et al. (2022). LABKIT: labeling and segmentation toolkit for big image data. Front. Comput. Sci 4:4. doi: 10.3389/fcomp.2022.777728

CrossRef Full Text | Google Scholar

Ashbrook, D. G., Arends, D., Prins, P., Mulligan, M. K., Roy, S., Williams, E. G., et al. (2021). A platform for experimental precision medicine: the extended BXD mouse family. Cell Syst 12, 235–247.e9. doi: 10.1016/j.cels.2020.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Avants, B. B., Tustison, N. J., Wu, J., Cook, P. A., and Gee, J. C. (2011). An open source multivariate framework for n-tissue segmentation with evaluation on public data. Neuroinformatics 9, 381–400. doi: 10.1007/s12021-011-9109-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Bandeira, F., Lent, R., and Herculano-Houzel, S. (2009). Changing numbers of neuronal and non-neuronal cells underlie postnatal brain growth in the rat. Proc. Natl. Acad. Sci. U. S. A. 106, 14108–14113. doi: 10.1073/pnas.0804650106

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonthius, D. J., McKim, R., Koele, L., Harb, H., Karacay, B., Mahoney, J., et al. (2004). Use of frozen sections to determine neuronal number in the murine hippocampus and neocortex using the optical disector and optical fractionator. Brain Res. Brain Res. Protoc. 14, 45–57. doi: 10.1016/j.brainresprot.2004.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Collins, C. E., Young, N. A., Flaherty, D. K., Airey, D. C., and Kaas, J. H. (2010). A rapid and reliable method of counting neurons and other cells in brain tissue: a comparison of flow cytometry and manual counting methods. Front. Neuroanat. 4:5. doi: 10.3389/neuro.05.005.2010

CrossRef Full Text | Google Scholar

Deniz, O. G., Altun, G., Kaplan, A. A., Yurt, K. K., von Bartheld, C. S., and Kaplan, S. (2018). A concise review of optical, physical and isotropic fractionator techniques in neuroscience studies, including recent developments. J. Neurosci. Methods 310, 45–53. doi: 10.1016/j.jneumeth.2018.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Ero, C., Gewaltig, M. O., Keller, D., and Markram, H. (2018). A cell atlas for the mouse brain. Front. Neuroinform. 12:84. doi: 10.3389/fninf.2018.00084

PubMed Abstract | CrossRef Full Text | Google Scholar

Herculano-Houzel, S., and Lent, R. (2005). Isotropic fractionator: a simple, rapid method for the quantification of total cell and neuron numbers in the brain. J. Neurosci. 25, 2518–2521. doi: 10.1523/JNEUROSCI.4526-04.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hillman, E. M. C., Voleti, V., Li, W., and Yu, H. (2019). Light-sheet microscopy in neuroscience. Annu. Rev. Neurosci. 42, 295–313. doi: 10.1146/annurev-neuro-070918-050357

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, G. A., Badea, A., Brandenburg, J., Cofer, G., Fubara, B., Liu, S., et al. (2010). Waxholm space: an image-based reference for coordinating mouse brain research. NeuroImage 53, 365–372. doi: 10.1016/j.neuroimage.2010.06.067

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, G. A., Benveniste, H., Black, R. D., Hedlund, L. W., Maronpot, R. R., and Smith, B. R. (1993). Histology by magnetic resonance microscopy. Magn. Reson. Q. 9, 1–30.

PubMed Abstract | Google Scholar

Johnson, G. A., Tian, Y., Ashbrook, D. G., Cofer, G. P., Cook, J. J., Gee, J. C., et al. (2023). Merged magnetic resonance and light sheet microscopy of the whole mouse brain. Proc. Natl. Acad. Sci. U. S. A. 120:e2218617120. doi: 10.1073/pnas.2218617120

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, G. A., Wang, N., Anderson, R. J., Chen, M., Cofer, G. P., Gee, J. C., et al. (2019). Whole mouse brain connectomics. J. Comp. Neurol. 527, 2146–2157. doi: 10.1002/cne.24560

PubMed Abstract | CrossRef Full Text | Google Scholar

Kasthuri, N., Hayworth, K. J., Berger, D. R., Schalek, R. L., Conchello, J. A., Knowles-Barley, S., et al. (2015). Saturated reconstruction of a volume of neocortex. Cells 162, 648–661. doi: 10.1016/j.cell.2015.06.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. Y., Cho, J. H., Murray, E., Bakh, N., Choi, H., Ohn, K., et al. (2015). Stochastic electrotransport selectively enhances the transport of highly electromobile molecules. Proc. Natl. Acad. Sci. U. S. A. 112, E6274–E6283. doi: 10.1073/pnas.1510133112

PubMed Abstract | CrossRef Full Text | Google Scholar

Krupa, O., Fragola, G., Hadden-Ford, E., Mory, J. T., Liu, T., Humphrey, Z., et al. (2021). NuMorph: tools for cortical cellular phenotyping in tissue-cleared whole-brain images. Cell Rep. 37:109802. doi: 10.1016/j.celrep.2021.109802

PubMed Abstract | CrossRef Full Text | Google Scholar

Lustig, M., Donoho, D., and Pauly, J. M. (2007). Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn. Reson. Med. 58, 1182–1195. doi: 10.1002/mrm.21391

PubMed Abstract | CrossRef Full Text | Google Scholar

Menegas, W., Bergan, J. F., Ogawa, S. K., Isogai, Y., Umadevi Venkataraju, K., Osten, P., et al. (2015). Dopamine neurons projecting to the posterior striatum form an anatomically distinct subclass. elife 4:e10032. doi: 10.7554/eLife.10032

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, E., Cho, J. H., Goodwin, D., Ku, T., Swaney, J., Kim, S. Y., et al. (2015). Simple, scalable proteomic imaging for high-dimensional profiling of intact systems. Cells 163, 1500–1514. doi: 10.1016/j.cell.2015.11.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, Y. G., Sohn, C. H., Chen, R., McCue, M., Yun, D. H., Drummond, G. T., et al. (2019). Protection of tissue physicochemical properties using polyfunctional crosslinkers. Nat. Biotechnol. 37:73-+. doi: 10.1038/nbt.4281

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavlov, Y. L. (2000). Random forests. Utrecht, Boston, VSP

Google Scholar

Peters, A., Josephson, K., and Vincent, S. L. (1991). Effects of aging on the neuroglial cells and pericytes within area 17 of the rhesus-monkey cerebral-cortex. Anat. Rec. 229, 384–398. doi: 10.1002/ar.1092290311

PubMed Abstract | CrossRef Full Text | Google Scholar

Pietzsch, T., Saalfeld, S., Preibisch, S., and Tomancak, P. (2015). BigDataViewer: visualization and processing for large image data sets. Nat. Methods 12, 481–483. doi: 10.1038/nmeth.3392

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, J. L., Ko, A. I., Wade, M. J., Tsou, S. K., McKeel, D. W., and Morris, J. C. (2001). Neuron number in the entorhinal cortex and CA1 in preclinical Alzheimer disease. Arch. Neurol. 58, 1395–1402. doi: 10.1001/archneur.58.9.1395

CrossRef Full Text | Google Scholar

Renier, N., Adams, E. L., Kirst, C., Wu, Z., Azevedo, R., Kohl, J., et al. (2016). Mapping of brain activity by automated volume analysis of immediate early genes. Cells 165, 1789–1802. doi: 10.1016/j.cell.2016.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosen, G. D., and Williams, R. W. (2001). Complex trait analysis of the mouse striatum: independent QTLs modulate volume and neuron number. BMC Neurosci. 2:5. doi: 10.1186/1471-2202-2-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmitz, C., and Hof, P. R. (2005). Design-based stereology in neuroscience. Neuroscience 130, 813–831. doi: 10.1016/j.neuroscience.2004.08.050

CrossRef Full Text | Google Scholar

Selemon, L. D., and Goldman-Rakic, P. S. (1999). The reduced neuropil hypothesis: a circuit based model of schizophrenia. Biol. Psychiatry 45, 17–25. doi: 10.1016/S0006-3223(98)00281-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Spocter, M. A., Hopkins, W. D., Barks, S. K., Bianchi, S., Hehmeyer, A. E., Anderson, S. M., et al. (2012). Neuropil distribution in the cerebral cortex differs between humans and chimpanzees. J. Comp. Neurol. 520, 2917–2929. doi: 10.1002/cne.23074

PubMed Abstract | CrossRef Full Text | Google Scholar

Susaki, E. A., Tainaka, K., Perrin, D., Kishino, F., Tawara, T., Watanabe, T. M., et al. (2014). Whole-brain imaging with single-cell resolution using chemical cocktails and computational analysis. Cells 157, 726–739. doi: 10.1016/j.cell.2014.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, Y., Cook, J. J., and Johnson, G. A. (2023). Restoring morphology of light sheet microscopy data based on magnetic resonance histology. Front. Neurosci. 16:1011895. doi: 10.3389/fnins.2022.1011895

CrossRef Full Text | Google Scholar

Ueda, H. R., Ertürk, A., Chung, K., Gradinaru, V., Chédotal, A., Tomancak, P., et al. (2020). Tissue clearing and its applications in neuroscience. Nat. Rev. Neurosci. 21, 61–79. doi: 10.1038/s41583-019-0250-1

PubMed Abstract | CrossRef Full Text | Google Scholar

von Bartheld, C. (2002). Counting particles in tissue sections: choices of methods and importance of calibration to minimize biases. Histol. Histopathol. 17, 639–648. doi: 10.14670/HH-17.639

CrossRef Full Text | Google Scholar

Wang, N., Anderson, R. J., Badea, A., Cofer, G., Dibb, R., Qi, Y., et al. (2018). Whole mouse brain structural connectomics using magnetic resonance histology. Brain Struct. Funct. 223, 4323–4335. doi: 10.1007/s00429-018-1750-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Ding, S. L., Li, Y., Royall, J., Feng, D., Lesnar, P., et al. (2020). The Allen mouse brain common coordinate framework: a 3D reference atlas. Cells 181, 936–953.e20. doi: 10.1016/j.cell.2020.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Wassie, A. T., Zhao, Y., and Boyden, E. S. (2019). Expansion microscopy: principles and uses in biological research. Nat. Methods 16, 33–41. doi: 10.1038/s41592-018-0219-4

CrossRef Full Text | Google Scholar

West, M. J. (1999). Stereological methods for estimating the total number of neurons and synapses: issues of precision and bias. Trends Neurosci. 22, 51–61. doi: 10.1016/S0166-2236(98)01362-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, R. W., and Rakic, P. (1988). Three-dimensional counting: an accurate and direct method to estimate numbers of cells in sectioned material. J. Comp. Neurol. 278, 344–352. doi: 10.1002/cne.902780305

CrossRef Full Text | Google Scholar

Wree, A., Schleicher, A., and Zilles, K. (1982). Estimation of volume fractions in nervous tissue with an image analyzer. J. Neurosci. Methods 6, 29–43. doi: 10.1016/0165-0270(82)90014-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, N. A., Flaherty, D. K., Airey, D. C., Varlan, P., Aworunse, F., Kaas, J. H., et al. (2012). Use of flow cytometry for high-throughput cell population estimates in brain tissue. Front. Neuroanat. 6:27. doi: 10.3389/fnana.2012.00027

CrossRef Full Text | Google Scholar

Zhang, C., Yan, C., Ren, M., Li, A., Quan, T., Gong, H., et al. (2017). A platform for stereological quantitative analysis of the brain-wide distribution of type-specific neurons. Sci. Rep. 7:14334. doi: 10.1038/s41598-017-14699-w

PubMed Abstract | CrossRef Full Text | Google Scholar