The variability of MR axon radii estimates in the human white matter

Abstract The noninvasive quantification of axonal morphology is an exciting avenue for gaining understanding of the function and structure of the central nervous system. Accurate non‐invasive mapping of micron‐sized axon radii using commonly applied neuroimaging techniques, that is, diffusion‐weighted MRI, has been bolstered by recent hardware developments, specifically MR gradient design. Here the whole brain characterization of the effective MR axon radius is presented and the inter‐ and intra‐scanner test–retest repeatability and reproducibility are evaluated to promote the further development of the effective MR axon radius as a neuroimaging biomarker. A coefficient‐of‐variability of approximately 10% in the voxelwise estimation of the effective MR radius is observed in the test–retest analysis, but it is shown that the performance can be improved fourfold using a customized along‐tract analysis.


| INTRODUCTION
The white matter of the central nervous system is an intricately organized system of neural pathways that link together anatomical areas to create functional circuits. These neural pathways are formed by bundles of densely packed micrometer-thin axons that are responsible for the transfer of information. The caliber of the axon and the presence of myelin are the most important morphological determinants that control the conduction velocity of action potentials (Drakesmith et al., 2019;Waxman, 1980).
In injury, axonal loss may occur depending on the extent of injury in affected white matter tracts (reviewed in Nashmi & Fehlings, 2001).
There is also evidence to suggest that altered features of axons, such as distributions of axon calibers or focal swellings, may contribute to various pathologies (Bartzokis, 2011) and neurodevelopmental disorders (Raven et al., 2020;Stassart et al., 2018). For example, in an animal model of Angelman syndrome, a rare genetic disorder linked to autism, widespread reductions in white matter volumes were linked to reduced numbers of axons with large radii (Judson et al., 2017). Similar observations were made in children with autism spectrum disorder (ASD), where electron microscopy identified a lower percentage of large-radii axons in the corpus callosum compared to age-matched typical developing children (Wegiel et al., 2018). Zikopoulos and Barbas (2010) reported a significantly lower relative density of extra-large axons in prefrontal white matter in brains of adults with ASD. These postmortem studies demonstrate the potential of the noninvasive quantification of axon radii (including the ability to perform longitudinal assessment in the same individual), for understanding neuropathology in clinical research and, potentially, diagnostics. Critically, and as discussed below, diffusion-weighted MRI methods for mapping axon diameter are more sensitive to larger axons than smaller axons, and so this preferential loss of axons with large radii in various disorders means that the technique holds great promise as a biomarker.
Novel insights in biophysical modeling and hardware developments improved the accuracy of axon radius mapping significantly (Jones et al., 2018;McNab et al., 2013;Veraart et al., 2020). However, MR axon radius mapping cannot replace in vivo microscopy. First, in the absence of strong a priori distributional assumptions (Assaf et al., 2008;, the information obtained from dMRI is typically limited to a single scalar representing the entire underlying axon distribution (e.g., Alexander et al., 2010;Fan et al., 2020;Veraart et al., 2020). Second, this scalar is strongly biased towards the largest axons, with nearly no sensitivity to the bulk of the axons (Burcaw et al., 2015;Nilsson et al., 2017). Note that we here adopt the term "effective MR radius" to refer to the MR-derived axon radius . In contrast, microscopy is a highly reproducible technique for the extraction of the bulk of the smaller axons, but suffers from a poor precision in the quantification of the underrepresented larger axons (Aboitiz, Scheibel, Fisher, & Zaidel, 1992).
Large-radii axons are nevertheless important in brain function, especially in mammals with increased brain size. First, axons of large radii are capable of more rapid conduction, which is advantageous for time-sensitive processes. Second, it has been hypothesized that the large-radius axons of long-range neurons are essential to maintaining neural synchrony (Buzsáki, Logothetis, & Singer, 2013). However, large-radii axons come at a disproportionate cost in terms of energy use and spatial constraints (Knowles, 2017;Perge et al., 2009). Histological studies have extensively reported axon radii to be in the range 0.25-1 μm for human brain (Aboitiz et al., 1992;Caminiti, Ghaziri, Galuske, Hof, & Innocenti, 2009;Liewald, Miller, Logothetis, Wagner, & Schüz, 2014;Tang, Nyengaard, Pakkenberg, & Gundersen, 1997), with only 1% of all axons having a radius greater than 1.5 μm (Caminiti et al., 2009). Indeed, despite different white matter tracts having very different lengths and interconnecting entirely different functional circuits, in common they share a skewed axon radius distribution, characterized by mostly thin axons (Aboitiz et al., 1992;Caminiti et al., 2009;Liewald et al., 2014;Perge et al., 2009;Tomasi, Caminiti, & Innocenti, 2012). The radii of the bulk of smaller axons do not vary significantly across mammals with varying brain size. However, it has been shown that larger brains have more large axons and an increased maximal radius (Olivares, Montiel, & Aboitiz, 2001;Schüz & Preiβl, 1996). Notably, this observation promotes MR axon diameter mapping in the human brain.
In recent work, the accuracy of effective MR radius estimation using dMRI was assessed through a comparison of microscopy and MRI in fixed white matter tissue of the rodent brain . In addition, the feasibility of the technique for in vivo human MRI has already been established by comparing MR-derived values in the human corpus callosum to values reported in the literature Veraart et al., 2020). Here, we will: (a) present the whole brain characterization of the effective MR radius; and (b) evaluate the inter-and intra-scanner test-retest reliability (repeatability and reproducibility) to promote the further development of the effective MR radius as a neuroimaging biomarker.

| Axon diameter mapping
In an experimental regime in which the extra-cellular signal can be assumed to be fully suppressed, the spherically-averaged signal S μ can be modeled as: with b = q 2 δ 2 (Δ − δ/3) . The b-value quantifies the diffusion-weighting strength for the monopolar Stejskal-Tanner pulse sequence with diffusion gradient duration δ and diffusion gradient separation Δ (Le Bihan et al., 1986;Stejskal, 1965). The prefactor p is a function of the intra-axonal signal fraction f and the parallel intra-axonal diffusivity D k a . The radial signal attenuation is modeled using the Gaussian phase approximation (Murday & Cotts, 1968) of the signal from protons trapped inside a cylinder with radius r: where q = γG is the diffusion-weighting wave vector with γ the gyromagnetic ratio for protons and G the gradient strength (Neuman, 1974; van Gelderen, DesPres, van Zijl, & Moonen, 1994). Furthermore, D 0 is the diffusivity of the axoplasm, α m is the m th root of dJ 1 (α)/dα = 0, and J 1 (α) is the Bessel function of the first kind. Here, t c = r 2 /D 0 is the diffusion time across the cylinder. The applicability of the Gaussian phase approximation in the context of axon diameter mapping has been studied in detail by Fan et al. (2020).
Under the assumption that the extra-cellar water is relatively mobile, that is, the extra-cellular diffusivity is nonzero in the radial direction, then its spherically-averaged signal decays exponentially fast, much faster than the intra-axonal signal that decays as 1= ffiffiffi b p .
In Veraart, Fieremans, and Novikov (2019), it has been observed that from b = 6 ms/μm 2 upwards the extra-cellular signal does not contribute significantly to the dMRI signal decay in the healthy human white matter. In comparison, earlier simulation studies reported the cut-off for extra-cellular signal to be as low as b = 3 ms/μm 2 (Raffelt et al., 2012). In this work, we adopt the higher threshold to minimize the impact of this potential confound.
The associated software is available for download .

| Effective MR radius
In the wide pulse limit (Neuman, 1974), r in Equations (1) and (2) denotes the effective MR radius, a scalar metric that represents the entire axon radius distribution captured within a single voxel, with minimal loss of accuracy . As explained in detail in Burcaw et al. (2015) and Veraart et al. (2020), r 4 equals the ratio between the sixth-order and second order moment of the axon radii distribution. The sixth-order in the numerator arises from the combination of biquadratic relation between ln S ⊥ c r ð Þ and r (Neuman, 1974), and of the subsequent volume-weighting that emphasizes the thickest axons by an extra quadratic factor (Alexander et al., 2010;Packer & Rees, 1972). Therefore, the effective MR radius is heavily weighted by the largest axons within the voxel or, more specifically, the largest Martin's radius in case of non-cylindrical axons (Andersson et al., 2020). Diffusion weighting was applied with b = 0.5,1,2.5,6, and 30 ms/μm 2 , for 30, 30, 30, 120, and 240 gradient directions that were isotropically distributed on a sphere, respectively (Jones, Horsfield, & Simmons, 1999). The diffusion gradients were characterized by Δ/δ = 30/15 ms and maximal gradient amplitude of 273 mT/m-see

| Diffusion MRI experiments
Supplementary material for considerations regarding this protocol design. The following scan parameters were kept constant: TR/ T E : 3500/66 ms, matrix: 88 × 88, and 54 slices with a spatial resolution of 2.5 × 2.5 × 2.5 mm 3 . Data were acquired with a multi-band blipped-CAIPI accelerated (SMS = 2) EPI sequence with additional GRAPPA acceleration (R = 2) (Setsompop et al., 2012). Partial Fourier encoding was turned off. In addition, non-diffusion-weighted images were acquired with the same (N = 23) and reversed phase encoding (N = 10) for susceptibility-induced geometrical distortion correction.
The axon diameter mapping pipeline employed here only uses the diffusion-weighted data with b = 6 and 30 ms/μm 2 , which were acquired in 24 min per session. The additional data with b ≤ 2.5 ms/μm 2 were only used for diffusion tensor imaging (DTI; Basser, Mattiello, & LeBihan, 1994) and diffusion kurtosis imaging (DKI; Jensen, Helpern, Ramani, Lu, & Kaczynski, 2005) analysis. The total scan time, covering both sessions, was 55 min.
The spherically averaged signal S μ b ð Þ was estimated as the zeroth order spherical harmonic coefficient for each b-value. The spherical harmonic coefficients, up to the sixth order, were estimated using a Maximum Likelihood estimator to account for the Rician data distribution (Sijbers, den Dekker, Scheunders, & Van Dyck, 1998). The spatially varying noise level was estimated prior to the fitting to boost the precision of the estimator (Veraart, Fieremans, & Novikov, 2016).
To minimize partial voluming effects, we retained the top 75% voxels from the tract density map. Moreover, we applied TractSeg for the automated segmentation of individual fiber tracts using the fiber orientation distribution functions as estimated using constrained spherical deconvolution (Wasserthal, Neher, & Maier-Hein, 2018  After computing the segment-averaged the effective MR radius can be estimated with a higher precision. Since the spherical mean of the signal is formally a rotationally invariant feature (Mirzaalian et al., 2016), the curvature of the underlying fiber within the segment is assumed not to have a significant impact on the signal averaging and therefore on the estimation.
This presented strategy is also suited to biophysical models that are derived from rotationally-invariant signal features (Novikov, Fieremans, Jespersen, & Kiselev, 2019;Raven et al., 2020) and that might suffer from poor robustness to noise or partial volume effects.

| Statistics
The voxel-wise test-retest reliability of each diffusion metric θ was evaluated per subject and per tract using the test-retest variability and the intraclass correlation coefficient (McGraw & Wong, 1996).
The test-retest variability (TRV) of estimates of parameter θ was computed across N voxels as: with Δ θ (x i ) and μ θ (x i ) the difference and the average of the test and retest estimates of parameter θ in the ith voxel x i , respectively. The TRV and ICC were also computed using the nodes derived from our tract-specific approach instead of the voxels to evaluate the increase in test-retest reliability of along-tract analysis instead of a voxel-wise analysis of the dMRI data.

| RESULTS
3.1 | Whole brain characterization of the effective MR radius In Figure 1, the voxel-wise maps of the effective MR radius are shown for 3 slices. The spatial variability of the effective MR radius is apparent, but, overall, the maps are broadly symmetrical. The inter-tract variability is further highlighted in Figures 2 and 3.
In Figure 2, the distribution of effective MR radius per tract is shown. All voxels of the left and right hemisphere, for all subjects, were considered. In addition, we show the tract-averaged effective MR radius per subject, and per tract, for the test and retest data to demonstrate that the inter-tract variability exceeds the inter-subject variability.
In Figure 3, the trend of the effective MR radius in the midsagittal cross section of the CC is shown. This cross section covers the various segments of the CC, including the rostrum, genu, body and splenium. The box plots show the median and 95% confidence interval of the effective MR radius across the 5 subjects.
In Figure 4, the along-tract analysis of the effective MR radius is shown. We show the average trend and its confidence interval, computed across all five subjects, including the test and retest data.
For completeness, the trends are also shown for each individual subject. The metric changes widely along and across the various tracts.

| Test-retest reliability
In Figure 6, Bland-Altman plots show the agreement between repeated measurements. We show the Bland-Altman plots for the five intrascanner test-retest experiments and the single inter-scanner experiment.
For the latter, only the "test" data set from each scanner was considered.
The absolute mean difference and its 95% confidence intervals are shown. The percentage differences were 1.51, 1.00, 0.65, −1.67, and 1.17% for the 5 subjects of the intra-scanner analysis. In comparison, the percentage difference for the inter-scanner repeatability was −1.19%.
The test-retest reproducibility in the voxel-wise estimation of the MR axon radius was quantified by the TRV. When including all segmented WM voxels, the TRV varied between 7.83 to 10.48% across the subjects on the first scanner. Hence, the TRV was slightly higher, yet of the same order of magnitude, as conventional DTI metrics (approximately 4% in this study; data not shown). The inter-scanner reproducibility performance (TRV =10.16%) was similar to the intra-scanner reproducibility (TRV = 10.48 and 9.04%, for the CUBRIC and MPI-CBS scanner, respectively).
In Table 1, we list all the TRV per subject and per tract. We observe small fluctuation across tracts, but overall, the test-retest reliability is fairly homogeneous across the brain. A dramatic reduction of TRV, by on average a factor of four is observed, when evaluating the TRV for along-tract segments instead of voxels. In Table 2 The box plots represent the average effective MR radius (μm) within segments of the human CC, including rostrum (dark blue) and genu to splenium (from left to right), for each of the 5 subjects. The median across the 5 subjects is shown by the red bar, while the boxes cover the 95% confidence intervals. The segmentation of the CC is shown in the inset mid-sagittal slice repeatability and reproducibility of the metric. In comparison to our previous work , we optimized the acquisition pro- human MRI is currently limited by the need for ultra-strong diffusionweighting gradients to achieve these scan settings (Jones et al., 2018;McNab et al., 2013).
The effective MR axon radius estimation is intrinsically sensitized towards large axons; therefore, it is not representative of the full distribution of axon radii present in a voxel (Alexander et al., 2010;Burcaw et al., 2015;Veraart et al., 2020). If it were, it would be much more straightforward to compare and validate with microscopy data. Sepehrband, Alexander, Clark, et al. (2016) studied the accuracy of various parametric distributions to describe the axon distribution in the mouse corpus callosum. All well-fitting distributions were described by at least two parameters (Sepehrband, Alexander, Clark, et al., 2016), so trying to reconstruct the parametric distribution from the effective MR radius alone is ill-posed. This problem is highlighted in Figure 7 where we show the relation between the average and effective radius of Gamma distributions with varying shape α and scale γ parameters. In order to obtain a unique mapping from the effective to the much lower average radii, at least one parameter, α or γ, must be known or chosen a priori. Unfortunately, distributions with different α and γ result in the same effective radius, with a widely varying average radius, while the corresponding shape of the distribution are all plausible candidates in describing realistic axon distributions. When fitting a Gamma distribution to the distributions shown in Wegiel et al. (2018), we conclude that the shape of the axon radius distribution, both in terms of α or γ, is significantly different between typically-developing children and children with ASD. This prevents us from fixing one of the distribution parameters and, as such, from mapping the effective to the average radius accurately.
On the bright side, the effective MR radius might be a very sensitive metric to distinguish between cohorts. Based on data reported in the postmortem study of Wegiel et al. (2018), we estimate that the percentage difference in effective MR radius is about 18% in the F I G U R E 4 The trend of the effective MR radius r (μm) along the tract (posterior to anterior or inferior to superior) for each individual measurements (5 subjects and 2 repetitions) are shown in shaded lines. In addition, we show the average (solid) and 95% confidence intervals (dashed) splenium of the corpus callosum when comparing typically developing children and children with ASD. With a voxelwise TRV of about 10%, one would only need a total sample size of N=14, with equal representation of both cohorts, to detect the difference in the effective MR radius with a statistical power of 0.96 (Faul, Erdfelder, Buchner, & Lang, 2009;Lakens, 2013). This preliminary power analysis highlights the feasibility of MR axon radius mapping in future research studies.
The good inter-site agreement also bodes well for the future inclusion of effective MR diameter mapping in studies of rare or difficult-torecruit disorders or diseases, where multi-site studies are needed to achieve the necessary statistical power to make a robust inference about a given pathology (cf. Jahanshad et al., 2013).
The test-retest reliability was good to excellent when adopting the along-tract analysis (Yeatman et al., 2012). This strategy is perfectly suited to the technique because the data are rotationally invariant along tract segments. As segments are averaged prior to model fitting, there is no loss of accuracy due to the underlying curvature of the tract (Mirzaalian et al., 2016). With a TRV of only a few percent, even in an inter-scanner evaluation, the sensitivity of the MR axon radius mapping was able to detect subtle changes across the brain, both inter-and along-tract within a single subject. Therefore, the technique is sensitive down to the individual segments. This is a potential limitation for very local effects; however, Yeatman, Wandell, and Mezer (2014) found sensitivity to lesions in patients with multiple sclerosis using along tract segments.
Moreover, "spherically-averaging" modeling approaches, such as the axon diameter mapping presented here and others (Kaden, Kruggel, & Alexander, 2016;Novikov, Veraart, Jelescu, & Fieremans, 2018;Reisert, Kellner, Dhital, Hennig, & Kiselev, 2017) entail an implicit assumption that dMRI data are perfectly shelled; that is, different gradient directions exist for a finite set of b-values. This assumption is usually unmet due to gradient nonlinearities and poses an unnecessary constraint on experimental design (Bammer et al., 2003). Although we corrected for the gradient nonlinearities by a spatially dependent scaling of the b-values, we could not account for the potential directional variability of the gradient nonlinearities due to the need for shelled data (see Afzali, Knutsson, Özarslan, and Jones (2020)). Alternatives to the spherical mean may well prove useful in the future, but is beyond the scope of the current work.
The consistent tract variability of the estimated effective MR radius within and between subjects is in agreement with histological data across various species. The observed "low-high-low" trend in axon radii across from the genu to splenium of the CC is in fair agreement with the variations in axon radius distributions previously reported in rat (Barazany, Basser, & Assaf, 2009;Sargon et al., 2003;Veraart et al., 2020), rhesus monkey (LaMantia & Rakic, 1990), and human (e.g., Aboitiz et al., 1992). The slightly larger radii in the rostrum of the corpus callosum are in qualitative agreement with Sargon et al. (2003). Overall, this result demonstrates that the effective MR radius detects differences in axon radius distribution across the CC. Furthermore, the large radii of the CST, in comparison to other tracts, have been reported extensively in the literature (Tomasi et al., 2012). However, the comparison of MR to histology is challenged by the need for the full axon radius distribution, which is usually not reported.
The consistent inter-tract variability was also observed in the recent work of Huang et al. (2020). However, despite some common findings, we did not observe an anterior-posterior gradient in effective MR radii. The difference in results might be rooted in modeling choices. In our approach, we aim to minimize modeling assumptions by using the acquisition itself, that is, high b-value, to filter out both extra-cellular signal and orientational dispersion prior to modeling the intra-axonal signal.
Even after eliminating orientational dispersion, it is not understood why we observe such strong along-tract variability of the MR axon radii mapping in certain tracts (e.g., CST). Various additional confounding factors have been discussed, for example, diameter variations due to curvature or undulations of the axons (Andersson et al., 2020;Brabec, Lasič, & Nilsson, 2020;Lee, Jespersen, et al., 2020;Lee, Paioannou, Kim, Novikov and Fieremans 2020;Nilsson, Lätt, Ståhlberg, van Westen, & Hagslätt, 2012). However, in agreement with our previous hypothesis, this observation might also be explained by the lack of specificity of the high b signal to axons only. The density and variability of cells in the human brain is considerable. It has been observed that the numerous and dispersed glial processes, from astrocytes in particular, are abundant in the white matter (Luse, 1956;Oberheim et al., 2009;Perge et al., 2009) and have radii that exceed even the largest axons (Oberheim et al., 2009). In comparison with other regions identified from the whole brain analysis, the MR axon radius is systematically very high in the most anterior interhemispheric connections of the white matter, see Figure 1. (LaMantia & Rakic, 1990) performed an in-depth cytological characterization in rhesus macaque of these tissues, identifying a functionally distinct sub-region of the anterior commissure -the basal telencephalic commissure. The axon radii distribution for these interhemispheric projections is comprised primarily of small axons; however, the basal telencephaplic commissure is encapsulated by glial processes remaining from neural migration in early fetal development (Lent, Uziel, Baudrimont, & Fallet, 2005 biophysical modeling using MRI in general, might be influenced by glial cells. In Palombo, Ligneul, and Valette (2017), the effective MR radius of glial processes was potentially larger than neuronal processes as shown with diffusion-weighted signal decay at high b-values for various metabolites. In along-tract analysis of metabolites, the ratio of Choline (Cho) to N-acetyl aspartate (NAA) had similar trends to MR axon radius in the CST (Govind et al., 2012). NAA is highly concentrated in neurons, whereas Cho have been shown to originate predominantly from glial cells. Further investigation is needed to disentangle these effects, which may be aided by tissue regions with high levels of naturally present glial cells, such as the basal telencephalic commissure.

| CONCLUSION
We demonstrated a good to excellent reliability in the quantification of micron-sized effective MR radii using human dMRI if ultra-strong diffusion-weighted gradients are employed. As a result, we were able to observe the subtle inter-and along-tract variability that has previously been reported in histological studies. However, our results foster the hypothesis that dMRI signals at high b-values might not be exclusively sensitive to neuronal processes or axons and that the contribution of glial processes to the dMRI signal needs to be better understood to allow for an unambiguous interpretation of morphological parameters such as the effective MR radius.

ACKNOWLEDGMENT
Research was performed as part of the Center of Advanced Imaging Evans for his support of the MR Lab at CUBRIC.

CONFLICT OF INTEREST
The Max Planck Institute for Human Cognitive and Brain Sciences has an institutional research agreement with Siemens Healthcare. NW was a speaker at an event organized by Siemens Healthcare and was reimbursed for the travel expenses. NW holds a patent on acquisition of MRI data during spoiler gradients (US 10,401,453 B2).

DATA AVAILABILITY STATEMENT
The data from this study are available from the corresponding author upon reasonable request.

Jelle Veraart
https://orcid.org/0000-0003-0781-0420 F I G U R E 7 A scatter plot between the average and effective radius of various Gamma distribution with varying shape (α) and scale (β) parameters are shown. Per color, the scale parameter is fixed. The full distribution is shown for each distribution with an effective radius of 3 (marked by the solid circles). Notably, the corresponding average radii vary widely