The variability of MR axon radii estimates in the human brain

The noninvasive quantification of axonal morphology provides an exciting avenue to gain 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, i.e., diffusion-weighted MRI, has been bolstered by recent hardware developments. In this work, we present the whole brain characterization of the effective MR axon radius and evaluate the interand intra-scanner test-retest repeatability and reproducibility to promote the further development of the effective MR axon radius as a neuroimaging biomarker. We observe a coefficient-of-variability of approximately 10% in the voxelwise estimation of the effective MR radius in the test-retest analysis, but demonstrate that the performance can be improved fourfold using a customized along-tract analyses.


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 micrometerthin 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 (Waxman, 1980;Drakesmith et al., 2019).
Axons are susceptible to a variety of insults, in part due to their unique morphology and energy requirements (Conforti et al., 2014;Perge et al., 2009;Wu et al., 2014). Axonal degeneration and/or dysfunction has been linked to physical trauma, oxygen and glucose deprivation, inflammation, and gene mutations (reviewed by Stassart et al., 2018). Axonal degeneration is an early feature of neurodegenerative diseases, such as Alzheimers disease (e.g. Blazquez-Llorca et al., 2017) and Multiple Sclerosis (e.g. Evangelou et al., 2001). In injury, axonal loss may occur depending on the extent of injury in affected white matter tracts (reviewed in Nashmi and 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 neurodevelopment disorders (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 of large radius (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 extralarge axons in prefrontal white matter in brains of adults with ASD. These postmortem studies demonstrate the huge 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.
Diffusion-weighted MRI (dMRI) is a particularly relevant neuroimaging modality to probe cellular features, far below the resolution of the imaging experiment (Assaf et al., 2008;Alexander et al., 2010;Romascano et al., 2020;McNab et al., 2013;Fan et al., 2020;Sepehrband et al., 2016b;Veraart et al., 2020;Huang et al., 2020). Indeed, dMRI is sensitive to a wide range of tissue microstructural parameters because the signal is sensitized to the micrometer length scale of the diffusion of water molecules (Tanner, 1979). The development of axon diameter mapping using dMRI has been challenged by various confounding factors that resulted in a significant and contested overestimation of the axon radius using MRI (Horowitz et al., 2015;Innocenti et al., 2015;Burcaw et al., 2015;Nilsson et al., 2017;Lee et al., 2018Lee et al., , 2020a. Novel insights in biophysical modeling and hardware developments improved the accuracy of axon radius mapping significantly (McNab et al., 2013;Jones et al., 2018;Veraart et al., 2020). However, MR axon radius mapping can not replace in vivo microscopy. First, in the absence of strong a priori distributional assumptions (Assaf et al., 2008;Sepehrband et al., 2016a), the information obtained from dMRI is typically limited to a single scalar representing the en-tire underlying axon distribution (e.g. Alexander et al., 2010;Veraart et al., 2020;Fan 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 will adopt the term "effective MR radius" to refer to the MR-derived axon radius (Veraart et al., 2020). 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 et al., 1992).
Large-radii axons are nevertheless important in brain function, especially in mammals with increased brain size. First, axons with larger 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 et al., 2013). However, large-radii axons come at a disproportionate price in terms of energy use and spatial constraints (Perge et al., 2009;Knowles, 2017). Histological studies have extensively reported axon radii to be in the range 0.25 − 1 µm for human brain (Aboitiz et al., 1992;Caminiti et al., 2009;Liewald et al., 2014;Tang et al., 1997), with only 1% of all axons having a radius larger 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;Tomasi et al., 2012;Liewald et al., 2014;Perge et al., 2009). 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 et al., 2001;Schüz and Preiβl, 1996). Notably, this observation promotes MR axon diameter mapping in the human brain.
In our recent work, the accuracy of the estimation of the effective MR radius using dMRI was assessed by a comparison of microscopy and MRI in fixed white matter tissue of the rodent brain (Veraart et al., 2020). 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. 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 powder-averaged signalS µ can be modelled as: with b = q 2 δ 2 (∆ − δ/3) (Veraart et al., 2020). The b-value quantifies the diffusion-weighting strength for the monopolar Stejskal-Tanner pulse sequence with diffusion gradient duration δ and diffusion gradient separation ∆ (Stejskal, 1965;Le Bihan et al., 1986). The prefactor β ≈ f √ D a is a function of the intraaxonal signal fraction f and the parallel intra-axonal diffusivity D a . The radial signal attenuation is modeled using the Gaussian phase approximation (Murday and Cotts, 1968) of the signal from protons trapped inside a cylinder radius r: where q = γG is the diffusion-weighting wave vector with γ the gyromagnetic ratio for protons and G the gradient strength van Gelderen et al. (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, i.e., 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/ √ b. In Veraart et al. (2019) we assessed that from b = 6ms/µm 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 = 3ms/µm (Raffelt et al., 2012). In this work, we adopt the higher threshold, i.e. b = 6ms/µm, to minimize the impact of this potential confound.

Effective MR radius
In the wide pulse limit (Neuman, 1974), r in Eqs. 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 Veraart et al. (2020). As explained in detail in Veraart et al. (2020) and Burcaw et al. (2015), r 4 equals the ratio between the 6 th -order and second order moment of the axon radii distribution. The 6 th -order in the numerator arises from the combination of biquadratic relation between lnS ⊥ c (r) and r (Neuman, 1974), and of the subsequent volume-weighting that emphasizes the thickest axons by an extra quadratic factor (Packer and C, 1972;Alexander et al., 2010). Therefore, the effective MR radius is heavily weighted by the largest axons within the voxel in the long-time limit (Neuman, 1974).

Diffusion MRI experiments
Five healthy adult volunteers were recruited and data were collected on two different scanning sessions with exactly the same imaging protocol on a Siemens Connectom 3T MR scanner using a 32-channel receiver coil and 300mT/m gradient coils at the Cardiff University Brain Research Imaging Centre (CUBRIC), in the UK. For each volunteer, the two scanning sessions ("test" and "retest") were performed one after the other, with only a short break (10 minutes) in between them. For both sessions, the subjects were positioned by the same operator. For one subject, the experiment was repeated 8 weeks after the initial experiments on the Siemens Connectom 3T MR scanner of the Max Planck Institute for Human Cognitive and Brain Science (MPI-CBS), in Germany using identical imaging protocols.
Data were collected under the approval of (a) the Cardiff University School of Psychology Ethics Committee (CUBRIC) and (b) the ethics committee of the Medical Faculty at Leipzig University. The participants gave written informed consent before participation in the study.
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 mins per session. The additional data with b ≤ 2.5ms/µm 2 were only used for diffusion tensor imaging (DTI; (Basser et al., 1994)) and diffusion kurtosis imaging (DKI; (Jensen et al., 2005))analysis. The total scan time, covering both sessions, was 55 minutes.
The acquisition protocol was optimized for the highest attainable precision in the estimation of the effective MR radius by evaluating the Cramér-Rao Lower bound on the variance of the estimator of the radius (Aitken and Silverstone, 1942). The MRI protocol was optimized for the number and amplitude of the b-values while considering the following constraints; (a) The diffusion gradient duration and separation was forced to be constant across the b-values to exclude the time-dependency of D a as a potential confounding factor (Lee et al., 2020b); (b) The minimal b-value was chosen to be 6 ms/µm 2 to suppress the extra-cellular signal (Veraart et al., 2019).
The data of this study are available from the corresponding author upon reasonable request.

Image processing
The diffusion-weighted images were corrected for Gibbs ringing (Kellner et al., 2016), geometric susceptibility-and eddy current distortions and subject motion , signal outliers . The bvalues were scaled to account for the gradient nonlinearities (Bammer et al., 2003;Rudrapatna et al., 2018). The data from both scan sessions were aligned to their common midway space using a rigid transformation prior (Maes et al., 1997).
The spherically-averaged signalS µ (b) was estimated as the zeroth order spherical harmonic coefficient for each b-value. The spherical harmonic coefficients, up to the 6 th order, were estimated using a Maximum Likelihood estimator to account for the Rician data distribution (Sijbers et al., 1998). The spatially varying noise level was estimated prior to the fitting to boost the precision of the estimator (Veraart et al., 2016).
The diffusion tensor and kurtosis tensor coefficents were estimated by fitting the DKI model to diffusion-weighted iamges with b ≤ 2.5ms/µm 2 using a weighted linear least squares estimator (Veraart et al., 2013).

Along-tract analysis
In addition to a voxelwise estimation of the effective MR radii, we estimated the effective MR radius along the length of individual tracts. Therefore, we introduce a strategy that is inspired by methods such as along fiber quantification (AFQ; Yeatman et al., 2012).
The spherically-averaged signals were compressed in twenty interspaced segments per white matter tract prior to computing the effective MR radius. For each tract, this compression included the following steps: (1) Generate 10,000 tract-specific streamlines; (2) For each b-value compute the interpolated value ofS µ (b) in each of the streamline nodes, which must be spaced along each streamline by a distance smaller than the voxel size; (3) Compute the center line of the fiber bundle (Klein et al., 2007); (4) Divide the center line in N segments with equal length L N , here N = 20. The tangent of the center line in the midpoint of the i th segment forms the axis of a cylinder with height L N . The radius of the cylinder equals the maximal distance between an individual streamline and the tangent line within the segment; (5) The signals of all streamline nodes within the segment-specific cylinder are averaged. To minimize partial voluming with neighboring tissue, the contribution of an individual streamline is weighted by the inverse of its distance to the centerline. Once this segment-averagedS µ (b) is computed for each bvalue, 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 novel strategy is suited to biophysical models that are derived from rotationally-invariant signal features and that 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 Coefficient of Variation and the Intraclass Correlation Coefficient (McGraw and Wong, 1996).
The test-retest coefficient of variation (CoV) 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 i th voxel x i , respectively.
The intraclass correlation coefficients (ICC) were calculated for two-way mixed effects, single measurement, with absolute agreement. ICC estimates were interpreted based on the following guidelines. Values less than 0.5 are indicative of poor reliability, values between 0.5 and 0.75 indicate moderate reliability, values between 0.75 and 0.9 indicate good reliability, and values greater than 0.90 indicate excellent reliability. ICC values were interpreted considering the 95% confidence interval.
The CoV 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 voxelwise analysis of the dMRI data.

Whole brain characterization of the effective MR radius
In Fig. 1, the voxelwise 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 Figs. 2 and 3.
In Fig. 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 Fig. 3, the trend of the effective MR radius in the midsagittal cross section of the CC is shown. This cross section covers the various bundles of the CC, including the rostrum, genu, body, and splenium of the CC. The box plots show the median and 95% confidence interval of the effective MR radius across the 5 subjects.
In Fig. 4, the along-tract analysis of the effective MR radius are shown. We show the average trend and its confidence interval, computed across all 5 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.

Correlation matrix
In Fig. 5, the correlation matrix is shown to display the Pearson's Correlation Coefficient ρ that was computed between all pairs of diffusion metrics: fractional anisotropy (FA), mean diffusivity (D), radial diffusivity (D ⊥ ), axial diffusivity (D ), mean kurtosis (K), radial kurtosis (K ⊥ ), axial kurtosis (K ), S µ (b = 1),S µ (b = 6),S µ (b = 30), β, and r. The calculation of ρ included segmented WM voxels for which all diffusion metric were within biophysically plausible bounds. The selected voxels are obtained from both the test and retest data of the five subjects, but are limited to the WM voxels with a tract density exceeding the subject-specific 75 th percentile to minimize partial voluming effects. In total 218,562 voxels were included in the correlation analysis.
In addition, scatter plots show the relation between r and FA,D ⊥ , andS µ (b = 30). The effective MR radius r shows no to small correlations with other diffusion metrics. No significant correlations were observed with radial kurtosis K ⊥ . A very weak linear relationship |ρ| < 0.1 was observed for FA, D ⊥ , K , S µ (b = 6), and β. The correlation coefficient between r and D,K, K ⊥ , andS µ (b = 6) was weak with ρ = 0.15, -0.18, -0.15, and -0.18, respectively. A strong negative correlation was observed between r andS µ (b = 30) with ρ = −0.61. This correlation analysis demonstrates that the effective MR radius provides additional information to the various metrics that are more routinely used in dMRI studies.

Test-retest reliability
In Fig. 6, Bland-Altman plots show the agreement between repeated measurements. We show the Bland-Altman plots for the five intra-scanner test-retest experiments and the single interscanner 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 voxelwise estimation of the MR axon radius was quantified by the CoV. When including all segmented WM voxels, the CoV varied between 8.84 to 11.83% across the subjects on the first scanner. Hence, the CoV 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 (CoV= 11.46%) was similar to the intra-scanner reproducibility (CoV = 11.83 and 10.20%, for the CUBRIC and MPI-CBS scanner, respectively).
In Table 1, we list all the CoV per subject and per tract. We observe small fluctuation across tracts, but overall, the testretest reliability is fairly homogeneous across the brain. A dramatic reduction of CoV, by on average a factor of four is observed, when evaluating the CoV for along-tract segments instead of voxels. In Table 2, we tabulate the voxelwise ICC for each subject and tract. The test/retest reliability ranged from moderate to good with ICC values ranging between 0.50 and 0.84 with a median value of 0.70. Evaluating the lower bound

Intra-scanner Scanner #1 vs Scanner #1
Inter-scanner Scanner #1 vs Scanner #2 Figure 6: The Bland-Altman shows high reproducibility between scan sessions on the same (left) and different scanner (right) and lack of systematic errors. The solid lines represent mean difference and ±1.96 × standard deviations of the difference. on the 95% confidence interval did not alter the classification of the reliability. For the approach in which the effective MR radius is computed in a tract-specific segment instead of voxels, the reliability is good to excellent with a median ICC of 0.84.

Discussion
The MR axon radius is a sensitive metric for the in vivo and noninvasive detection and quantification of large axonal (and possibly glial) projections in the white matter tracts. We demonstrated this by evaluating the along-and inter-tract variability in comparison to repeatibility and reproducibility of the metric. In comparison to our previous work (Veraart et al., 2020), we optimized the acquisition protocol in consideration of participant comfort and limitations imposed by in vivo research studies by 1) reducing the total acquisition time to 24 mins.
The acquisition of two b-shells with strong diffusion-weighting is required for the accurate and precise estimation of the effective MR axon diameter. The precision of the estimator benefits from a wide range of b-values with a minimal value of b = 6 ms/µm 2 to filter out the extra-cellular signal. We found that b = 30 ms/µm 2 provided a good compromise between T E and diffusion time. It is crucial that the gradient duration δ, separation ∆, and magnitude G are tuned to a minimal sensitivity to the axon diameter at b = 6 ms/µm 2 and a maximal sensitivity at the second b > 6 ms/µm 2 -shell. To realise these conditions in practice, for a given T E , we first maximize ∆ (and δ) to achieve b = 6 ms/µm 2 with minimal G; (with our experimental setup G = 122 mT/m for b = 6 ms/µm 2 ). Next, we maximize G while keeping the ∆ and δ constant to achieve the highest attainable b-value. This strategy will result in the maximal contrast of dMRI signal inside narrow axons for a given T E . Widespread deployment of this technique in human MRI is currently limited by the availability of ultra-strong diffusion-weighting gradients (McNab et al., 2013;Jones et al., 2018).
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 et al. (2016a) studied the accuracy of various parametric distribution to describe the axon distribution in the mouse corpus callosum. All well-fitting distributions are described by at least two parameters (Sepehrband et al., 2016a), so trying to reconstruct the parametric distribution from the effective MR radius alone is ill-posed. This problem is highlighted in Fig. 7.
In Fig. 7 we show a scatter plot 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 distributions 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 the effective MR radius is about 18% in the splenium of the corpus callosum when comparing cohorts of typically-developing children and children with ASD. With a voxelwise CoV of about 10%, one would only need a total sample size of 14, with equal representation of both cohorts, to detect the difference in the effective MR radius with a statistical power of 0.96 (Lakens, 2013;Faul et al., 2009). 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-to-recruit 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 adopt-    ing 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 CoV 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, in patients with multiple sclerosis, Yeatman et al. (2014) found sensitivity to lesions along tracts using a similar analysis approach. Moreover, powder averaging modeling approaches such as the presented axon diameter mapping and others (Kaden et al., 2016;Novikov et al., 2018;Reisert et al., 2017) entail an implicit assumption that diffusion MRI data are perfectly shelled, i.e. different gradient directions for a finite set of b-values. However, 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 the bvalues, we could not account for the potential directional variability of the gradient nonlinearities due to the need for shelled data.
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 "lowhigh-low" trend in axon radii across the CC is in fair agreement with the variations in axon radius distributions previously reported in rat (Sargon et al., 2003;Barazany et al., 2009;Veraart et al., 2020), rhesus monkey (LaMantia and Rakic, 1990), and human (e.g. Aboitiz et al., 1992). The slightly increased 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 the differences of 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 most often not reported.
The consistent inter-tract variability was also observed in a 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, i.e. 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, including, but not limited to, the curvature or undulations of the axons (Nilsson et al., 2012;Brabec et al., 2020;Lee et al., 2020a). 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). Even a relatively small fraction of these large glial processes might bias the MR axon diameter mapping due the nature of how signals are encoded.
In comparison with other regions identified from the whole brain analysis, the MR axon radius is systematically very high in the most anterior inter-hemispheric connections of the white matter, see Fig. 1. (LaMantia and Rakic, 1990) performed an indepth cytological characterization in rhesus macaque of these tissues, identifying a sub-region of the anterior commissure, the so-called basal telencephalic commissure, as having a distinct functional target. The axon radii distribution for both commissures is comprised primarily of small axons; however, the basal telencephaplic commissure is encapsulated by tightly wound glial processes that provided a bridge for neural migration in early development (Lent et al., 2005). The glial cells responsible were identified as GFAP-positive fibrous astrocytes. The high occurrence of astrocytic processes in this area are a likely contribution to the high axon radii estimates in that region.
MRI spectroscopy may offer some insight towards the hypothesis that MR axon diameter mapping, and the interpretation of biophysical modeling using MRI in general, might be influenced by glial cells. In Palombo et al. (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 telecephalic commissure.  Figure 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.

Conclusion
We demonstrated a good to excellent reliability in the quantification of micron-sized effective MR radii using human dMRI if 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 might foster the hypothesis that dMRI signals at high b-values might not be exclusivily 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.