Automated 3D Axonal Morphometry of White Matter

Abstract

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.

Introduction

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)⁶. SBEM combines scanning electron microscopy (SEM) with back-scattered electron detection and low beam energies¹⁰. 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¹², TrakEM2¹³, Microscopy Image Browser¹⁴, and CATMAID¹⁵), or interactive processing of data by combining automated analysis and proof-reading capabilities (e.g., rhoANA¹⁶, ilastik¹⁷, and SegEM¹⁸). 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²³, 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 matter components from 3D-SBEM datasets. The previous segmentation methods cannot be used to address the segmentation of white matter for several reasons. First, using manual or semi-automated segmentation software tools (e.g.²⁷, TrackEM2¹³ and ilastik¹⁷) or pipelines (Chklovskii et al.²⁰), would be prohibitively time consuming. Providing annotated data, i.e., training data, for supervised learning-based segmentation methods, such as SegEM¹⁸ and SyConn²⁸, is also very time consuming, or requires several annotators. Second, our interest lies in segmentation of several white matter constituents, such as myelin, myelinated and unmyelinated axons, cell bodies, mitochondria, and vacuoles. The methods utilizing binary pixel classifiers^23,24, which assign a pixel as either a cell boundary or a cell-interior, are inappropriate for this type of multiclass segmentation problem. Especially if subcellular structures, such as mitochondria, are not labeled separately, the clustering step of these methods fails to correctly merge regions within a complete cell²⁶. Some studies have addressed the multiclass segmentation^22,25,26, and segmented mitochondria as a subcellular structure. These approaches, however, are only valid for grey matter, which does not contain myelin. In the SBEM images of white matter, the segmentation of mitochondria in the presence of myelin is extremely difficult because the signal intensity and textural features of mitochondria and myelin are highly similar. Interestingly, a previous study tracked axons in a SBEM volume of the optic tract using Kalman-snakes²⁹ initialized either manually, or automatically using watershed filtering³⁰. However, this approach fails in tracking full length of axons throughout the SBEM volume. Therefore, the automated segmentation of SBEM images of white matter requires a specifically developed method to address these problems. In conclusion, currently there exist no automated methods to quantify the axonal morphometry in 3D EM images.

We developed a novel pipeline for AutomatiC 3D Segmentation and morphometry Of axoNs (ACSON) in mesoscale SBEM volumes of white matter. The automated pipeline eliminates the need for time-consuming manual segmentation of 3D datasets and enables full 3D analysis of the white matter ultrastructure. ACSON segments the main cellular and subcellular components of the corpus callosum. To confirm the accuracy of the automated segmentation, the automated segmentation was evaluated against manual annotation by an expert. We quantified the actual cross-sections of the segmented myelinated axons based on their diameter and eccentricity. We analyzed the morphological features of SBEM datasets from the ipsilateral and contralateral sides of the corpus callosum in two sham-operated and three TBI rats.

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³¹ 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.

Figure 1

White matter ACSON segmentation pipeline. (a) A 2D representative image of the SBEM dataset from the contralateral corpus callosum of Sham-1 dataset. (b) The same image denoised with BM4D. (c) Boundary mask B: union of myelin obtained with BVG and dilated Canny edges. (d) The Euclidean distance transform was calculated for every slice of B individually to define the location of the seeds for the BVG algorithm. (e) The primary result of BVG segmentation. As the volume of cells/cell process exceeded ϑ, voxels corresponding to cells/cell process remained unlabeled in the primary segmentation (black regions). (f) The SLIC supervoxel technique refined the segmentation. (g) Segmentation of cells and mitochondria refined by the SLIC supervoxels. (h) A 4 labels map of myelin, intra-axonal space of myelinated and unmyelinated axons, oligodendrocyte cell body and its processes. (i) 3D rendering of myelinated axons in the contralateral corpus callosum of a Sham-1 dataset. The 3D rendering depicts myelinated axons with different thicknesses running along the dataset and organizing bundles with different orientations.

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.

Figure 2

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.

Table 1 We evaluated ACSON segmentation by comparing it to the manual expert segmentation using precision (Prec.), recall (Rec.), weighted Jaccard index (WJI), and weighted Dice coefficient (WDC) metrics. The precision and recall metrics evaluate the segmentation in the tissue-type level and weighted Jaccard index and weighted Dice coefficients evaluate each axon as its own region. Therefore, the latter metrics are very sensitive to split-and-merge type errors.

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, Supplementary video shows the 3D rendering of myelin, cell, and myelinated axons in the contralateral corpus callosum of sham-1 dataset.

ACSON morphometry pipeline automatically quantifies the segmented myelinated axons

The ACSON morphometry pipeline quantified cross-sections of the intra-axonal space of myelinated axons. A cross-section is the intersection of a segmented myelinated axon and a perpendicular plane to the axonal skeleton. To detect the axonal skeleton, we defined three points in the myelinated axon domain: one with the largest distance from the axon surface and two endpoints of the axons, i.e., the tips of the axon. The minimum-cost path connecting these three points was defined as the axonal skeleton (see Fig. 3a, which shows the cross-sectional planes at three randomly selected positions along an axon). The orientation of the cross-sectional planes at each skeleton point was determined by the unit tangent vector at that point. For each cross-section, we measured the length of minor and major axes of the fitted ellipses, equivalent diameter³³ and eccentricity. The equivalent 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³⁴. 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.

Figure 3

ACSON morphometry. (a) 3D reconstruction of a representative intra-axonal space of one myelinated axon and its mitochondria. Three intersecting planes at randomly selected positions show the cross-section of the myelinated axon. (b) The cross-sectional diameter of a myelinated axon varies along its length.

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.

Figure 4

A comparison between the 2D and 3D morphological analyses. (a) 3D reconstruction of intra-axonal space of a representative myelinated axon with an axis that does not align with image cross-sections. The myelinated axon was quantified at three planes parallel to the x-y plane simulating the 2D morphometry and three cross-sections for the 3D morphometry. The quantified parameters were minor axis (MinAxis), major axis (MajAxis), equivalent diameter (EqDiameter), and eccentricity (Ecc). When the angle between the normal planes to the axonal axis and the image cross-sections grew, the relative difference (RD) between the 2D and 3D quantifications increased (p2 and p3). (b-d) Percentage of the relative difference between the 2D and 3D measurements of (b) the minor axis [median: 16.64%, median absolute deviation (MAD): 10.66%], (c) the major axis [median: 18.25%, MAD: 16.68%], (d) the equivalent diameter [median: 16.23%, MAD: 9.60%], and (d) the eccentricity [median: 11.34%, MAD: 6.30%] for all myelinated axons in contralateral corpus callosum of Sham-1 dataset.

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³⁵. 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.

Figure 5

Box-plots of equivalent diameter (EqDiameter), minor axis (MinAxis), major axis (MajAxis), and axonal eccentricity measured for myelinated axons in sham-operated and TBI animals. The measurements are medians of the cross-sectional measurements of each myelinated axon. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. 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.

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 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.

Table 2 Variance components analysis. The variance was partitioned into cross-sectional, axonal, and animal levels.

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³⁶ with no measurement of the myelin thickness³⁷, denoted as $aggregate\,g-ratio=\sqrt{1-Myeli{n}^{\ast }}$. We defined Myelin* as the ratio of the myelin volume to the myelin volume plus the intra-axonal volume of all myelinated axons. In addition, we calculated the density for myelinated axons and mitochondria defined as ${\rho }_{N}=\frac{N}{{V}_{nc}}$, 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⁴¹ and in the hippocampus and cerebellum⁴² indicated highly irregular axial shapes with periodic varicosities containing organelles⁴². 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.⁴⁷ and Fieremans et al.⁴⁸ have shown the effect of structural disorder along axons and micron-level sample architecture while modeling dMRI. Thus, realistic modeling of cells⁴⁹ and myelinated axons⁵⁰ 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.³² 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)⁵¹. 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⁵² measuring axon diameter in the rat corpus callosum. Kim et al.⁵² 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 mitochondria, increasing the cross-sectional diameter locally³⁴. Wake et al.⁵³ 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³⁵. Also, averaging the lowest levels of hierarchy has sub-optimal statistical power⁵⁴. 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⁵⁵ 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⁵⁸ and for up to 1 year in rats⁵⁹. 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⁶⁰ or regeneration⁵⁵. 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⁶¹. 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)⁶². 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.

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⁶³ (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² base and an approximately 600 × 600 μm² 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¹⁴ 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⁶⁴. To suppress these complex noise patterns, we applied a non-local BM4D algorithm⁶⁵. 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 $\hat{E}\subset X$ as the set of edges of z, using a Canny edge detector⁶⁶. We set the parameter values of the Canny edge detector as follows: the SD of the Gaussian filter was $\,\sqrt{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 $\hat{E}$ 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₁, V₂, …, V_n, E ⊂ X, in which V_i ∩ V_j = $\varnothing$, V_i ∩ E = $\varnothing$, ∀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 $\delta (x)=|z(x)-\frac{1}{|{S}_{k}|}\sum _{s\in {S}_{k}}z(s)|$ and set the similarity threshold δ_T to 0.1. An iteration terminates, if A = $\varnothing$, 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₁. BVG was initiated with a random low-intensity voxel S₁ with z(S₁) ≤ 0.4. This one seed was sufficient to segment V₁ 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₁). To segment other structures V₂, …, 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₁ (Fig. 1c). Denoting each 2D-slice of B as B_i, we computed the Euclidean distance transform⁶⁷ 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⁶, which equals 12.5 μm³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³ 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⁶⁸ 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}_{s}={d}_{int}+\frac{c}{\rho }{d}_{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³ 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 ${V}_{i}^{^{\prime} }=\mathop{\cup }\limits_{q\in I}S{V}_{q}$ where $I=\{q^{\prime} |\frac{|S{V}_{q^{\prime} }\cap {V}_{i}|}{|S{V}_{q^{\prime} }|}\ge \mathrm{0.8,}\forall j < i:\frac{|S{V}_{q^{\prime} }\cap {V}_{j}|}{|S{V}_{q^{\prime} }|}\mathrm{ < 0.8}\}$. 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}_{i}^{^{\prime} }{\mathrm{| > 10}}^{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₁ and all mitochondria, denoted as MY. In our SBEM-datasets, vacuoles appeared brighter than all of the other ultrastructures. Thus, if $|{V}_{i}^{^{\prime} }| < 2\times {10}^{4}$ and $\frac{1}{|{V}_{i}^{^{\prime} }|}\sum _{\forall s\in {V}_{i}^{\text{'}}}z(s)\ge 0.85$, we labeled ${V}_{i}^{^{\prime} }$ as a vacuole. We defined the remaining volumes, not mitochondria nor vacuole, as axons denoted as AX_i, i = 1, …, m. To distinguish if an axon AX_i was a myelinated or an unmyelinated axon, we studied a thick hollow cylinder enclosing the 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 $\frac{|(\mathop{\cup }\limits_{q\in {{\rm{\Lambda }}}_{i}}S{V}_{q})\,\cap MY|}{|{(\mathop{\cup }\limits_{q\in {{\rm{\Lambda }}}_{i}}SV)}_{q}|}\ge 0.7$, 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⁶ for the segmentation of V₂, …, 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 γ⁶⁹. To detect the myelinated axon skeleton γ with sub-voxel precision, we adapted a method from Van Uitert & Bitter⁷⁰. 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,

$$\gamma ={\rm{\arg }}\mathop{{\rm{\min }}}\limits_{P}{\int }_{{x}_{e1}}^{{x}_{e2}}H(P(\varsigma ))d\varsigma ,$$

(1)

where ς traces the path P, and H is the cost function. To enforce the minimum-cost path to run at the middle of the object, the cost function H should be higher if the path moves away from the center. Points x*, x_e1, and x_e2 and solving equation (1) was defined by solving an eikonal equation on the axonal domain Ω. The eikonal equation is a non-linear partial differential equation defined as a special case of wave propagation in which the front Γ advances monotonically with speed F(x) > 0. The eikonal equation can be formulated as |∇T(x)|F(x) = 1, where T|_Γ = 0. The solution, T(x), is the shortest time needed to travel from Γ to any point x ∈ Ω, with the speed F(x) > 0. Although the eikonal equation can be solved with the FMM⁷¹, we used 3D MSFM⁵¹. MSFM combines multiple stencils and second-order approximation of the directional derivatives over the FMM to improve the accuracy of solving the eikonal equation on Cartesian domains. To find x*, we computed the time-crossing map T₁(x) from the myelinated axon interface Γ with constant speed F₁(x) = 1, x ∈ Ω. The global maximum of T₁, where T₁(x) ≤ T₁(x*), ∀x ∈ Ω was defined as x*. To find x_e1, x_e2, and γ, we calculated a new time-crossing map T₂(x), starting at x* to every voxel in Ω, with a non-constant speed ${F}_{2}(x)={(\frac{{T}_{1}(x)}{{T}_{1}({x}^{\ast })})}^{2}$ for x ∈ Ω. T₁(x)|_Γ = 0. Using H(x) = 1 − F₂(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₂ map, i.e., the global maximum of T₂, 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₃(x), starting from x_e1 to every voxel in Ω, with speed F₂(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 $-\frac{\nabla {T}_{2}(x)}{|\nabla {T}_{2}(x)|}$ until x* was reached. x* is the global minimum of T₂(x), so that we were guaranteed to find it with backtracking. The backtracking procedure can be described by the ordinary differential equation $\frac{d{\gamma }_{i}}{d\varsigma }=-\frac{\nabla {T}_{2}}{|\nabla {T}_{2}|},\,{\gamma }_{i}(\varsigma =0)\mathrm{=}{x}_{ei}$, 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 γ = γ₁ ∪ γ₂. Note that computing the skeleton in this way prevented the skeleton from cutting corners⁷². 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₂(x) and T₃(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₁, x₂) was foreground iff x ∈ C. We defined the center point of C as $c=\frac{1}{|C|}\sum _{x\in C}x$. 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⁷³, 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 as in the automatic one. It also annotated each segmented region as myelin, myelinated or unmyelinated axon. In the annotated images, mitochondria and vacuoles were included into intra-axonal space.

Precision and recall

For a tissue-type level evaluation, we used the precision and recall as in the previous studies^28,32. Let A and B be the sets of voxels of a particular tissue-type (myelin, myelinated axon, unmyelinated axon) in the manual and automated segmentations, respectively. We defined $Precision=\frac{|A\cap B|}{|B|}$, and $Recall=\frac{|A\cap B|}{|A|}$. 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.

Weighted Jaccard index and weighted Dice coefficient

To further evaluate the automated segmentation, we used Jaccard index and Dice coefficients in the region level. The Jaccard index⁷⁴ and Dice coefficient⁷⁵ is defined by $J(A,B)=\frac{|A\cap B|}{|A\cup B|}$ and $Dice\,Coef(A,B)=\frac{\mathrm{2|}A\cap B|}{|A|+|B|}$, 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 $\sum _{i\mathrm{=1}}^{a^{\prime} }{w}_{i}J({A}_{i},{B}_{b(i)})$, and $\sum _{i\mathrm{=1}}^{a^{\prime} }{w}_{i}Dice\,Coef({A}_{i},{B}_{b(i)})$, respectively, where b(i) is the index of the region best matching A_i and the weight ${w}_{i}=\frac{|{A}_{i}|}{|{\sum }_{i\mathrm{=1}}^{a^{\prime} }{A}_{i}|}$.

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 $relative\,difference=\frac{|{q}_{2D}-{q}_{3D}|}{{\rm{\max }}({q}_{2D},{q}_{3D})}$, 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³⁵. 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³⁵.