Mapping tissue microstructure of brain white matter in vivo in health and disease using diffusion MRI

Abstract Diffusion magnetic resonance imaging offers unique in vivo sensitivity to tissue microstructure in brain white matter, which undergoes significant changes during development and is compromised in virtually every neurological disorder. Yet, the challenge is to develop biomarkers that are specific to micrometer-scale cellular features in a human MRI scan of a few minutes. Here, we quantify the sensitivity and specificity of a multicompartment diffusion modeling framework to the density, orientation, and integrity of axons. We demonstrate that using a machine learning-based estimator, our biophysical model captures the morphological changes of axons in early development, acute ischemia, and multiple sclerosis (total N = 821). The methodology of microstructure mapping is widely applicable in clinical settings and in large imaging consortium data to study development, aging, and pathology.


Introduction
Diffusion magnetic resonance imaging (dMRI) is a non-invasive technique that maps the probability density S(t, x) of water molecules' displacements x(t) in each imaging voxel (Jones, 2010).With typical displacements ⟨x 2 (t)⟩ ∼ 10 µm during diffusion times t ∼ 50 ms used in the clinic, the dMRI signal becomes uniquely sensitive to how tissue structure on the micrometer scale restricts the diffusion of water molecules, opening a window into cellular-level details such as cell density, shape, orientation, and permeability of cell membranes (Novikov et al., 2019;Alexander et al., 2019).Thanks to this unique in vivo contrast, dMRI is particularly promising in detecting mi-crostructural changes related to developmental and pathological processes in the brain white matter (WM), including myelination and demyelination, axonal growth and axonal loss, pruning, beading, and inflammation (Horsfield and Jones, 2002).
The greatest technical challenge of clinical dMRI is to uncover the exact relationship between cellular-level features and the dMRI signal -i.e., to make dMRI not just sensitive, but specific to tissue microstructure.This would turn an empirical diagnostic technique into a quantitative and reproducible scientific measurement paradigm enabling improved understanding of the changes that underlie development, aging and disease, and tracking of its progression.So far, widely adopted dMRI techniques, such as diffusion tensor imaging (DTI) 1 arXiv: 2307.16386v2 [physics.med-ph]19 Dec 2023 (Basser et al., 1994) and diffusion kurtosis imaging (DKI) (Jensen et al., 2005), offer ways to represent the dMRI signal S(t, q) (the Fourier transform of the displacements probability density S(t, x)) as expansions up to q 2 and q 4 , correspondingly.However, these empirical signal representations inherently lack specificity to cellular-level phenomena, as they do not rely on any assumptions about tissue microgeometry.
In the pursuit of specificity, there has been a growing interest (Novikov et al., 2018a) in biophysical models that directly parameterize relevant tissue microgeometry and thus offer ways to quantify its changes in health and disease (Novikov et al., 2019;Alexander et al., 2019;Jelescu and Budde, 2017).For WM, the Standard Model (SM) has been proposed as an overarching framework (Novikov et al., 2018b(Novikov et al., , 2019;;Reisert et al., 2017), unifying multicompartment model-based strategies over the past two decades (Kroenke et al., 2004;Jespersen et al., 2007;Assaf et al., 2004;Alexander et al., 2010;Fieremans et al., 2011;Zhang et al., 2012;Kaden et al., 2016;Reisert et al., 2017;Novikov et al., 2018b), see Fig. 1A.In the SM, an elementary fiber fascicle is comprised of two non-exchanging compartments, the intra-and extra-axonal spaces (IAS and EAS).The SM offers an exciting potential of specificity to cellular-level biological phenomena, as its scalar parameters f , D a , D ∥ e , D ⊥ e , as defined in Fig. 1A and described in detail below, are by design more specific to micrometer-level manifestations of pathological processes.
In the intracellular space, the axonal water fraction (f ) characterizes the relative contributions of IAS and EAS water.A decrease in f typically indicates axonal loss, suggesting a lower density of axons within the sampled voxel-a potential hallmark of neurodegeneration (Jelescu et al., 2016b).The intra-axonal diffusivity (D a ) assesses axonal integrity.For instance, axonal injury, such as beading, disrupts the uniform diffusion pathway by introducing variations in axonal caliber, leading to the possible slowing or transient restriction of water molecules within axons (Budde and Frank, 2010).Furthermore, within the EAS, the radial diffusivity (D ⊥ e ) reflects the condition of the myelin sheath.The process of demyelination will reduce the complexity of pathways available to water molecules moving perpendicularly to the axons, resulting in increased D ⊥ e (Jelescu et al., 2016b).The axial diffusivity within the EAS (D ∥ e ) detects a wider array of extra-axonal alterations.Decreases in D ∥ e could indicate pathological events such as astrogliosis or microglial activation, both associated with neuroinflammatory responses (Xie et al., 2010) The practical relevance of the SM is rooted in its assumptions that make it compatible with clinically feasible diffusion acquisitions as well as presence of a large number of publicly available dMRI datasets with multiple diffusion weightings b = q 2 t up to 2-3 ms/µm 2 , such as UK Biobank (Miller et al., 2016), Human Connectome Project (Glasser et al., 2016), Alzheimer's Disease Neuroimaging Initiative (Jack Jr et al., 2008) and Adolescent Brain Cognitive Development (Casey et al., 2018), together comprising hundreds of thousands of subjects.The application of the SM to these datasets presents unparalleled opportunities for large-scale in vivo studies of tissue microstructure in health and disease.
However, SM parameter estimation has proven to be challenging due to "shallow" (almost degenerate) directions in the likelihood landscape (Jelescu et al., 2016a;Novikov et al., 2018b).
Conventional maximum likelihood estimation (MLE) suffers from low precision in this degenerate likelihood landscape.Recently, machine-learning (ML) based approaches have emerged as a promising tool for increased precision and speed.In the field of quantitative magnetic resonance imaging, ML have been used to estimate T 1 and T 2 (Cohen et al., 2018), myelin water fraction (Liu et al., 2020), susceptibility (Yoon et al., 2018), and dMRI model parameters (Golkov et al., 2016;Reisert et al., 2017;Palombo et al., 2020).
Here we quantify the sensitivity and specificity of ML-based SM estimation for relatively short (∼6 min), clinically feasible dMRI protocols.Using dMRI acquired on patients during routine brain scans, we demonstrate how SM parameters are able to capture specific cellular-level changes in early development, stroke and multiple sclerosis.Our results open the way to apply this modern methodology in clinical settings and to large imaging consortium data (Miller et al., 2016;Glasser et al., 2016;Jack Jr et al., 2008;Casey et al., 2018) for investigating development, aging and pathology.

Assumptions of the Standard Model
According to the SM, the dMRI signal originates from a collection of identical fiber fascicles in a WM voxel, that are oriented based on an arbitrary orientation density function (ODF) P(n), Fig. 1A.The following SM assumptions specify the physics of diffusion inside an elementary fascicle.
• First, the signal from a fascicle is a sum of contributions from non-exchanging spin populations in the IAS and EAS.Water exchange can be neglected since myelin layers form a practically impermeable boundary for the relevant diffusion times.
• Second, the fascicle's IAS is represented as a collection of aligned zero-radius cylinders ("sticks"), such that diffusion occurs only along the stick, while the transverse diffusion is negligible since axonal diameters of ∼ 1 µm are much smaller than typical diffusion displacements in a clinical MRI measurement.The bulk along-stick diffusion coefficient D a is reduced relative to that of free water D 0 = 3 µm 2 /ms due to intra-axonal organelles and micro-variations of axonal shape, such as beads (Lee et al., 2020).(Diffusivities and b-values are in the units of µm 2 /ms and ms/µm 2 throughout this work.) • Third, diffusion in the fascicle's EAS is assumed to be anisotropic and Gaussian, characterized by the axially-symmetric tensor with parallel and perpendicular eigenvalues D ∥ e and D ⊥ e .This means that diffusion at clinical diffusion times is assumed to be in the long-time limit, and any residual diffusion time-dependence (Novikov et al., 2014a) is neglected.
This "impermeable stick" assumption has been verified in vivo in human WM by observing the distinct functional form at strong diffusion weightings that is inherent to such stick compartment (McKinnon et al., 2017;Veraart et al., 2019).The two SM compartments, IAS and EAS, define a response kernel for a fiber fascicle, which is a local bundle of aligned sticks with the extra-neurite space sur- e .The IAS is modeled as sticks with zero radial diffusivity, as axons are much narrower than the diffusion length of the dMRI measurement.Within a macroscopic imaging voxel, elementary fascicles contribute to the directional dMRI signal according to their orientation distribution function (ODF) P(n).(B) The Sensitivity-Specificity Matrix (SSM) S ij is defined in Eq. 8.An ideal SSM is an identity matrix.The diagonal elements measuring sensitivity are in bold (some diagonal elements may not be on the diagonal line because the parameters that are not estimated are left out).The nonzero off-diagonal elements reveal spurious correlations between model parameters, and are a hallmark of decreased specificity.The SSM is color-coded to highlight the elements that are greater in absolute values.
rounding them.The kernel's signal is: (1) where ξ=ĝ•n is the scalar product between the symmetry axis n of the kernel and the gradient direction ĝ.Further compartments, such as isotropic cerebrospinal fluid (CSF), ∼ f iso e −bD0 , can in principle be added.However, given typical multi-shell protocols with moderate b, the CSF fraction f iso , which is a lot smaller than the intra axonal water fraction f in WM, is very difficult to estimate.Fig. S1 shows that multi-shell protocols are not sensitive enough to estimate f iso at realistic SNR.Introducing the free-water compartment will further increase the difficulty of estimating other SM parameters, especially EAS diffusivities (as the EAS signal is similar in its functional form to that of the CSF).Therefore, in this work, we use a two-compartment kernel (1) without CSF.
Such multicompartment fascicles are distributed in a voxel based on the fiber ODF.All fascicles are assumed to have the same compartment fractions and diffusivities, and differ from each other only by their orientation (Christiaens et al., 2020).Thus, the SM signal, measured along gradient direction ĝ, is a convolution between fiber response kernel K (b, ĝ•n) and the ODF P (n) on a unit sphere: where the ODF P (n) is normalized to |n|=1 P (n)dn = 1.

Rotational invariants in the spherical harmonic basis
We factorize the kernel from the ODF parameters in Eq. (2) using the spherical harmonic (SH) basis (Reisert et al., 2017;Novikov et al., 2018b;Tournier et al., 2007): where S lm and p lm are the SH coefficients of the signal S ĝ (b) and of the ODF up to order l max which practically depends on the dMRI sampling and signal-to-noise ratio (SNR).
The functions K l (b) are projections of the kernel response onto the Legendre polynomials P l (ξ): To factor out the dependence on the choice of the physical basis in three-dimensional space (via m = −l . . .l), the rotational invariants are defined as follows (Reisert et al., 2017;Novikov et al., 2018b): From the relationship between rotational invariants and SH coefficients, one can relate the rotational invariants to the kernel parameters (Reisert et al., 2017;Novikov et al., 2018b): This enables a compression of raw directional dMRI measurements S ĝ (b) to a small number of data features S l without loss of information.Here, p 0 ≡ 1 under the ODF normalization; the remaining ODF invariants, one for each l, characterize its anisotropy, with the normalization factor chosen so that 0 < p l < 1.Among these p l , l = 2, 4, 6..., p 2 has the lowest order and thus highest SNR.
Combining p l of the ODF with the kernel parameters, the SM parameters of interest are defined as θ = f, D a , D ∥ e , D ⊥ e ; p 2 , p 4 , . . . .We will focus on p 2 , as the most easily interpretable ODF alignment metric.

Degeneracy of the estimation landscape
For any diffusion direction, the SM signal is a sum of decaying exponentials.Parameter estimation for models of such kind is generally illposed.Specific near-degenerate dimensions in the MLE landscape, have been established for SM numerically (Jelescu et al., 2016a) and analytically (Novikov et al., 2018b).In such a shallow MLE landscape, different combinations of model parameters may become indistinguishable in the presence of realistic noise, causing unstable estimation results.

Feature count of the rotational invariants
Fundamentally, the number of independent parameters one can in principle determine is tied to the information content of the dMRI signal (the number of independent features accessible from data at a given noise level).The estimation of the fascicle and ODF parameters factorizes in the spherical harmonics basis (Reisert et al., 2017;Novikov et al., 2018b), Eq. ( 7), such that the number of independent scalar signal features N b N l is a product of the number N b of the b-shells in the q-space, and the number N l of the independent rotational invariants S l (b) of the signal (constructed from its spherical harmonics S lm (b), Eq. ( 6)) accessible at a given noise level.
In this work, as well as in publicly available dMRI datasets (Miller et al., 2016;Glasser et al., 2016;Jack Jr et al., 2008;Casey et al., 2018), all our acquisitions have N b = 2 shells, and we use S l (b) with l = 0, 2, 4, such that N l = 3, yielding overall 2 × 3 = 6 independent measurements for the fascicle response.This exactly matches the number of independent SM parameters θ = {f, D a , D ∥ e , D ⊥ e ; p 2 , p 4 } contributing at these l, Eq. ( 7).Note that employing the invariant S 4 (b) is crucial, since the system ( 7 2.5.Unconstrained parameter estimation with machine learning As a faster alternative to MLE, a machine learning (ML)-based estimator was proposed to directly map rotational invariants S l (b) to the SM parameters θ i (Reisert et al., 2017).The ML-based estimator uses a "soft" prior, as the prior distribution (training sets) implicitly regularizes the estimation in the training process, instead of imposing hard constraints on the model.Here we use an extended version of this method, dubbed Standard Model Imaging (SMI) (Coelho et al., 2022), applicable to multi-dimensional dMRI.SMI uses third-order polynomial regression to map S l (b) with l = 0, 2, 4 to a set e ; p 2 is used for training the ML-based estimator throughout this study with mean [0.5, 2, 2, 0.7, 0.45] and variance [0.06, 1, 1, 0.1, 0.06].
The five aforementioned estimators (SMI, NODDI, SMT, WMTI and W-WMTI) effectively measure the same set of parameters under the SM framework, but adopt different constraints (as summarized in Table S1), therefore resulting in different outcomes and trends.It has become a crucial need to determine the most reliable (sensitive and specific) diffusion model estimator for routinely used multi-shell protocols, which will enable leveraging the enormous publicly available datasets.To evaluate the performance of SM estimators, below we propose a metric to quantify the sensitivity and specificity of parameter estimation, similar to the concept of the confusion matrix in classification.We then apply these estimators in various in vivo datasets, including early development, acute ischemia and multiple sclerosis (total N = 821), and compare them in light of the current knowledge of relevant (patho)physiological processes in the WM.

Subjects and dMRI acquisition
We studied various datasets that included twoshell dMRI scans ranging from 5 to 7 minutes long, that were acquired on patients referred for routine clinical brain MRI in the department of Radiology at New York University (NYU) and Medical University of South Carolina, indicating the potential for clinical translation of the proposed methods.Institutional Internal Review Board approval with waiver of consent was obtained for retrospective study.

Early development
For assessment of human development, brain MRIs were studied of 59 pediatric subjects (30 females) who underwent DKI imaging as part of a routine MRI exam under sedation at NYU School of Medicine from June 2009 to October 2010 (Paydar et al., 2014;Jelescu et al., 2015).The subjects ranged from 1 day to 2 years and 9 months in age, and all underwent brain MR imaging for nonneurological indications.All the included exams were interpreted as normal by fellowship-trained board-certified neuroradiologists, and were reevaluated by a board-certified pediatric neuroradiologist for normalcy prior to inclusion.

Ischemia
For assessment of (sub)acute stroke (Hui et al., 2012), clinical and MRI data was reviewed for consecutive patients at the Medical University of South Carolina who were admitted due to acute onset of neurological symptoms and were subsequently diagnosed with acute or subacute stroke in the middle cerebral artery territory as the cause for neurological impairments.A total of 28 patients admitted to this institution between August 2011 and February 2012 were included.These patients underwent MRI 7 hours to 3 weeks after symptom onset (82% of the stroke patients were scanned within the first week of symptom onset).Patients with a history of brain neoplasm or intracranial hemorrhages were excluded from study.

Multiple sclerosis and healthy controls
For assessment of multiple sclerosis (MS), we studied 177 subjects (age 48.47 ± 9.78 years old, 119 females) identified with a clinical diagnosis of MS using the McDonald criteria (Polman et al., 2011) who were referred for MRI of the head at NYU Langone Health between November 2014 and June 2020.Within one year of the MRI, disability status was assessed using the Patient Determined Disease Steps (PDDS) questionnaire, a validated nine-point patient-reported metric of disease severity (Kister et al., 2013).MS patients were separated into mild MS (0 ⩽ PDDS ⩽ 3) and severe MS (4 ⩽ PDDS ⩽ 7) based on the need of canes for walking.
In total, 557 healthy controls (age 45.29 ± 13.94 years old, 388 females) were selected from headache patients with normal brain MRI, and no history of neurological disorder.The subjects were referred for MRI of the head at NYU Langone Health.To compare with the MS patients, 177 of healthy subjects (age 48.47 ± 9.76 years old, 119 females) matched by age and sex were chosen as controls.Moreover, 177 adults aged between 25 and 35 years old (age 30.28 ± 2.91 years old, 94 females) were selected out of the cohort to establish the normative values of SM parameters to compare with the pediatric population.
To mitigate the partial volume effects, we shrink each WM region by 1 voxel relative to the 1mm template.The genu of corpus callosum (GCC) was chosen as the region of interest (ROI) for the MS study because it is a large homogeneous region in the corpus callosum with relatively few outliers.
MS lesions, identified using icometrix (Rakić et al., 2021), are notably heterogeneous (Rovaris et al., 2005) with potential exchange between the IAS and EAS in the case of unmyelinated axons and possibly additional compartments due to increased inflammation.Characterizing MS lesions is beyond the scope of this study.Hence, we focus on comparing the normal appearing white matter (NAWM) in the MS subjects with that of healthy controls.For the stroke patients, the WM mask was determined by fractional anisotropy greater than 0.2 to include more voxels to the ROI for small ischemic lesions.
SM parameters were estimated using the five WM estimators described, i.e.SMI, NODDI, SMT, WMTI, W-WMTI.The mean of an ROI was extracted for further analysis after excluding outliers.The voxels with unphysical SM parameter values were first excluded, then parameter values ±2σ away from the ROI mean were considered as outliers, where σ is the standard deviation within an ROI.Typically, fewer than 5% of the voxels are discarded.

Statistical analysis
An exponential function (A•exp(−t/τ )+B) was fitted to the dataset of early development.The absolute value of τ was used to quantify the pace of the exponential growth or decay.For the stroke patients, the mean of ischemic lesions and their contralateral hemisphere in the WM were compared pairwise.The relative changes of ischemic lesions from their contralateral regions were calculated to quantify the degree of change for each SM parameter.In the MS study, MS patients were separated into two group: mild MS patients, severe MS patients based on the PDDS score.ANCOVA was used to study the group difference between every two groups covarying for age.

Sensitivity-Specificity Matrix
To quantify sensitivity and spurious correlations of parameter estimation, we consider the Sensitivity-Specificity Matrix (SSM) in noise propagations whose elements quantify relative changes of an estimated parameter θj with respect to the actual change of parameter θ i .Here, the angular brackets denote averaging over the distribution of ground truths (the test set) of all N θ parameters.The normalization by the mean values ⟨θ i ⟩ is introduced for convenience, to make the off-diagonal elements dimensionless (and are redundant for the diagonal elements).Practically, we evaluate the SSM from a linear regression of the estimates θj with respect to ground truths θ i of all N θ parameters in a test set.Likewise, a linear regression of the estimated parameters was applied to the prior mean of SMI to demonstrate the dependency of ML-based estimation on the prior.We can define a matrix quantifying such impact as where µ θi are the mean values of the prior distribution for each model parameter.While fixing the variance of the prior distribution of SM parameters at [0.06, 1, 1, 0.1, 0.06], the prior mean was varied from 90% to 110% of the reference [0.5, 2, 2, 0.7, 0.45] at the step of 2.5% for each parameter separately.
The synthetic data was generated based on the two-compartment SM using a two-shell protocol (same as in vivo controls).The ODF was simulated by spherical harmonics up to l = 4.The ratio of p 4 to p 2 was set between 0.75 and 0.85 according to histology results (Lee et al., 2019).Gaussian noise was added to the signal at SNR = 25 with respect to b 0 images.To evaluate the SSM, the ground truth of

Sensitivity and specificity of SM parameters in simulations
Results in Fig. 1B show that SMI provides the most trustworthy estimates of SM parameters.In particular, it estimates p 2 , f and D a almost free of spurious correlations.On the other hand, the SSM (Fig. 1B) reveals that NODDI parameters have non-negligible spurious correlations with other SM parameters, notably SSM(D a , p2 ) = 0.51 and SSM(D a , f ) = −0.48.These spurious correlations are particularly concerning in the case of a significant change in D a , e.g., in pathology, that would translate into apparent changes of f and p 2 .The estimation of D a by SMT has a combination of contrasts from D a , D ∥ e and f , and the estimation of D a by WMTI has a combination of contrasts of D a and p 2 .These spurious correlations revealed by the SSM are caused by limited information obtainable from multi-shell dMRI scans and imposing hard constraints on the SM.
We show in Fig. S3 the relationship between SM estimates and the prior mean given the same variance, which suggests the lower the sensitivity and specificity are, the more strongly the estimates are influenced by the prior distribution (training set).Their relationship is close to linear at the realistic SNR of dMRI (Fig. S3).Furthermore, the scatter plots of estimated parameters against ground truth for numerical noise propagations are shown in Fig. S4A for all five estimators.Specificity of the SM parameters has also been validated through histological studies.Coronado-Leija et al. ( 2023) recently performed in a rat model of chronic traumatic brain injury an extensive comparison of diffusion MRI derived SM parameters against histology, evaluating different models including SMI and NODDI.They has demonstrated that the SM parameters for fiber dispersion are in excellent agreement with those derived from 3d electron microscopic images.Furthermore, the intra-axonal diffusivity agrees with the estimate from histology (based on the variation in the axon diameter).This work provides robust validation for SM parameters and demonstrate their specificity to geometric microscopic properties of WM tissues.
SM parameters have been demonstrated to be sensitive to various WM processes.Notably, Jelescu et al. (2016b) reported that demyelination leads to a decrease in f and an increase in D ⊥ e , in a mouse model of de-and remyelination, suggesting there is no one-to-one correspondence between these two SM parameters, and both myelination and demyelination can influence f and D ⊥ e simultaneously.On the other hand, pruning will predominantly reduce the anisotropy of axon fibers (p 2 ) and to a lesser extent also might result in axonal loss (f ).To maintain clarity in interpretation of the following SM findings, we focus on discussing the primary effects of specific processes.

Disentanglement of sequential processes in early development
Pediatric subjects (N=59) aged between 0 and 3 years old were selected for the study of early development (Paydar et al., 2014;Jelescu et al., 2015).As a reference for the pediatric data, 177 healthy controls aged between 25 and 35 were also selected.More detailed parameter distribution for healthy controls in the WM are presented in Fig. S4B.Fig. 2 shows the changes in the SM parameters of p 2 , D ⊥ e and f during early development in WM.
The results by NODDI for newborns are consistent with previous studies using NODDI (Dean III et al., 2017;Kunz et al., 2014).DiPiero et al. (2023) offers a comprehensive review of dMRI studies focusing on early development.The overall trends of early human brain development are largely consistent across the five estimators for p 2 , D ⊥ e and f (Fig. 2A).Yet, they differ in the pace of developmental processes, which can be quantified by the time constant τ of an exponential functional form A exp(−t/τ ) + B for its dynamics.
On the other hand, Fig. S5 shows that the compartmental diffusivities are remarkably stable across the early development up to adulthood.Potentially higher fluid content in early development would result in an overestimation in D ∥ e and an under-estimation of D a by SMI, which is not the case here.
Pruning is generally completed earlier than when the bulk of myelination occurs, as evidenced by animal histology (e.g., rats (Gorgels, 1990), cats (Remahl and Hildebrand, 1982;Berbel and Innocenti, 1988), rhesus monkeys (LaMantia and Rakic, 1990)) and human imaging studies (Natu et al., 2019;Cafiero et al., 2019;Fletcher et al., 2021).Pruning removes redundant axon collaterals and synaptic connections, increasing the anisotropy of ODF (p 2 ).SMI stands out as the only estimator that captures the rapid pruning after birth, as it shows p 2 has the smallest time constant among three SM parameters p 2 , f and D ⊥ e .Both axon diameter growth and myelination span extensively, persisting into adulthood (Miller et al., 2012).Results from SMI indicate that D ⊥ e and f are characterized by different time constants and carry distinct information, as opposed to NODDI linking these two parameters by the tortuosity relation.Moreover, the "impermeable stick" assumption of the SM might not hold for axons that are not yet myelinated.Thus, the markedly low axonal water fraction at birth likely indicates the fraction of impermeable, myelinated axons rather than all axons.Therefore, f is more indicative of myelination, with SMI exhibiting a time constant of 1.03±0.35years.Among all SM estimators, this is closet to the 1∼2 years of time constant measured from myelin water fraction (Deoni et al., 2012).Furthermore, regional variances of diffusion in early development have been explored in various studies.Paydar et al. (2014) and Jelescu et al. (2015) discussed the order of development in DKI, WMTI and NODDI, respectively.Their findings are generally consistent with the estimations of SMI.Regarding the intra-axonal fraction (f ) estimated by SMI, we observed the following time constants: corticospinal tract (0.51 ± 0.38), splenium of corpus callosum (0.53 ± 0.16), posterior limb of internal capsule (1.06 ± 0.34), and genu of corpus callosum (1.21 ± 0.37).These early development trends in individual WM regions conform to the neuroscience principle that WM matures in a posterior-to-anterior and inferior-to-superior manner (Colby et al., 2011).

Detection of axonal beading in ischemic lesions
Subjects (N=28) suffering from stroke and imaged with MRI from 6 hours to 2 weeks after ischemic onset were selected for the study of (sub)acute ischemia (Hui et al., 2012).Fig. 3A shows the relative change of SM parameters in ischemic WM lesions compared to the same regions in the contralateral hemisphere.It is exemplified by the parametric maps of a stroke patient scanned 26 hours after the onset of ischemia in Fig. 3B.Remarkably, the D a map of SMI reveals that D a drops below 0.5 µm 2 /ms in a large portion of the ischemic region while it is around 1.5 µm 2 /ms in the contralateral hemisphere.In line with the parametric maps, the bar plots of SMI show D a and D ⊥ e experience the largest decrease of ∼ 30% averaged over all stroke patients, while f increases by ∼ 20% and D ∥ e decreases by ∼ 15%.
Cytotoxic edema, commonly observed after a stroke, arises from the movement of water from the interstitial space into cells (Heo et al., 2005).This phenomenon occurs when extracellular Na + and other cations accumulate inside neurons and astrocytes, partly due to the breakdown of energydependent extrusion mechanisms (Liang et al., 2007).dMRI is acutely sensitive to cerebral ischemia; the ADC from a dMRI scan drops dramatically within the infarcted region shortly after an ischemic event.From a microstructure perspective, Budde and Frank (2010) proposed that axonal beading could lead to a significant reduction in ADC, a hypothesis supported by simulations and in vitro experiments.Furthermore, Novikov et al. (2014b) demonstrated that the structural disorder in the brain post-stroke is one-dimensional, based on diffusion time-dependence, a finding that corroborates the axonal beading model.In this context of axonal beading and cytotoxic edema, we anticipate the most profound change in intra-axonal diffusion being reduced due to beading (D a ), to a lesser extent an increase in axonal space (f ) from the beading, and a reduction in the extra-cellular diffusion, axially (D ∥ e ) and radially (D ⊥ e ), due to extracellular constriction.SMI results are consistent with anticipated changes from axonal beading and cytotoxic edema, always falling within the SM parameter range.On the other hand, NODDI and SMT produce unphysical f and p 2 values, indicating the assumptions D a = D ∥ e = 1.7 µm 2 /ms or D a = D ∥ e might not apply during ischemia.

Detection of axonal loss and demyelination in MS
MS patients (N=177, 134 mild, 43 severe) were compared with 177 age and sex-matched controls in the GCC.MS lesions have been masked out from the GCC before comparison.The severity of MS is determined based on Patient Determined Disease Steps (PDDS) questionnaire (mild: 0 ⩽ PDDS ⩽ 3, severe: 4 ⩽ PDDS ⩽ 7) (Kister et al., 2013), where the clinical distinction between mild and severe MS is determined based on the ability to walk without a cane.The mean values of GCC normalized for age are shown in Fig. 4.
MS patients show lower f and higher D ⊥ e compared to controls, which is consistent with known MS pathology, i.e., axonal loss and demyelination (Trapp and Nave, 2008).The decrease of p 2 in MS patients could be induced by the activation of microglia in the neuroinflammatory response (Voet et al., 2019).An increasing number of microglia, which are morphologically plastic and considerably larger in size than axon diameters (Lawson et al., 1990), may reduce the apparent anisotropy of ODF.While the five estimators detect largely consistent changes between controls and MS patients in p 2 , f and D ⊥ e , SMI is the only estimator that captures the continuous trend of MS severity from mild to severe cases in all three parameters.Nevertheless, the changes in D a and D ∥ e are inconsistent among the five estimators (discussed below).

Source of discrepancies between SM estimators explained by the SSM
Cellular and pathological specificity is the primary motivation for biophysical modeling of dMRI signals in brain WM.Yet, due to the limited information available in multi-shell dMRI protocols, biophysical models commonly employ model constraints to stabilize parameter estimation.These constraints tend to introduce unknown systematic biases and lead to discrepancies in parameter quantification and group comparisons.Despite efforts to compare the results of different WM estimators for the same dataset (Jelescu et al., 2015;Beck et al., 2021), the source of discrepancies so far has not been fully explained.We address this gap by proposing the SSM as a means to quantify the sensitivity and specificity of parameter estimation and thereby provide an explanation for the source of discrepancies in different estimators.The usefulness of SSM is illustrated here in various clinical datasets.
In early development, the SSM suggests SMI has the highest specificity in estimating p 2 .Indeed, SMI is the only SM estimator that captures the rapid axon elimination shortly after birth.
In ischemic lesions, the SSM predicts that SMI can capture the drop in D a , while both SMT, and to an even larger extent NODDI, will instead cause the estimated f to increase, since NODDI completely fixes D a and SMT still allows D a to be estimated.
In MS, the primary pathological processes, axonal loss, demyelination and inflammation, affect f ,

ML-based estimator vs MLE
For the multi-shell dMRI protocols employed in clinical and large-scale studies (Miller et al., 2016;Glasser et al., 2016;Jack Jr et al., 2008;Casey et al., 2018), the limited information necessitates employing constraints to stabilize MLE in a "shallow" (almost degenerate) optimization landscape (Jelescu et al., 2016a;Novikov et al., 2018b).Constraints like D a = D ∥ e in SMT, or further D a = D ∥ e = 1.7 µm 2 /ms in NODDI significantly narrow down the shallow likelihood landscape and increase precision for f and p 2 , which is why these constrained estimators have been embraced.However, such constrained estimators introduce biases to the estimation and result in spurious findings (Jelescu et al., 2015(Jelescu et al., , 2016a;;Novikov et al., 2018b;Lampinen et al., 2017).WMTI takes a different approach by relating SM parameters to DKI metrics.To establish this analytically, WMTI assumes fiber alignment, essentially forcing p 2 close to 1.It is no surprise to see all of parameters estimated by WMTI are correlated with p 2 (Fig. 1B).However, all SM parameters are not fixed throughout development, aging, and pathology.Hence, while an SM parameter can be fixed during the estimation process, the influence of a parameter may be reflected as changes in the other parameters, leading to spurious correlations.
ML-based methods provide a promising alternative and hold some unique advantages over MLE, which usually relies on applying hard constraints to achieve robustness.First, they optimize the training error of all samples as a whole rather than one sample at a time as in MLE, which allows the ML-based estimator to learn the trade-off between bias and variance.Second, ML-based estimators can learn an individual mapping from signals to every parameter, whereas MLE can only search for the most probable combination of parameters.The learning approach enhances the sen-sitivity of estimators to parameters that are confounded by others.Third, ML-based estimators use a "soft" prior distribution as training sets to regularize the estimation instead of applying hard constraints.
The estimation of p 2 or f by SMI does not depend on the diffusivity as evidenced by the lack of spurious correlations between them.This is because the ML model learns to estimate p 2 and f through the training samples with varying diffusivities and is able to resolve them.In comparison, imposing hard constraints on diffusivities like in NODDI and SMT to stabilize the estimation process will inevitably lead to spurious correlations between the estimated f and the ground truth of diffusivities, as shown in their SSM (Fig. 1).Fourth, SMI relies on rotational invariants up to S 4 (b), which is more informative than, say, SMT relying only on S 0 (b).

Limitations of the ML-based estimator
By releasing hard constraints and using a "soft" prior in the training, SMI achieves higher sensitivity and specificity than conventional MLE algorithms.However, there are still several caveats.First, despite a large number of datasets available, only very limited information regarding the tissue microstructure of WM can be obtained from multishell protocols, which results in the spurious correlations revealed by the SSM.To further improve the parameter estimation of the SM, incorporating extra "orthogonal" measurements is key in extracting complimentary information about the tissue microstructure for the ML algorithms to learn from.Planar and spherical diffusion encodings, and measuring dMRI signals at multiple echo times have been shown to significantly improve the estimation precision and accuracy (Topgaard, 2017;Veraart et al., 2018;Lampinen et al., 2020;Coelho et al., 2022;Reisert et al., 2019;Coelho et al., 2019).
Second, the SM may not be fully valid during early developmental stages or under certain pathological conditions.
During very early development, there is minimal myelin around axons, potentially leading to water exchange between the IAS and EAS at the diffusion times employed in clinical diffusion MR protocols.We hypothesize that the SM parameter estimates might only measure myelinated axonal volume fractions when applied to developing WM.This hypothesis gains support from observations made from 0 to 3 years of age.During this period, there is a notable trend: the SM parameters-including f , D ⊥ e , p 2 , gradually approach adult levels (Fig. 2).The consistent trajectories of SM parameters like f , D ⊥ e and p 2 underscore the ongoing myelination process and other developmental activities.These findings suggest that SM parameters can effectively capture developmental changes.
Third, since SMI is an ML-based estimator, it is inherently influenced by the choice of training sets.Figure S3 illustrates that alterations to the prior mean induce a shift in the overall distribution.The amount of shift for a given parameter is the complement of its sensitivity (diagonal elements of the SSM).Thus, the impact of the prior distribution is more pronounced for parameters that are more difficult to estimate, such as compartmental diffusivities.As a result, the absolute values of estimates by SMI are only comparable using the same prior distribution.The ML-based estimator aims to minimize the mean squared error, which is a composite of both bias and variance.In this minimization process, a reduction in variance is attained at the expense of increased bias.This dynamic essentially nudges the estimates closer to the prior mean.The bias introduced in this manner is systematic, meaning it will not affect the outcome in group comparisons.To better control the prior mean, we use a Gaussian prior centered at a predefined value, with a cutoff at the physical bounds.The variance of the Gaussian prior mediates the bias-variance tradeoff inherent to the estimation process.A reduction in the variance of prior enhances the precision of estimates located at the center of the prior distribution, albeit potentially increasing the bias of estimates at the outer range of the prior.
It's crucial to understand that the bias observed for the ML-based estimation is not inherently a product of the prior itself.Instead, it emerges due to the insufficient information provided by the acquisition protocol and SNR, which are imperative for accurately resolving the parameters.As detailed in the supplement, the prior distributions adopted in this study do cover the entire physical range and are general enough to process various datasets.The mean of the prior distribution is substantiated by results from an extensive protocol with maximum b-value up to 10 ms/µm 2 (Novikov et al., 2018b).The variance of the prior distribution is selected to optimize the SSM for the multi-shell dMRI protocol at realistic SNR.This prior is recommended to be used in future studies with multi-shell dMRI protocols for comparison purposes.

Conclusion
In conclusion, because of the unique advantages enabled by ML-based estimators, as well as by relying upon the complementary rotational invariants S l (b) with l = 0, 2, 4, SMI captures the biologically sensible morphological trends across three different clinical datasets, including early development, acute ischemia and MS.The SSM was proposed as a novel metric to measure the sensitivity and specificity of SM estimation, which largely explains the source of discrepancies between SM estimators and clearly demonstrates the highest sensitivity and specificity of SMI among five SM estimators of interest.Hence, SMI can serve as a powerful tool in clinical settings and for large imaging consortium data to study WM microstructure in a wide range of neurological diseases.
Figure S1: Free water estimation with two-shell protocols using NODDI and SMI.(A) NODDI estimate of free water volume fraction fiso is plotted against the ground truth f iso of synthetic data for SNR=20 and 200, and for NODDI diffusivity constraints followed and released.Correlation ρ between estimates and ground truth is indicated on each plot.(B) SMI estimate fiso against the ground truth with SNR=20 and constraints released, same as the last plot in (A).The prior distribution of f iso is uniform between 0 and 0.2.As a result, the SMI estimate fiso is roughly near the prior mean 0.1 due to lack of sensitivity to the CSF compartment.(C) Exemplary parametric maps of a 43-year-old female control.Fractional anisotropy (FA) and free water volume fraction (FW) estimated by NODDI and SMI are presented.NODDI and SMI can distinguish fiber tracts and ventricles, but within a fiber tract, NODDI exhibits noisy estimates while SMI lacks contrast.
1 Figure S3: Sensitivity-to-prior matrix.Sensitivity-to-prior matrix was evaluated by applying linear regression to the prior mean of SMI to demonstrate the dependency of ML-based estimation on the prior.As is numerically validated in this figure, it is also possible to prove analytically that S ij + P ij ≈ I for a linear model estimated by a linear regressor, and this relationship becomes an approximation for a nonlinear model or regressor.This result suggests for ML-based estimators, deviations from the identity matrix in the SSM are related to the bias that is introduced to the estimator by the prior distribution (training set).

Figure 1 :
Figure 1: Sensitivity and specificity matrix of SM estimators.(A) Elementary fiber fascicles of the SM, consisting of the IAS and EAS, are described by at least 4 independent parameters: f , Da, D ∥ e , D ⊥e .The IAS is modeled as sticks with zero radial diffusivity, as axons are much narrower than the diffusion length of the dMRI measurement.Within a macroscopic imaging voxel, elementary fascicles contribute to the directional dMRI signal according to their orientation distribution function (ODF) P(n).(B) The Sensitivity-Specificity Matrix (SSM) S ij is defined in Eq. 8.An ideal SSM is an identity matrix.The diagonal elements measuring sensitivity are in bold (some diagonal elements may not be on the diagonal line because the parameters that are not estimated are left out).The nonzero off-diagonal elements reveal spurious correlations between model parameters, and are a hallmark of decreased specificity.The SSM is color-coded to highlight the elements that are greater in absolute values.
both assume D a = D ∥ e , with NODDI further fixing D a = D ∥ e = 1.7µm 2 /ms.Both NODDI and SMT use a tortuosity model to derive D ⊥ e : D ⊥ e = D ∥ e • (1 − f ) imply specific fiber ODF shapes and impose a square-root branch choice, D a ⩽ D ∥ e for WMTI and D a ⩾ D ∥ e ) with S 0 (b) and S 2 (b) for two b-shells yields only N b N l = 4 independent measurements for 5 nonlinearly-coupled parameters θ = {f, D a , D ∥ e , D ⊥ e ; p 2 }.Practically, the signal invariants S l (b) decrease quickly with l; fortunately, maps of S l (b) up to l = 4, as shown in Fig. S2A, display clear WM structure.Fig. S2B-C further demonstrates that S 4 (b) remains informative for the acquisitions discussed in this work.

Figure 2 :
Figure 2: Development trends of diffusion parameters in WM.Parameters are ordered by the pace of development processes (based on histology studies): pruning (p 2 ) occurs faster than myelination (D ⊥ e ), followed by axonal growth (f ).Dots represent the WM mean of the pediatric subjects and error bars indicate their 95% confidence interval.The exponential fit of the development data is plotted as a solid line within the 95% confidence interval, while the 95% confidence interval of its time constant τ is indicated on each plot.As a reference for the pediatric data, the dashed line and its neighboring shaded area represent the mean and standard deviation of the WM mean from 177 controls aged between 25 and 35.

Figure 3 :
Figure 3: White matter microstructure parameter changes in ischemic lesions as compared to the contralateral hemisphere.(A) Mean relative percent changes of 28 subjects in diffusion parameters from normal (contralateral hemisphere) to (sub)acute ischemic tissues are shown in the bar graph (error bars indicate the 95% confidence interval).(B) Parametric maps of WM are overlaid on the b=0 images of an exemplary stroke patient (scanned 26 hours after the onset of ischemia).In the middle of the top row (within the white box) is the mean DWI (mDWI) signal averaged over different directions of the b = 1 µm 2 /ms shell, where the ischemic lesion is clearly shown.

Figure 4 :
Figure 4: Box plots of white matter microstructure parameters in MS patients and controls.In comparison to controls (blue), MS patients are separated into mild MS (yellow) and severe MS (red).The mean of GCC is shown in the box plot after correcting for age.Please note that all diffusivities estimated by NODDI are constrained, i.e.Da = D ∥ e = 1.7µm 2 /ms, D ⊥ e = D ∥ e • (1 − f ), while SMT assumes Da and D ∥ e are identical but not fixed.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 '+' marker symbol.The significance levels of the statistical ANCOVA test are displayed in asterisks on top of every two groups (*: p<.05; **: p<.01; ***: p<.001).

Figure S2 :
Figure S2: Rotational invariants of the dMRI signals.(A) Parametric maps of rotational invariants from a 3-month-old female infant and a 28-year-old female adult.b-value is in the unit of ms/µm 2 .(B) All rotational invariants of GCC mean excluding MS lesions (corrected for age) exhibit a significant decline (p-value indicated on the top right corner) as the severity of MS disability increases.On the leftmost position of the x-axis, 'C' stands for controls.(C) To account for the SNR differences across subjects, the magenta dots represent the simulated rotational invariants of typical SM parameter combinations for newborn infants at the mean SNR level detected in each pediatric subject.In the trajectory of S 4 (b = 1) alone, which has a lower SNR than the rest of rotational invariants as shown in (A), the effect of varying SNR overtakes the change caused by development.

Figure S4 :
Figure S4: Histograms of the SM parameter estimation.(A) Simulations of realistic noise propagation for all estimators.Here, the estimates of each available SM parameter are plotted against the ground truth in scatter plots.Brighter color suggests higher data point density.Synthetic data are generated based on the SM for a two-shell protocol (b = 1, 2) and realistic signal-to-noise ratio SNR = 25 for b = 0.The Pearson correlation coefficient ρ between estimates and ground truth is indicated on each plot, where stronger correlations imply higher sensitivity.(B) Probability distributions of SM parameters from in vivo data.Distributions are made up of over 200,000 WM voxels from 177 young controls aged between 25 and 35. 4

Figure S5 :
Figure S5: Development trends of compartment axial diffusivities.Early development trends between age 0 and 3.The dots represent the mean value of the entire WM for the pediatric subjects and the error bars indicate its 95% confidence interval.As a reference for the pediatric data, the dashed line and its neighboring shaded area represent the mean and standard deviation of the WM mean for 177 controls aged between 25 and 35.

Table 1 :
Conceptual comparison between WM estimators