Automated 3D Axonal Morphometry of White Matter

Axonal structure underlies white matter functionality and plays a major role in brain connectivity. The current literature on the axonal structure is based on the analysis of two-dimensional (2D) cross-sections, which, as we demonstrate, is precarious. To be able to quantify three-dimensional (3D) axonal morphology, we developed a novel pipeline, called ACSON (AutomatiC 3D Segmentation and morphometry Of axoNs), for automated 3D segmentation and morphometric analysis of the white matter ultrastructure. The automated pipeline eliminates the need for time-consuming manual segmentation of 3D datasets. ACSON segments myelin, myelinated and unmyelinated axons, mitochondria, cells and vacuoles, and analyzes the morphology of myelinated axons. We applied the pipeline to serial block-face scanning electron microscopy images of the corpus callosum of sham-operated (n = 2) and brain injured (n = 3) rats 5 months after the injury. The 3D morphometry showed that cross-sections of myelinated axons were elliptic rather than circular, and their diameter varied substantially along their longitudinal axis. It also showed a significant reduction in the myelinated axon diameter of the ipsilateral corpus callosum of rats 5 months after brain injury, indicating ongoing axonal alterations even at this chronic time-point.

Electron microscopy (EM) techniques are used extensively to assess brain tissue ultrastructure. Studies have reported the morphology, distribution, and interactions of different cellular components in both healthy and pathological brain using transmission electron microscopy (TEM) [1][2][3][4] . The ultra-thin sections prepared for TEM can only provide 2-dimensional (2D) information, limiting the full characterization of 3-dimensional (3D) cellular structures. Recent advanced EM techniques allow for new possibilities to study the ultrastructure of the brain in 3D [5][6][7][8][9] . One of these techniques is serial block-face scanning electron microscopy (SBEM) 6 . SBEM combines scanning electron microscopy (SEM) with back-scattered electron detection and low beam energies 10 . Images are acquired from the block-face of a sample each time an ultra-microtome inside the vacuum chamber removes the top section from a block-face to expose a new surface for imaging. The result is a stack of high-resolution, high-contrast images of tissue. Compared to other 3D EM techniques, such as focused ion beam (FIB), serial section TEM, or 3D-tomography, SBEM enables imaging of up to several hundreds of micrometers of tissue at nanoscopic resolution without manual tissue sectioning 5,11 . Thus, SBEM is the method of choice for mesoscale imaging of brain tissue ultrastructure.
Despite substantial progress in 3D image acquisition techniques, segmentation and quantification of SBEM data remain challenging. To date, several software tools have been developed that focus on either manual annotation (e.g., KNOSSOS 12 , TrakEM2 13 , Microscopy Image Browser 14 , and CATMAID 15 ), or interactive processing of data by combining automated analysis and proof-reading capabilities (e.g., rhoANA 16 , ilastik 17 , and SegEM 18 ). In addition to these software tools, a variety of studies have also proposed segmentation pipelines for analyzing large amounts of TEM data. Recent studies [19][20][21][22][23][24][25][26] initially identified cellular boundaries using pixel-wise classification methods, followed by over-segmentation of the intracellular regions in each 2D image. This procedure requires merging the results within and between consecutive images using different strategies (e.g., watershed merge tree 23 , agglomerative or hierarchical clustering [19][20][21]25,26 , and joint segmentation of several images in anisotropic datasets 22,24 ).
Although the EM segmentation approaches cited above have yielded impressive results, they have focused on the neuronal reconstruction of grey matter. In this study, we address quantification of white matter ultrastructure and particularly the morphometry of myelinated axons in sham-operated and animals after traumatic brain injury (TBI). Characterization of the white matter ultrastructure requires the segmentation of the white Results ACSON segmentation pipeline automatically annotates the white matter ultrastructure. We devised the ACSON segmentation pipeline to annotate the ultrastructure in SBEM volumes of white matter. The segmentation procedure of ACSON labeled white matter voxels as myelin, myelinated and unmyelinated axons, cells, mitochondria, and vacuoles. In addition, ACSON provided a separate label for each individual axon. The ACSON segmentation pipeline, illustrated in Fig. 1, comprised the following steps: (1) denoising the SBEM volume; (2) segmenting the volume using the bounded volume growing (BVG) technique, which integrates seeded region growing 31 and edge detection algorithms in 3D; (3) refining the segmentation with supervoxel techniques; (4) segmenting the subcellular structure and cells, and annotating myelinated and unmyelinated axons.
Evaluation of the ACSON segmentation pipeline. We quantified the accuracy of the ACSON segmentation pipeline against manual segmentation by an expert. The expert (A.S.) manually segmented three 2D images (images 50, 55, and 60) from the contralateral dataset of one of the sham-operated rats (Sham-1). The images were selected to be 0.2 μm apart. The expert had no access to the automated segmentations of the dataset. Figure 2b shows the manual segmentations and the corresponding images produced by the automated segmentation. The segmentation accuracy was quantified using the precision (positive predictive value) and recall (sensitivity) in the tissue-type level similar to the previous studies 28,32 , and weighted Jaccard index and weighted Dice coefficients in the region level (see Materials and Methods). The precision and recall obtain their maximum value, one, if the automated segmentation correctly assigned voxels to myelin, myelinated or unmyelinated axon labels. The evaluation metrics in the tissue-type level, however, fail to account for topological differences between the ground truth and the automated segmentation. We used weighted Jaccard index and weighted Dice coefficients to evaluate the segmentation accuracy in the region level. Thus, each axon was considered to be its own region, which is a much more stringent evaluation criterion than considering all axons as a single region. The maximum value for these metrics is one, which occurs when a segmented region by ACSON perfectly matches a region segmented manually. If no overlap occurs, the metric is equal to zero. Table 1 reports the precision, recall, weighted Jaccard index, and weighted Dice coefficient values of the three slices shown in Fig. 2 and Supplementary Fig. S1. The results showed an excellent agreement between the automated and manual segmentations for myelin (precision ≥ 0.86, recall ≥ 0.88, weighted Jaccard index ≥ 0.78, and weighted Dice coefficients ≥ 0.87) and myelinated axons (precision ≥ 0.84, recall ≥ 0.88, weighted Jaccard index ≥ 0.80, and weighted Dice coefficients ≥ 0.88). For unmyelinated axons, in the tissue-type level the agreement was good (precision ≥ 0.72 and recall ≥ 0.68). The weighted Jaccard index and weighted Dice coefficients of unmyelinated axons showed approximately 0.37 and 0.50 agreement, respectively, which indicated topological differences and differences in perceiving ultrastructures such as myelin pockets and delamination between the automated and manual segmentations. These metrics are sensitive to minor displacements in the location of the boundaries as demonstrated in Supplementary Fig. S2.
We further evaluated the segmentation of myelinated axons for the split (over-segmentation) and merge (under-segmentation) errors. Supplementary Fig. S3 illustrates the distribution of the length of myelinated axons in Sham-1 dataset. Myelinated axons with lengths smaller than 5 μm appeared in the corners of the SBEM volumes. These myelinated axons have partially traversed the SBEM volume. We did not observe split and merge errors in myelinated axons visually, confirming the excellent Jaccard index and Dice coefficient values for myelinated axons. The 3D rendering of myelinated axons for all datasets is shown in Supplementary Fig. S4. Also,  Manual expert segmentation and automated segmentation of three images from the contralateral corpus callosum of Sham-1 dataset. We used the Hungarian algorithm to match the color of segmented regions between automated and manual segmentation panels. The red arrowheads indicate delamination in the myelin sheaths. These substructures were annotated as myelin in the manual annotation, while the automated segmentation excluded the myelin delamination from the myelin-labeled structures. Occasionally, the automated segmentation erroneously merged neighboring unmyelinated axons when the space between the membranes was poorly resolved (white arrowheads). The automated segmentation of myelinated axons was in excellent agreement with the manual segmentation, as shown in Table 1 www.nature.com/scientificreports www.nature.com/scientificreports/ diameter is the diameter of a circle with the same area as the cross-section, and the eccentricity is a measurement of how much a conic section deviates from being circular. For example, the eccentricity of a circle is zero and the eccentricity of an ellipse is greater than zero but less than one. Figure 3b shows the measurements of minor axes along a myelinated axon. It is apparent that the myelinated axon diameter was not constant. We present more examples of cross-sectional variation along myelinated axons in Supplementary Fig. S5 and show that the histogram of cross-sectional diameters within a myelinated axon was bimodal, which can partially be related to the location of mitochondria 34 . Note that mitochondria was included in the intra-axonal space of myelinated axons, while vacuoles excluded, because they appeared between myelinated axon membrane and myelin. Substantial differences between 2D and 3D morphological analyses. We compared the traditional 2D morphometry and the proposed 3D morphometry pipelines for myelinated axons in one dataset of the sham-operated rats (Sham-1). In Fig. 4a, a representative myelinated axon is elongated in the z direction. We sampled the myelinated axon parallel to the x-y plane at three points denoted as p1, p2, and p3 in the figure. Figure 4a shows that when the axonal axis was nearly perpendicular to the sampling plane (point p1), the relative difference between the 2D and 3D quantifications was small. However, when the axonal axis was not aligned with three main orientations, the relative difference between the 2D and 3D quantifications increased and was substantial (points p2 and p3). In addition, as the relative difference varied along a myelinated axon, a single 2D measurement was noisy. We compared the 2D and 3D measurements for all myelinated axons in contralateral corpus callosum of Sham-1 dataset as shown in Fig. 4b-e. The comparisons showed that the median of the relative difference, in percentage, was 16.64% for the minor axis, 18.25% for the major axis, 16.23% for the equivalent diameter, and 11.34% for the eccentricity. This indicates substantial differences between the 2D and 3D based measurements of morphometry of myelinated axons. As the cross-sectional measurements of a myelinated axon varied along it, sampling a single cross-section for measurement, as in 2D morphometry, gave an incomplete account of the parameters to be quantified. The 3D analysis, instead, quantified the morphology of a myelinated axon while taking into account the morphological variations along its axis.
3D morphometry of the ultrastructure of the corpus callosum. We quantified the morphological and volumetric aspects of the white matter ultrastructure in our SBEM datasets. For the morphological analysis, we thresholded myelinated axons based on their length, and preserved those myelinated axons whose length was long enough to run approximately one third of the SBEM volumes, which was equal to 5 μm. Supplementary  Fig. S3a shows the thresholded myelinated axons which partially traversed the SBEM volume. In addition, the thresholding can eliminate myelinated axons which were split erroneously or subcellular structures that were mistakenly labeled as myelinated axons.
For each myelinated axon, we considered the median of the cross-sectional measurements (minor and major axes, equivalent diameter and eccentricity) as our measurement variable. The box plots of these measurements are shown in Fig. 5. We subjected these quantities to nested (hierarchical) 1-way analysis of variance (ANOVA) separately for each hemisphere 35 . The nested ANOVA tests whether there exists a significant variation in means among groups, or among subgroups within groups. The myelinated axons were nested under rats and the identities of rats were nested under the groups (sham-operated and TBI). We set the alpha-threshold defining the statistical significance as 0.05 for all analyses. The ANOVA was performed using the anovan function of MATLAB R2017b. Myelinated axon and rat identity were treated as random effects and the group was treated as fixed effect.
We observed a significant reduction of the equivalent diameter (F = 14.4, p = 0.029) and the major axis (F = 26.4, p = 0.012) in the ipsilateral corpus callosum of TBI rats as compared to the ipsilateral corpus callosum www.nature.com/scientificreports www.nature.com/scientificreports/ of sham-operated rats (see Supplementary Table S1 and Fig. 5). The contralateral corpus callosum, however, did not show a significant difference between sham-operated and TBI rats for any of the measures (see Supplementary  Table S1 and Fig. 5). We measured that the equivalent diameter was about 15% greater than the minor axis of the fitted ellipses. The eccentricity in all datasets were markedly different from zero indicating that the axonal cross-sections were not circular, but elliptical.
We also partitioned the variability of the measurements in cross-sectional, axonal, and animal levels, known as variance component analysis. The variance components describe what percentage of the total variance is attributable to each level. The design matrix of a 4-levels nested ANOVA for all measurements was very large (2 groups, i.e., sham-operated and TBI, 5 animals, approximately 250 myelinated axons per animal, and approximately 250 cross-sectional measurements per myelinated axon) making it impossible to perform the computations with the complete dataset. Therefore, we sampled 10000 of the measurements (without repetition), computed the variance components and repeated the procedure for 10 times. The variance components (Table 2) indicated that the greatest part of the variance was attributed to variance of cross-sections, then between myelinated axons, and the least amount of variance was attributed to between sham-operated and TBI animals.
We also quantified the volumetric aspects of myelin in our 3D-SBEM datasets. The ultrastructure volumetry was dataset-dependent, preventing a direct comparison between datasets. For example, the volume of cell body/ processes varied among datasets [1.3-22.1%], influencing the volume of the other ultrastructures (see the results of volumetric analysis in Supplementary Table S2). Therefore, we calculated the implicit representation of the g-ratio 36 with no measurement of the myelin thickness 37 , denoted as − = − ⁎ aggregate g ratio M yelin 1 . We defined Myelin* as the ratio of the myelin volume to the myelin volume plus the intra-axonal volume of all The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using the 'o' symbol. Nested ANOVA showed a significant reduction in the diameter of myelinated axons (the equivalent diameter and the length of major axis) in the ipsilateral corpus callosum of rats after TBI.  Table 2. Variance components analysis. The variance was partitioned into cross-sectional, axonal, and animal levels.

Contralateral
www.nature.com/scientificreports www.nature.com/scientificreports/ myelinated axons. In addition, we calculated the density for myelinated axons and mitochondria defined as where N is the number of myelinated axons or mitochondria counted once in a volume V nc which is the volume of non-cellular ultrastructure. Supplementary Table S2 shows the aggregate g-ratio and the density of myelinated axons and mitochondria for all datasets. We did not run statistical hypothesis tests for comparing these metrics between sham-operated and TBI rats. However, visually, we observed myelin delamination and pockets in myelin sheath more frequent in rats after TBI compared to sham-operated rats.
Computation time. On a 4-core Intel CPU 3.41 GHz machine with 64 GiB RAM using MATLAB R2017b, the computation times for one dataset were approximately as follows: block-matching and 4D (BM4D) filtering 5-7 h, segmentation process 1 day, and skeletonization and cross-sectional analysis 6 h. Correcting the segmentation of mitochondria and proofreading of myelinated axons versus unmyelinated axons was accomplished in 7 h for the entire SBEM-volume. The manual expert annotation of a single slice of the SBEM datasets required 10 h.

Discussion
Previous studies that quantified white matter were limited to 2D morphometry, which simplifies the assumptions about axonal morphology. In this paper, we reported an extensive 3D morphological analysis of SBEM volumes. For this, we devised a novel pipeline, termed ACSON, for automated segmentation and morphometry of the cellular and subcellular components of the corpus callosum in SBEM datasets. ACSON segmented white matter into myelin, the intra-axonal space of myelinated and unmyelinated axons, cell bodies and their processes, and subcellular compartments, such as mitochondria and vacuoles. The segmentation accuracy evaluations revealed excellent agreement between the automated and manual segmentation of myelin and myelinated axons. ACSON quantified the morphology of the segmented myelinated axons. The 3D morphometry demonstrated a substantial variation in the diameter of myelinated axons along their longitudinal axis. The results indicated that the cross-sections of a myelinated axon are more likely to be elliptical than circular. To compare sham-operated and TBI animals, we used nested ANOVA and variance component analysis. After TBI, we found a significant reduction in the diameter of myelinated axons in the ipsilateral of corpus callosum, indicating that the alterations persisted for several months after the injury.
Traditionally, studies have modeled axons as straight cylinders, and utilized 2D-EM sections to assess the axonal morphology [38][39][40] . However, studies of unmyelinated axons in peripheral nerves 41 and in the hippocampus and cerebellum 42 indicated highly irregular axial shapes with periodic varicosities containing organelles 42 . In addition to studying the anatomical alterations of brain ultrastructure in disease, models of diffusion magnetic resonance imaging (dMRI) 43,44 or electrophysiology 45,46 in many studies assume simplified geometries for axons. However, Novikov et al. 47 and Fieremans et al. 48 have shown the effect of structural disorder along axons and micron-level sample architecture while modeling dMRI. Thus, realistic modeling of cells 49 and myelinated axons 50 have raised attention for numerical simulation of dMRI. Segmentation of high-resolution EM volumes provides realistic properties such as axonal diameter, eccentricity of cross-sections and orientation dispersion that can substitute simplified models of axons.
Our segmentation and morphometry pipeline is automated. The ACSON segmentation pipeline requires tuning several parameters, such as the thresholds for measuring similarity or vacuole intensity. These parameters, however, are easy to set. Annotating the mitochondria was the only step that remained a challenge and required human intervention. Compared with manual segmentation, which required, on average, 10 hours for a single slice, annotating the mitochondria required approximately six hours for the entire dataset comprising 285 SBEM images. In an earlier study, Lucchi et al. 32 specifically targeted segmentation of mitochondria in high resolution FIB-SEM datasets of 5-6 nm × 5-6 nm planar resolution, which is a much finer resolution than in our datasets. They reported that mitochondria and myelin were difficult to discriminate whenever they were in close proximity. The high resolution of their FIB-SEM datasets enabled them to outline the prominent shape-features of mitochondria and they segmented mitochondria almost in the absence of myelin. Unfortunately, their technique is not applicable to our data due to our larger tissue samples, lower resolution, and proximity of mitochondria to myelin.
The ACSON morphometry pipeline required no user input parameters and extracted a sub-voxel precise and naturally smooth axis for each individual myelinated axon. We assumed a myelinated axon skeleton with no branches and only two endpoints. This allowed us to optimize the computation time for skeleton extraction. The computational efficiency was crucial because we solved the eikonal equation for several thousands of axons with multi-stencils fast marching (MSFM), which is more accurate but also more time-consuming than the fast marching method (FMM) 51 . We quantified myelinated axons for their minor and major axes, equivalent diameter and eccentricity. The minor and major axes and the equivalent diameter, all can describe the cross-sectional diameter of myelinated axons. Estimating the minor and major axes, however, requires fitting an ellipse to the cross-sections. Ellipse fitting step regularizes the measurements, which reduces the cross-sectional variance along myelinated axons. The equivalent diameter is measured directly from the area of cross-sections and it can capture the cross-sectional variance along myelinated axons, however, it measured the cross-sectional diameter 15% greater than the minor axis. Our morphological analyses of myelinated axons in sham-operated animals were in line with a previous study 52 measuring axon diameter in the rat corpus callosum. Kim et al. 52 measured an average myelinated axon diameter of 0.35 μm in the splenium of the corpus callosum in male and female rats at postnatal day 60 using 2D electron microscope. The authors measured the myelinated axon diameter as the minimum diameter of each axon. For a comparison to this study, we measured an average myelinated axon diameter as equivalent diameter (Sham-1: 0.44 μm and Sham-2: 0.5 μm) and minor axis (Sham-1: 0.37 μm and Sham-2: 0.42 μm) for the contralateral corpus callosum. Note that our tissue samples were acquired from the body of the corpus callosum at approximately −3.80 mm from bregma from male adult rats. Moreover, the variability of the myelinated axon diameter along its longitudinal axis can be related to the accumulation of organelles, such as www.nature.com/scientificreports www.nature.com/scientificreports/ mitochondria, increasing the cross-sectional diameter locally 34 . Wake et al. 53 also indicated that accumulation of vesicles containing neurotransmitter along an axon induces a local rise in the form of varicosities.
The segmentation accuracy evaluations demonstrated substantial agreement between the automated and manual segmentations as shown by the precision and recall metrics. These evaluations, however, did not provide a realistic estimation of the segmentation quality when the goal is to separate distinct axons. The weighted Jaccard index and weighted Dice coefficients computed in the region level were much more stringent evaluation measurements, and demonstrated excellent results for the segmentation of myelin and myelinated axons. However, the weighted Jaccard index and weighted Dice coefficients indicated 37% and 50% agreement in the segmentation of unmyelinated axons, respectively. There were several possible reasons for the decreased accuracy in the segmentation of unmyelinated axons. First, pockets in the myelin sheaths were included in the myelin label, while ACSON annotated these volumes individually. This potentially introduced false positives reducing the weighted mean of the Dice coefficients. Second, faintly resolved unmyelinated axons might have been included into the cell body/ process annotation. Finally, the cellular boundaries of unmyelinated axons were often difficult to detect, which resulted in the erroneous merging of neighboring unmyelinated axons. Overall, as the precision and recall scores of unmyelinated axons were good, we investigated the volumetry of unmyelinated axons.
We have provided a framework for the statistical inference of data with nested structure. Simply pooling the measurements from all animals is not a valid approach 35 . Also, averaging the lowest levels of hierarchy has sub-optimal statistical power 54 . The nested ANOVA analysis, however, can capture the variability of the measurements in all levels. In addition to testing the equality of the means at each level through nested ANOVA, we partitioned the variance into different levels. The variance component analysis showed that variation among animals was relatively small compared to the variation among cross-sections and myelinated axons.
It is becoming increasingly clear that white matter pathology plays a major role in brain disorders. After TBI, the white matter pathology is extremely complex. The initial response of an axon to a brain injury 55 can be either degeneration or regeneration. Studies of the white matter ultrastructure indicate axonal swelling in the early stages of TBI 56,57 . Axonal damage persists for years after injury in humans 58 and for up to 1 year in rats 59 . In the present study, we observed morphological changes in the corpus callosum 5 months after severe TBI in rats. We found that the axonal diameter of myelinated axons in the ipsilateral corpus callosum of TBI rats had reduced significantly. The reduction of the diameter of myelinated axons might indicate a prolonged axonal degeneration 60 or regeneration 55 . In analysis of myelin, the automated pipeline has excluded the pockets and delamination from the myelin label. Visually, we observed more frequent delamination in TBI animals, which might indicate active myelin processes still ongoing 5 months after the injury. Myelin delamination was also observed in sham-operated rats, which may be part of the natural dynamics of healthy myelin. Note that, we observed pockets in the myelin sheaths more frequently in animals after injury. The delamination of the myelin sheaths in healthy and injured brains at this chronic time point is currently unclear, and warrants further studies.
When characterizing the ultrastructural morphometry, the tissue shrinkage caused by fixation, staining, and sectioning might affect quantification of the axonal diameter and volumes 61 . In addition, the locations from which the specimens were obtained might influence the quantifications. The SBEM datasets in this study were consistently imaged at a specific location in the corpus callosum in both sham-operated and TBI animals, and in both hemispheres. Due to the small tissue size, however, the environment might change from one sample to another. For example, one of our datasets from the contralateral corpus callosum in a TBI rat contained more cell body/ process volume than the other datasets. The small sample size can be reflected as well in the orientation of the myelinated axons. In a coronal view of the rat brain, the corpus callosum mainly contains in-plane parallel axons oriented in the latero-medial orientation, however, there are also bundles of axons oriented dorso-ventrally, and even rostro-caudally. For that reason, our datasets can contain different axonal populations with different orientations. Thus, studies including more SBEM datasets from more subjects and/or locations in the corpus callosum, as well as bigger sample size are necessary to increase our understanding of the 3D ultrastructure of the corpus callosum.

Materials and Methods
Animal model, tissue preparation, and SBEM imaging. Animals. Five adult male Sprague-Dawley rats (10-weeks old, weight 320 and 380 g, Harlan Netherlands B.V., Horst, Netherlands) were used in the study. The animals were singly housed in a room (22 ± 1°C, 50-60% humidity) with 12 h light/dark cycle and free access to food and water. All animal procedures were approved by the Animal Care and Use Committee of the Provincial Government of Southern Finland and performed according to the guidelines set by the European Community Council Directive 86/609/EEC. Traumatic brain injury model. TBI was induced by lateral fluid percussion injury in three rats (TBI-1, TBI-2, TBI-3) 62 . Rats were anesthetized with a single intraperitoneal injection (6 mL/kg) of a mixture of sodium pentobarbital (58 mg/kg), magnesium sulphate (127.2 mg/kg), propylene glycol (42.8%), and absolute ethanol (11.6%). A craniectomy (5 mm in diameter) was performed between bregma and lambda on the left convexity (anterior edge 2.0 mm posterior to bregma; lateral edge adjacent to the left lateral ridge). Lateral fluid percussion injury was induced in one rat by a transient fluid pulse impact (21 ms) against the exposed intact dura using a fluid-percussion device (Amscien Instruments, Richmond, VA, USA). The impact pressure was adjusted to 3.1 atm to induce a severe injury. Two sham-operated rats (Sham-1, Sham-2) underwent all the same surgical procedures except for the impact.
Tissue processing. Five months after TBI or sham operation, the rats were transcardially perfused using 0.9% NaCl (30 mL/min) for 2 min followed by 4% paraformaldehyde (30 mL/min) at 4 °C for 25 min. The brains were removed from the skull and post-fixed in 4% paraformaldehyde /1% glutaraldehyde overnight at 4 °C. (2019) 9:6084 | https://doi.org/10.1038/s41598-019-42648-2 www.nature.com/scientificreports www.nature.com/scientificreports/ Tissue preparation for SBEM. The brains were sectioned into 1-mm thick coronal sections with a vibrating blade microtome (VT1000s, Leica Instruments, Germany). From each brain, a section at approximately 3.80 mm from bregma was selected and two samples from the ipsilateral and the contralateral corpus callosum were further dissected, as shown in Supplementary Fig. S6a. The samples were stained using an enhanced staining protocol 63 (see Supplementary Fig. S6b). First, the samples were immersed in 2% paraformaldehyde in 0.15 M cacodylate buffer containing 2 mM calcium chloride (pH = 7.4), and then washed five times for 3 min in cold 0.15 Μ cacodylate buffer containing 2 mM calcium chloride (pH = 7.4). After washing, the samples were incubated for 1 h on ice in a solution containing 3% potassium ferrocyanide in 0.3 M cacodylate buffer with 4 mM calcium chloride combined with an equal volume of 4% aqueous osmium tetroxide. They were then washed in double distilled H2O (ddH2O) at room temperature (5 × 3 min). Thereafter, the samples were placed in a solution of 0.01 mg/mL thiocarbohydrazide solution at room temperature for 20 min. The samples were then rinsed again in ddH2O (5 × 3 min), and placed in 2% osmium tetroxide in ddH2O at room temperature. Following the second exposure to osmium, the samples were washed in ddH2O (5 × 3 min), and then incubated in 1% uranyl acetate overnight at 4 °C. The following day, the samples were washed in ddH2O (5 × 3 min) and en bloc Walton's lead aspartate staining was performed. In this step, the samples were incubated in 0.0066 mg/mL lead nitrate in 0.03 M aspartic acid (pH = 5.5) at 60 °C for 30 min, after which the samples were washed in ddH2O at room temperature (5 × 3 min), and dehydrated using ice-cold solutions of freshly prepared 20%, 50%, 70%, 90%, 100%, and 100% (anhydrous) ethanol for 5 min each, and finally placed in ice-cold anhydrous acetone at room temperature for 10 min. Embedding was performed in Durcupan ACM resin (Electron Microscopy Sciences, Hatfield, PA, USA). First, the samples were placed into 25% Durcupan #1 (without component C):acetone, then into 50% Durcupan #1 :acetone, and after into 75% Durcupan #1 :acetone overnight. The following day, they were placed in 100% Durcupan #1 for 2 in a 50 oven (2 times), and into 100% Durcupan #2 (4-component mixture) for 2 h in a 50 °C oven. Finally, the samples were embedded in 100% Durcupan #2 in Beem embedding capsules (Electron Microscopy Sciences) and baked in a 60 °C oven for 48 h.
After selecting the area within the samples, as shown in Supplementary Fig. S6c, the blocks were further trimmed into a pyramidal shape with a 1 × 1 mm 2 base and an approximately 600 × 600 μm 2 top (face), which assured the stability of the block while being cut in the SBEM microscope. The tissue was exposed on all four sides, bottom, and top of the pyramid. The blocks were then mounted on aluminum specimen pins using conductive silver epoxy (CircuitWorks CW2400). Silver paint (Ted Pella, Redding, CA, USA) was used to electrically ground the exposed block edges to the aluminum pins, except for the block face or the edges of the embedded tissue. The entire surface of the specimen was then sputtered with a thin layer of platinum coating to improve conductivity and reduce charging during the sectioning process.
SBEM data acquisition. All SBEM data were acquired on an SEM microscope (Quanta 250 Field Emission Gun; FEI Co., Hillsboro, OR, USA), equipped with the 3View system (Gatan Inc., Pleasanton, CA, USA) using a backscattered electron detector (Gatan Inc.). The top of the mounted block or face was the x-y plane, and the z direction was the direction of the cutting.
All the samples were imaged with a beam voltage of 2.5 kV and a pressure of 0.15 Torr. The datasets were acquired with a resolution of 13-18.3 nm × 13-18.3 nm × 50 nm amounting to an area of 13.3-18.7 μm × 13.3-18.7 μm × 14.25 μm in the x, y, and z directions, respectively. After imaging, Microscopy Image Browser 14 was used to apply lateral registration to the slices. We quantified the registration using cross correlation analysis between successive slices and subtracting the running average (window size = 25) to preserve the directionality of axons while registration. Supplementary Fig. S6d shows a representative SBEM volume of the contralateral corpus callosum of the sham-operated rat. We also show two representative images cropped from the sham-operated and TBI volumes in Supplementary Fig. S6e and f, respectively. ACSON segmentation pipeline. The ACSON segmentation pipeline annotates the ultrastructure in SBEM volumes of white matter. The pipeline began with denoising of the SBEM volumes, and proceeded by segmenting the volumes using BVG. The segmented volumes were refined using supervoxel techniques, and, finally, the subcellular structures, cells, and myelinated and unmyelinated axons were annotated.
Denoising. SBEM images are degraded by noise from different sources, such as noise in the primary beam, secondary emission noise, and noise in the final detection system 64 . To suppress these complex noise patterns, we applied a non-local BM4D algorithm 65 . Unlike local averaging filters, which smooth an image by averaging values in the neighborhood of a target voxel, non-local filtering considers all the voxels in the image, weighted by how similar these voxels are to the target voxel. BM4D, in particular, enhances a sparse representation in the transform-domain by grouping similar 3D image patches (i.e., continuous 3D blocks of voxels) into 4D data arrays called, groups. The steps to realize the filtering are the 4D transformation of 4D groups, shrinkage of the transform spectrum, and inverse 4D transformation. While BM4D has been extensively used for denoising datasets from diverse imaging modalities, its application in 3D-SBEM datasets is novel. We applied the default parameter values of BM4D for denoising, which automatically estimates the noise type and variance. Figure 1a shows a slice of SBEM volume before filtering, and Fig. 1b illustrates the BM4D output in which the noise has been strongly attenuated.
Segmentation. For the segmentation of a 3D-SBEM image, we devised a hybrid technique, named BVG, that integrates seeded region growing and edge detection. To elaborate BVG, we denote a 3D-SBEM image after denoising with BM4D as z(x):X → [0, 1], where x ∈ X is a 3D spatial coordinate. Note that the intensity can range from 0 to 1. We defined ⊂ E X as the set of edges of z, using a Canny edge detector 66 . We set the parameter values of the Canny edge detector as follows: the SD of the Gaussian filter was 2 and thresholds for weak and strong edges were 0.25 and 0.6 times the maximum gradient magnitude. In SBEM volumes with resolution anisotropy and a coarser resolution in the z direction, regions in successive slices did not appear continuously, and areas close to the structure boundaries overlapped. Therefore, we dilated the set of edge coordinates Ê in-plane with a 3 × 3 square structuring element. The dilated edges are denoted as E. The edge dilation was proportional to the resolution anisotropy. We then used BVG to segment z into n + 1 distinct volumes V 1 , V 2 , …, V n , E ⊂ X, in which V i ∩ V j = ∅, V i ∩ E = ∅, ∀i, j = 1, …, n, i ≠ j. BVG is a serial segmentation algorithm, meaning that segmentation of V i starts only when V i−1 is segmented. To segment V i , BVG begins with one voxel called the seed, denoted as S k ⊂ V i , which iteratively grows-k is the iteration number-and finally results in the volume V i . N(S k ) is in the neighborhood of S k defined as N(S k ) = {r|r ∉ S k , ∃s ∈ S k :r ∈ N(s), r ∉ E, r ∉ V 1,…,i−1 }, where N(s) is the 3D-neighborhood of voxel s. In each iteration, BVG appends a set of voxels A to S k , where A = {x|x ∈ N(S k ), δ(x) ≤ δ T } and δ(x) measures the similarity of voxel x to the set S k . We defined the measure of similarity as and set the similarity threshold δ T to 0.1. An iteration terminates, if A = ∅, or |S k | ≥ ϑ, where ϑ is a volume constraint. If V i grew larger than ϑ, the results were discarded and the voxels within V i were freed for other regions to compete for them. The segmentation was initiated by annotating the low-intensity structures, i.e., myelin and mitochondria, which were considered together as V 1 . BVG was initiated with a random low-intensity voxel S 1 with z(S 1 ) ≤ 0.4. This one seed was sufficient to segment V 1 because myelin is a connected structure in a consecutive SBEM image. We defined N(s) using 26 neighbors, and set ϑ = ∞. Figure 1c shows a slice of canny edges together with the segmented myelin and mitochondria (V 1 ). To segment other structures V 2 , …, V n , we needed a more advanced seeding mechanism. Referring to Fig. 1c, we noticed that other structures are surrounded by myelin and edges. Therefore, we first generated a binary mask, B, of the union of the dilated edges E and the myelin-mitochondria segment V 1 (Fig. 1c). Denoting each 2D-slice of B as B i , we computed the Euclidean distance transform 67 for every B i individually, defined as DT i and shown in Fig. 1d. The pixel value in the distance transform DT i is the shortest distance from that pixel to a set of pixels B i . We defined the location of the seeds by extracting the regional maxima of each DT i (Fig. 1d). To segment V i , BVG was initiated with a seed from the set of extracted regional maxima not belonging to any previously segmented V j , j = 1, …, i − 1. We defined N(s) using 6 neighbors and set ϑ = 10 6 , which equals 12.5 μm 3 of tissue or 1.5 times the volume of the largest axons in the dataset. Figure 1e shows one slice of the primary segmentation of the white matter ultrastructure, not belonging to B. The seed extraction overestimates the number of segmented volumes. This does not pose a problem, however, as the serial nature of BVG does not permit repetitive segmentation of an already segmented volume.
Segmentation post processing. The segmentation with BVG may result in small volumes, e.g., smaller than 5 × 10 3 voxels, which actually belong to larger segments. As well, dilated Canny edges E should be assigned with a label as Fig. 1f shows. Therefore, we refined the segmentation by utilizing the SLIC supervoxel 68 technique to relabel the small volumes and attach them to larger ones. Supervoxels group nearby voxels with similar intensity values into clusters. Particularly, SLIC clusters voxels based on a distance measure defined by = + ρ D d d s i nt c sp , where d int ensures intensity similarity and d sp enforces voxel proximity to the supervoxels. In SLIC, the initial supervoxel centers are defined at regular grid steps ρ, and their compactness is controlled by c. Also, the spacing parameter s allows accounting for resolution anisotropy in x, y, and z directions. We assigned the SLIC arguments to produce compact and large supervoxels, while accounting for the resolution anisotropy. Thus, we set c = 23, ρ = 11 and s = [1, 1, vx/vz], where vx and vz are the voxel size in x and z directions, respectively. Note that, voxel size in x and y directions was equal. Supplementary Fig. S7 shows the effect of altering ρ, c and s while generating supervoxels. We refined the large volumes V i with more than 5 × 10 3 voxels by the SLIC supervoxels. In more detail, suppose that we have generated Q supervoxels SV q , q = 1, …, Q. Then, we re-defined V i as Refining the segmentation with SLIC technique eliminated most of the small volumes by relabeling and attaching them to the larger segments. Note that as the edges were included in the supervoxels, the supervoxel-based refinement also labeled voxels belonging to the edges (Fig. 1f).
Annotating subcellular structures and cells. The segmentation of myelin sometimes included mitochondria because the boundaries between these two structures were not clearly resolved. We wanted to label mitochondria separately, however, and include them as part of the myelinated axons. Not including mitochondria in the myelinated axon domain produces cavities, as shown in Fig. 1e,f. The cavities can be used to label the mitochondria. To detect the cavities in the myelinated axons, on each large volume, | |> ′ V 10 i 4 , using a 3D distance transform, we propagated the surface of the volume for 1 μm. The surface of the enlarged volume was then propagated for −1 μm shrinking of the volume. Applying this procedure to each large volume altered the morphology of the volume, and closed those cavities smaller than 1 μm. The difference between the altered volume and V′ was considered a potential mitochondrion, M i . We refined M i with SLIC supervoxels with the same parameter values and techniques mentioned in the segmentation post-processing section. Note that because some of the cavities were due to myelin, annotating the mitochondria was finalized using human supervision to check for myelin. Figure 1g shows the final result of the mitochondria segmentation. The myelin segment was then re-defined as the set difference of V 1 and all mitochondria, denoted as MY. www.nature.com/scientificreports www.nature.com/scientificreports/ axon. The enclosing cylinder was formed by those supervoxels having a common face with the axon AX i . If the enclosing cylinder contained myelin above a threshold, the axon has been considered as a myelinated axon. In more detail, let Λ i be the indexes of supervoxels enclosing an axon AX i . If ≥ .
, we considered AX i to be a myelinated axon. Note that because unmyelinated axons can be surrounded by several myelinated axons, they can be miss-classified as myelinated axons, thus requiring a proof reading after the classification.
To label cells and cell-processes, we considered a straightforward approach as the volumes of cells were expected to be larger than the volumes of any other structure, excluding myelin. Recall that we set the volume threshold ϑ = 10 6 for the segmentation of V 2 , …, V n , which leaves some voxels unlabeled. These unlabeled voxels X′ comprised cells and cell processes. We segmented X′ into n′ cells using connected component analysis. In general, we detected 1-4 cell bodies/process in each SBEM-volume. Figure 1h demonstrates the final segmentation results of myelin, myelinated axons, unmyelinated axons, oligodendrocyte cell body and its processes. Mitochondria and vacuoles belonging to myelinated axons were colored the same as their corresponding myelinated axons. Figure 1i and Supplementary Fig. S4 show the 3D rendering of myelinated axons in contralateral corpus callosum of Sham-1 dataset.
ACSON morphometry pipeline. We defined a cross-section as the intersection of a segmented myelinated axon Ω and a perpendicular plane to the axonal skeleton γ 69 . To detect the myelinated axon skeleton γ with sub-voxel precision, we adapted a method from Van Uitert & Bitter 70 . First, we defined three points in the myelinated axon domain Ω: x* with the largest distance from the myelinated axon surface Γ, and x e1 and x e2 as the endpoints of the axons, i.e., the tips of the axon. The minimum-cost path connecting x e1 to x e2 through x* was the axon skeleton γ. The path was found in two steps, first from x e1 to x*, and then from x e2 to x*. Mathematically, for x ∈ Ω. T 1 (x)| Γ = 0. Using H(x) = 1 − F 2 (x) to define the cost ensured that voxels in the middle of the myelinated axon were reached prior to the voxels close to Γ. We defined the furthest point from x* on the T 2 map, i.e., the global maximum of T 2 , as x e1 . Similarly, x e2 was defined as the furthest point from x e1 , at the global maximum of the time-crossing map T 3 (x), starting from x e1 to every voxel in Ω, with speed F 2 (x). For both endpoints, we determined the minimum-cost path between x ei and x*, γ i , i = 1, 2, by backtracking, starting from x ei and progressing along the negative gradient − until x* was reached. x* is the global minimum of T 2 (x), so that we were guaranteed to find it with backtracking. The backtracking procedure can be described by the ordinary differential equation , where ς traces γ i . We used a 4th order Runge-Kutta scheme, with a 25 nm step size, to solve the ordinary differential equation with sub-voxel accuracy. The myelinated axon skeleton was formed as γ = γ 1 ∪ γ 2 . Note that computing the skeleton in this way prevented the skeleton from cutting corners 72 . Figure 3a shows a 3D reconstruction of a myelinated axon, its mitochondria, and the extracted skeleton (axonal axis). Note that x e1 and x e2 defined as the global maxima of T 2 (x) and T 3 (x), lie on the myelinated axon surface Γ, and not in the center of the myelinated axon. The cost function H, however, forces the skeleton to immediately move away from the surface Γ toward the center. Therefore, we dropped the first 1 μm at both ends of γ in our later calculations.
To determine the cross-sectional planes perpendicular to γ, we formed a moving reference frame of the size 8 μm × 8 μm with 50 nm resolution. At each skeleton point ς, the unit tangent vector to γ was used to define the orientation of the reference frame. The intersection of the reference frame with the myelinated axon defined the cross-section of the myelinated axon.
The intensity values of a cross-section ranged between 0 and 1. Each cross-section was thresholded at 0.5, resulting in a 2D binary image C denoted as C : X → {0, 1}, where the point x = (x 1 , x 2 ) was foreground iff x ∈ C. We defined the center point of C as = ∑ | | ∈ c x C x C 1 . By translating the binary 2D cross-section C to the center of 2D Cartesian coordinate C t = {y ∈ X:y = x − c}, we found an ellipse that had the same normalized second central moment as C t . The cross-sectional morphology of myelinated axons were quantified by computing the minor and major axes and the eccentricity of the fitted ellipse 73 , and the diameter of a circle with the same area as the cross-section, called equivalent diameter.
Evaluation of segmentation accuracy. Manual segmentation. The manual segmentation by A.S. defined each ultrastructure as its own region, i.e., different axons had distinct labels in the manual segmentation . The maximum for the precision and recall is equal to one when the automated segmentation perfectly matches the manual segmentation. These metrics do not describe topological differences between the manual and automated segmentations. For example, these metrics do not penalize the automatic segmentation for incorrectly dividing a single axon into two axons. , where A and B are the regions segmented manually and automatically, respectively. The maximum for these metrics is equal to one occurring when A perfectly matches B. If no overlap occurs between A and B, these metrics are equal to zero. Let A i , i = 1, …, a′, and B j , j = 1, …, b′ be the regions in the manual and automated segmentation, respectively. To assign A i and the best matching B j , we formed a similarity matrix based on Dice coefficients for any possible pair of A i and B j , where the element (i, j) of the similarity matrix was Dice Coef(A i , B j ). We used the Hungarian algorithm 76,77 to match the regions. We defined the weighted mean of the Jaccard index and the weighted mean of the Dice coefficients as ∑ = ′ w J A B ( , )

Comparison of 2D and 3D morphological analyses.
To simulate 2D morphometry, we assumed that the myelinated axon morphometry is quantified on a single image of the 3D image stack. For each myelinated axon, we determined the best orientation of the image stack for the 2D quantifications by extracting the Euler angles of a fitted ellipsoid to the segmented myelinated axon. We randomly selected a single image in that orientation to present the myelinated axon. For example, if a myelinated axon was elongated parallel to the z-axis, we selected a random image parallel to the x − y plane. Each myelinated axon was quantified separately for its minor and major axes, equivalent diameter, and eccentricity. The relative difference between the 2D and 3D quantifications for each myelinated axon was defined as , where q is the quantity of interest, i.e., minor and major axes, equivalent diameter, and the eccentricity, measured by the 2D or 3D procedures. Note that, for the 3D morphometry, median of the measurements along the axonal axis were quantified Fig. 4a.
Statistical analysis. Nested ANOVA. Nested (hierarchical) ANOVA is a parametric hypothesis testing and an extension of 1-way ANOVA. A nested ANOVA is used when there is one measurement variable and more than one nominal variable, and the nominal variables are nested 35 . The nominal variables being nested means that each value of one nominal variable (the subgroups) is found in combination with only one value of the higher-level nominal variable (the groups). We considered lower-level variables (cross-section, axon, animal) as random effects variables and the top level (group, either sham or TBI) as a fixed effect variable. The null hypotheses were whether there existed significant variation in means among groups at each level. The analysis was performed using the anovan function of MATLAB R2017b, with type II sum of squares. Because, the design matrix of anovan grows quickly when the nesting levels increases, we measured the median of the cross-sectional quantities and assigned them to the lowest level of nested ANOVA. At the cross-sectional level, the distributions of equivalent diameter, minor and major axes and eccentricity were multimodal ( Supplementary Fig. S5), thus median of cross-sections was preferred to mean. The analysis was performed separately for two hemispheres.
Variance components. Nested ANOVA partitions the variability of the measurements into different levels. The variance components describe what percentage of the total variance is attributable to each level 35 .

Data Availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. The source code of ACSON is available at https://github.com/aAbdz/ACSON.