DeepACSON: Automated Segmentation of White Matter in 3D Electron Microscopy

To trace the entirety of ultrastructures in large-scale three-dimensional electron microscopy (3D-EM) images of brain tissue, automated segmentation techniques are essential. The current automated techniques use deep convolutional neural networks (DCNNs), approaching human-level segmentation, but designed for tracing small numbers of neurites and for high-resolution datasets. Therefore, we developed DeepACSON to segment 3D-EM datasets of white matter at low-resolutions and for tens of thousands of myelinated axons. DeepACSON performs DCNN-based semantic-segmentation and shape decomposition-based instance-segmentation: a top-down design to account for severe membrane discontinuities inherited with low-resolution imaging. We applied DeepACSON to ten serial block-face scanning electron microscopy datasets of rats after sham-operation or traumatic brain injury, segmenting hundreds of thousands of long-span myelinated axons, thousands of cell nuclei, and millions of mitochondria with excellent evaluation scores. DeepACSON quantified the morphological and spatial aspects of white matter ultrastructures.


Introduction
Recent advances in automated serial-sectioning electron microscopy (EM) techniques allow for acquiring 3-dimensional (3D) image datasets of brain ultrastructures, with high-resolutions and from hundreds of micrometers of tissue [1][2][3] . The acquired image datasets are of big sizes, ranging from gigabytes to hundreds of terabytes or even petabytes of data. Converting such big datasets into quantitative information, such as mapping neuronal connections known as connectome or morphological data of the ultrastructures, requires annotation of individual components. But manual annotation of even a small section of a 3D-EM dataset is a tedious task, consuming thousands of hours of experts' time. For example, Berning et al. 3 reported that manual annotation of 215 neurites in 5 × 10 8 voxels required 1 500 h of time of trained students, or we estimated that the segmentation of 15 × 15 × 15 µm 3 white matter ultrastructures requires 2 400 h of an expert's time 4 . Semi-automated segmentation methods based on machine learning approaches 3,5 have improved the rate of segmentation. However, these methods still require a considerable amount of manual interaction as the segmentation is driven on the manually extracted skeletons of neurites, proofreading, and correction of errors.
Deep convolutional neural networks (DCNNs) have emerged as the key to the state-of-the-art automated segmentation techniques [6][7][8][9][10][11] . DCNN-based segmentation methods generally comprise two steps: generating a membrane probability map from 3D-EM datasets referred to as semantic-segmentation, followed by instance-segmentation that turns probability maps into individual object instances. DCNNs are typically used for semantic-segmentation, and more traditional image analysis techniques favoring a bottom-up design, i.e., over-segmentation and subsequent merge, are used for instance-segmentation.
An example is DeepEM3D 7 , applying a 3D watershed transform on the membrane probability maps. To avoid instancesegmentation of a two-step design, Funke et al. 8 proposed learning voxel affinities. Their method, however, continued with a 3D watershed segmentation and iterative region agglomeration. Recently, Januszewski et al. 9 suggested flood-filling networks, a single-object tracking technique, and Meirovitch et al. 11 introduced cross-classification clustering, a multi-object tracking technique, merging the semantic-and instance-segmentation into recurrent neural networks. These techniques have a bottom-up design, offering a sophisticated training-based alternative for over-segmentation and subsequent merge.
The bottom-up design of the state-of-the-art EM segmentation techniques [6][7][8][9][10][11] is subjected to local optimization. Hence, segmentation requires distinctive image features such as high contrast between the cellular membrane and other structures.
Therefore, common to these techniques is that they are operating on 3D-EM images with high in-plane resolutions of 4-10 nm  . Another limitation of the bottom-up design is that it is computationally inefficient for segmenting tens or hundreds of thousands of neuronal processes on over a million cubic micrometers of tissue. In the segmentation of low-resolution 3D-EM datasets of white matter, we encountered both of the above mentioned difficulties: first, the membrane of myelinated axons is not resolved, meaning that there is no distinctive feature to differentiate the intra-and extra-axonal space of a myelinated axon at locations such as nodes of Ranvier. Second, tens of thousands of myelinated axons pass through a volume of about a million cubic micrometers. Therefore, we developed DeepACSON, a Deep learning-based AutomatiC Segmentation of axONs. DeepACSON is also a two-step segmentation technique, but it approaches the segmentation problem from a top-down perspective, i.e., under-segmentation and subsequent split. This approach is made possible by the use of a-priori knowledge of the topology of myelinated axons and cell nuclei. DeepACSON performs DCNN-based semantic-segmentation and 3D shape analysis based instance-segmentation.
We obtained our 3D-EM datasets using serial block-face scanning electron microscopy (SBEM) from the corpus callosum and cingulum of rats after sham-operation or traumatic brain injury (TBI). The SBEM images were acquired at two resolutions simultaneously: high-resolution images, 15 × 15 × 50 nm 3 voxel size in a small field-of-view 15 × 15 × 15 µm 3 , and lowresolution images, 50 × 50 × 50 nm 3 voxel size in a large field-of-view 200 × 100 × 65 µm 3 . We segmented the high-resolution datasets using our earlier automated ACSON pipeline 4 . We used the ACSON segmentation of the high-resolution datasets to train DeepACSON, eliminating the need for manually annotated training sets. DeepACSON segmented ten large field-of-view SBEM datasets that sum to 1.09 × 10 7 µm 3 of white matter tissue into myelin, myelinated axons, mitochondria, and cell nuclei. From these ten datasets, we segmented about 288 000 myelinated axons span to 9 m, 2 600 cell nuclei, and millions of mitochondria with excellent evaluation scores. Our segmentation of white matter enabled quantifying axonal morphology, such as the axonal diameter, eccentricity and tortuosity, and also spatial organization of ultrastructures, such as the spatial distribution of mitochondria, and the density of cells and myelinated axons.

Results
DeepACSON segmented and analyzed ten low-resolution large field-of-view SBEM datasets. We acquired the images from the ipsi-and contralateral of white matter from two sham-operated rats and three TBI rats. Figure 1 illustrates different steps of DeepACSON, where DCNNs perform semantic-segmentation, and 3D shape-analysis techniques perform instancesegmentation.

Datasets
We acquired SBEM images of the white matter simultaneously at low-and high-resolution (Fig. 2a). The images were from the ipsi-and contralateral corpus callosum of the sham-operated and TBI animals. The low-resolution datasets were 200 × 100 × 65 µm 3 in size, with a voxel size of 50 × 50 × 50 nm 3 . Two-thirds of the low-resolution images were from the corpus callosum and one-third from the cingulum (Supplementary Table S1). The high-resolution datasets were 15 × 15 × 15 µm 3 in size and imaged with an anisotropic voxel size of 15 × 15 × 50 nm 3 . The high-resolution images were acquired from the corpus callosum. Figure 2a shows the contralateral corpus callosum and cingulum of a sham-operated rat in low-and high-resolution.
In Fig. 2b, we show images of the low-and high-resolution datasets from the same location. Ultrastructural components such as myelin, myelinated axons, mitochondria, and cell nuclei were resolved in both settings. However, other components such as unmyelinated axons (2b fuchsia panel, asterisks) and exposed axonal membrane, e.g., at nodes of Ranvier (2b cyan panel, arrowheads), were only resolved in the high-resolution datasets. Figure 2b (purple panel) shows a cell nucleus from the low-resolution dataset, the membrane of which is resolved, but not continuously.

Semantic-segmentation of white matter ultrastructures
We segmented the high-resolution SBEM datasets using our ACSON pipeline 4 . First, we down-sampled the high-resolution datasets and their corresponding segmentation to the resolution of low-resolution datasets to train the DCNNs (Fig. 1, step 1). Figure 2c shows a 3D rendering of myelinated axons of a high-resolution SBEM dataset segmented by the ACSON pipeline.
Only four of the total ten datasets were used as training data, while the remaining six datasets were reserved for testing. We trained two DCNNs denoted as DCNN-mAx and DCNN-cN. DCNN-mAx generated the probability maps of myelin, myelinated axons, and mitochondria, and DCNN-cN generated the probability maps of cell nuclei and the membrane of cell nuclei (Fig. 1, step 2). The receptive field of DCNN-mAx was (45,45,15) voxels in (x, y, z) directions with 2 881 978 trainable parameters, and the receptive field of DCNN-cN was (45,45,15) voxels in (x, y, z) directions with 2,881,777 trainable parameters. We trained the DCNNs on a single NVIDIA Tesla P100-16 GB graphics processing unit (GPU) for one day. From 85 GB of aligned Step 1: The ACSON segmentation of the high-resolution SBEM datasets, down-sampled to the resolution of the low-resolution datasets, was used for training DeepACSON. We trained two DCNNs denoted as DCNN-mAx and DCNN-cN.
Step 2: DCNN-mAx returned the probability maps of myelin, myelinated axons, and mitochondria, and DCNN-cN returned the probability maps of cell nuclei and the membrane of cell nuclei.
Step 3: The segmentation of myelin was finalized by thresholding the myelin probability map. We performed the initial segmentation of myelinated axons by the binarization and connected component analysis. Subsequently, the geometry of the segmented components was rectified using our newly developed cylindrical shape decomposition 12 (CSD) technique. We performed the segmentation of cell nuclei in a geometric deformable model (GDM) framework by applying elastic deformations to the initial segmentation of cell nuclei.
Step 4: The segmentation of myelinated axons and cell nuclei was finalized by eliminating non-axonal and non-nucleus structures with support vector machines (SVMs). and denoised low-resolution SBEM datasets, we acquired 1.

Instance-segmentation of myelin, myelinated axons, and mitochondria
In low-resolution SBEM datasets, myelin was a single connected ultrastructural component. The segmentation of myelin was accomplished by binarizing its probability map (see Fig. 1, step 3 and Fig. 3b). To segment the intra-axonal space of myelinated axons, we binarized their probability map and applied connected component analysis. This preliminary segmentation was under-segmented (divided into too few segments) because the exposed axonal membrane was not fully resolved, e.g., at nodes (a) We acquired SBEM images of the white matter, corpus callosum (cc) and cingulum (cg), simultaneously at high-and lowresolution. The field-of-view of the low-resolution dataset is 204.80 × 102.20 × 65.30 µm 3 equivalent to 4 096 × 2 044 × 1 306 voxels in x, y and z directions respectively, which is about 400 times larger than the field-of-view of the high-resolution datasets. (b) Images of the low-and high-resolution datasets acquired from the same location (the orange-rendered volume in (a)). The visualization of the high-and low-resolution images shows that myelin, myelinated axons, mitochondria, and cell nuclei were resolved in both settings, while the axonal membrane at nodes of Ranvier (cyan panel, arrowheads) and unmyelinated axons (fuchsia panel, asterisks) were only resolved in the high-resolution images. The purple panel shows a cell nucleus from the low-resolution dataset (a), which membrane is resolved, but not continuously. (c) A 3D rendering of myelinated axons in the high-resolution SBEM dataset (contralateral sham #25) segmented by the automated ACSON pipeline.
of Ranvier (Fig. 3c). Therefore, the preliminary segmentation of the intra-axonal space of myelinated axons required a further shape analysis step ( Fig. 1, step 3). We devised an algorithm for the decomposition of cylindrical objects, called cylindrical shape decomposition 12 (CSD) ( Supplementary Fig. S5). The technique decomposes under-segmented myelinated axons and subsequently rectifies their shape using generalized cylinders (Fig. 3c, Fig. 3c 1 , and Supplementary Fig. S6). To finalize the segmentation, we excluded non-axonal instances from the segmentation using a support vector machine (SVM) with a quadratic kernel (as implemented in MATLAB, Statistics and Machine Learning Toolbox, fitcsvm function) (Fig. 1, step 4).

Instance-segmentation of cell nuclei
The membrane of cell nuclei was not fully resolved at low-resolution (50 × 50 × 50 nm 3 voxel size), and hence its probability map exhibited discontinuity. We segmented cell nuclei in a geometric deformable model (GDM) framework, where the initial segmentation of a cell nucleus was rectified for its topological errors (see Fig. 1, step 3 and Fig. 4c and d). Non-nucleus instances were excluded from the segmentation using an SVM with a quadratic kernel (Fig. 1, step 4). Figure 4e shows a 3D rendering of cell nuclei in the contralateral corpus callosum and cingulum of a sham-operated rat dataset.

Evaluations
We evaluated DCNN-mAx on a tissue-type level against a test set: six unseen high-resolution SBEM datasets labeled by the automated ACSON pipeline, and down-sampled to the resolution of low-resolution datasets. We calculated the precision (positive predictive value), recall (sensitivity), and F1 score (harmonic mean of precision and recall) for the segmentation  Table S2). The average F1 score for myelin was 0.886 ± 0.049, and for myelinated axons was 0.861 ± 0.037.
We evaluated the performance of the two SVMs by a leave-one-group-out (LOGO) cross-validation, where the classifier was trained excluding the data from one group of animals (sham-operated or TBI) from training and evaluated against it (Supplementary Table S2). F1 scores were 0.988 and 0.955 for eliminating non-axonal structures while excluding sham datasets and TBI datasets, respectively. For eliminating non-nucleus structures, F1 scores were 0.969 and 0.981, excluding sham datasets and TBI datasets, respectively. The LOGO cross-validation showed that the performance of the SVMs was equally robust in all datasets, regardless of the experimental condition, i.e., sham-operated or TBI.
An expert (A.S.) evaluated DeepACSON segmentation of myelinated axons and mitochondria in an object-level using a GUI-based software tool we developed for this purpose ( Supplementary Fig. S1). For this, we randomly sampled every dataset for five non-overlapping windows of size 300 × 300 voxels (10 datasets, 50 samples). The expert had no access to the dataset ID nor the sampling location. The expert evaluated the sampled images of the final segmentation by counting the number of true-positives, false-positives, and false-negatives. This resulted in following evaluation scores: myelinated axons (precision: 0.965 ± 0.027, recall: 0.877 ± 0.061, and F1 score: 0.918 ± 0.038) and mitochondria (precision: 0.856 ± 0.100, recall: 0.804 ± 0.091, and F1 score: 0.823 ± 0.067).

White matter 3D morphology analysis
We quantified several morphological and volumetric aspects of the segmented ultrastructures to demonstrate the applications of the white matter segmentation (Fig. 5). For every myelinated axon, we extracted cross-sections along its axonal skeleton, i.e., central axis, with a plane perpendicular to the skeleton. The cross-sections were quantified by the equivalent diameter, i.e., the diameter of a circle with the same area as the cross-section, and the eccentricity, i.e., how much a conic section deviates from being circular, of the fitted ellipse. (Fig. 5a and b). We measured the tortuosity of myelinated axons as τ = l Gd l Ed , where l Gd is the geodesic distance and l Ed is the Euclidean distance between the two endpoints of the axonal skeleton. To quantify the spatial distribution of mitochondria, we measured the inter-mitochondrial distances along each myelinated axon. For this, we projected the centroid of mitochondria on the axonal skeleton and measured the geodesic distance between the consecutive projected centroids (Fig. 5a and c).
The analysis indicated that the diameter of myelinated axons varies substantially along the axonal skeleton. Moreover, the distribution of diameters was bimodal, which can partially be related to the location of mitochondria 13 (Fig. 5a, b, and d). The cross-sections of myelinated axons were elliptic rather than circular (Fig. 5d). Myelinated axons were almost straight as the mean of τ in the corpus callosum was 1.068 and in the cingulum was 1.072, but there was a big variation in the tortuosity as the standard deviation of τ was 0.731 and 0.396 in the corpus callosum and cingulum, respectively (see Fig. 5a and d).
The distribution of the inter-mitochondrial distance along a myelinated axon was bimodal because mitochondria were either accumulated or appeared distant from each other (Fig. 5a, c, and d).
For the statistical hypothesis testing between the sham-operated and TBI animals, myelinated axons were represented by the median of equivalent diameters, the median of the eccentricities, tortuosity, and the mean of inter-mitochondrial distances ( Fig. 5d). We subjected the measurements to the nested (hierarchical) 1-way analysis of variance (ANOVA) separately for each hemisphere 14 . 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 the rats' ID, and the rats' ID was nested under the groups (sham-operated and TBI). We set the alpha-threshold defining the statistical significance as 0.05 for all analyses. The rats' IDs were treated as random effects, and the group was treated as a fixed effect. The equivalent diameter was significantly smaller in the ipsi-(F = 16.27, p = 0.027) and contralateral cingulum (F = 29.28, p = 0.011) and in the ipsilateral corpus callosum (F = 15.75, p = 0.047) of TBI rats as compared to sham-operated rats. Also, the tortuosity in the ipsilateral cingulum of TBI rats was significantly greater than sham-operated rats (F = 25.23, p = 0.018). The equivalent diameter of the contralateral corpus callosum of TBI rats was slightly smaller than sham-operated rats, but not significantly (F = 5.78, p = 0.095). In the ipsilateral cingulum, the inter-mitochondrial distance was slightly smaller in TBI rats, but not significantly (F = 6.27, p = 0.086). We did  Figure 5g shows examples of 3D-rendered myelinated axons in the cingulum, demonstrating the high variability of axonal diameters. Fig. 5h demonstrates the sparsity of myelinated axons due to the injury in the cingulum and corpus callosum of TBI rats.

Computation time
The pipeline required about five days to segment one raw low-resolution SBEM dataset of size 4 000 × 2 000 × 1 200 voxels into its final segmentation. The details are presented in Supplementary Table S4.

Discussion
We developed DeepACSON applying DCNN-based semantic-segmentation and introducing 3D-shape analysis as an instance- proposing a 3D watershed transform and subsequent merge and continuing with detecting X-shape objects to find under-segmentation errors. Unlike this method, which makes restrictive assumptions about the geometry of under-segmented objects, we address a generic geometry for the under-segmented objects as X-shape objects are a special case in our CSD algorithm. Also, we reconstruct/rectify an under-segmented object on its own by applying generalized cylinders, as opposed to MaskExtend, that requires accessing the watershed segments and merging them at a high merge-threshold 10 . At a more general level, Falk et al. 18 developed a plugin in ImageJ software 19 and addressed under-segmentation errors in cells by inserting an artificial one-pixel wide background ridge between touching instances in the training set. This approach is similar to our technique for segmenting cell nuclei, where we inserted an artificial membrane around cell nuclei using morphological operations and assigned high-weights while training. Another approach is CDeep3M 6 , which uses a cloud-based residual inception-network architecture for the segmentation of image datasets from diverse modalities. This implementation assumes to obtain highly accurate semantic-segmentation from the network, hence continues to instance-segmentation with thresholding.
This assumption ignores the topological errors of the segmentation. The simultaneous high-and low-resolution imaging enables to use the high-resolution images as the training data to segment low-resolution datasets. The high-resolution images can be segmented using automated methods; hence, we provided a human-annotation-free training set. In our case, we used the ACSON pipeline to segment the high-resolution SBEM datasets 4 .
ACSON is fully automated, requires no training data, and segments one of our high-resolution datasets, containing about 400 myelinated axons, in approximately 24 hours. Furthermore, high-resolution imaging resolves the membrane of tissue ultrastructures in several voxels wide, as opposed to low-resolution imaging that does not fully resolve the membrane, resulting in a more accurate segmentation of ultrastructures. Accurate segmentation of ultrastructures yields generating an accurate training set, which technically means that we have smaller label-noise, i.e., correlating confusing labels with underlying tissue ultrastructures, in the training set.
The validation of DeepACSON demonstrated an equally robust segmentation performance in datasets of both sham-operated and TBI rats, remarking that abnormalities in the shape of myelinated axons, such as myelin delamination, are more frequent in TBI datasets 4 . As well, the expert evaluation showed that the sensitivity of DeepACSON was smaller than its precision; we had a bigger number of false-negatives, i.e., myelinated axons that were not segmented, than false-positives. We speculate that very thin myelinated axons have been over-segmented into small-components. And because we filtered out small-components before the shape decomposition step, it increased the number of false-negatives. For example, an axon with a diameter smaller than 0.2 µm was resolved with less than 13 cross-sectional voxels at 50 × 50 × 50 nm 3 resolution, which was prone to over-segmentation. Segmentation of myelinated axons at low-resolutions can be difficult for an expert as well. In Mikula et al.
study 21 , an expert made an error in every 80 nodes of Ranvier in white matter, where the resolution was 40 × 40 × 40 nm 3 and for the axonal diameter greater than 0.5 µm.
Large field-of-view imaging makes possible quantifying parameters whose measurement in a small field-of-view is not reliable, because the measurement may reflect a very local characteristic of the underlying ultrastructure. Therefore, we quantified the tortuosity of myelinated axons, inter-mitochondrial distance, and cell density only in the large field-of-view datasets. We also evaluated DeepACSON quantification. For that, we compared DeepACSON measurements of the axonal diameter and the eccentricity in the low-resolution datasets with measurements from the high-resolution datasets in our earlier work 4 . This comparison indicated that the quantification in the low-resolution datasets was consistent with the high-resolution We quantified morphological changes in myelinated axons in white matter five months after a severe TBI. The white matter pathology after TBI is extremely complex. The initial response of an axon to a brain injury can be either degeneration or regeneration 22,23 , however, morphological alterations in axons can persist for years after injury in humans 24,25 and for up to one year in rats 26 . We found that the axonal diameter of myelinated axons in the ipsilateral corpus callosum, ipsilateral cingulum, and contralateral cingulum of TBI rats has significantly reduced. We further measured that the density of myelinated axons has reduced significantly in TBI, which reflects degeneration of axons after injury 25 . We also found that the tortuosity of myelinated axons has increased in the ipsilateral cingulum of TBI animals. We speculate that prolonged damage in microtubules might underlie the increase in axonal tortuosity 27 .
Ultrastructural tissue modeling is an active research field aiming to bridge the gap between macroscopic measurements and cellular-and sub-cellular tissue level. Currently, tissue modeling are based on simplistic representations and assumptions of the ultrastructural features. These assumptions typically assume axons to be perfect cylinders, or the variation of the axonal diameter along axons are neglected 28,29 . Segmentation of brain tissue ultrastructures in 3D-EM datasets can substitute simplistic biophysical models and provide realistic models. Such realistic 3D tissue models open the possibility to investigate underlying reasons for the diffusion magnetic resonance imaging contrast and its macroscopic changes in brain diseases [30][31][32][33] or to investigate conduction velocity in myelinated and unmyelinated axons in electrophysiology 34,35 .
The animals were singly housed in a room ( The brains were sectioned into 1-mm thick coronal sections with a vibrating blade microtome. Sections from −3.80 mm from bregma from each brain were selected and further dissected into smaller samples containing the areas of interest. We collected two samples from each brain: the ipsi-and contralateral of the cingulum and corpus callosum. The samples were osmium stained using an enhanced staining protocol 37 . After selecting the area within the samples, 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 microscope. Silver paint 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. The details of the animal model and tissue preparation are described in ACSON study 4 .
The face of the blocks was in the x-y plane, and the cutting was in z direction. All blocks were imaged with a beam voltage of 2.5 kV and a pressure of 0.15 Torr. We acquired the high-and low-resolution datasets consistently at one specific location in the white matter in both sham-operated and TBI animals and in both hemispheres. The images were collected with an unsigned 16-bits per voxel. Supplementary Table S1 shows the volume size of the low-resolution datasets.

DeepACSON segmentation pipeline
We devised DeepACSON: a DCNN-based segmentation pipeline to automatically annotate white matter ultrastructures in the large field-of-view low-resolution 3D-EM datasets (Fig. 1). DeepACSON pipeline annotated white matter such that myelin and each myelinated axon, mitochondrion, and cell nucleus carries its label. Figure 1 shows the four steps of the DeepACSON pipeline.
Step 1: Training set generation for semantic-segmentation. DeepACSON requires labeled EM datasets as the training material to learn to perform the segmentation task. To automatically provide a training set for DeepACSON, we utilized our previously developed ACSON segmentation pipeline 4  The images were converted to the unsigned 8-bits per voxel format. The size of low-resolution SBEM datasets ranged from 10-15 GB, summing up to 85 GB (8-bits per voxel). To denoise SBEM datasets, we used a non-local filtering algorithm named block-matching and 4D filtering (BM4D) 39 . The algorithm was run in its low-complexity mode to reduce the computation time.
We processed the low-resolution SBEM datasets in non-overlapping patches with an approximate size of 1000 × 1000 × 300 to account for RAM requirements and enabling parallel BM4D denoising on CPU clusters (see Supplementary Fig. S2).
Step 1: training sets for DCNNs. We used high-resolution SBEM datasets of the corpus callosum segmented by the automated ACSON pipeline to generate a training set for DeepACSON. The automated ACSON pipeline was designed for segmenting high-resolution small field-of-view SBEM datasets into myelin, myelinated and unmyelinated axons, mitochondria, and cells.
The ACSON pipeline segments white matter based on a bounded volume growing algorithm in which seeds are defined automatically using regional maxima of the distance transform of myelin maps 4 . We generated two training sets: one for training DCNN-mAx and one for training DCNN-cN. The segmented high-resolution datasets were modified before using in the training sets. The labels of unmyelinated axons and cells were set equal to the background label to generate the DCNN-mAx training set. The labels of all ultrastructures except for the cell nuclei were set equal to the background label to generate the DCNN-cN training set. We created an additional label for the membrane of cell nuclei by a morphological gradient with a flat 11 × 11 square structuring element applied to the segmented cell nuclei. The membrane of cell nuclei was relatively thick: 22 voxels or 0.33 µm wide. We over-emphasized the membrane of cell nuclei because the number of voxels representing it was very small as compared to the nucleus and background labels. After modifying labels, we uniformly down-sampled the high-resolution datasets in both training sets by a window size of 3 × 3 in x − y plane to conform to the resolution of the low-resolution datasets. We generated nine down-sampled datasets out of each high-resolution dataset by shifting the sampling window for one voxel at a time in x and y directions. The DCNN-mAx training set included datasets of the ipsi-and contralateral of sham-operated and TBI rats. It is important to include both intact and injured axons to enrich the training set.
The DCNN-cN training set included only ten cell nuclei.
Step 2: semantic-segmentation of the white matter ultrastructure, architecture and implementation of DCNNs. We implemented the DCNNs and data augmentation using ElektroNN2 library 5 based on the Python Theano framework 40 .
ElektroNN was optimized for short model training times and fast inference on large datasets by eliminating the redundant computations of sliding window approaches 41 . The training set was augmented using randomized histogram distortions to become invariant against varying contrast and brightness gradients. In addition, half of the training set underwent image warping consisting of affine and random perspective transformations. We determined the optimal architectural parameters and hyperparameters of the DCNNs experimentally in the training set. We used the same architectural parameters for DCNN-mAx and DCNN-cN, as shown in Supplementary Fig. S3; receptive field was (45,45,15) voxels in (x, y, z) directions with 25 205 prediction neurons and 2 881 777 trainable parameters. The DCNNs were trained using the softmax function in tandem with the multinomial negative log-likelihood as a loss function. We set the batch size equal to 1 because there was a sufficient number of prediction neurons for updating the gradient. For the optimization, we used Adam optimizer 42 and set its initial learning rate α = 5 × 10 −5 , exponential decay rate for the first moment β 1 = 0.9, exponential decay rate for the second-moment β 2 = 0.999, and the weight decay = 5 × 10 −6 (we use the same notation as in 42 ). Also, after every 1 000 training steps, the learning rate α was decayed by a factor of 0.995. For DCNN-mAx, we set the class weights to 0.2, 0.5, 1, and 1 for the background, myelin, myelinated axons, and mitochondria, respectively. The class weights for DCNN-cN were set to 0.1, 1 and, 4 for the background, cell nuclei, and the membrane of cell nuclei, respectively. The weights were proportional to the inverse of the volume of ultrastructure in the training datasets.
Step 3.1: myelin, myelinated axons, and mitochondria segmentation. DCNN-mAx generated probability maps of myelin, myelinated axons, and mitochondria. We binarized the probability maps of myelin, myelinated axons, and mitochondria with a threshold of θ myelin = 0.5, θ mAxon = 0.8, and θ mitochondria = 0.5. The binarized volume of myelin required no further processing, and it was considered as the final myelin segmentation. The binarized volume of mitochondria was dilated with a flat cubic structuring element of size three voxels, and unified with the binarized volume of myelinated axons to form the intra-axonal space of myelinated axons, I IAS . Initial segmentation of individual myelinated axons was obtained using connected component analysis on I IAS . The segmented volumes underwent a morphological volume closing with a flat cubic structuring element of size three voxels to close the holes. We excluded volumes whose size was smaller than 4 712 voxels (the volume of a circular cylinder with a base radius of 0.25 µm and height of 3 µm, resolved with an isotropic 50 × 50 × 50 nm 3 voxel size).
In the low-resolution datasets, the axonal membrane was not resolved, e.g., at nodes of Ranvier. Thus, the initial segmentation of myelinated axons might exhibit under-segmentation error (Supplementary Fig. S6). We developed a specific shape analysis algorithm, CSD, to decompose an object into its semantic-components. As shown in Supplementary Fig. S5a, b, we defined the object as the union of several tubular components; similar to an axon with under-segmentation errors. The decomposition of such an object requires, first, to identify its semantic components, and second, to correct the ambiguity at intersections of branches, i.e., to decompose the shape at intersecting branches. The mathematical details of the algorithm and its validation are described in a different manuscript 12 , but we summarize the algorithm here: Step 3.1.1: skeleton partitioning and critical points. We propose to decompose an object into its semantic-components using its curve skeleton and cross-sectional analysis. We determined the curve skeleton of the object Ω by applying the method from Hassouna & Farag 43 . This method also returns the collection of skeleton branches Γ = {γ 1 , . . . , γ n } (see supplementary  Fig. S5e). We cut the object at critical points to obtain object-parts. We assigned the same label to the object-parts along the same sub-skeleton to acquire a semantic-component ( Supplementary Fig. S5f).
Step 3.1.2: object reconstruction. To reconstruct a semantic component at an intersection, we used generalized cylinders. We constructed a generalized cylinder Φ : R 2 → R 3 as a linear homotopy between two simple closed curves C c 1 and C c 2 elongated on a simple curve ζ ⊂ R 3 . To define C c 1 and C c 2 , let x c 1 , x c 2 ∈ ψ be two critical points in the proximity of the joint j on a sub-skeleton ψ. We cut the surface of the object ∂ Ω at these two critical points. The cross-sectional curves of ∂ Ω at x c 1 and x c 2 are C c 1 and C c 2 . We defined ζ to be the spline interpolation of ψ between x c 1 and x c 2 . The acquired generalized cylinder was used to reconstruct the object at the intersection which interior includes j (Supplementary Fig. S5f).
Step where i = 1, . . . , n and n ≥ n . We used the volumes in which the number of voxels was greater than 1 000 as initial surfaces for a secondary deformation. The evolution of a surface started at ∂U and stopped at the desired nucleus membrane because we defined the speed function equal to PR m . Each initial surface ∂U was deformed for 300 iterations, where the smoothing term of GDMs preserved the surface sphericity from taking irregularities at sites where membrane probability map PR m was discontinuous.
Step To generate a training set, we randomly picked five datasets out of ten segmented large field-of-view datasets. From the selected datasets, we manually selected representative myelinated axons and non-axonal instances from both the cingulum and corpus callosum (1 097 myelinated axons and 189 non-axonal instances; less than 0.2% of the total number of components segmented as myelinated axons in all ten datasets). To determine the optimal SVM hyperparameters, the regularization parameter C and kernel parameter σ , we selected the pair that minimized the error sum of 5-fold cross-validation on the training set. Using Bayesian optimization algorithm 47,48 to search for the optimal C and σ pair constrained in [10 −6 , 10 6 ], we set C = 1.12 and σ = 8.48. The trained SVM was applied to eliminate non-axonal structures from all the datasets. We applied the same approach to eliminate false-positives from cell nuclei. For each component segmented as a cell nucleus, we calculated the same features (as mentioned for myelinated axons) before and after the geometric deformation. To generate a training set, out of ten segmented large field-of-view datasets, we picked four random datasets.
From the selected datasets, we manually selected representative cell nuclei and non-nucleus instances (263 cell nuclei and 113 non-nucleus instances; 7.3 % of the total number of components segmented as cell nuclei in all ten datasets) to train an SVM with a quadratic kernel. We set the optimal hyperparameters of the SVM as C = 6.51 and σ = 4.17 using aforementioned methods. Figure S1. DeepACSON evaluation. We developed a GUI-based software tool, gACSON, in Matlab to load and visualize the segmentation for proofreading. gACSON is designed for the visualization of the large-scale image-datasets/segmentation. Using gACSON, an expert evaluated the DeepACSON segmentation of myelinated axons and mitochondria at the object-level. We randomly sampled each low-resolution dataset and its corresponding segmentation by non-overlapping images of size 300 × 300 voxels. The sampled images were quantified for the number of true-positives (TP), false-positives (FP), and false-negatives (FN) to calculate the precision, recall, and F1 score. The expert had no access to the dataset ID nor the sampling location. Figure S2. BM4D filtering of the low-resolution SBEM datasets. BM4D recognized the noise as Gaussian-distributed with the standard deviation in the range [17,22] in our low-resolution SBEM datasets. Shown images were acquired from the cingulum and corpus callosum of sham #25 low-resolution dataset. Figure S3. The architecture of DCNNs used in DeepACSON. We used the same architecture for DCNN-mAx and DCNN-cN. The size of the convolutional kernels is denoted as Conv(x, y, z). The number of channels/feature maps created from a layer of convolutional kernels is denoted by @n. The size of the max pooling operation is denoted as Max pool(x, y, z).  (e) On a sub-skeleton ψ and in the proximity of a junction-point j ∈ ψ, we defined two decomposition intervals. The boundaries of decomposition intervals are shown with red filled-circles. In each interval, the cross-section of the object was swept along ψ and towards the joint j to find a critical point. At a critical point, the normalized Hausdorff distance H ρ between a cross-sectional contour and the mean of visited cross-sectional contours exceeds θ H . Sweeping directions are shown with arrows. (f) We cut the object at critical points to obtain object-parts. The object-parts along the same sub-skeleton were assigned the same label to construct a semantic-component. The semantic-components were further reconstructed between their comprising object-parts using generalized cylinders, magnified in (f 1 -f 4 ). The synthetic object in (a) comprised of seven object-parts, and our algorithm decomposed it into three semantic-components.

Supplementary Information
v Figure S6. Decomposition of under-segmented myelinated axons into their semantic axonal components using the CSD algorithm.
vi  x