Noninvasive quantification of axon radii using diffusion MRI
Abstract
Axon caliber plays a crucial role in determining conduction velocity and, consequently, in the timing and synchronization of neural activation. Noninvasive measurement of axon radii could have significant impact on the understanding of healthy and diseased neural processes. Until now, accurate axon radius mapping has eluded in vivo neuroimaging, mainly due to a lack of sensitivity of the MRI signal to micronsized axons. Here, we show how – when confounding factors such as extraaxonal water and axonal orientation dispersion are eliminated – heavily diffusionweighted MRI signals become sensitive to axon radii. However, diffusion MRI is only capable of estimating a single metric, the effective radius, representing the entire axon radius distribution within a voxel that emphasizes the larger axons. Our findings, both in rodents and humans, enable noninvasive mapping of critical information on axon radii, as well as resolve the longstanding debate on whether axon radii can be quantified.
Introduction
Axons facilitate connectivity between distant neurons. Along with myelination, the axon radius determines the conduction velocity, thereby shaping the timing of neuronal computations and communication (Waxman, 1980). Using a model of action potential neurophysiology (Rushton, 1951), it has been shown that the axon radius explains the largest proportion of variance in conduction velocity (Drakesmith et al., 2019). Histological studies demonstrated that axon sizes vary widely within the human brain, ranging from 0.1$\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ to more than 3$\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ (Aboitiz et al., 1992; Innocenti et al., 2015; Liewald et al., 2014), and across species (Olivares et al., 2001; Schüz and Preiβl, 1996; Liewald et al., 2014). Moreover, axon radii have been shown to be altered in various disease processes. For example, direct axon counting in postmortem tissue has suggested that smaller axons may be preferentially susceptible to axonal injury in multiple sclerosis (Evangelou et al., 2001) due to inflammation (Campbell et al., 2014). Electron microscopy has revealed a higher percentage of smallradius axons and a lower percentage of largeradius axons in several anatomically and functionally distinct segments of the corpus callosum in autistic subjects compared to healthy controls (Wegiel et al., 2018). From the animal literature, morphometric analysis of adult rat brains showed reduced axonal radii without axonal loss after chronic alcohol feeding (Kjellström and Conradi, 1993). Such studies indicate that noninvasive metrics capable of reporting on features of the axon radius distribution could provide important neuroimaging biomarkers for basic research and clinical applications.
A particularly relevant neuroimaging modality attuned to the microarchitecture of living brain tissue is diffusionweighted MRI (dMRI). dMRI is sensitive to the thermal motion of water molecules and their interference with microscopic boundaries, such as imparted by cells and subcellular structures in the brain (Tanner, 1979; Le Bihan, 2003; Le Bihan et al., 1986; Callaghan et al., 1988; Basser et al., 1994; Jones, 2010; Beaulieu, 2002; Novikov et al., 2019). Applications of dMRI specialize in revealing macroscopic brain connections (Jbabdi et al., 2015) and in the interpretation of contrast differences in diffusionweighted images (Moseley et al., 1990; Baron et al., 2015). However, reproducible and specific biomarkers for studying disease onset and progression noninvasively and quantitatively in the entire brain, in particular visavis axonal properties, would confer clear advantages. Several studies have used various methods to report on axon radius parameters; still, despite many attempts, axon radius mapping using dMRI remains highly contested (Assaf et al., 2013; Horowitz et al., 2015; Alexander et al., 2010; Innocenti et al., 2015; Xu et al., 2014; Burcaw et al., 2015; Ong et al., 2008; Ong and Wehrli, 2010). Discrepancies between histology and dMRIderived axon radii uncovered various confounding factors, for example orientation dispersion (Drobnjak et al., 2016; Nilsson et al., 2012), timedependent extraaxonal diffusion overshadowing the intraaxonal signal at low diffusion weighting (Burcaw et al., 2015; Fieremans et al., 2016; Lee et al., 2018), weak signal attenuation for typically very narrow axons, especially in the realistic experimental regime of long diffusion gradient duration (van Gelderen et al., 1994; Neuman, 1974), and/or putative shrinkage during tissue preparation (Barazany et al., 2009; Innocenti et al., 2015; Aboitiz et al., 1992).
Recent advances in biophysical modeling and hardware prompted a revival of MR axon radius mapping (McNab et al., 2013; Huang et al., 2015; Jones et al., 2018). First, several of the most crucial confounding factors have been removed using powderaveraging concepts (Callaghan et al., 1979; Jespersen et al., 2013; Kaden et al., 2016). Averaging diffusionweighted signals that are isotropically distributed on a sphere with constant diffusionweighting strength $b$ has been shown to factor out the orientation dispersion (Jespersen et al., 2013; Kaden et al., 2016; Mollink et al., 2017), thereby eliminating one of the most important confounding factors in axon radius mapping (Nilsson et al., 2012). Second, gradient systems capable of producing relatively strong gradient pulses have been introduced in human scanners (Jones et al., 2018). Third, it has been shown that dMRI can be made specific to a particular water population restricted into long, yet micrometerthin cylindrical objects by imparting high diffusionweighting regimes (McKinnon et al., 2017; Veraart et al., 2019). Often, an axon is too narrow to yield a measurable diffusionweighted MR signal decay, hence the popular use of ‘sticks’ (Behrens et al., 2003; Kroenke et al., 2004) when referring to axons (and possibly glial cell processes) within the context of biophysical modeling of white matter using dMRI.
The intuition behind promoting specificity to intraaxonal water comes from Callaghan’s model (Callaghan et al., 1979) of diffusion inside infinitely narrow onedimensional randomlyoriented cylinders, as applied to intraneurite diffusion by Kroenke et al. (2004). The spatial Fourier transform ${e}^{{D}_{a}^{\parallel}{(\mathbf{\mathbf{q}}\widehat{\mathbf{\mathbf{n}}})}^{2}t}$ of the diffusion propagator (with respect to the diffusion wave vector $\mathbf{\mathbf{q}}$) for a single stick as measured with MRI (Callaghan, 1991), averaged over the orientations $\widehat{\mathbf{\mathbf{n}}}$ of the sticks, yields the asymptotic scaleinvariant power law $\overline{S}=\int \phantom{\rule{negativethinmathspace}{0ex}}\mathrm{d}\mathrm{cos}\theta \phantom{\rule{thinmathspace}{0ex}}{e}^{{D}_{a}^{\parallel}{q}^{2}t{\mathrm{cos}}^{2}\theta}\sim 1/{b}^{\alpha}$ as a function of the diffusion weighting parameter $b={q}^{2}t$ (Le Bihan et al., 1986), with the scaling exponent $\alpha =1/2$. Evidently, this power law scaling should be only approximate, for $q\ll 1/r$, where $r$ is the cylinder radius. Its observation (McKinnon et al., 2017; Veraart et al., 2019) in the range $6\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}/\mu {\mathrm{m}}^{2}\le b\phantom{\rule{thinmathspace}{0ex}}\le 10\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}/\mu {\mathrm{m}}^{2}$ is a manifestation of the insensitivity of dMRI to the transverse axonal dimensions on clinical scanners. However, for sufficiently strong diffusion weighting, the power law scaling eventually breaks down, and the dMRI measurement becomes sensitive to the axonal diameter.
Technically, this work addresses the detection and the interpretation of the deviation from the radiusinsensitive $\alpha =\mathrm{\hspace{0.17em}1}/2$ power law signal behavior at the largest possible $b$ (by varying $q$ at fixed diffusion time $t$), in rat and human white matter. Indeed, either sensitivity of MR to a finite axonal radius, or a notable exchange rate between intra and extraaxonal water at the clinical dMRI time scales $t$∼100 $\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}$, would alter the very particular power law scaling (Kroenke et al., 2004; Jensen et al., 2016; McKinnon et al., 2017; Veraart et al., 2019).
Following theoretical considerations, we demonstrate the breaking of the power law scaling at very high $b$values in ex vivo rodent brains, from which metrics associated with the axon radius distribution can be mapped quantitatively. Confocal microscopy of the rat corpus callosum (CC) validated that (i) the signal arises mainly from the intraaxonal space, and (ii) the MRderived axon radius metrics are in good quantitative agreement with those derived from histology. We then observe the same signal signatures in living human brain on the Connectom 3T scanner, that is, a high performance research scanner with a maximal gradient amplitude of 300 mT/m – a fourfold increase compared to stateofthe art clinical scanners (Glasser et al., 2016). Our findings both validate the mechanism with which axon radii are weighted in dMRI (Burcaw et al., 2015), and demonstrate the accuracy of which properties of the radius distributions can be estimated. After validating and evaluating our methodology in rat and human brain, we further discuss the impact of axon radius measurements in health and disease.
Theory
Power law scaling
In most biophysical models for diffusion in white matter, axons (and possibly glial cell processes) are represented by zeroradius impermeable ‘sticks’, characterized by locally onedimensional diffusion, that is radial intraaxonal diffusivity ${D}_{a}^{\u27c2}\equiv \mathrm{\hspace{0.17em}0}$ (Kroenke et al., 2004; Behrens et al., 2003; Jespersen et al., 2007; Jespersen et al., 2010; Fieremans et al., 2011; Sotiropoulos et al., 2012; Zhang et al., 2012; Novikov et al., 2014; Novikov et al., 2018; Novikov et al., 2019; Jensen et al., 2016; Reisert et al., 2017; McKinnon et al., 2017; Veraart et al., 2019). The stick model then yields an asymptotic intraaxonal orientationally averaged signal decay,
with an intercept $\gamma $ (discussed below), the power law exponent $\alpha =1/2$, and the coefficient $\beta =\sqrt{\pi /4}\cdot f/{\left({D}_{a}^{\parallel}\right)}^{1/2}$ where $f$ is the ${T}_{2}$weighted axonal water fraction (Veraart et al., 2018; Lampinen et al., 2019) and ${D}_{a}^{\parallel}$ the parallel intraaxonal diffusivity. This particular signal decay only holds in the absence of extraaxonal signal, which is assumed to decay exponentially fast and, as such, to be fully suppressed at sufficiently high $b$values (McKinnon et al., 2017; Veraart et al., 2019). Therefore, we restrict our in vivo analysis to $b>6\mathrm{ms}/\mu {\mathrm{m}}^{2}$ (Veraart et al., 2019). Our lower bound on the $b$value is significantly higher than previous predictions from Monte Carlo simulations (Raffelt et al., 2012), thereby minimizing the likelihood of residual extraaxonal signal contributions. For the ex vivo analysis, we increase this lower bound to $b=20\mathrm{ms}/\mu {\mathrm{m}}^{2}$ to compensate for the reduced diffusivities in fixed tissue (Shepherd et al., 2009).
Breaking of the power law
The following computations always assume that the signal is normalized to ${S}_{b=0}\equiv 1$. Sensitivity of MR to either finite axon radius or notable exchange rate between intra and extracellular water would break the ${b}^{1/2}$scaling at large $b$ as follows:
Finite axon radius: A finite ${D}_{a}^{\u27c2}>\mathrm{\hspace{0.17em}0}$ results in a truncated power law:
with ${f}_{\mathrm{im}}\equiv {S}_{b\to \mathrm{\infty}}\ge 0$ the signal fraction of a fully restricted immobile water compartment (the socalled ‘dot’ compartment) (Stanisz et al., 1997; Dhital et al., 2018; Tax et al., 2019; Veraart et al., 2019). If ${D}_{a}^{\u27c2}=0$, the power law offset (found from extrapolating the signal to $b\to \mathrm{\infty}$), should give the fraction of the dot compartment, $\gamma \equiv {f}_{\mathrm{im}}$ (Veraart et al., 2019). However, for nonzero ${D}_{a}^{\u27c2}$ and finite $b$, the Taylor expansion of,
around any finite point ${\xi}_{0}$ predicts the $\xi \to 0$ intercept $\gamma <{f}_{\mathrm{im}}$, Figure 1. The always negative difference $\u03f5=\gamma {f}_{\mathrm{im}}<0$ depends on $\beta $, ${D}_{a}$, and ${\xi}_{0}$; its maximal magnitude ${\u03f5}_{\mathrm{max}}=\beta \sqrt{2{D}_{a}^{\u27c2}/e}=f\cdot \sqrt{\frac{\pi}{2e}\cdot {D}_{a}^{\u27c2}/{D}_{a}^{\parallel}}$ is achieved at the curve’s inflection point ${\xi}_{*}^{2}=2{D}_{a}^{\u27c2}$. Hence, the lower bound ${f}_{\mathrm{im}}{\u03f5}_{\mathrm{max}}$ for the $\xi \to 0$ intercept $\gamma $ may be negative. A negative $\gamma $ is biophysically implausible if the stick model holds, ${D}_{a}^{\u27c2}\equiv 0$; however, $\gamma <0$ becomes a natural consequence of a finite ${D}_{a}^{\u27c2}>0$ (and hence, of a finite axonal diameter), Figure 1, in the case when the extrapolated negative intercept overcomes the positive immobile fraction ${f}_{\mathrm{im}}$. Recently, ${f}_{\mathrm{im}}$ was shown to be negligible in healthy human white matter (Dhital et al., 2018; Tax et al., 2019; Veraart et al., 2019). Therefore, a negative intercept is a novel hallmark of MR sensitivity to the inner axon diameter, even if the signal scaling might appear linear as a function of ${b}^{1/2}$ for $b$ranges accessible on human MR scanners. Importantly, the finite axon radius model, Equation 3, is poorly conditioned as a result of which the simultaneous estimation of ${f}_{\mathrm{im}}$ and ${D}_{a}^{\u27c2}$ is practically impossible, especially for human MR experiments, see Figure 1. An accurate and precise measurement of ${D}_{a}^{\u27c2}$ depends on the prior knowledge of ${f}_{\mathrm{im}}$ and requires a dedicated measurement (Dhital et al., 2018; Tax et al., 2019).
Exchange: The spherical integration of the twocompartment ‘Kärger’ model (Kärger, 1985) with a finite exchange rate $\mathcal{R}>0$ yields approximately the following signal decay:
(4) $\overline{S}(b)\approx \beta \left({b}^{1/2}+c\cdot {b}^{3/2}\right)+{f}_{\mathrm{im}},c\propto \mathcal{R}{T}_{\mathrm{E}}/{D}_{e}^{\u27c2}>0,$
with ${D}_{e}^{\u27c2}$ the radial diffusivity in the interstitial space, and ${T}_{\mathrm{E}}$, the echo time, during which exchange can happen. Importantly, Equation 4 is convex as a function of $\xi ={b}^{1/2}$.
The relative fit quality of the models (i.e., Equations 1, 2, and 4) to the dMRI signal decays can be evaluated qualitatively (convex versus concave shape) or statistically by means of the corrected Akaike information criterion (AICc) (Burnham and Anderson, 2002).
From diffusivity to effective MR radius
The radial signal attenuation ${S}_{c}^{\u27c2}(r)$ inside the cylinder of radius $r$ in the Gaussian phase approximation (van Gelderen et al., 1994):
with $b={g}^{2}{\delta}^{2}\left(\mathrm{\Delta}\delta /3\right)$ and ${t}_{c}={r}^{2}/{D}_{0}$ defines the connection between the intraaxonal radial diffusivity ${D}_{a}^{\u27c2}$ and the radius $r$. Here, ${D}_{0}$ is the diffusivity of the axoplasm, $g$ the gradient of the Larmor frequency, ${\alpha}_{m}$ is the $m$^{th} root of $\mathrm{d}{J}_{1}(\alpha )/\mathrm{d}\alpha =0$, where ${J}_{1}(\alpha )$ is the Bessel function of the first kind, and $\delta $ and $\mathrm{\Delta}$ are the gradient duration and separation, respectively (Stejskal, 1965).
In the longpulse limit, that is when $\delta \gg {t}_{c}$, the dependence on $\mathrm{\Delta}$ drops out (Neuman, 1974), and Equation 5 approaches the Neuman’s limit
This limit practically applies to the majority of axons. Importantly, the attenuation is proportional to the fourth power of the radius $r$ and, as such, it is very weak for narrow axons. Hence the low sensitivity of dMRI to the inner axon diameter.
For an unknown distribution $h(r)$ of axons with radii $r$, the total intraaxonal signal attenuation becomes a volumeaverage over the histogram bins ${r}_{i}$ (Packer and Rees, 1972):
such that the signal contribution of an axon scales quadratically with its radius $r$. The Taylor expansion of the net signal attenuation ${S}^{\u27c2}$ demonstrates the sensitivity of the dMRI signal to the distribution’s higher order moments:
where the effective axon radius:
captures the contribution from the whole distribution $h(r)$ in a single metric (Burcaw et al., 2015). The ability to represent the whole distribution by the ratio of its 6^{th} and 2^{nd} moments relies on almost all axons falling into the Neuman’s limit, Equation 6. Representing ${S}_{c}^{\u27c2}({r}_{\mathrm{eff}})\equiv {e}^{b{D}_{a}^{\u27c2}}$, we can calculate
the MRI estimate of ${r}_{\text{eff}}$ after estimating ${D}_{a}^{\u27c2}$ from the orientationaveraged signal using Equation 2.
Note that the effective radius, Equation 9, is heavily weighted by the tail of $h(r)$. Physically, this happens due to the combination of the weak NMR signal attenuation by small radii, $\mathrm{l}\mathrm{n}\phantom{\rule{thinmathspace}{0ex}}S\sim {r}^{4}$, in the diffusionnarrowing (Neuman’s) regime (Neuman, 1974), and of the subsequent volumeweighting that emphasizes the thickest axons by an extra factor of ${r}^{2}$ (Packer and Rees, 1972; Alexander et al., 2010). The error associated with these modeling assumptions is discussed in the Results section.
Results
Simulations
Accuracy
First, we evaluate the accuracy of axon radius mapping as a function of $r$ for axon radius distributions extracted from histology; Figure 2 (left and middle panels). We used a simulation framework based on the matrix formalism for diffusion signal attenuation within fully restricted cylinders (Callaghan, 1997), as implemented in the MISST toolbox (Drobnjak et al., 2010), while mimicking the entire experimental setup, for both the human and preclinical experiments.
In the case of diffusion restricted in a single cylinder with radius $r$, the error in the estimated radius $\widehat{r}$ increases with $r$. Indeed, the missing higherorder $\mathcal{O}({g}^{4})$ corrections to Equations (5)(6) set an upper bound on the achievable accuracy for large axons, as estimated recently (Lee et al., 2018).
The combined error in the estimation of ${r}_{\mathrm{eff}}$ associated with the approximations made in Equation 6 and Equation 8 is only 5% for the human setup when considering the axon radius distribution provided by Aboitiz et al. (1992); the distributions of Caminiti et al. (2009) result in a subpercent error. Additionally, we show the errors in the estimation of ${r}_{\mathrm{eff}}$ for the axon caliber distributions that were observed in our different histological sections while considering the scan parameters from our fixed tissue experiments. The shorter diffusion timings increase the approximation errors, leading to an underestimation up to 9%.
Feasibility and precision
Figure 2 (right panel) shows a theoretical lower bound on the 95% confidence interval in the voxelwise estimation of ${D}_{a}^{\u27c2}$ from Equation 2, as predicted using a CramérRao lower bound analysis (Kay, 1993). Using the dependence ${D}_{a}^{\u27c2}\approx {D}_{a}^{\u27c2}({r}_{\mathrm{eff}})$, Equation 8, that approximately identifies ${r}_{\mathrm{eff}}$ with the single cylinder radius in van Gelderen’s model, can be used to translate the lower bound on ${D}_{a}^{\u27c2}$ to that on ${r}_{\mathrm{eff}}$.
Notably, it follows from Figure 2, that an estimate of ${D}_{a}^{\u27c2}$ exceeds zero with a statistical threshold of $p>0.05$, if the corresponding ${r}_{\mathrm{eff}}>1.41\mu \mathrm{m}$ and ${r}_{\mathrm{eff}}>0.76\mu \mathrm{m}$, when mimicking the diffusion acquisition and SNR on the Siemens Connectom (${G}_{\mathrm{max}}=300\mathrm{mT}/\mathrm{m}$) and Bruker Aeon (${G}_{\mathrm{max}}=1500\mathrm{mT}/\mathrm{m}$) MR scanners, respectively. In comparison, for a typical acquisition on a modern clinical scanner with ${G}_{\mathrm{max}}=80\mathrm{mT}/\mathrm{m}$, this lower bound is 3.2 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ (Veraart et al., 2019).
Preclinical data
Dot compartment
Because ${f}_{\mathrm{im}}$ has been reported to be significant in fixed tissue by Stanisz et al. (1997), we first have to estimate the signal fraction of the immobile water compartment ${f}_{\mathrm{im}}$ in the three fixed brain samples from a dedicated MRI acquisition (see Materials and methods). The estimate ${\widehat{f}}_{\mathrm{im}}$ is in the range of 8 to 17% with a median value of 13%. The range is defined here by the 5^{th} and 95^{th} percentile of the distribution of the estimated dot fractions in all CC voxels, across the three samples.
Breaking the power law
Figure 3 shows the signal decay, averaged across all CC voxels, based on diffusion measurements in the three rat brain samples with $b$ up to 100$\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}/\mu {\mathrm{m}}^{2}$. We notice that an extrapolation to infinite $b$, that is $1/\sqrt{b}\to 0$, yields a small but significant negative offset $\gamma $, of the order of a few per cent of the nonattenuated ${S}_{b=0}$ signal, in all three samples after subtracting the ${\widehat{f}}_{\mathrm{im}}$ from the diffusion measurements.
We reevaluate the validity of a perfect stick assumption in the high $b$regimes using a AIC analysis. To study fit robustness with respect to the number of degrees of freedom by considering the full, nested, and extended models to Equation 1. Specifically, we evaluated the following models:
${f}_{\mathrm{im}}+\beta {b}^{\alpha}$
${f}_{\mathrm{i}\mathrm{m}}+\beta {b}^{1/2};$
${f}_{\mathrm{im}}+\beta {e}^{b{D}_{a}^{\u27c2}}{b}^{1/2}$
${f}_{\mathrm{im}}+\beta \left({b}^{1/2}+c\cdot {b}^{3/2}\right)$
$\beta {b}^{\alpha}$
$\beta {b}^{1/2}$
$\beta {e}^{b{D}_{a}^{\u27c2}}{b}^{1/2}$
$\beta \left({b}^{1/2}+c\cdot {b}^{3/2}\right)$
Our analysis shows that a truncated power law (vii), which explicitly accounts for ${D}_{a}^{\u27c2}>0$ (and hence does not require a negative intercept parameter), yet sets ${f}_{\mathrm{im}}=0$, fits the experimental data significantly better than pure power law forms (models (i), (ii), (v), and (vi)) (the difference in $\mathrm{A}\mathrm{I}\mathrm{C}\mathrm{c}<2$), with or without an offset $\gamma $, if the immobile (dot) compartment is corrected for, that is when using ${\overline{S}}_{\ast}(b)=\overline{S}(b){\widehat{f}}_{\mathrm{im}}$. Without dotcorrecting the ex vivo data, the power law form (ii) with an intercept outperforms the other models. In that case, the intercept $\gamma $ is negative, while ${f}_{\mathrm{im}}$ is defined to be positive. Hence, the intercept encodes both the still water fraction and a negative offset to the intercept associated with the sensitivity to the axon diameter, such that the overall $\gamma <0$.
Axon radius estimation and histological validation
Axon radii were estimated from the diffusion MRI data for the different CC ROIs (Figure 4) along with the axon radius distributions extracted histologically, Figure 4. The errors between the associated tailweighted effective radii and MRderived ${r}_{\mathrm{MR}}$ vary between 5 and 21% in the different ROIs. Notably, a consistent residual overestimation was observed, whereas the previous simulations (Figure 2) predicted an underestimation.
To further examine the correlations between the MRderived parameters and underlying microstructure, we analyzed 16 patches with the same size as an imaging voxel, that is 100 × 100 $\phantom{\rule{thinmathspace}{0ex}}\mu {\mathrm{m}}^{2}$ within the genu of the CC of the second sample, CC#2. The confocal light microscopy images of two of those patches are shown for the various stainings, Figure 5. Two notable biological components other than axons were highlighted, namely, astrocytic processes and (neuronal or glial) cell bodies, which were found to have volume fractions of about 5%. Although the radius of the astrocytic processes cannot be measured due to their random orientations w.r.t. the image plane, it is clear that some of the processes have a large diameter, for example up to 7 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ in the first patch. The average cell body radius is 2.6 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$, with an effective radius of 4.3 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. It is worth highlighting that the ${T}_{2}$weighted signal fraction of both cell types remains unknown since the corresponding relaxation times are unknown. This unknown difference between compartmental volume (histology) and signal (MRI) fractions remains the Achilles’s heel of comparisons between MRI measurements and histological evaluations of tissue microstructure.
Within each of the 16 patches, we extracted the axon radii distribution and derived the average $\overline{r}$ and effective radius ${r}_{\mathrm{eff}}$. The box plots of those metrics are shown in Figure 5. The median $\overline{r}$ and median effective radius ${r}_{\mathrm{eff}}$ across all patches are 0.61 and 1.06 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ , respectively. In comparison, the median ${r}_{\mathrm{MR}}$, derived from dMRI in 16 voxels within the genu of the CC, is 1.16, 1.10, and 1.19 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ for the three rat samples, respectively. The median MRderived effective axon radius is between 81 and 97% larger than the median $\overline{r}$, whereas the error to the median ${r}_{\mathrm{eff}}$, as derived from histology, is only 3.4 to 12.8%.
Parameter maps
ROI measurements provided robust estimation, but a remaining question is whether dMRI could be used to map the effective MR radius in a voxelwise manner. Figure 6 shows the maps of the MRderived effective axon radii for all three rat CC’s. The maps appear smooth with very few outlier voxels, suggesting that the estimation is robust even when voxelwise data is used. Furthermore, the qualitative trends are in good agreement with previously reported observations of larger axons in the body of the CC in comparison to the genu and splenium (Barazany et al., 2009). Intersubject variability is not very large and can be attributed to slightly different slice positions. When computing the effective radius of the CCaveraged signal ${\overline{r}}_{\mathrm{MR}}$, the intersubject variability nearly nullifies. Indeed, we estimate ${\overline{r}}_{\mathrm{MR}}=1.22,1.25$ and 1.25 µm in the three samples, respectively.
Towards human applicability
Breaking the power law
To assess the applicability of this approach in more realistic settings available for human imaging, experiments were performed in human subjects on the Connectom scanner, which is capable of producing 300 mT/m gradient amplitudes. The dMRI signal decay curves, averaged across all WM voxels, with $b$values up to 25 ms/µm^{2} for the four human subjects are shown in Figure 3. Importantly, we find that — in excellent correspondence with the previous preclinical data — the linear extrapolation of the signal decay as a function of $1/\sqrt{b}$ to $1/\sqrt{b}\to 0$ produces a significant negative offset $\gamma $ in all subjects.
Note that the dot compartment was not measured directly, because previous dedicated studies revealed a negligible dot compartment, that is ${f}_{\mathrm{im}}=0$, in the healthy white matter (Dhital et al., 2018; Veraart et al., 2019; Tax et al., 2019); see Discussion.
The AICc analysis of various models demonstrated that also for the human white matter, the truncated power law (vii) with ${D}_{a}^{\u27c2}>0$ and negligible dot compartment ${f}_{\mathrm{im}}=0$ fits the data significantly better than pure power law forms, with or without intercept. However, this statistical analysis cannot be interpreted as a datadriven justification for ${f}_{\mathrm{im}}=0$ because of the degeneracy of Equation 3, as highlighted in Figure 1.
Comparison with histology
Since direct histological evaluation in volunteers is unfeasible, we turn to validate the MRderived metrics in humans with previously reported literature of human corpus callosum microstructure (Aboitiz et al., 1992; Innocenti et al., 2015). In Figure 7, the MRderived metrics were directly compared with axon radius distributions of multiple histological studies (Aboitiz et al., 1992; Innocenti et al., 2015). Various reports and histological studies show a good correspondence for the bulk of the distributions, represented by the average radius $\overline{r}$, that is the average radius $\overline{r}$ only ranges between 0.54 and 0.69 µm. In histological samples, the corresponding effective radius ${r}_{\text{eff}}$ dominated by large axons, shows strong variability. Indeed, compared to $\overline{r}$, ${r}_{\text{eff}}$ varies more than threefold, from 0.91 to 2.9$\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$.
The four human subjects show good correspondence in terms of ${r}_{\text{MR}}$. In Figure 7, we show the individual and combined distributions describing ${r}_{\text{MR}}$ for all voxels in the midbody of the CC for all four subjects. It is apparent that the combined distribution falls almost entirely within the range spanned by ${r}_{\text{eff}}$values as predicted from histology – even without introducing a putative axonal shrinkage factor (maximally 35% [Aboitiz et al., 1992], and typically within 10% [Tang et al., 1997]).
Parameter distribution and maps
In Figure 8, the distribution and map of ${D}_{a}^{\u27c2}$ for WM voxels in all human subjects, estimated using the ODFindependent model, Equation 2 with ${f}_{\mathrm{im}}=0$, are shown. Considering the statistical bound from Figure 2, it is to be expected that the estimated ${D}_{a}^{\u27c2}$ is biophysically meaningful for the vast majority of WM voxels for the Connectom scanner (Figure 8 shows a representative map of ${D}_{a}^{\u27c2}$ and associated ${r}_{\mathrm{MR}}$ for a single subject of the Connectom cohort), whereas similar measurements on a modern clinical scanner result in a biophysically implausible negative ${D}_{a}^{\u27c2}$ in approximately 35% of all WM voxels. Note that the data from a clinical scanner (Siemens Prisma with 80 mT/m gradients) are adopted from our recent work (Veraart et al., 2019). This suggests that estimating ${D}_{a}^{\u27c2}$ and the associated effective axonal radius ${r}_{\mathrm{MR}}$ is only possible on MR systems with ultrastrong gradients (Jones et al., 2018; Huang et al., 2015). The spatial variability as well as the observed asymmetry between the hemispheres in, for example, the occipital lobes was noted for all subjects. However, our cohort is too small and not sufficiently characterized to study the whole brain characterization or the role of lateralization in large axons of human brain (Eichert et al., 2019; Liewald et al., 2014; Lebel and Beaulieu, 2009).
Gray matter
It is worth examining the power law scaling also in areas outside the white mater. Therefore, Figure 9 shows the diffusionweighted signal decay, averaged over all cortical gray matter (GM) voxels as a function of $b$ in the human subjects. The signal scaling in the WM is shown for qualitative comparison. The nonlinear scaling of the isotropicallyaveraged signal as a function of $1/\sqrt{b}$ of all human subjects indicates strong deviations from the ‘stick’ model in the cortical GM, (McKinnon et al., 2017; Palombo et al., 2019). Accounting for a finite neurite radius, Equation 2, does not describe the data well either. Instead, the convex signal decay as a function of $1/\sqrt{b}$ at high $b$values is in good agreement with the anisotropic exchange model that we derived from the expansion of the anisotropic Kärger model in the powers of inverse $b$, Equation 4. Both the finite radius and exchange model predictions are shown in the absence of an immobile water fraction. The exchange model fits the data better than all other evaluated models in all subjects according to an AIC analysis (data not shown). The residence time within the neurites $1/\mathcal{R}$ varies from approximately 10 to 15 $\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}$ or 20 to 30 $\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}$ if we assume ${D}_{e}^{\u27c2}=1\mu {\mathrm{m}}^{2}/\mathrm{ms}$ or ${D}_{e}^{\u27c2}=0.5\mu {\mathrm{m}}^{2}/\mathrm{ms}$, respectively. Dedicated experiments are required for a more precise measurement of the exchange rate.
Discussion
What do we measure with dMRI?
Noninvasively estimating metrics associated with axon radius distributions is a formidable task, yet it could have a strong impact for numerous areas of research including neuroscience, biomedicine and even for clinical research and applications. Histological studies have extensively reported axon diameters $2r$ to be in the range 0.5 − 2 µm for human WM (Aboitiz et al., 1992; Caminiti et al., 2009; Liewald et al., 2014; Tang et al., 1997), with only 1% of all axons having a diameter larger than 3 µm (Caminiti et al., 2009). A vigorous debate has emerged in the MRI and neuroanatomy communities as in vivo, MRIderived axon diameters are reported to fall within the range 3.5 −15 µm (Alexander et al., 2010; Horowitz et al., 2015; Huang et al., 2015). On the MRI side, the discrepancy has been attributed to the long diffusion pulses that strongly reduce the signal attenuation of protons restricted in a narrow cylinder (van Gelderen et al., 1994; Burcaw et al., 2015). Therefore, the timedependence of extraaxonal diffusion ${D}_{e}^{\u27c2}(t)$ (Burcaw et al., 2015; Fieremans et al., 2016; Lee et al., 2018), and the undulation or alongaxon caliber variation (Nilsson et al., 2012; Brabec, 2019; Özarslan et al., 2018; Lee et al., 2019), potentially overshadow the relatively small ${D}_{a}^{\u27c2}$.
On the other hand, shrinkage during tissue fixation has been suggested as a potential shortcoming of histology (Barazany et al., 2009; Horowitz et al., 2015), implying that in vivo axons are thicker than their histologically reported values.
This study aimed at investigating what the dMRI signal can measure in terms of axon radius, as well as provide insight into the aforementioned debate. Our wide range of diffusion weightings in both human and preclinical dMRI enables a suppression of the extraaxonal contribution (that otherwise biases the radii Burcaw et al., 2015; Fieremans et al., 2016), thereby allowing us to shed light on this controversy. We claim that the effective MR radii measured in this study (${r}_{\mathrm{MR}}$) quantitatively agree with those derived from histology — to the extent that histology correctly captures the tail of the axonal radii distribution $h(r)$. That is, ${r}_{\mathrm{MR}}$ obtained from dMRI appears to be a selfaveraging quantity in each imaging voxel, as large MRI voxels ensure that the moments of $h(r)$ sampled from a voxel represent well the ‘true’ underlying $h(r)$ in that WM region, so that the spatial variations in ${r}_{\mathrm{MR}}$ stem mainly from genuine biological variations of the tails of axon distributions across the brain.
Mesoscopic fluctuations
Histologyderived ${r}_{\mathrm{eff}}$ are prone to mesoscopic fluctuations due to small sampling sizes, Figure 7. Despite a good correspondence of the bulk of axon radii distributions obtained from different histological studies and samples (Aboitiz et al., 1992; Caminiti et al., 2009; Liewald et al., 2014; Tang et al., 1997), the tail of the distribution is typically coarsely sampled, with only a few spikes representing the occasional observation of large axons within the relatively small histological sections, Figure 7. It is precisely for the detectable large $r$, that the relative fluctuations for the bin counts ${N}_{i}$ are observed for bin values of ${N}_{i}\sim 1$ (Table 2 of Aboitiz et al. (1992) and Figure 5), according to the Poissonian statistics governing each ${N}_{i}$. Not surprisingly, ${r}_{\text{eff}}$ derived from discrete histological histograms exhibits strong fluctuations, as depicted by dotted vertical lines in Figure 7a and the error bars in Figure 5.
Humans
Although the average radius $\overline{r}$, as reported in human literature (Aboitiz et al., 1992; Innocenti et al., 2015), only ranges between 0.54 and 0.69$\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$, the corresponding ${r}_{\text{eff}}$ varies from 0.91 to 2.9$\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. In comparison, the average dMRIderived ${r}_{\text{MR}}$ estimated from the four Connectom data sets within the same regionofinterest, the midbody of the CC, only varied from 2.48 to 2.82$\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ (Figure 7b).
Rodents
The average radius, as measured in this study, varies between 0.54 and 0.68 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ across the 16 patches of the genu of the CC, while the associated ${r}_{\text{eff}}$ varies from 0.81 to 1.30 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. The MRderived effective axon radii ${r}_{\text{MR}}$ vary similarly, that is 0.94 to 1.4 $\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ across several voxels within the genu of the CC for all three scanned samples.
For dMRI, the variability in the estimation of ${r}_{\text{MR}}$ is determined by thermal MRI noise, and genuine anatomical – intervoxel and intersubject – differences. For human MRI, the mesoscopic fluctuations are much weaker, due to the large MRI voxels in comparison to the histological samples. Indeed, the variance in the estimation of ${r}_{\text{MR}}$ is expected to decay inversely with the number of axons within a field of view. However, for rodent MRI, in which the MRI voxels have the same surface area as the histological patches, the precision in the estimation of the effective radius is similar for both modalities.
Overall, dMRI provides a precise measurement of the largest axons, which are captured within an MRI voxel. In contrast, histology, so far, mainly probes the bulk that consists of smaller axons with high precision. Therefore, both modalities are complementary, especially in human MRI for which the voxels are significantly larger than a typical histological sample.
Measuring the bulk of the axon distribution using MRI
As the signal attenuation inside axons, Equations (5)(6), scales as $\mathrm{ln}S\sim {g}^{2}{r}_{\mathrm{e}\mathrm{f}\mathrm{f}}^{4}$, getting to twotimes smaller ${r}_{\mathrm{eff}}$ would require another fourfold increase in gradients. However, even with stronger gradient systems, the main bottleneck might be the missing prior knowledge about the shape of the expected axon radius distribution $h(r)$. Even when assuming a particular functional form of $h(r)$, one is limited to estimating a single parameter to describe the axon radius distribution, whereas realistic distributions such as the generalized extreme value distribution (Sepehrband et al., 2016) are parameterized by at least two variables. Hence, the reconstruction of $h(r)$ from only ${r}_{\mathrm{eff}}$ is technically illposed and, as such, prone to mis– or over–interpretation due to biases towards userdefined distribution shapes and parameters, even more so if confounding factors such as dispersion or fixed diffusivities are ignored (Assaf et al., 2008; Barazany et al., 2009; Horowitz et al., 2015; McNab et al., 2013).
With unknown $h(r)$, only a single metric representing the entire distribution, that is ${r}_{\mathrm{eff}}$, for which the strength of tailweighting is determined by the gradient pulse width, can be estimated reliably. In the best case, that is the narrowpulse limit $\delta \ll {t}_{c}$, see 'From ${D}_{a}^{\u27c2}$ to effective MR radius’, ${r}_{\mathrm{eff}}$ will depend on the fourth rather than the sixth order moment of $h(r)$, that is ${r}_{\text{eff}}\equiv \sqrt{\u27e8{r}^{4}\u27e9/\u27e8{r}^{2}\u27e9}$ (Burcaw et al., 2015), thereby reducing, but not eliminating the difference between ${r}_{\mathrm{eff}}$ and $\overline{r}$. Other methods, such as oscillating gradient diffusion weighting or double diffusion encoding, may provide other sources of contrast encoding different aspects of the size distribution (Jiang et al., 2016), although the lowfrequency limit of the oscillatinggradient attenuation has been shown to be equivalent to the Neuman’s limit, not providing any extra information (Novikov et al., 2019). It can be hypothesized that the combination of methods could perhaps recover more accurate information on the underlying $h(r)$.
Dot fraction
The presence of isotropic immobile water ${f}_{\mathrm{im}}$ has been conjectured by Stanisz et al. (1997) as water possibly restricted inside the soma of various cell types, such as neurons or oligodendrocytes. Several previous studies, for example Veraart et al. (2019), Tax et al. (2019), and Dhital et al. (2019), demonstrated with various diffusion encoding strategies that in vivo dMRI is practically insensitive, that is < 0.2%, to such signal contributions in the healthy white matter of the living human brain, excluding the cerebellum (Tax et al., 2019). The lack of sensitivity of dMRI to immobile water might be explained by a small volume fraction, a short T_{2} relaxation time, and/or a fast water exchange rate on the scale of our diffusion time $\mathrm{\Delta}=30\mathrm{ms}$ for treating them as coming from separate compartments. In contrast, the dot compartment has been observed in fixed brain samples in various studies (Stanisz et al., 1997; Alexander et al., 2010). The origin of this signal contribution is not well understood yet, but the still water compartment needs to be considered when validating or studying biophysical models in fixed tissue.
In this work, for the human experiments, we build upon the previously reported observations to fix ${f}_{\mathrm{im}}=0$ in the healthy white matter to avoid fitting degeneracies that are associated with the poor conditioning of model (iii). However, any underestimation of the dot compartment, for example due to fixing ${f}_{\mathrm{im}}=0$, leads directly to an underestimation of the effective MR radius, see Figure 1. Therefore, in future studies, especially those that focus on the developing, aging, or pathological brain, we encourage the independent measurement of the dot compartment to complement the axon radius acquisitions. The fast measurement of the dot compartment is promoted by the availability of spherical diffusionencoding, as demonstrated by Dhital et al. (2019), and Tax et al. (2019).
In our ex vivo experiments, the measurement of the dot compartment is based on the diffusionweighing in the direction parallel to the average fiber direction in the CC at the maximal $b$value of 100 ms/µm^{2}. The measured signal fraction of such a still water compartment in our fixed brain samples was in the range of 8 to 17%, in line with Stanisz et al. (1997). Applying a radial or planar diffusionweighting filter prior to this measurement would suppress any contribution of anisotropic signal compartments, such as crossing or dispersed axons, to the isotropically restricted dot compartment. Although we aimed to minimize this confounding factor by using a very high $b$value (Dhital et al., 2019), the dot compartment fraction, and as such the effective MR radius, might be slightly overestimated because of various complex fiber configurations. Additional confounding factors are listed in the following section.
Confounding factors
The apparent discrepancy between histology and dMRI, when confounding factors such as extraaxonal water (Burcaw et al., 2015; Fieremans et al., 2016; Lee et al., 2018) and orientational dispersion (Drobnjak et al., 2016; Nilsson et al., 2012) are addressed, is mainly a result of the difference between $\overline{r}$ and ${r}_{\mathrm{eff}}$ – that is between the bulk and the tail of axonal distribution. This already provides an important insight into the discussion on why the radii reported in the literature vary so much between the methods. When comparing applestoapples, despite the excellent agreement observed in this study between ${r}_{\mathrm{MR}}$ and its histological counterpart ${r}_{\mathrm{eff}}$, in our own histological validation we still observed a small, yet consistent overestimation of between 5 and 20% in axon radius ${r}_{\mathrm{MR}}$ using dMRI. Aside from the previously discussed dot compartment, various other factors might contribute to this discrepancy.
First, an underlying assumption of all studies targeting the measurement of the axon radius is specificity: that the signal observed at these high $b$values could be attributed exclusively to the intraaxonal space. However, this assumption is not established nor in our opinion is it justified given that water resides in all cellular compartments of the central nervous system. We cannot exclude that water trapped in other ‘stick’like features such as the radiating processes of astrocytes contribute observable signals; it has been previously reported, but also observed in our histological sample, that such glial processes can have large diameters, up to 7$\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ in our sample. In the future, this contribution could be investigated using the increased cellular specificity of (diffusionweighted) spectroscopy (Palombo et al., 2016; Shemesh et al., 2017; Ligneul et al., 2019).
Second, Stanisz et al. (1997) and, more recently, Palombo et al. (2019) demonstrated that at shorter diffusion times, the signal contribution from cell bodies might be characterized by a specific $b$value dependent signature (Neuman, 1974 and Murday and Cotts, 1968) that might enable the extraction of MR effective cell body sizes in both the white and gray matter (Palombo et al., 2019). In this study, the potential $b$value dependent signal contribution of cell bodies was unaccounted for, and, given our and other (e.g. SampaioBaptista et al., 2019) observations of a finite cell body volume fraction, the axon radius measurements could be biased. However, deviations to the power law scaling due to the presence of cell bodies are more likely to be expected in the GM because of larger volume density of large somas in comparison to the WM (Palombo et al., 2019). In our histological sample, we observed a significant volume fraction of cell bodies in the genu of the CC, that is 5%, but due to unknown compartmental relaxation times, the associated, yet more important, signal fraction is unknown (Lampinen et al., 2019). Regardless, a biophysical model parameterized by various volume fractions, axon radii, soma sizes, and compartmental diffusivities may be poorly conditioned and degenerate.
Third, alongaxon undulations (Nilsson et al., 2012) and curvature (Özarslan et al., 2018) might result in an increased apparent radial diffusivity and, as such, contribute to an overestimation of the axon radius using dMRI, especially for increasingly long diffusion times (Lee et al., 2019; Brabec, 2019).
Finally, the estimation of the MR effective axon radius depends on the unknown intrinsic diffusivity ${D}_{0}$ of the axoplasm. In ex vivo samples of a wellaligned WM bundle, one could estimate ${D}_{0}$ directly by exploring the time dependence of the apparent diffusivity at very short diffusion times, (Mitra et al., 1993). In this study, we were not able to achieve a reproducible and precise estimate of ${D}_{0}$ and opted to use the longitudinal diffusivity $D}_{a}^{\parallel$ as a proxy for ${D}_{0}$, with ${D}_{a}^{\parallel}\le {D}_{0}$. Therefore, we might actually underestimate the positive bias in the estimation of ${\widehat{r}}_{\mathrm{eff}}\sim {\left({D}_{0}\right)}^{1/4}$ (Equation 10). However, the propagation of the error in the estimation of ${D}_{0}$ to ${\widehat{r}}_{\mathrm{eff}}$ is strongly reduced by the fourth root relation between both metrics.
Interspecies variability
In our study, the effective MR radius in humans was significantly higher than in rats when comparing similar regions of interest, for example the midbody of the CC. This difference is in agreement with several studies that compared the callosal fiber composition as a function of the brain size of various mammals and concluded that large brains have more large axons and an increased maximal radius, whereas the bulk of axons is not altered (Olivares et al., 2001; Schüz and Preiβl, 1996; Liewald et al., 2014). Since the effective MR radius is predominantly sensitive to the larger axons, observed differences between humans and rats will be amplified when comparing effective MR radii. Overall, this observation favors future application of MR axon radius mapping in species with relatively large brain sizes.
Gray matter
Although this work mainly focuses on the WM, we do report significantly different signal scaling for the cortical GM. We suggest that the proton exchange between dendrites and interstitial water might explain this scaling behavior, in particular due to the convex scaling with ${b}^{1/2}$. However, the abundance of cell bodies in the gray matter might confound this analysis (Palombo et al., 2019). Moreover, the study of the cortical GM is challenged by its low SNR and susceptibility to partial voluming. Nonetheless, we conclude that the stick assumption does not hold in the cortical GM and that biophysical models building upon that assumption must be interpreted with caution if applied to tissue regions outside of WM.
Conclusion
In summary, we provide a realistic perspective on MR axon radius mapping by showing MRderived effective radii that have good quantitative agreement with histology. First, we compared the MRderived axon radii directly to confocal microscopy of the same rat brain samples. Second, the distribution of dMRIderived ${r}_{\text{MR}}$ of the living human brain falls almost entirely within the range spanned by histologyderived ${r}_{\text{eff}}$ that has been reported in the literature — even without introducing a putative axonal shrinkage factor. This estimation is inherently bound to a single scalar ${r}_{\text{eff}}$ that encodes moments of the axon distribution, which is – by virtue of the signal encoding – dominated by the largest axons. Therefore, the average axon radius $\overline{r}$ and ${r}_{\text{eff}}$ can be practically considered as two complementary metrics probing the underlying axon caliber distribution: histology, so far, mainly probes its bulk, that is $\overline{r}$, while dMRI probes ${r}_{\mathrm{MR}}={\widehat{r}}_{\mathrm{eff}}$, that is its tail. Due to the intrinsic bias of MRderived axon radii to larger axons, clinical applications should focus on pathologies that specifically target those larger axons, until other methods are developed that probe the smaller axon diameter.
Materials and methods
MRI of fixed rat brain tissue
Ethics
Animals used in this study were handled in agreement with the European FELASA guidelines and all procedures were approved by the Champalimaud Animal Welfare Body and by the national authorities, Direção Geral de Alimentação e Veterinária, Lisbon, Portugal, under the approved protocol number 0421/000/000/2016. All animal care procedures were conducted in agreement with the European Directive 2010/63, at the vivarium of the Champalimaud Foundation, a research facility part of CONGENTO, project number Lisboa01–0145FEDER022170.
Sample preparation
Request a detailed protocolThree Long Evans rats (Female, 12weeksold) were transcardially perfused using 4% paraformaldehyde. The extracted brains were kept for 24 hr in 4% paraformaldehyde and washed using PBS over two days (changed daily). Given our focus on the CC of the rat brains, we will refer to the samples as CC#1, CC#2, and CC#3.
MRI scanning
Request a detailed protocolMultishell dMRI data: The three samples were scanned on an 16.4T MR scanner (Bruker BioSpin) at room temperature with $\mathrm{\Delta}/\delta =\mathrm{\hspace{0.17em}20}/7.1\mathrm{ms}$ interfaced with an AVANCE IIIHD console and a micro2.5 imaging probe with maximal gradient amplitude ${G}_{\mathrm{max}}=\mathrm{\hspace{0.17em}1500}\mathrm{mT}/\mathrm{m}$. Diffusionweighting was applied using a RARE sequence in the midsagittal plane along 60 gradient directions for a densely sampled spectrum of 18 different $b$values up to 100 ms/µm^{2}. Furthermore, $\mathrm{TR}/{T}_{\mathrm{E}}=2400/30.4\mathrm{ms}$ and the spatial resolution was 100 × 100 × 850µm^{3}.
‘dot fraction’ ${f}_{\mathrm{im}}$: we acquired 60 repeated measurements of diffusionweighing applied in the direction parallel to the average fiber direction in the corpus callosum (CC) at the maximal $b$value of 100 ms/µm^{2}. The average fiber orientation was defined as the first eigenvector of the dyadic tensor (Jones, 2003) that was computed from the voxelwise first eigenvectors of the diffusion tensors that were estimated by fitting the DTI model to the lowest $b$values, that is $b<5\mathrm{ms}/\mu {\mathrm{m}}^{2}$, of the multishell data in each voxel within the manually segmented CC (Basser et al., 1994).
The average SNR for ${S}_{b=0}$ was 195 and some examples of the acquired images at various low and high $b$values are shown in Figure 10, where the quality of the raw data can be evaluated. Notably, since images are spherically averaged over many directions, the signal is characterized by high SNR even at high $b$values.
Data analysis
Request a detailed protocolFrom the multishell data, the sphericallyaveraged signals $\overline{S}(b)$ are estimated per $b$value as the zero^{th} order spherical harmonic using a Rician maximum likelihood estimator of the even order spherical harmonic coefficients up to the 6^{th} order (Sijbers et al., 1998).
The spatially localized dot fraction ${f}_{\mathrm{im}}$ is computed as the signal estimated from the repeated ($N=60$) measurements in the direction parallel to the principal fibre direction using a Rician maximum likelihood estimator with precomputed noise level (Veraart et al., 2016), normalized by the respective nondiffusion weighted signal. We compute the dotfree signal in each voxel as follows ${\overline{S}}_{\ast}(b)=\overline{S}(b){f}_{\mathrm{im}}$. In the remainder of the work, analyses were done on both ‘dot contaminated’ $\overline{S}(b)$ and ‘dot free’ ${\overline{S}}_{\ast}(b)$ signals.
The intraaxonal radial diffusivity ${\widehat{D}}_{a}^{\u27c2}$ is estimated voxelwise by fitting Equation 2 with ${f}_{\mathrm{im}}=0$ to ${\overline{S}}^{\ast}(b\ge 20)$ using a nonlinear least squares estimator (code is available for download on GitHub [Veraart and Novikov, 2019; https://github.com/NYUDiffusionMRI/AxonRadiusMapping; copy archived at https://github.com/elifesciencespublications/AxonRadiusMapping]).
Thereafter, the estimated effective axon radius ${r}_{\mathrm{MR}}$ is derived from ${\widehat{D}}_{a}^{\u27c2}$ using Equation 6. The alternative approach, that is the simultaneous estimation of ${\widehat{D}}_{a}^{\u27c2}$ and ${f}_{\mathrm{im}}$ from $\overline{S}(b\ge 20)$ is very poorly conditioned. Hence, disentangling both parameters from only the linearlyencoded multishell data is impossible, even at unrealistically high SNR.
Histology of fixed rat brain tissue
Full details of the immunohistochemistry for sample preparation, confocal microscopy, and image analysis are provided in Nunes et al. (2017). Studyspecific elements are described below.
Sample preparation
Request a detailed protocolAfter MRI scanning, freefloating horizontal sections 50 $\mu \mathrm{m}$thick were collected from the medial lateral center of two rat brains, CC#1 and CC#2, corresponding to the MR imaged volume. For CC#2, we used antibodies against neurofilaments 160/200 (axonal marker; SigmaAldrich, cat.# N2912) and GFAP (astrocytes marker; ThermoFisher Scientific, cat.# PA110019), as well as a staining for cell bodies using DAPI (SigmaAdrich, cat.# D9542). For CC#1, the staining was limited to the neurofilaments to focus on the axon radius count.
Confocal microscopy
Request a detailed protocolA Zeiss LSM 710 laser scanning confocal microscope was used for immunohistochemistry image acquisition. A tile scan using a 10× objective (EC Plan Neofluar, numerical aperture = 0.3, Zeiss, Germany) was used to cover the entire CC. Various ROIs were imaged using a 63× immersion objective (Plan Apochromat, numerical aperture = 1.4, Zeiss, Germany) in confocal mode, with pixel resolution of 65 × 65 × 150 nm^{3} and fieldofview of 135 × 135 µm^{2} (Figure 11). The placement of the ROIs is shown in Figures 4 and 5 for CC#1. and CC#2, respectively.
Confocal microscopy data analysis
Request a detailed protocolImages were processed using the ImageJ software. Noise suppression of the confocal single frames was done using a subsequent application of a 2D anisotropic diffusion filter and bandpass filtering in the frequency domain, (Nunes et al., 2017). Thereafter, axons were identified as particles with a minimum area size of 0.2 µm^{2} and a circularity larger than 0.4 in the confocal images that were stained for neurofilaments. The number of extracted axons varied from about 500 to 2000, depending on the placement of the ROI. The long axes of fitted ellipsoids served as proxies for the respective axon radii. For each ROI, we obtain a distribution of axon radii $h(r)$ from which we compute the associate effective radius ${r}_{\mathrm{eff}}$ using Equation 9.
In vivo MRI of human brain
Ethics
Data were acquired after obtaining written informed consent and consent to publish. The project was approved by the Cardiff University School of Psychology Ethics Committee (approval number EC.06.05.02.891).
Subjects
Four healthy volunteers (3 males and 1 female between 22 and 45 years) were recruited for this study. We will refer to the human subjects as 22M, 25F, 32M, and 45M to encode both the age and gender. The data were collected under the approval of the Cardiff University School of Psychology Ethics Committee.
MRI scanning
Request a detailed protocolAll four subjects were scanned on a Siemens Connectom 3T MR scanner using a 32channel receiver coil. The 300 mT/m gradient system was used to achieve $b$values up to 25$\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}/\mu {\mathrm{m}}^{2}$. The diffusion gradients were characterized by $\mathrm{\Delta}/\delta =30/13\mathrm{ms}$ and maximal gradient amplitude of 289 mT/m. Diffusion weighting was applied along 60 isotropically distributed gradient directions (Jones et al., 1999) for $b=$ 1, 3, 5, 7, 9, 11, 12.1, 13.5, 15, 16.9, 19.1, 21.7, and 25$\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}/\mu {\mathrm{m}}^{2}$, with TR/T_{E }: 3500/62$\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{s}$, matrix: 74 × 74, and 42 slices with a spatial resolution of 3 × 3 × 3$\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}{\mathrm{m}}^{3}$. The average SNR for ${S}_{b=0}$ was 52. See Figure 11 for the image quality and contrast at various $b$values.
The dot compartment was not measured directly (see Dhital et al., 2018 and Tax et al., 2019).
Data analysis
Request a detailed protocolImage processing was done according to the DESIGNER pipeline (AdesAron et al., 2018) using the FSL (Smith et al., 2004) and MRtrix (Tournier et al., 2019) software packages. In particular, MPPCA noise estimation and denoising (Veraart et al., 2016) were used for estimating noise maps $\sigma (x)$ by exploiting the inherent redundancy in dMRI data. The positive signal bias, inherent to lowSNR magnitude MR data, was removed by using the method of moments (Koay and Basser, 2006), where the denoised signal was used as a proxy for the Rician expectation value. Denoised and Ricefloorcorrected images were subsequently corrected for Gibbs ringing (Kellner et al., 2016), geometric eddy current distortions and subject motion (Andersson and Sotiropoulos, 2016). The pipeline is available on https://github.com/NYUDiffusionMRI/DESIGNER (AdesAron and Veraart, 2018). We used tractdensity imaging (Calamante et al., 2010) based on wholebrain probabilistic fibertracking (Tournier et al., 2019) of the $b=5\mathrm{ms}/\mu {\mathrm{m}}^{2}$shell for identifying all WM voxels. To avoid voxels affected by partial voluming with the gray matter (GM), an additional, more conservative, segmentation was obtained by omitting all voxels with a fractional anisotropy smaller than 0.6. In addition, the cortical GM was segmented using FreeSurfer (Dale et al., 1999).
Data availability
All source data files generated or analysed during this study have been deposited in Dryad Digital Repository (http://doi.org/10.5061/dryad.4qrfj6q66).

Dryad Digital RepositoryData from: Noninvasive quantication of axon radii using diffusion MRI.https://doi.org/10.5061/dryad.4qrfj6q66
References

Fiber composition of the human corpus callosumBrain Research 598:143–153.https://doi.org/10.1016/00068993(92)90178C

Axcaliber: a method for measuring axon diameter distribution from diffusion MRIMagnetic Resonance in Medicine 59:1347–1354.https://doi.org/10.1002/mrm.21577

MR diffusion tensor spectroscopy and imagingBiophysical Journal 66:259–267.https://doi.org/10.1016/S00063495(94)807751

The basis of anisotropic water diffusion in the nervous system  a technical reviewNMR in Biomedicine 15:435–455.https://doi.org/10.1002/nbm.782

Characterization and propagation of uncertainty in diffusionweighted MR imagingMagnetic Resonance in Medicine 50:1077–1088.https://doi.org/10.1002/mrm.10609

BookInformation and Likelihood Theory: A Basis for Model Selection and InferenceIn: Burnham K. P, Anderson D. R, editors. Model Selection and Multimodel Inference. Springer. pp. 49–97.

NMR microscopy of dynamic displacements: kspace and qspace imagingJournal of Physics E: Scientific Instruments 21:820–822.https://doi.org/10.1088/00223735/21/8/017

A simple matrix formalism for spin Echo analysis of restricted diffusion under generalized gradient waveformsJournal of Magnetic Resonance 129:74–84.https://doi.org/10.1006/jmre.1997.1233

The central role of mitochondria in axonal degeneration in multiple sclerosisMultiple Sclerosis Journal 20:1806–1813.https://doi.org/10.1177/1352458514544537

Optimizing gradient waveforms for microstructure sensitivity in diffusionweighted MRJournal of Magnetic Resonance 206:41–51.https://doi.org/10.1016/j.jmr.2010.05.017

PGSE, OGSE, and sensitivity to axon diameter in diffusion MRI: insight from a simulation studyMagnetic Resonance in Medicine 75:688–700.https://doi.org/10.1002/mrm.25631

The human connectome project's neuroimaging approachNature Neuroscience 19:1175–1187.https://doi.org/10.1038/nn.4361

In vivo correlation between axon diameter and conduction velocity in the human brainBrain Structure and Function 220:1777–1788.https://doi.org/10.1007/s0042901408710

Comments on the paper by Horowitz et al. (2014)Brain Structure and Function 220:1789–1790.https://doi.org/10.1007/s0042901409747

Measuring macroscopic brain connections in vivoNature Neuroscience 18:1546–1555.https://doi.org/10.1038/nn.4134

Quantification of cell size using temporal diffusion spectroscopyMagnetic Resonance in Medicine 75:1076–1085.https://doi.org/10.1002/mrm.25684

Determining and visualizing uncertainty in estimates of fiber orientation from diffusion tensor MRIMagnetic Resonance in Medicine 49:7–12.https://doi.org/10.1002/mrm.10331

Quantitative mapping of the peraxon diffusion coefficients in brain white matterMagnetic Resonance in Medicine 75:1752–1763.https://doi.org/10.1002/mrm.25734

NMR selfdiffusion studies in heterogeneous systemsAdvances in Colloid and Interface Science 23:129–148.https://doi.org/10.1016/00018686(85)80018X

Gibbsringing artifact removal based on local subvoxelshiftsMagnetic Resonance in Medicine 76:1574–1581.https://doi.org/10.1002/mrm.26054

Analytically exact correction scheme for signal extraction from noisy magnitude MR signalsJournal of Magnetic Resonance 179:317–322.https://doi.org/10.1016/j.jmr.2006.01.016

On the nature of the NAA diffusion attenuated MR signal in the central nervous systemMagnetic Resonance in Medicine 52:1052–1059.https://doi.org/10.1002/mrm.20260

Searching for the neurite density with diffusion MRI: challenges for biophysical modelingHuman Brain Mapping 40:2529–2545.https://doi.org/10.1002/hbm.24542

Looking into the functional architecture of the brain with diffusion MRINature Reviews Neuroscience 4:469–480.https://doi.org/10.1038/nrn1119

Dependence on bvalue of the directionaveraged diffusionweighted imaging signal in brainMagnetic Resonance Imaging 36:121–127.https://doi.org/10.1016/j.mri.2016.10.026

Early detection of regional cerebral ischemia in cats: comparison of diffusion and T2weighted MRI and spectroscopyMagnetic Resonance in Medicine 14:330–346.https://doi.org/10.1002/mrm.1910140218

Self‐diffusion coefficient of liquid lithiumThe Journal of Chemical Physics 48:4938–4945.https://doi.org/10.1063/1.1668160

Spin Echo of spins diffusing in a bounded mediumThe Journal of Chemical Physics 60:4508–4511.https://doi.org/10.1063/1.1680931

Mapping axonal density and average diameter using nonmonotonic timedependent gradientecho MRIJournal of Magnetic Resonance 277:117–130.https://doi.org/10.1016/j.jmr.2017.02.017

Species differences and similarities in the fine structure of the mammalian corpus callosumBrain, Behavior and Evolution 57:98–105.https://doi.org/10.1159/000047229

Pulsed NMR studies of restricted diffusion. I. droplet size distributions in emulsionsJournal of Colloid and Interface Science 40:206–218.https://doi.org/10.1016/00219797(72)900100

A theory of the effects of fibre size in medullated nerveThe Journal of Physiology 115:101–122.https://doi.org/10.1113/jphysiol.1951.sp004655

BookOligodendrocytesIn: Lyons D. A, Kegel L, editors. Magnetic Resonance Techniques for Imaging White Matter. Springer. pp. 397–407.

Basic connectivity of the cerebral cortex and some considerations on the corpus callosumNeuroscience & Biobehavioral Reviews 20:567–570.https://doi.org/10.1016/01497634(95)000690

Parametric probability distribution functions for axon diameters of corpus callosumFrontiers in Neuroanatomy 10:1–9.https://doi.org/10.3389/fnana.2016.00059

Aldehyde fixative solutions alter the water relaxation and diffusion properties of nervous tissueMagnetic Resonance in Medicine 62:26–34.https://doi.org/10.1002/mrm.21977

Maximumlikelihood estimation of rician distribution parametersIEEE Transactions on Medical Imaging 17:357–361.https://doi.org/10.1109/42.712125

An analytical model of restricted diffusion in bovine optic nerveMagnetic Resonance in Medicine 37:103–111.https://doi.org/10.1002/mrm.1910370115

Use of spin echoes in a pulsed magnetic‐field gradient to study anisotropic, restricted diffusion and flowThe Journal of Chemical Physics 43:3597–3603.https://doi.org/10.1063/1.1696526

AgeInduced white matter changes in the human brain: a stereological investigationNeurobiology of Aging 18:609–615.https://doi.org/10.1016/S01974580(97)001553

Self diffusion of water in frog muscleBiophysical Journal 28:107–116.https://doi.org/10.1016/S00063495(79)851620

Evaluation of restricted diffusion in cylinders. Phosphocreatine in rabbit leg muscleJournal of Magnetic Resonance, Series B 103:255–260.https://doi.org/10.1006/jmrb.1994.1038

Diffusion MRI noise mapping using random matrix theoryMagnetic Resonance in Medicine 76:1582–1593.https://doi.org/10.1002/mrm.26059

Axon radius mappingAxon radius mapping, https://github.com/NYUDiffusionMRI/AxonRadiusMapping.
Decision letter

Floris P de LangeSenior Editor; Radboud University, Netherlands

Birte ForstmannReviewing Editor; University of Amsterdam, Netherlands

Birte ForstmannReviewer; University of Amsterdam, Netherlands

Saad JbabdiReviewer; University of Oxford, United Kingdom

Robert MulkernReviewer
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This study focuses on a longstanding and important question in the field of diffusion MRI, namely the measurement of axon diameters. Importantly, until now, the accurate estimation of axon diameter mapping with noninvasive techniques such as diffusion MRI has been elusive due to a lack of sensitivity in the signal. The authors provide compelling evidence using sophisticated modeling that axon diameters can be estimated for largest axons when eliminating confounding factors such as extraaxonal water and axonal orientation dispersion. Data of fixed rat brains and optical microscopy of the same specimen are presented showing good quantitative agreement for MRderived axon diameters. Finally, in vivo data from Connectom 3T scanning is presented which shows the feasibility of mapping axon diameters in healthy subjects. The work is therefore of interest to a broad scientific audience ranging from physicists to cognitive neuroscientists.
Decision letter after peer review:
Thank you for submitting your article "Noninvasive quantification of axon diameter using diffusion MRI' for consideration by eLife. Your article has been reviewed by three peer reviewers, including Birte Forstmann as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen Floris de Lange as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
Overall, this manuscript is well written, interesting, timely and will help resolve the debate in the field. We therefore suggest revising the manuscript to address the points raised by the reviewers which are outline below.
Reviewer 1:
This paper introduces an interesting and surprisingly simple method for estimating axon diameter (potentially in vivo). Their approach relies on three key ingredients:
– ensure that the diffusion sensitisation (b value) is high enough to eliminate extraaxonal water
– use powder averaging across diffusion gradient orientations to eliminate confounds due to orientation of axons
– estimate intraaxonal perpendicular diffusivity from a simple powerlaw formula predicted by the theory
The intraaxonal perpendicular diffusivity is then used to estimate an "effective axonal radius" which describes the tail of the axon radius distribution under certain assumptions.
In general, this is a potentially important contribution. I have the following comments which I am sure the authors will be able to address:
1) Presentation
I found the presentation of the theory unnecessarily dense and difficult to follow. A simple diagram might help. Something like Figure 1 but with a single line showing the prediction from a model with f_{im} as well as D_{perp}>0 (so with nonzero intercept of the nonlinear part and a negative intercept of the linear part). It would also be helpful to indicate on the same diagram that the slope on the righthand side depends on D_{para}, and that the intercept is f_{im}, but the 'dashed line' intercept is a function of f_{im} and D_{perp} (and D_{para}?). The xaxis could be doublelabelled with both $1/\sqrt{b}$ and b. And also it would help to have vertical sections in the graph indicating the b value regimes (e.g. clinicallyfeasible, vs. Connectom vs. smallbore scanners, vs. lowb regime).
2) Modelling
Having played with the models (i)(viii) an little, I see lots of degeneracy between f_{im} and D_{perp} over a wide range of parameter values which is not surprising: unless some curvature is visible in the data (above the noise level), it is difficult to disentangle the contributions of f_{im} and D_{perp} to the negative intercept. In the data that is shown (e.g. Figure 5), the points fall in a straight line and so there is no curvature to help disentangle f_{im} and D_{perp}. The authors assumed that f_{im}=0 for the in vivo data, but is this really justified (could there not still be very slow diffusing water that is unmodelled) and what about ex vivo?
In general, some sort of analysis of when the degeneracy breaks down (as a function of the max b value attainable and the other params like D_{para}) would be helpful here. For example looking at the full posterior distribution and not just point estimates of the parameters (I don't find AIC values very helpful compared to looking at the full posterior distribution).
3) Axon diameter estimation
It would be helpful if the authors could unpack how they get to r_{eff} as a function of D_{perp}. I can see that the diffusion in a cylinder formula gives rise to a r^{4} dependence of the log(S) in the regime that the acquisition are made in. But then to go from there to r_{eff} (which is the ratio of 6th to 2nd moment) is a stretch for me. Is it simply by doing a Taylor expansion of S=exp(a*r^{4}) around zero inside the integral in Equation 7? If that is the case then perhaps an appendix would not hurt. Also, it is not clear how accurate it is to use the Taylor expansion.
Also on axon diameters, the authors make it quite clear that they don't like methods that make explicit distributional assumptions on the axon diameter (e.g. AxCaliber) – but I think it would be interesting to compare them with the author's approach. Looking at the histology data that is presented, one wonders how accurate a gamma distribution would be. With a distributional assumption there would be no need for the Taylor expansion above, and everything can be done keeping the exponentials and directly inverting the equations to get the parameters of the gamma distribution. How does this compare to the r_{eff} proposed by the authors?
4) Are the results biologically sound?
Generally, I found that there was not enough in terms of showing results that indicate the technique actually works well, e.g. in Figure 8: is there a way to show a similar map from histology?. Or in general show that interareas variation in radius from histo correlate with interarea variation from MR in the ex vivo data. The only comparison between the two modalities is done in Figure 7 in a single region.
Similarly for the in vivo data: Figure 10: is there any evidence that the intervoxel variation is meaningful? There are zones of reduced radius in the lateral frontal lobes near the cortex – are those meaningful?
5) Exchange
The authors present an interesting extra source of information, in that an exchange model makes a different prediction at high b values and they found evidence for exchange in GM. Can you convince the reader that this is not just partial volume effects (e.g. multiple pools of water with different diffusion coefficients and no exchange?) would that for example induce a curve with the opposite convexity? As GM is likely to have more partial volume issues I think this is a valid question that needs addressing.
Also on exchange: My understanding of the Karger model is that it assumes that exchange happens in situ (a molecule would change its behaviour from e.g. slow to fast diffusion with some probability instantaneously). But in reality exchange happens at the membrane. Does that invalidate the equations? Can the equations be derived here in an appendix?
6) Presentation of the data
Single voxel data is never shown and so it is difficult to tell how noisy the signal vs. 1/sqrt(b) curves actually are.
7) Data sharing
The authors are to be commended on sharing their data. However the way they have done it is not optimal in that they only provide raw data with no particular documentation or curation. The shared data set would strongly benefit if you would add the following:
– include preprocessed data not just raw data (including the D_{perp} and r_{eff} maps) – or at least provide code to generate the maps and do the preprocessing
– match data format between human and rodent
– include documentation
– avoid lsm format as it is proprietary – maybe use tiff instead?
– include processed histo data?
Reviewer 2:
This study focuses on a longstanding and important question in the field of diffusion MRI, namely the measurement of axon diameters. Importantly, until now, the accurate estimation of axon diameter mapping with noninvasive techniques such as diffusion MRI has been elusive due to a lack of sensitivity in the signal. The authors provide compelling evidence using sophisticated modeling that axon diameters can be estimated for largest axons when eliminating confounding factors such as extraaxonal water and axonal orientation dispersion. Data of fixed rat brains and optical microscopy of the same specimen are presented showing good quantitative agreement for MRderived axon diameters. Finally, in vivo data from Connectom 3T scanning is presented which shows the feasibility of mapping axon diameters in healthy subjects.
Major comments:
Generalizability of the data:
My main concern is that the MRIbased axon diameter modeling was only evaluated in the corpus callosum. It would be important to see whether the modeling also holds in other fiber tracts, e.g., frontooccipital fasciculus.
This is something that the authors should ideally address, but in case this is not feasible, at least comment on.
Reviewer 3:
This is an impressive work combining wellthought out theory with experimental data only recently available, particularly for the human studies, using the Connectom system to provide gradient strengths some 4 times larger than available on commercial scanners. The mix of preclinical data with rat CC for which histological distributions of axon diameters was measured, with human data (using somewhat less gradient strength than available on the animal system) and only literature histology is justified and actually adds strength to the comparison of experimental with theoretical considerations. The lack of any attempt to measure the "dot" component in humans is less justifiable in my view though that might have significantly added to the scan time and further comments on this might be appreciated. The authors recognize the limitations of their assessment in having to rely upon a rather "weighted" version of the distribution which gives an output index well into the tail of the distribution, the larger axons, but at least the measures are getting closer to the actual size of the median axon values than those reported in the past with more standard gradient strengths and perhaps dubious modeling. It also would be helpful to perhaps add to Figure 1 or another figure the curves that would be anticipated theoretically from the exchange model of Equation 4 at such high b values, emphasizing the difference between concave and convex theoretical curves that the authors, I assume, deem to eliminate the latter model given the experimental data. Finally as a major point, in the Data Analysis section the authors explain r_{eff} or r_{MR} from the data but this description is difficult to follow. For example, in Equation 2, how are the O(b^{2}) taken into account if they are. Then, assume we now have Da(perpendicular) how does one use that with Equations 5 and 9 to get r_{eff}. This must be clarified. People should be able to replicate this calculation from what is in this text.
https://doi.org/10.7554/eLife.49855.sa1Author response
Reviewer 1:
[…]
In general, this is a potentially important contribution. I have the following comments which I am sure the authors will be able to address:
1) Presentation
I found the presentation of the theory unnecessarily dense and difficult to follow. A simple diagram might help. Something like Figure 1 but with a single line showing the prediction from a model with f_{im} as well as D_{perp}>0 (so with nonzero intercept of the nonlinear part and a negative intercept of the linear part). It would also be helpful to indicate on the same diagram that the slope on the righthand side depends on D_{para}, and that the intercept is f_{im}, but the 'dashed line' intercept is a function of f_{im} and D_{perp} (and D_{para}?). The xaxis could be doublelabelled with both 1/sqrt(b) and b. And also it would help to have vertical sections in the graph indicating the b value regimes (e.g. clinicallyfeasible, vs. Connectom vs. smallbore scanners, vs. lowb regime).
The theoretical sections, mainly “Breaking the power law” has been revised to improve readability. Following the suggestion of the reviewer, we amended Figure 1 to visualize the difference between f_{im}and the intercept. While f_{im}is a signal fraction of full restricted isotropic diffusion, the intercept is a convoluted metric that combines f_{im}with a b value dependent offset that relates to the sensitivity of MR to the radius.
2) Modelling
Having played with the models (i)(viii) an little, I see lots of degeneracy between f_{im} and D_{perp} over a wide range of parameter values which is not surprising: unless some curvature is visible in the data (above the noise level), it is difficult to disentangle the contributions of f_{im} and D_{perp} to the negative intercept. In the data that is shown (e.g. Figure 5), the points fall in a straight line and so there is no curvature to help disentangle f_{im} and D_{perp}. The authors assumed that f_{im}=0 for the in vivo data, but is this really justified (could there not still be very slow diffusing water that is unmodelled) and what about ex vivo?
In general, some sort of analysis of when the degeneracy breaks down (as a function of the max b value attainable and the other params like D_{para}) would be helpful here. For example looking at the full posterior distribution and not just point estimates of the parameters (I don't find AIC values very helpful compared to looking at the full posterior distribution).
The degeneracy is an intrinsic limitation of multicompartmental modeling of diffusion MRI data. Various strategies to resolve this problem have been attempted during the past decade, ranging from imposing constraints or priors, to complement the data with orthogonal measurements. The former strategy has been contested because of the potential biases that might arise from inaccurate priors, whereas the exploring novel acquisition strategies to resolve the degeneracies is currently widely studied. A promising avenue is complementing the classical StjeskalTanner diffusionweighting – as used in this study – with planar and spherical diffusionweighted strategies. Whereas linear encoding is bestsuited to probe elongated cellular structures, such as axons, the spherical encoding is most sensitive to spherical objects such as cell bodies. Tax et al., 2019, and Dhital et al., 2018, used this spherical encoding to quantify the signal fraction f_{im}in healthy white matter. In previous work (Veraart et al., 2019) we made a similar observation, i.e. f_{im}<0.2%, using linear diffusion encoding at high b when only considering the parallel diffusion directions. Here, we adopt the conclusive result of f_{im}not being significant for the human data in the major white matter structures.
Tissue fixation and possibly temperatureinduced alterations of the microstructure has resulted in a nonzero f_{im}in ex vivo studies. To avoid the need for the simultaneous estimation of f_{im}and D_{perp}, we included a dedicated experiment to quantify f_{im}prior to the axon diameter mapping. We recommend a similar strategy when studying clinical cohorts in future studies.
The optimization landscape of Equation 2 (shown in updated Figure 1) reveals that disentangling f_{im}and D_{perp}is impossible, even at unrealistically high SNR. Indeed, the contrast in error function along a valley that runs through the landscape is minimal in comparison to the noise floor, see Author response image 1. In the revised manuscript, we highlight the intrinsic degeneracy of axon diameter mapping with unknown dot compartment explicitly. We dedicate a more extensive section to the dot compartment in which we discuss the degeneracy, overinterpretation of AIC analysis, and need for measurement of the dot compartment in atypical cohorts.
3) Axon diameter estimation
It would be helpful if the authors could unpack how they get to reff as a function of Dperp. I can see that the diffusion in a cylinder formula gives rise to a r4 dependence of the log(S) in the regime that the acquisition are made in. But then to go from there to reff (which is the ratio of 6th to 2nd moment) is a stretch for me. Is it simply by doing a Taylor expansion of S=exp(a*r4) around zero inside the integral in Equation 7? If that is the case then perhaps an appendix would not hurt. Also, it is not clear how accurate it is to use the Taylor expansion.
r_{eff} is derived from the Taylor expansion of <r^{2} exp(a*r^{4})>/<r^{2}>. We revised section “From D_{perp} to effective MR radius” to guide the reader through the mathematics behind the effective MR radius more rigorously. The accuracy of the Taylor expansion, and the model in general, for realistic settings has been addressed in a simulation section, Figure 4.
Also on axon diameters, the authors make it quite clear that they don't like methods that make explicit distributional assumptions on the axon diameter (e.g. AxCaliber) – but I think it would be interesting to compare them with the author's approach. Looking at the histology data that is presented, one wonders how accurate a Gamma distribution would be. With a distributional assumption there would be no need for the Taylor expansion above, and everything can be done keeping the exponentials and directly inverting the equations to get the parameters of the Gamma distribution. How does this compare to the reff proposed by the authors?
Sepehrband et al., 2016, published a comprehensive study on the accuracy of various parametric distribution to describe the axon distribution in the mouse corpus callosum. The generalized extreme valuedistribution consistently fitted the observed distributions better than other distribution functions, including the Gamma distribution. Most importantly, wellfitting distributions are parametrized by two or more parameters, so trying reconstructing the parametric distribution from the MR effective radius is illposed if these parameters are to be estimated from the same data. We simply don’t have enough information to estimate all parameters of such a wellfitting distribution. Instead of making additional assumption that will further bias the axon diameter estimates, we encourage future research directions in which the current approach is complemented with oscillating or short gradient pulse experiments to decode additional, lowerorder, moments of the distribution. This might provide an avenue to reconstruct the parametric distribution based on MRI, but current hardware limits prevent us from conducting this experiment.
4) Are the results biologically sound?
Generally, I found that there was not enough in terms of showing results that indicate the technique actually works well, e.g. in Figure 8: is there a way to show a similar map from histology?. Or in general show that interareas variation in radius from histo correlate with interarea variation from MR in the ex vivo data. The only comparison between the two modalities is done in Figure 7 in a single region.
Similarly for the in vivo data: Figure 10: is there any evidence that the intervoxel variation is meaningful? There are zones of reduced radius in the lateral frontal lobes near the cortex – are those meaningful?
The direct comparison between MRI and histology is done in 20 different patches, covering 4 different locations within the corpus callosum and two samples. For all those samples, we reported the accuracy, demonstrating good quantitative agreement, despite a small remaining overestimation. The potential confounding factors that might have contributed to this observation, yet we were unable to eliminate, are discussed in depth in the Discussion. Aside from this direct comparison, we do rely on a qualitative and indirect comparison for other brain structures, including the human corpus callosum. The trend of axon sizes within the corpus callosum, as shown in Figure 8, has been reported in various species, including the rat (Barazany et al., 2009 and Sargon et al., 2003), rhesus monkey (Lamantia et al., 1990), and human (e.g. Aboitiz et al., 1992). Although we don’t emphasize it in the manuscript, one might even argue that the slightly increased radii in the rostrum of the corpus callosum are in qualitative agreement with Sargon et al., 2003.
A thorough comparison of our MR results with the histological literature is challenged by the dependency on the axon diameter distribution to calculate the corresponding effective radius, as discussed in issue #4 above. The corpus callosum is best characterized in that respect and, as such, our validation component is limited to that commissural fiber. Please consider Author response image 2 in which we show the same trend in our human data.
The lack of such data in various regions across the human brain, especially in regions such as the frontal lobes, prevents us from making strong claims on meaning. Instead, we would claim that, now, we can start relying on axon diameter mapping using MRI to explore inter and alongtract variability in the entire living brain.
5) Exchange
The authors present an interesting extra source of information, in that an exchange model makes a different prediction at high b values and they found evidence for exchange in GM. Can you convince the reader that this is not just partial volume effects (e.g. multiple pools of water with different diffusion coefficients and no exchange?) would that for example induce a curve with the opposite convexity? As GM is likely to have more partial volume issues I think this is a valid question that needs addressing.
We don’t want to make the claim that exchange is the only possible explanation for the different signal decay, nor can we convince the reader that it is the most likely explanation. Palombo et el., 2019, recently demonstrated that the presence of a (non)exchanging compartments such as somas might confound our interpretation. However, the functional form of the exchange with a “stick” compartment is novel, and we believe that bringing it up as a viable possibility has value. Especially because the convex scaling is not expected from a partialvolume contribution of a nonexchanging nonstick compartment, such as a spherical cell. By reporting the signal decay in the gray matter and posing the hypothesis of exchange, we would like to thereby trigger further research to explore different relevant biophysical processes in the gray matter, in order to develop the most adequate models of diffusion in GM. Hence, as stated in the Discussion, we acknowledge that exchange is just one out of a few different avenues to be explored and investigated.
Also on exchange: My understanding of the Karger model is that it assumes that exchange happens in situ (a molecule would change its behaviour from e.g. slow to fast diffusion with some probability instantaneously). But in reality exchange happens at the membrane. Does that invalidate the equations? Can the equations be derived here in an appendix?
As clarified by Fieremans et al., 2010, the Karger model (KM) is applicable for diffusion times t long enough so that the medium, coarsegrainedover the corresponding diffusion length L(t), can be viewed as one where exchange happens everywhere, as the reviewer has indeed noted.
For the intraneurite water, this means that the exchange should be slow enough, so that it is “barrierlimited” rather than “diffusionlimited” (per Karger’s original terminology). Our estimates of the residence times exceeding 10ms ensure that this is the case, as most neurite diameters are under 1 micron and the corresponding diffusion time is therefore an order of magnitude slower. These estimates are in line with other measurements, e.g., by Yang et al., 2018, whose exchange times exceeded 100 ms. For the extraneurite water, technically, the KM applicability requires L(t) ≫ l_{c},or, equivalently, t ≫ t_{c},where l_{c}is the correlation length of the structure. Practically, the same requirements ensure that the tortuosity limit has been reached in each compartment separately, i.e., diffusion in each compartment has become Gaussian (neglecting the exchange); this Gaussian diffusion is what is assumed by the Karger model’s equations (our only difference is that we made this Gaussian diffusion anisotropic, as compared to isotropic diffusion originally considered by Karger). For our diffusion time, the diffusion length in the plane transverse to the neurites, Lt=2∙2Det≈2∙2∙130≈11μmexceeds the mean distance between them; it also exceeds the distance between beads along neurites (37 microns). The exact tissue parameters remain controversial while the EM segmentation and tissue quantification is underway. However, the beauty of employing very strong diffusion weighting is that the extraneurite compartment is suppressed, therefore the above issues become irrelevant.
All in all, we believe KM may be asymptotically applicable here, in the view of suppressing the extraneurite signal at strong b, as well as because the estimated exchange time a posteriorijustifies the barrierlimited exchange for the intraneurite water. However, as KM may still be too simplistic, while giving an approximate range, we are not making strong claims about the numerical value of the exchange time in gray matter, because it is an entirely different active area of investigation, and the relevance of different biophysical effects there is still under intense debate. Our relatively short estimated exchange time (on the scale of clinically used diffusion times) may indeed explain the difference in the functional form between high b signal in GM from a “stick” seen in WM, and adds value to the current debate about the need to include exchange in GM modeling.
6) Presentation of the data
Single voxel data is never shown and so it is difficult to tell how noisy the signal vs 1/sqrt(b) curves actually are.
Please consider Author response image 3 in which we show the sphericallyaveraged signal decay as a function of 1/sqrt(b) for all individual voxels of the WM of one human subject (faded colored lines) and for the average across all WM voxels (blue line with markers, cf. Figure 5B). In addition, to highlight the precision, we show the signal decay in 10 arbitrarily chosen individual WM voxels. We would like to highlight that Figures 8 and 10 in the manuscript show voxelbyvoxel results that demonstrate the precision qualitatively.
7) Data sharing
The authors are to be commended on sharing their data. However the way they have done it is not optimal in that they only provide raw data with no particular documentation or curation. The shared data set would strongly benefit if you would add the following:
– include preprocessed data not just raw data (including the D_{perp} and r_{eff} maps) – or at least provide code to generate the maps and do the preprocessing
– match data format between human and rodent
– include documentation
– avoid lsm format as it is proprietary – maybe use tiff instead?
– include processed histo data?
Since there is no consensus within the neuroimaging community on how to process dMRI data, especially with a strong diffusionweighting, we opt to share the raw data, for both MRI and microscopy, and not limit the users to the considerations and choices that we made. We would like to highlight that all image processing was done using publicly available software tools, i.e. FSL, Freesurfer, MRtrix, and ImageJ – as listed in the Key resources table. The code for the axon diameter fitting is also made available on our GitHub.
We now match the data format, NiFTi, for humans and rodents; provide the microscopic images in LSM and TIFF format, and include a more comprehensive documentation of the data in the Dryad Digital Repository.
Reviewer 2:
[…]
Major comments:
Generalizability of the data:
My main concern is that the MRIbased axon diameter modeling was only evaluated in the corpus callosum. It would be important to see whether the modeling also holds in other fiber tracts, e.g., frontooccipital fasciculus.
This is something that the authors should ideally address, but in case this is not feasible, at least comment on.
Our work focuses on theory, validation, and human feasibility. To demonstrate the human feasibility, we opted to limit ourselves to the corpus callosum because it’s the tract which has been characterized most thoroughly across multiple histological studies. Nonetheless, in Author response image 4, we show that the intertract variability of the effective MR axon diameter is higher than the intersubject variability to provide on optimistic outlook on whole brain characterization in future studies. A comment is made in the revised manuscript.
Reviewer 3:
[…]
The lack of any attempt to measure the "dot" component in humans is less justifiable in my view though that might have significantly added to the scan time and further comments on this might be appreciated.
The justification of a negligible dot compartment in the healthy white matter was found in various recent publications, Dhital et al., 2018, Veraart et al., 2019, and Tax et al., 2019, which concluded independently from each other, using different techniques, that such a compartment is not significant in the major white matter pathways in the healthy adult brain. Note that the study of Tax et al. was performed on exactly the same scanner. That being said, when moving towards clinical or preclinical applications, we encourage the independent measurement of the dot compartment to complement to axon diameter acquisitions. The fast measurement of the dot compartment is promoted by the availability of spherical diffusionencoding. We do comment more extensively on this discussion point in the revised manuscript.
It also would be helpful to perhaps add to Figure 1 or another figure the curves that would be anticipated theoretically from the exchange model of Equation 4 at such high b values, emphasizing the difference between concave and convex theoretical curves that the authors, I assume, deem to eliminate the latter model given the experimental data.
The curve is now shown in Figure 1 to showcase the significantly different signal decay.
Finally as a major point, in the Data Analysis section the authors explain r_{eff} or rMR from the data but this description is difficult to follow. For example, in Equation 2, how are the O(b^{2}) taken into account if they are. Then, assume we now have Da(perpendicular) how does one use that with Equations 5 and 9 to get r_{eff}. This must be clarified. People should be able to replicate this calculation from what is in this text.
This comment is in line with comments 1 and 3 from reviewer #1. We include the derivation of r_{eff} using Equations 5 to 8 explicitly. The error associated with modelling approximations, e.g. omitting O(b^{2}) or the Taylor expansion in Equation 9, is estimated using a realistic simulation framework, as shown in Figure 4.
https://doi.org/10.7554/eLife.49855.sa2Article and author information
Author details
Funding
Fonds Wetenschappelijk Onderzoek (12S1615N)
 Jelle Veraart
National Institute of Neurological Disorders and Stroke (R01NS088040)
 Els Fieremans
 Dmitry S Novikov
National Institute of Biomedical Imaging and Bioengineering (P41 EB017183)
 Els Fieremans
 Dmitry S Novikov
H2020 European Research Council (679058)
 Noam Shemesh
 Daniel Nunes
Engineering and Physical Sciences Research Council (EP/M029778/1)
 Derek K Jones
Wellcome (096646/Z/11/Z)
 Derek K Jones
Wellcome (104943/Z/14/Z)
 Derek K Jones
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
All authors would like to thank Prof. Mark D Does (Vanderbilt University) for the remmiRARE sequence used in this study, supported by grant number NIH EB019980, and Dr. Erika Raven for discussions and comments on the manuscript. JV is a Postdoctoral Fellow of the Research Foundation  Flanders (FWO; grant number 12S1615N). DN and NS are supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Starting Grant, agreement No. 679058). EF and DSN were supported by the NIH/NINDS award R01NS088040 and research was performed as part of the Center of Advanced Imaging Innovation and Research (CAI2R, www.cai2r.net), an NIBIB Biomedical Technology Resource Center (NIH P41 EB017183). The Connectom data were acquired at the UK National Facility for in vivo MR Imaging of Human Tissue Microstructure funded by the EPSRC (grant EP/M029778/1), and The Wolfson Foundation. DKJ is supported by a Wellcome Trust Investigator Award (096646/Z/11/Z) and a Wellcome Trust Strategic Award (104943/Z/14/Z).
Ethics
Human subjects: Human studies were carried out under a protocol (EC.06.05.02.891) approved by Cardiff University School of Psychology Ethics Committee. Written informed consent and consent to publish was obtained from all participants.
Animal experimentation: Animals used in this study were handled in agreement with the European FELASA guidelines and all procedures were approved by the Champalimaud Animal Welfare Body and by the national authorities, Direção Geral de Alimentação e Veterinária, Lisbon, Portugal, under the approved protocol number 0421/000/000/2016. All animal care procedures were conducted in agreement with the European Directive 2010/63, at the vivarium of the Champalimaud Foundation, a research facility part of CONGENTO, project number Lisboa010145FEDER022170.
Senior Editor
 Floris P de Lange, Radboud University, Netherlands
Reviewing Editor
 Birte Forstmann, University of Amsterdam, Netherlands
Reviewers
 Birte Forstmann, University of Amsterdam, Netherlands
 Saad Jbabdi, University of Oxford, United Kingdom
 Robert Mulkern
Publication history
 Received: July 2, 2019
 Accepted: January 7, 2020
 Version of Record published: February 12, 2020 (version 1)
 Version of Record updated: February 13, 2020 (version 2)
Copyright
© 2020, Veraart et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 4,430
 Page views

 576
 Downloads

 70
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
The automatic initiation of actions can be highly functional. But occasionally these actions cannot be withheld and are released at inappropriate times, impulsively. Striatal activity has been shown to participate in the timing of action sequence initiation and it has been linked to impulsivity. Using a selfinitiated task, we trained adult male rats to withhold a rewarded action sequence until a waiting time interval has elapsed. By analyzing neuronal activity we show that the striatal response preceding the initiation of the learned sequence is strongly modulated by the time subjects wait before eliciting the sequence. Interestingly, the modulation is steeper in adolescent rats, which show a strong prevalence of impulsive responses compared to adults. We hypothesize this anticipatory striatal activity reflects the animals’ subjective reward expectation, based on the elapsed waiting time, while the steeper waiting modulation in adolescence reflects agerelated differences in temporal discounting, internal urgency states, or explore–exploit balance.

 Computational and Systems Biology
 Neuroscience
How dynamic interactions between nervous system regions in mammals performs online motor control remains an unsolved problem. In this paper, we show that feedback control is a simple, yet powerful way to understand the neural dynamics of sensorimotor control. We make our case using a minimal model comprising spinal cord, sensory and motor cortex, coupled by long connections that are plastic. It succeeds in learning how to perform reaching movements of a planar arm with 6 muscles in several directions from scratch. The model satisfies biological plausibility constraints, like neural implementation, transmission delays, local synaptic learning and continuous online learning. Using differential Hebbian plasticity the model can go from motor babbling to reaching arbitrary targets in less than 10 min of in silico time. Moreover, independently of the learning mechanism, properly configured feedback control has many emergent properties: neural populations in motor cortex show directional tuning and oscillatory dynamics, the spinal cord creates convergent force fields that add linearly, and movements are ataxic (as in a motor system without a cerebellum).