Volume electron microscopy in injured rat brain validates white matter microstructure metrics from diffusion MRI

Abstract Biophysical modeling of diffusion MRI (dMRI) offers the exciting potential of bridging the gap between the macroscopic MRI resolution and microscopic cellular features, effectively turning the MRI scanner into a noninvasive in vivo microscope. In brain white matter, the Standard Model (SM) interprets the dMRI signal in terms of axon dispersion, intra- and extra-axonal water fractions, and diffusivities. However, for SM to be fully applicable and correctly interpreted, it needs to be carefully evaluated using histology. Here, we perform a comprehensive histological validation of the SM parameters, by characterizing white matter (WM) microstructure in sham and injured rat brains using volume electron microscopy and ex vivo dMRI. Sensitivity is evaluated by how well each SM metric correlates with its histological counterpart, and specificity by the lack of correlation with other, non-corresponding histological features. Compared to previously developed SM estimators with constraints, our results show that SMI is the most sensitive and specific. Furthermore, we derive the functional form of the fiber orientation distribution based on its exponentially decreasing rotational invariants. This comprehensive comparison with histology may facilitate the clinical adoption of in vivo dMRI-derived SM parameters as biomarkers for neurological disorders.


INTRODUCTION
Diffusion MRI (dMRI) has the remarkable ability [1][2][3] to probe the Brownian motion of water molecules on length scales comparable to cellular structures (∼ 10 µm), two orders of magnitude smaller than the standard clinical MRI resolution (∼ 1 mm).In this way, dMRI carries information about the cellular-level architecture restricting the diffusion of water molecules, and offers the exciting potential to characterize tissue microstructure in vivo and non-invasively in clinical settings 4 .However, dMRI signals are indirect measures of tissue structural properties, and the task of identifying and extracting the relevant information remains an active field of research [2][3][4][5][6] .
Biophysical modeling of the dMRI signal strives to provide metrics that are not only sensitive but also specific to the underlying tissue microstructure 4,5 .Biophysical models rely on assumptions (they "draw a picture") about the geometry of the tissue that restricts the diffusion of the water.If these assumptions accurately capture the relevant properties that are encoded in the dMRI signal, such biophysical models can bridge the gap between the millimeter macroscopic MRI voxel resolution and the micrometer-level tissue architecture.
In brain white matter (WM), at sufficiently long diffusion times, the dMRI signal of an elementary WM fiber segment, or fascicle, can be modeled by two nonexchanging anisotropic Gaussian compartments 2 : an intra-axonal compartment in which axons are repre- * ricardo.coronadoleija@nyulangone.org, rleija@cimat.mxsented by impermeable cylinders with zero radius 7 ; and an extra-axonal compartment represented by an axially symmetric tensor aligned with the fascicle.Given sufficient information in the dMRI measurements, a third, isotropic compartment could be considered.A macroscopic MRI voxel usually contains multiple orientationally dispersed fiber fascicles that contribute to the directional dMRI signal according to their orientation distribution.This general picture, illustrated in Fig. 1d and mathematically defined in Methods section, is known as the Standard Model 2 (SM) of diffusion in WM, as it unifies many previously proposed biophysical models relying on similar assumptions [8][9][10][11][12][13][14][15][16][17][18][19] .
In this work, we perform a comprehensive histological validation of the SM parameters, by characterizing brain WM microstructure on sham-operated and traumatic brain injury (TBI) rats, using volume (3d) electron microscopy (EM) and ex vivo dMRI.The 3d nature of our histological analysis, the variety of the EM samples (corpus callosum, cingulum, in the ipsi-and contralateral sides, for TBI and sham-operated animals), and the large number of segmented 3d axons (in the order of ten thousand per sample) -all together enable an unprecedented characterization of tissue microstructure.For the dMRI analysis, we estimate the SM parameters using four publicly available approaches: white matter tract integrity (WMTI) 10 , neurite orientation dispersion and density imaging (NODDI) 11 , spherical mean technique (SMT) 13 , and standard model imaging (SMI) 14,16,19 .We refer to these as SM estimators, as they target the same biophysical SM parameters, yet differ in the additional constraints (see Table I) and ways of performing model fitting (max- e .One voxel contains multiple orientationally dispersed fascicles that contribute to the directional dMRI signal according to their FOD P(n).(e) EM images transverse to the axons of the cingulum for the rat brains in (a), with segmented IAS, EAS, and myelin.Representative segmented 3d axons are also shown.Using skeletons from thousands of segmented axon, FODs for each axon population in the volumes were constructed.Reduction in axon density, changes in axon morphology, and increase in axon dispersion can be observed for the TBI animal with respect to the shamoperated animal.(f) dMRI-derived SM parametric maps reveal specific differences between sham and TBI that are observed in corresponding 3d EM histology.likelihood versus machine learning).While all these estimators provide SM parameters that correlate significantly (p < 0.05) with their histological counterparts, indicating sensitivity, we find that SMI shows the highest specificity by presenting smallest cross-correlations with other, non-corresponding histological features.This comparison of SM parameters against 3d EM validates SM metrics as specific tissue biomarkers valuable for clinical applications and neuroscience research, as well as reveals the strengths and limitations of the chosen estimator.
The parameters provided by SM describe specific cellular features and could be a powerful tool for the study of pathological conditions.The intra-axonal water fraction f is a potential marker for axon density that could quan-tify axonal loss [21][22][23][24] , the intra-axonal diffusivity D a is a potential marker for axon morphology that could quantify axonal injury such as beading [25][26][27] , the extra-axonal axial diffusivity D ∥ e is a potential marker for changes in the extra-axonal space such as inflammation, and the extra-axonal radial diffusivity D ⊥ e is a potential marker for demyelination 23,28,29 .
Aside from the SM kernel parameters related to properties of fiber fascicles, SM includes the fiber orientation distribution P(n) (FOD), which can be described using its spherical harmonics decomposition without assuming any particular functional form.P(n) properties can be described using FOD rotational invariants p l , which are the energies of the FOD spherical harmonics coefficients p lm in each degree-l sector (Eq. 5 in Methods section).
Here, the second-order rotational invariant p 2 is of special importance as it is inversely related 16 to the average angle of dispersion of the axons θ via Eq. 6. P(n) and θ p2 are informative for structural brain connectivity, surgical planning, as well as to assess pruning in neurodevelopment 12 .

METHODS
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
Five adult male Sprague-Dawley rats (10 weeks old, weights between 320 and 380 g, Harlan Netherlands B.V., Horst, Netherlands) were used in the study.The animals were housed in a room (22 ± 1 • C, 50-60% humidity) with 12 h light/dark cycle and free access to food and water.
TBI was induced by lateral fluid percussion injury in three rats 49 .Rats were anesthetized with a single intraperitoneal injection, and a craniotomy (5 mm 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 by a transient fluid pulse impact (21-23 ms) against the exposed intact dura using a fluid-percussion device.The impact pressure was adjusted to 3.2-3.4atm to induce a severe injury.The sham operation on the other two rats included all the surgical procedures except the impact.See 50 and 51 for detailed information about the animal model.
Five months after TBI or sham operation, the rats were transcardially perfused using 0.9% NaCl for 2 min followed by 4% paraformaldehyde (PFA).The brains were removed from the skull, post-fixed in 4% PFA / 1% glutaraldehyde overnight at 4 o C, and then placed in 0.9% NaCl for at least 12 h to remove excess PFA.

Ex vivo dMRI imaging
The extracted rat brains were scanned ex vivo at room temperature (21 o C) in a vertical 9.4 T/89 mm magnet (Oxford Instruments PLC, Abingdon, UK) interfaced with a DirectDrive console (Varian Inc., Palo Alto, CA, USA) using a quadrature volume RF-coil (Ø = 20 mm; Rapid Biomedical GmbH, Rimpar, Germany) as both transmitter and receiver.During imaging, the brains were immersed in perfluoropolyether (Solexis Galden ®, Solvay, Houston, TX, USA) to avoid signals from the surrounding area.

3d EM imaging
After the ex vivo dMRI acquisitions, brains were placed in 0.9% NaCl for at least 4 h to remove excess perfluoropolyether and then sectioned into 1 mm thick coronal sections with a vibrating blade microtome.For each brain, sections at 3.80 from bregma were selected and further dissected into smaller samples containing the areas of interest (Figure 1a and Fig. 1b).Ten WM samples, two from each brain were collected: the ipsi-and contralateral of the cingulum and corpus callosum.The samples were stained using an enhanced protocol with heavy metals, dehydrated and embedded in Durcupan ACM resin.Before mounting the specimens, the excess resin in the hardened tissue blocks was trimmed.After selecting the area of interest for imaging 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.The tissue was exposed on all four sides, bottom, and top of the pyramid.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 blocks were imaged in a scanning electron 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 face of the blocks was in the x-y plane, and the cutting was in z direction.All blocks were imaged using a beam voltage of 2.5 kV and a pressure of 0.15 Torr.After imaging, Microscopy Image Browser (MIB; http://mib.helsinki.fi) was used to process and align the EM image stacks.
Two datasets were acquired simultaneously at low-and high-resolution consistently at one specific location in the white matter in both the ipsi-and contralateral hemispheres for all sham-operated and TBI animals as shown in Fig. 1b.The low-resolution datasets were imaged from large tissue volumes of 200×100×60 µm 3 with resolution 50×50×50 nm 3 .Two-thirds of the acquired volumes correspond to the corpus callosum and one-third to the cingulum.Only the corpus callosum region in the ipsilateral side for one rat is missing.The high-resolution datasets were imaged from small tissue volumes of 15 × 15 × 15 µm 3 on the corpus callosum with resolution 15 × 15 × 50 nm 3 .For detailed information about 3d EM tissue preparation, acquisition, and pre-processing, see 50 and 51 .
The labor-intensive sample preparation for EM studies limits the number of samples that can be collected, especially for the large 3d EM samples used in this work 50 .Yet, thanks to the features of these large 3d EM samples (comparable to the MRI voxel size) and the automated axon segmentation 51,52 , this study was able to include a large number of segmented 3d axons per population (in the order of ten thousand) in a diverse set of ROIs, for sham and TBI conditions, allowing for unprecedented microstructure characterization in comparison with previous histological studies.

3d EM Segmentation
Myelin, myelinated axons, and cell nuclei were automatically segmented from the low-resolution datasets using the deep-learning-based pipeline DeepACSON 52,53 (https://github.com/aAbdz/DeepACSON),which is a method suited for the challenging task of segmenting large datasets with limited visualization of the cell membranes caused by low resolution.
First, the high-resolution datasets (small tissue volumes) were segmented using the ACSON pipeline 51 , which is an automated method that integrates edge detection and seeded region growing algorithms in 3d, refining the segmentation with the SLIC (Simple Linear Iterative Clustering) supervoxel technique.These highresolution images and their corresponding segmentations were down-sampled to match the low-resolution images (large tissue volumes) and used as human-annotation-free training sets for the DeepACSON pipeline 52 .DeepAC-SON uses two convolutional neural networks to compute probability maps of myelin, myelinated axons, mitochondria, cell nuclei, and cell nuclei membranes.Then, these probability maps are binarized, and the segmentations are refined.Because of the low resolution of the images, the myelinated axons are prone to under-segmentation; for this reason, the segmented axons are rectified using a cylindrical shape decomposition method 53 , and finally a support vector machine is used to exclude non-axonal segmented objects.
Further proofreading was performed by considering only the myelin attached to the segmented axons, thus removing spurious segmented objects in the extra-axonal space that were labeled as myelin.For the computation of volume fractions all segmented axons were considered, but for the quantifications of properties such as: axon diameter, cross-sectional area variation, and the fiber orientation distribution; only axons larger than 10 µm, were used.In order to avoid biases due to oblique cross sections at the extremes caused by the inclination of the axons, 1 µm was removed at the two ends of each axon.

Microstructural metrics derived from 3d EM
Using the segmented 3d EM samples, the intra-axonal, extra-axonal, and myelin volumes were computed as the number of voxels covering each compartment.The dMRI related intra-axonal volume fraction metric f was computed as the intra-axonal volume divided by the total volume of the sample, excluding the myelin compartment due to its short T 2 compared with the TE of the dMRI acquisition 54 , which makes it dMRI invisible.
Morphological properties of the segmented axons were obtained following the approach in 55 (https://github.com/NYU-DiffusionMRI/RaW-seg).To compute axon diameters and cross-sectional area variations, each axon's main direction was first aligned to the z-axis.Then, the axon skeleton was created as the line connecting the center of mass of each slice.Assuming a circular cylinder axon, the cross-sectional area A(z) and the diameter 2r ≡ 2 A(z)/π of a circle with the same area were computed for each slice perpendicular to the axon skeleton.
Recently, the exact axial tortuosity (Λ ∥ ) relation between the varying cross-sectional areas A(z) along the main axis of the axon and the intra-axonal diffusivity was derived 56 from the Fick-Jacob's equation: where A is the mean cross-sectional area, and D 0 is the diffusivity for a perfect cylinder axon.
Dispersion angle for an axon population is defined 16 as: where the individual axon segment's dispersion angle θ i is the angle between the axon segment's direction with respect to the main direction of the axon population.Dispersion angle θ u caused only by undulations 27 can be computed using Eq. 2 on the axons aligned to the z-axis.Assuming free diffusivity for perfect cylinders (D 0 = D w = 2.0 µm 2 /ms for dMRI acquired at room temperature), the value of SM intra-axonal diffusivity D a can be predicted from histology using Eq. 1, and θ u computed with Eq. 2 on each axon, aggregating the contribution for all of them: where the weights w k are proportional to the axon volumes, with k w k ≡ 1. Fiber orientation distributions on each WM region for each animal were computed using the corresponding skeletons of the aligned axons by rotating them back to the original main directions and projecting their tangent vectors on a 3d triangulated spherical surface 55 ; creating in this way a histogram P(n) that was normalized such that dnP(n) ≡ 1.These FODs were decomposed on spherical harmonic (SH) basis: where Y lm are the SH basis, and p lm are the SH coefficients.Because of FOD antipodal symmetry, only even harmonic orders l are considered up to a maximum order of l max = 16.Also, for a normalized FOD p 0 = 1.To simulate coarse-grained effects caused by diffusion in the calculation of the FOD, the axon skeletons were smoothed by a Gaussian kernel with variance 55 , where L = √ 2Dt is the diffusion length, D = 2µm 2 /ms, and t is the diffusion time of the dMRI acquisition.
From the SH decomposition of P(n) in Eq. 4, FODs rotational invariants p l were be computed 14,16 : where p 2 is of special interest as it is related with the average dispersion angle of the axons 16,55 in Eq. 2 by Some models represent the FOD P(n) as a Watson distribution 11,15 , which is defined by its concentration parameter κ.A relation between κ and p 2 is given in 15 : It has been suggested 14 that the rotational invariants p l of the FODs follow an exponential decay with respect to the harmonic order l, given by the Poisson kernel p l = p 0 •λ l with p 0 = 1 and λ ≤ 1.In 55 , by using segmented 3d axons in the normal mouse corpus callosum, it was found that this functional form was not properly normalized and an extra term was required where C ≥ 1, for l = 2, 4, ....In our work, given the large number of segmented 3d axons, and the variety of samples with and without pathology, this hypothesis was further validated.

Standard Model
The Standard Model 2 of the dMRI signal in the measured direction û for brain WM (Figure 1d), can be written as: This equation describes the dMRI signal S(b, ξ, û) as a convolution between the FOD P(n) and the SM kernel K(b, ξ, û • n), which represents the fiber segment, and is the sum of two, non-exchanging, intra-and extra-axonal Gaussian compartments: , where b is the diffusion weighting, û • n is the cosine of the angle between the symmetry axes of the SM kernel n and the measured direction û, and ξ are the parameters of the SM kernel: In SH basis, the convolution between the FOD and the SM kernel in Eq. 8 is a product Orientation information can be factored out by computing the signal, FOD and SM kernel rotational invariants 14,16 : with for the diffusion MRI signal, and Although all these approaches aim to estimate essentially the same SM parameters, the main differences between them are the number of independent parameters (second column), the constraints/interdependencies imposed on parameters (third column), the functional form of the FOD (fourth column), as well as different strategies for estimation (WMTI: linear diffusion and kurtosis tensors fit followed by explicit analytical formulas based on kurtosis tensor parameters; NODDI and SMT: max-likelihood; SMI: machine learning approach based on the polynomial regression from signal's rotational invariants to SM paramteres).All these differences lead to different values in parametric maps, Fig. 2, and correlations with histological counterparts, Figs.3-4.
for the SM kernel, where P l are the Legendre polynomials.
The estimation of the SM parameters from conventional Pulsed-Gradient Spin-Echo (PGSE) dMRI is challenging because the problem is ill-posed 16,57 .It has been shown that "advanced" dMRI acquisitions are needed for an accurate, precise, and robust estimation of all the SM metrics 16,18,19 .Nonetheless, there is still a need to correctly estimate the SM parameters on the typical PGSE data used in this study, as the results can be extrapolated to the majority of clinical human studies that use similar dMRI protocols including the large multicenter dMRI datasets currently being collected [30][31][32][33] .Several strategies have been proposed to solve the degeneracy of the SM parameter estimation in PGSE data.The principal characteristics of these strategies, regarding to estimated parameters and constraints, are summarized in Table I.
White matter tract integrity 10 uses the diffusion and kurtosis tensors to derive formulas for the SM parameters, however, its mostly valid for aligned axon populations.Neurite orientation dispersion and density imaging (NODDI) 11 is the most widely used estimator for f and the FOD (as Watson distribution), however, it fixes the axial diffusivities, and constrains the extraaxonal radial diffusivity using the first order tortuosity approximation 58 .Spherical mean technique (SMT) 13 uses the spherical mean (0 th order rotational invariant) on each shell of the dMRI signal to factor out the FOD and focus only in the estimation of the SM kernel parameters, however, it uses two of NODDIs strong assumptions.Standard model imaging (SMI) 14,19 uses a machine learning approach to solve the maximum-likelihood estimation of the system 16 in Eq. 11, here, the rotational invariants of the dMRI signal (Eq.12) are mapped to the SM kernel parameters (Eq.9) via cubic polynomial regression, however, when the information in the data is insufficient, the estimation of the parameters could be biased by the training data.
SM parameters also include the FOD P(n), which can be described without assuming any functional form.Once the SM kernel parameters ξ are obtained, the FODs SH coefficients p lm can be estimated from Eq. 10.Also, the FODs SH rotational invariants p l can be directly estimated from Eq. 11.

dMRI Processing
In order to improve the precision and accuracy of the estimated SM maps, dMRI data was processed following the steps in the DESIGNER pipeline (https://github.com/NYU-DiffusionMRI/DESIGNER) 59 adapted for ex vivo dMRI rat brains.Magnitude and phase volumes were reconstructed from the data to perform MP-PCA complex-denoising 60,61 , which has shown to highly reduce the Rician noise floor.Images were also corrected for Gibbs-ringing 62 , and motion distortions 63 .WMTI maps were computed using tools provided by the DE-SIGNER pipeline.NODDI, SMT and SMI maps were estimated with their respective toolboxes: SMT (https: //github.com/ekaden/smt),NODDI (http://mig.cs.ucl.ac.uk/index.php?n=Tutorial.NODDImatlab), and SMI (https://github.com/NYU-DiffusionMRI/SMI).
For NODDI, the axial diffusivity was fixed to the suggested value for ex vivo, D ∥ = D a = D ∥ e = 0.6 µm 2 /ms.(inversely related to dispersion) from SMI, NODDI, WMTI, and SMT for a sham-operated and a TBI animal.Changes in TBI with respect to sham, such as the reduction in f and p2, agree with histology for all methods.However, evident quantitative differences between their maps are observed.
For SMT and SMI, the upper bound for the diffusivities was set to the free diffusion of water at room temperature ∼ 2 µm 2 /ms.For SMI, rotational invariants up to l = 6 were used for parameter mapping.SM free water fraction f w metric can be estimated by SMI and NODDI, however, imaging was performed using b = 2, 3, 4 ms/µm 2 at room temperature (21 o C), for which free water signal is negligible (≈ 1.5% of S 0 for the lower b-vale).For this reason, for SMI we set f w = 0.However, for NODDI f w was estimated as it was shown to improve the correlations of its other parameters.Manually delineated ROIs were drawn by identifying, in the diffusion tensor ellipsoids, the corpus callosum, cingulum and their crossing region in the ipsi-and contralateral side at the location where the 3d EM samples were collected (Fig. S5).
FODs second-order rotational invariant p 2 is sensitive to axon dispersion, as shown in Fig. 5a, and its affected by complex fiber configuration such as crossing fibers in the voxel.In the dMRI volumes, we observed partial volume in the cingulum ROI which contained a small fraction of corpus callosum axons in the perpendicular direction.For this reason, the SM maps in Fig. 1f and Fig. 2 show increased dispersion (higher θ p2 and lower p 2 ) in the cingulum voxels with respect to the corpus callosum voxels, contrary to what it is observed in the 3d EM volumes.SMI estimates the full non-parametric FOD, from which we segmented the lobes corresponding to the different axon populations 64 , then from the segmented lobes we computed per-bundle rotational invariants p l .

Statistics
Pearson correlations (ρ), p-values (p), and 95% confidence intervals (CI) between each dMRI and 3d EM derived metric pairs were computed using MATLABs command corrcoef, which uses a t-test with N −2 degrees of freedom for p-values (N is the number of samples), and for the confidence intervals uses an asymptotic Gaussian distribution of 1 2 ln 1+ρ 1−ρ , with an approximate variance equal to 1/(N − 3).Here, the total number of samples used in the analysis is 28 instead of 30, as the corpus callosum region in the 3d EM is missing, hence no cc and no cc+cg, for one sham-operated animal.We corrected for multiple comparisons using the False Discovery Rate (FDR) computed with MATLABs command mafdr that uses the Benjamini-Hochberg procedure.Concordance correlation coefficients ρ c were also computed between SM dispersion angle θ p2 and its histological counterpart θ, as agreement in their values is expected for these metrics in particular.It is important to mention that in this study a correlation does not necessarily imply causality, as a pathological condition such as TBI is very complex, with many microstructural changes happening simultaneously, and it is challenging to get the full picture of how a particular SM parameter is affected by them.Also, statistical significance, in our case p < 0.05, could be affected by the number of samples and the assumptions made by the statistical test, such as Gaussianity.

RESULTS
The TBI rat brains show visible damage due to the impact on 1 mm-thick brain sections (Fig. 1a) and corresponding fractional anisotropy (FA) maps (Fig. 1c), as opposed to sham-operated rat brains showing no visible damage.To reveal the corresponding changes caused by TBI at the microstructural scale, we automatically segmented the myelin and myelinated axons on 3d EM, as illustrated in Fig. 1b.Fig. 1e shows the reduction of axonal density, changes in axon morphology (diameter 2r, axial tortuosity Λ ∥ and undulations, as defined in Meth- ods section), as well as an increase in axonal dispersion in the TBI sample.Remarkably, these microstructural changes are also detected non-invasively by ex vivo MRI, as shown on the SM maps in Fig. 1f for the axon water fraction f , intra-axonal diffusivity D a , and axon dispersion angle θ p2 .
Parametric maps obtained from the same ex vivo dMRI measurement but using the different estimators (SMI, WMTI, NODDI, SMT) are compared in Fig. 2 for f and p 2 (its relation to dispersion angle via Eq.6 is demonstrated in Fig. 5a).Consistent changes are observed between sham-operated and TBI animals that are in qualitative agreement with the histology shown in Fig. 1e.Yet, the maps also demonstrate quantitative differences among the distinct estimators.
Sensitivity and specificity are evaluated by comparing SM parameters f , θ p2 and D a with their 3d EM counterparts on six ROIs: cc, cg, cc+cg for the ipsi-and contralateral hemispheres for each animal, as shown in Fig. S5.While histological f and θ are computed directly from the segmented axons, Da is predicted from histology (Eq. 3) by considering the axons axial tortuosity 56 Λ ∥ and the dispersion angle caused only by axon undulations 27 θ u , assuming free diffusivity D w = 2 µm 2 /ms (at room temperature) for perfect cylinder axons.The influence of axon diameter 2r in the SM parameters is also evaluated.
The diagonals of the dashed submatrices in Fig. 3 indicate that all SM estimators provide metrics that correlate significantly (p < 0.05) with their corresponding tissue properties derived from 3d EM histology, except for f from NODDI.This suggests good overall sensitivity of all SM estimators.Here, SMI achieves the highest correlations for the three SM metrics, with only SMT providing a similar correlation for f .
Regarding specificity, the off-diagonal elements of the dashed submatrices in Fig. 3 show that only SMI has no significant correlations with non-corresponding histolog-ical features, while the other estimators reveal multiple spurious cross-correlations.Looking at each SM parameter individually (columns), dMRI f is specific to its histological counterpart for SMI, WMTI, and SMT, while for NODDI, it correlates with 3d EM θ. dMRI θ p2 is specific to histological θ for SMI, SMT and NODDI, while for WMTI, it correlates also with 3d EM f .dMRI D a is specific to its predicted value from histology Da only for SMI, while for SMT, it correlates also with 3d EM f , and for WMTI, it shows correlations also with 3d EM f and θ.The highest specificity of SMI parameters can be understood by the lack of hard constraints and assumptions as opposed to the other estimators (see Table I).
In the scatter plots for dMRI f against EM f of Fig. 4a, agreement with the identity line is not expected.For one part, the resolution of the EM samples does not allow to segment unmyelinated axons.Furthermore, dMRI data was acquired using a single TE, so SM f is T 2 -weighted, while EM f is not.For these reasons, and assuming intraaxonal T 2 longer than extra-axonal T 2 , 17,65 an overestimation of dMRI f is expected.
On the other hand, agreement with the identity line is expected when comparing dMRI θ p2 against EM θ in the scatter plots of Fig. 4, as unmyelinated axons are expected to have similar dispersion than myelinated axons, and the fiber orientation distribution is expected to be unaffected by relaxation.In Fig. 4a, better agreement with the identity line is observed for θ p2 from SMI than from the other methods, this is also indicated by its highest Pearson correlation (ρ) values shown in Fig. 3 and its highest concordance correlation coefficients ρ c , with respect to the other methods: SMI ρ c = 0.850, WMTI ρ c = 0.332, SMT ρ c = 0.191, and NODDI ρ c = 0.569.
In the scatter plots of Fig. 4a that compare dMRI D a against its predicted value from histology Da (Eq.3), we observe that D a from SMI agrees better with the identity line.For the other methods, the estimation of D a could  .NODDI keeps D a fixed and therefore the scatter plots for this metric are not shown.Correlation between D a and axon diameter is observed for SMI, but not for the other methods, however, it presents the lower value among the significant correlations of SMI.

Figure 4b scatter plots reveal positive correlations for SMI and WMTI D ⊥
e with histological f , while negative correlations are observed for SMT and NODDI.The negative correlations follow from the extra-axonal tortuosity approximation used by NODDI and SMT 11,13,58 , that relates D ⊥ e with f and D ∥ e as shown in Table I, and as also observed in Monte Carlo simulation studies that mimic axonal loss 28,58 .SMI and WMTI, on the other hand, don't impose this constraint, and show a somewhat surprising opposite sign of correlation.
Using the segmented 3d axons, we validate in Fig. 5a that the dispersion angle θ p2 , computed from the histological FOD p 2 (Eq.6), is a good approximation for the true angle of dispersion θ of the axons.While dMRI p 2 , and hence θ p2 , can be straightforwardly obtained from all SM estimators (Eq.11), SMI also estimates nonparametric FODs in SH basis.The rotational invariants p l computed from these dMRI FODs (Eq.5) correlate with their corresponding histological p l , as shown in Fig. 5b.These higher order p l allow to further validate their exponential decay 14,55 with respect to the increasing order l, for l ≥ 2 in Eq. 7, for both 3d EM and dMRI, as shown in Fig. 5c for the ipsilateral cingulum of all animals.Furthermore, by comparing the exponential decay parameters C and λ, from SMI and histology, we observe in Fig. 5d that C shows better agreement.It should be noted that dMRI parameters C and λ are directly computed from the rotational invariants p 2 , p 4 and p 6 estimated with SMI, without imposing the exponential decay functional form as constraint during the fitting.

DISCUSSION
Biophysical modeling bridges the resolution gap between the macroscopic dMRI voxel and microscopic tissue features, and provides a powerful tool to study development, aging, and diseases at the cellular level where these processes originate.Numerous human in vivo studies in brain WM 12,26,[34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]66 have been conducted using SM parameters derived from different estimators: WMTI, NODDI, SMT and SMI. Hoever, the interpretation of these and future studies significantly relies on the biological accuracy and promised microstructural specificity of the metrics provided by these estimators 4,5 , which highly depends on how correctly they describe the relevant features of the cellular environment, and how accurately and precisely they are estimated.
Previous validation efforts so far have been limited to 2d histology or incomplete.Indeed, the intra-axonal water fraction f has been evaluated in rodents by comparing against axon density or axon volume fraction metrics derived from 2d optical and electron microscopy (EM) for normal tissue 21 , and models for traumatic brain injury (TBI) 67 , demyelination 23 , and axonal degeneration 24 ; f also has been evaluated in postmortem human tissue for multiple sclerosis 68 .However, the use of 2d microscopic images limited the scope of these studies, since 3d tissue properties, such as the FOD and axon morphology metrics (that reveal injuries along the axons), could not be quantified, and specificity of f was not evaluated.Other studies have compared the FOD and/or dispersion parameters against related histological metrics derived, not from segmented axons, but from structure tensor analysis of 2d optical microscopy [68][69][70] , 3d confocal 71 and 3d electron microscopy 72 .However, the previous 3d studies focused on healthy tissue and did not evaluate any other metric besides dispersion, preventing the validation of the specificity of SM metrics.Furthermore, previous histological works did not assess the differences between the available SM estimators, such as WMTI, NODDI, SMT, and SMI, and focused the analysis on only one of them, biasing the results to the particular estimator used in the study.
In this work, we characterized WM microstructure of sham-operated and TBI rat brains using 3d EM and ex vivo dMRI to validate the sensitivity and specificity of the SM parameters.The large number of segmented 3d axons per sample (in the order of the ten thousand), and the variety and large size of these samples (from corpus callosum and cingulum, in the ipsi-and contralateral hemispheres of the brain, for sham-operated and TBI rats), allowed an unprecedented histological microstructure characterization (Figure 1e) revealing the following WM microstructure changes due to TBI: reduction in axon density, injuries along the axon (observed by the reduction in axon diameter, and by the increase in beading and undulation), and increase in axon dispersion.The simultaneous changes in several microstructural properties indicate the complexity of the damage in chronic TBI, and the need for sensitive and specific dMRI biomarkers in order to properly disentangle between them.
The direct comparison of 3d EM derived f vs dMRI intra-axonal water fraction f (Figure 3) reveals the strongest sensitivity for f estimated with SMT and SMI, followed by WMTI, as well as strong specificity for these estimators, as there are no correlations between f and θ nor D a .On the other hand, f from NODDI correlates with θ but not with histological f , which may suggest that imposing constraints biases NODDI parameter estimation, especially in pathological scenarios 73 .Nonetheless, the otherwise overall agreement between SM estimators in terms of sensitivity and specificity is in line with the consistent findings in many human in vivo studies that report changes in f with development, aging and disease.Indeed, f has been shown to increase in devel-opment and decrease in aging using NODDI 12,[34][35][36][37]39 , WMTI 12,36,38,39 , SMT 36,39 and SMI 36 .In neurological disorders, f is observed to increase in stroke ischemic lesions using NODDI 40 , WMTI 26 , and SMI 36 , while decrease in TBI 47,48 using NODDI, Alzheimer's disease using WMTI 41,42 and NODDI 43 , and in Parkinson's disease 66 using NODDI. Furthermore f has been suggesting as surrogate marker to track disease as it correlates with disability scores in multiple sclerosis, using WMTI 44 , NODDI, SMT 46 , and SMI 36 . Theconsistent correspondence between f estimated by WMTI, SMT, and SMI, with 3d EM f shows that changes in this metric can be interpreted as mainly driven by changes in axon volume fractions.
For dispersion, Fig. 3 shows that the dMRI dispersion angle θ p2 is sensitive and specific to the 3d EM true angle of dispersion of the axons θ for SMI, SMT, and NODDI, while for WMTI it shows to be sensitive but not specific, correlating also with histology f .This histological validation may help interpret human in vivo studies that showed changes in dMRI derived axon dispersion, namely reduced dispersion in Alzheimer's disease 43 , and chronic TBI 48 , while increased dispersion in stroke 36,40 , multiple sclerosis 36,45 , and Parkinson's disease, 66 .The quantitative differences obtained by distinct SM estimators, as shown in Fig. 2 and Fig. 4a, are also observed in human in vivo studies using different estimators.In particular, dispersion has shown to decrease during neurodevelopment using WMTI 12 and SMI 36 in the genu corpus callosum, which has been associated with axonal pruning, while weak or no significant changes have been observed using NODDI 12,[34][35][36] .Our histological validation helps understand these discrepancies and highlight the importance of SM estimators and their underlying constraints.Interestingly, SMI provides a good one-to-one correspondence with histological dispersion angle θ p2 , which confirms that the lack of hard constraints (Table I) in SM estimation results in the best overall accuracy for SMI 36 .
Using microscopy data to validate compartmental specific diffusivities D a , D ∥ e , and D ⊥ e is challenging, as there is no straightforward way to directly match them with histological metrics 4 .However, insight can be gained from studying changes in these diffusivities under certain pathological conditions 23,[25][26][27] , or by performing Monte Carlo simulations 25,27,56 in realistic substrates mimicking brain tissue in different (pathological) states.Here, we compared SM D a against the expected diffusivity from histology Da (Figure 4a and Eq. 3) based on the axon tortuosity Λ ∥ 56 and the dispersion from axon undulation θ u 27 .Correlations between D a and Da are observed for SMI, SMT and WMTI.However, specificity is observed only for D a from SMI, as D a from SMT and WMTI both correlates with other non-corresponding histological features.Our result underscore the previously observed dependency of D a on axon diameter variation and axon beading, obtained using Monte Carlo simulations 25,27,56 and human in vivo studies 26,27,36 .
For the extra-axonal perpendicular diffusivity D ⊥ e , Fig. 4b shows different dependencies on EM-derived volume fraction between the different SM estimators.The positive correlation between histological f and SMI/WMTI D ⊥ e , contradicts the extra-axonal tortuosity approximation used as a hard constraint in NODDI and SMT.While intuitively, D ⊥ e is expected to decrease with increasing f , Monte Carlo simulations of axon loss have also shown opposite trends 28,58 .The behavior of D ∥ e and D ⊥ e during pathological conditions is still an ongoing topic of investigation, as the complex processes happening in the extra-axonal space are not fully understood.Our results suggest that using functional relations between D ⊥ e and f is unjustified, and both parameters should rather be estimated independently, particularly in pathological conditions.Interestingly, for SMI and WMTI that do not employ the tortuosity constraint, several human in vivo studies have found changes in their D ⊥ e for neurodevelopment 12,36 , aging 38,39 , stroke 26,36 , Alzheimer's disease 41,42 , and multiple sclerosis 36,44 .
As an outlook, this work bridges the gap between micrometer-scale cellular architecture and the macroscopic MRI resolution, by comprehensively validating a key white matter diffusion MRI model and its four publicly available estimators, WMTI, NODDI, SMT, and SMI.The uniqueness of the present approach is in the 3-dimensional EM, the use of pathology (TBI) as well as sham-operated animals, the variety of the samples, and the large number of segmented 3d axons.The specific correspondence with histology indicates that the validated open-source Standard Model estimator will become a powerful tool for neuroscience research, providing valuable non-invasive imaging markers to study brain development, aging, disease, disorders and injuries at the cellular level where these processes originate.

DATA AND CODE AVAILABILITY
The 3d EM volumes and their segmentation are publicly available at https://etsin.fairdata.fi/dataset/f8ccc23a-1f1a-4c98-86b7-b63652a809c3.The dMRI volumes will also be made publicly available.Several publicly available toolboxes were used in this study as indicated in the methods section.Extra in-house codes created for the analysis of the 3d EM and dMRI data will be made publicly available by the authors.

FIG. 1 .
FIG. 1. SM parameter maps reveal specific microstructure properties that correspond to 3d EM histology.(a) Brain sections (1 mm-thick) of sham-operated and TBI rats.(b) Location of the 3d EM samples within sham-operated rat brain and automatically segmented myelinated axons.(c) FA maps of rat brains in (a), indicating the corpus callosum (cc) and cingulum (cg) at the location where the 3d EM samples were collected, brain sagittal contours adapted from 20 .(d) In the Standard Model, a fiber fascicle is modeled by an intra-axonal compartment, in which axons are represented by cylinders with zero-radius, and an extra-axonal compartment represented by an axially symmetric tensor (see Methods).These compartments are described by parameters f , Da, D ∥ e and D ⊥e .One voxel contains multiple orientationally dispersed fascicles that contribute to the directional dMRI signal according to their FOD P(n).(e) EM images transverse to the axons of the cingulum for the rat brains in (a), with segmented IAS, EAS, and myelin.Representative segmented 3d axons are also shown.Using skeletons from thousands of segmented axon, FODs for each axon population in the volumes were constructed.Reduction in axon density, changes in axon morphology, and increase in axon dispersion can be observed for the TBI animal with respect to the shamoperated animal.(f) dMRI-derived SM parametric maps reveal specific differences between sham and TBI that are observed in corresponding 3d EM histology.

2 FIG. 2 .
FIG.2.Standard Model parameter maps from different estimators.Maps for f (intra-axonal water fraction) and p2 (inversely related to dispersion) from SMI, NODDI, WMTI, and SMT for a sham-operated and a TBI animal.Changes in TBI with respect to sham, such as the reduction in f and p2, agree with histology for all methods.However, evident quantitative differences between their maps are observed.

FIG. 3 .
FIG.3.Comparison between dMRI and 3d EM derived metrics.For each estimator, matrices show the Pearson correlation coefficients comparing 3d EM histology derived metrics (rows) against dMRI estimated SM parameters (columns), as plotted in Fig.4.Significant correlations (p < 0.05) are highlighted in bold, where * indicates that significance remains after adjusting for multiple comparisons using the false discovery rate.Supplementary Figs.S1-S4 lists also the p-values and 95% confidence intervals.Submatrices for corresponding dMRI and 3d EM metrics, in which significant correlations are expected only in the diagonal, are highlighted by dashed lines.

FIG. 4 .
FIG. 4. Scatter plots for dMRI derived parameters against histology metrics.Linear regressions are shown using all (black lines and their 95% confidence intervals), TBI-only (red lines) and sham-only (blue lines) samples.Thicker lines indicate a correlation with p < 0.05.The numerical values (using all samples) for the Pearson correlations coefficients are shown in Fig. 3, their p-values and 95% confidence intervals are listed in Fig. S1-S4.(a) SM parameters f , Da, and θp 2 from SMI, WMTI, SMT, and NODDI compared against their corresponding 3d EM metrics.(b) SM extra-axonal diffusivities D ⊥ e

FIG. 5 .
FIG. 5. Dispersion angle θ, FODs rotational invariants p l , and their exponential decay.(a) For 3d EM histology, dispersion angle θp 2 (computed from p2) corresponds to the true average angle of deviation θ of the axons with respect to the main direction of the axon population.The agreement is better for single bundle FODs (• markers) than for crossing fibers (x markers).(b) FODs rotational invariants p2, p4, p6 correlate with their histological counterparts.(c) FODs p l follow an exponential decay p l = Cλ l for l = 2, 4, .... Examples from dMRI and 3d EM in the ipsilateral cingulum for all animals.(d) Comparison between dMRI (SMI) and 3d EM derived parameters C and λ.
Research was supported by the National Institute of Neurological Disorders and Stroke of the NIH under awards R01 NS088040 and R21 NS081230, and by the Hirschl foundation, and was performed at the Center of Advanced Imaging Innovation and Research (CAI2R, www.cai2r.net),a Biomedical Technology Resource Center supported by NIBIB with the award P41 EB017183.A.S. was funded by the Academy of Finland (grant #323385) and the Erkko Foundation.H.H.L. was funded by the Office of the Director of the NIH under award DP5 OD031854.7. AUTHOR CONTRIBUTIONS R.C.L, E.F., and D.S.N. conceived the project and designed the study.R.C.L. performed the formal analysis of the 3d EM and dMRI data, extracted the microstructural metrics and compared them.A.A. segmented the 3d EM data.A.S. provided animal models, dMRI and EM imaging.H.H.L., S.C., B.A.A., and Y.L. provided software for the processing and analysis of the EM and dMRI data.A.A., J.T, and A.S. contributed with the processing and segmentation of the 3d EM data.R.A.S. contributed with the analysis of the dMRI data.R.C.L. and E.F. wrote the manuscript.E.F. and D.S.N. supervised the project.All authors commented, edited and approved the final manuscript.
8. COMPETING INTERESTS E.F. and D.S.N. are co-inventors in technology related to this research with US patents US10360472B2 and US10698065B2.
FIG.S5.Six ROIs were analyzed on each animal: cc, cg and cc+cg for the ipsi-and contralateral hemispheres.