Graphical Abstract. Surface-based relaxation model to predict axon radius.

Highlights

- Diffusion-relaxation MRI data were acquired using a high b-value acquisition.

- A diffusion-relaxation model to estimate the intra-axonal T2 and T1 was proposed.

- The histological inner axon radius modulated the estimated relaxation times.

- A surface-based relaxation model predicted the axon radius in the corpus callosum.

- The predicted axon radii agreed with the mean effective histological radius.

1. Introduction

The speed of action potentials along axons is partly determined by their radii (Goldstein and Rall, 1974). Axon radius explains the biggest variance in conduction speed, as demonstrated by previous studies (Hursh, 1939), with larger axons conducting faster than those with smaller radii (Waxman and Bennett, 1972; Costa et al., 2018; Drakesmith et al., 2019). Therefore, accurately measuring axon radii in vivo is essential for better understanding the neural mechanisms underlying brain function and their impact on diseases.

The diffusion Magnetic Resonance Imaging (dMRI) signal is sensitive to axon radii if strong diffusion encoding gradients (i.e., up to 300 mT/m in Connectom scanners (Jones et al., 2018) and 1,500 mT/m in animal preclinical scanners) are used (Assaf et al., 2004, 2008; Assaf and Basser, 2005; Alexander, 2008; Dyrby et al., 2013; Duval et al., 2015; De Santis et al., 2016; Veraart et al., 2020; Barakovic et al., 2021a). However, the main limitation of this approach is that the dMRI signals from axons with radii smaller than ∼1–2 μm are practically indistinguishable from each other, even when the most advanced human Connectom scanners with ultra-strong (300 mT/m) gradients are employed in the data acquisition (Nilsson et al., 2017). Today, the challenge is that the peak of the axon radius distribution per voxel is below one micrometer in most brain regions, as observed in histology. Hence, most axon radii are below the lower bound for detection (Edgar and Griffiths, 2014; Dyrby et al., 2018). For an overview of the different strategies that have been employed to measure axon radius with dMRI, the reader is referred to Assaf and Basser (2005), Assaf et al. (2008, 2013), Alexander et al. (2010, 2019) Dyrby et al. (2013, 2018), Novikov et al. (2019), Fan et al. (2020), Jelescu et al. (2020), Veraart et al. (2020), Barakovic et al. (2021a), Pizzolato et al. (2023).

Theoretical reasons explain the lower sensitivity of dMRI to the inner radius of smaller axons. The commonly employed model [i.e., Gaussian phase approximation in the long-pulse limit (van Gelderen et al., 1994)] predicts an intra-axonal dMRI signal attenuation that depends on the fourth power of the radius r. Moreover, since the measured intra-axonal signal per voxel is the sum of all the individual intra-axonal signals weighted by each axon’s contribution to the signal (scaling by an extra-factor r2), larger axons contribute more than smaller axons to the measured signal. After considering these two factors together, an approximate expression for the mean “effective” dMRI-based radius reff per voxel can be derived, which depends on the higher-order moments of the unknown axon radius distribution. The resulting analytical expression re⁢f⁢f≈(⟨r6⟩⁢/⁢⟨r2⟩)14 (where ⟨⟩ denotes the average over the distribution) demonstrates that the estimate is heavily weighted by the right-hand tail of the axon radius distribution (Burcaw et al., 2015; Veraart et al., 2020). Consequently, the estimated mean axon radius is mainly affected by the bigger axons from the fractions of axons larger than the lower bound. This explains why estimations may appear overestimated compared to histology (Alexander et al., 2010; Dyrby et al., 2018).

Finding another source of MRI contrast sensitive to the size of axons smaller than the diffusion resolution limit is essential. Various studies in porous media have demonstrated that the interaction between the water molecules and the confining pore surface reduces the observed transverse T2 relaxation time (Brownstein and Tarr, 1977). This surface-based relaxation mechanism allows pore size to be estimated (Hurlimann et al., 1994; Slijkerman and Hofman, 1998; Sørland et al., 2007; Mohnke and Hughes, 2014; Müller-Petke et al., 2015). Notably, a similar T2 relaxation model to predict the size of cells was proposed previously (Brownstein and Tarr, 1979), and the idea of applying it to estimate the axon radius was suggested by Kaden and Alexander (2013). However, there is a lack of validation studies to demonstrate whether the inner axon radius modulates the intra-axonal relaxation times. This might be explained by the fact that approaches to estimating the intra-axonal relaxation times have only been developed recently (Veraart et al., 2018; McKinnon and Jensen, 2019; Barakovic et al., 2021b; Tax et al., 2021; Pizzolato et al., 2022). Furthermore, to our knowledge, no dataset is available that offers the combined histological information and relaxometry MRI data from the same sample, which are necessary for the estimation and comparison of these parameters.

The dMRI signals arising from the intra-axonal space can be isolated if a sufficiently high b-value is employed (i.e., b > 4,000 s/mm2 for in vivo data), which significantly attenuates the signal from spins experiencing large displacements (Jensen et al., 2016; McKinnon and Jensen, 2019). As the confining axonal geometry restricts the self-diffusion motion of spins inside axons (assuming a slow exchange between the intra- and extra-axonal spaces), the strongly diffusion-weighted MRI signal should come from the intra-axonal spins. Thus, it is possible to fit a diffusion-relaxation model of intra-axonal relaxation to strongly diffusion-weighted MRI data collected at multiple diffusion gradient directions and different echo times. This approach, combined with taking the spherical mean (orientational average), was employed previously to estimate the mean intra-axonal T2 time per voxel (McKinnon and Jensen, 2019) and bundle (Barakovic et al., 2021b), unconfounded by fiber orientation effects.

This proof of concept study investigates whether the intra-axonal T2 and T1 relaxation times are related to the inner axon radius and whether they can be employed to predict the mean effective radius. To do this, (1) we implemented two acquisition protocols and measured diffusion-T1-T2 and diffusion-T2 weighted MRI data from three healthy volunteers, one of them scanned using both sequences; (2) we employed a diffusion-relaxation model to enable the estimation of both intra-axonal T2 and T1 relaxation times by using the spherical mean signals from the acquired data; (3) we fitted the estimated relaxation times to a surface-based relaxation model that depends on the histological axon radius; (4) using histology from some brain regions we calibrated the surface-based relaxation model to enable predicting axon radius in other brain regions, and (5) we compared the MRI-based estimated axonal radii with those obtained from two postmortem histological human brain datasets in several regions in the midsagittal Corpus Callosum (CC) cross-section. Additional details are provided at the end of the next section.

2. Theory 2.1. Surface-based relaxation model

Inspired by the standard surface-based relaxation model used in porous media (Zimmerman and Brittin, 1957; Brownstein and Tarr, 1979), we propose the following model described in Figure 1 and Eqs. (1)–(2). We assume that in the intra-axonal space, there are two distinct water pools in fast exchange (Zimmerman and Brittin, 1957): the surface water immediately adjacent to the axonal membrane, e.g., see Le Bihan (2007), and the cytoplasmic water (i.e., axoplasm). The T2 and T1 relaxation times of the surface water are shorter because this water layer is in a more ordered state (both spatially and orientationally) than pure water (Halle, 1999; Finney et al., 2004) and the cytoplasmic water, due to the strong water-tissue interactions (Levy and Onuchic, 2004; Zhang et al., 2007). Moreover, the relaxation times of the cytoplasmic water are expected to be smaller than those of pure water and Cerebrospinal fluid because the water molecules in this pool could interact with cytoskeletal elements and a higher number of macromolecules (Beaulieu, 2002). The fast exchange assumption is reasonable if we consider that water molecules, on average, travel distances much larger than the axon radius for typical diffusion and echo times, as employed in this study.

Figure 1. Transmission electron micrograph of a myelinated axon (adapted) illustrating the employed relaxation model for the intra-axonal space, composed of two pools (arbitrarily colored in green and red for illustrative purposes) in fast exchange (Zimmerman and Brittin, 1957). This model is equivalent to the Brownstein and Tarr (1977) model in the fast diffusion limit. The structured water (Le Bihan, 2007) adjacent to the inner axon surface (red) has a shorter T2 than the cytoplasmic water (green). As the cytoplasmic water (i.e., axoplasm) interacts with large proteins, organelles, and cytoskeletal elements (LoPachin et al., 1991; Beaulieu, 2002), its T2 is shorter than pure water. An equivalent model was assumed for the T1 relaxation. [This transmission electron micrograph was deposited into the public domain by the Electron Microscopy Facility at Trinity College]. This is a file from the Wikimedia Commons, a collection of freely usable media files, under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation (Source: https://en.wikipedia.org/wiki/File:Myelinated_neuron.jpg). This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license (CC BY-SA 3.0). Any copy and remix of the original file must be distributed under the same or compatible license as the original.

According to the general model provided by Zimmerman and Brittin (1957), the inverse of the observed intra-axonal T2 can be modeled by the linear combination of the inverse relaxation times of the surface water and the cytoplasmic water pools, weighted by their volume fractions. Although the volume of the surface water layer is much smaller than the total intra-axonal volume, its T2 time (T2s) is much shorter than that of the cytoplasmic water (T2c). It thus could have a non-negligible impact on the observed intra-axonal (T2a) time. These assumptions are summarized in the following model:

1T2a=V-S⁢εV⁢1T2c+S⁢εV⁢1T2s(1)

≈1T2c+2⁢εr⁢1T2s 

=1T2c+2⁢ρ2r, 

where ρ2=ε⁢/⁢T2s is the T2 surface relaxivity; S is the surface area of the axonal membrane; V is the intra-axonal volume; ε is the thickness of the water layer. Note that when assuming a cylindrical axonal geometry, as commonly done in dMRI, the surface-to-volume ratio depends on the inner axon radius, S/V = 2/r. An equivalent expression was obtained for the intra-axonal T1 time.

1T1a≈1T1c+2⁢ρ1r,(2)

Where ρ1=ε⁢/⁢T1s is the longitudinal surface relaxivity.

2.2. Axon radius estimation from intra-axonal relaxation times

By inverting Eqs. (1) or (2) it is possible to predict the inner axon radius from the estimated intra-axonal T2a and T1a relaxation times, respectively.

r≈2⁢ρ21T2a-1T2c,r≈2⁢ρ11T1a-1T1c.(3)

However, this approach requires knowing T2c and ρ2 or T1c and ρ1 in advance. As these parameters are unknown and cannot be estimated for each brain voxel without additional data, here we propose a histologically-informed calibration approach to calculate them.

The calibration is based on assuming that any dependence of T2c or T1c on the axon radius, owing to potential changes in the intra-axonal structure with the axon size (e.g., density of proteins, organelles, and cytoskeletal elements), is weak and can be neglected. That is, the dependence of T2a and T1a on the axon radius is dominated by the surface-to-volume ratio terms in Eqs. (1) and (2). Therefore, we assume that T2c, ρ2, T1c and ρ1 are constant across axons with different sizes. Nevertheless, we noted that the calibration process is equally helpful in predicting axon radius when T2c or T1c linearly varies with the radius. In that case, the linear models [Eqs. (1) and (2)] can be rewritten in terms of two alternative parameters. For more details, see the discussion subsection “Is the cytoplasmic T2 constant?”

In this study, we collected in vivo diffusion-T1-T2 MRI data in a human brain to estimate T2a and T1a. We employed a reduced diffusion-T2 relaxation sequence to validate our model further by scanning the same subject and two additional healthy volunteers, which allowed us to estimate T2a. Subsequently, we used histological information from four regions of interest (ROIs) located in the CC of a postmortem human brain to measure the mean histological axon radii. The mean intra-axonal relaxation times and histological axonal radii estimated in the four ROIs were combined to estimate T2c and ρ2, and T1c and ρ1 via linear regression (i.e., calibration step) from Eqs. (1) and (2). Then, using the calibrated parameters, we predicted axon radius in another eleven CC ROIs for each scanned subject via Eq. (3). Finally, we employed a second histological dataset containing data from nine postmortem human brains to further validate our results. All the details are provided in the “Methods” section.

3. Methods 3.1. Intra-axonal diffusion-relaxation models

As in McKinnon and Jensen (2019), we assume that for b = 6,000 s/mm2 the in vivo dMRI signal comes from the intra-axonal space. Thus, the diffusion-T1-T2 relaxation model for the measured signal M for a given b, diffusion gradient unit vector g^, echo time (TE), repetition time (TR), and inversion time (TI) is

M⁢(b,g^,T⁢E,T⁢I)=k⁢P⁢D⁢fa⁢Ma⁢(b,g^)⁢exp⁡(-T⁢ET2a)(4)

 |1-2⁢exp⁡(-T⁢IT1a)+exp⁡(-T⁢RT1a)|+η,

where k is a scalar that depends on the MRI machine, pulse sequence, image-reconstruction algorithm, digital converter, etc.; PD is the proton density; fa is the intra-axonal water volume fraction; Ma⁢(b,g^) denotes the orientation-dependent diffusion-weighted signal from the intra-axonal compartment; η is the experimental noise, assumed to be additive; |x| denotes the absolute value of x; T2a and T1a are the intra-axonal relaxation times.

Following the approach of Edén (2003), Lasiè et al. (2014), Kaden et al. (2016a,b), Eq. (4) can be simplified by computing the orientation-averaged spherical mean signal M¯ as:

M¯⁢(b,T⁢E,T⁢I)≈K⁢exp⁡(-T⁢ET2a)(5)

 |1-2⁢exp⁡(-T⁢IT1a)+exp⁡(-T⁢RT1a)|,

where T2a and T1a are the parameters to be estimated, along with the constant K (per voxel) that is proportional to the intra-axonal water volume fraction (i.e., K=k⁢P⁢D⁢fa⁢M¯a⁢(b,g^)); it also depends on the intra-axonal diffusivities via M¯a.

It is important to note that in Eqs. (4)-(5), the T1 relaxation terms follow the standard relaxation model (Bydder et al., 1998), which assumes an ideal inversion pulse (Pykett et al., 1983; Barral et al., 2010). Other acquisition sequences may require different models. For a comprehensive review of alternative relaxometry sequences and models, please refer to Stikov et al. (2015).

The diffusion-relaxation model in Eq. (5) is a more general version of the model proposed by McKinnon and Jensen (2019) for an inversion recovery sequence incorporating T1 relaxation. The diffusion-T2 model for dMRI data collected at multiple TEs (McKinnon and Jensen, 2019) without considering T1 effects is,

M¯⁢(b,T⁢E,T⁢I)≈K⁢exp⁡(-T⁢ET2a).(6)

3.2. MRI data acquisition and preprocessing

Human brain MRI data were acquired from three healthy volunteers, and one of them was scanned twice on a Siemens Connectom 3T system with 300 mT/m diffusion gradients (Cardiff University Brain Research Centre, Wales, UK). The ethics committee approved the study, and the participant provided written informed consent.

Two diffusion-relaxation protocols were implemented. A longer diffusion-T1-T2-weighted imaging sequence was designed to obtain independent estimates of the axon radius from the first subject’s intra-axonal T1 and T2 times (male, 28 years old). A reduced diffusion-T2 protocol was employed to scan three subjects (age-range = 28–39 years, mean-age = 32.3 ± 4.8 years, males), including the first subject that was also scanned with the longer sequence. Accordingly, for the second sequence, the axon radii were estimated from the intra-axonal T2 times.

The diffusion-T1-T2 relaxation sequence comprised four images with b = 0 s/mm2 and 48 diffusion directions at b = 6,000 s/mm2 (diffusion gradient, 275 mT/m; diffusion times Δ/δ = 22/8 ms) for each of the following nine (TE, TI) combinations (in ms): (80, 200), (110, 200), (110, 331), (150, 200), (80, 906), (110, 906), (110, 1,500), (150, 906), (150, 1,500). The TIs were chosen empirically from relatively small to large values to obtain maps with different visual contrasts without nullifying the WM signal. The lowest TE was set to minimize the contribution of the myelin water (Mackay et al., 1994) to the measured signal, and the largest TE was chosen as a trade-off between image contrast and noise. For each (TE, TI) pair, one additional image with b = 0 s/mm2 and opposite phase encoding direction was acquired to correct susceptibility distortions (Andersson et al., 2003; Andersson and Sotiropoulos, 2016). Figure 2 shows the nine pairs of TEs and TIs. The TR was 5,000 ms, and the voxel size was 2.5 × 2.5 × 3.5 mm3. Ten slices were acquired with matrix size and field of view of 88 × 88 and 220 × 220 mm2, respectively. The acceleration factor was 2, and the total acquisition time was 42 min.

Figure 2. Orientation-averaged spherical mean signals for each pair of TE and TI (TE, TI) in ms. These images were used to fit the diffusion-relaxation model in Eq. (5).

The diffusion-T2 protocol employed a dMRI sequence that was repeated by changing the TE, using the following four values TEs = (73, 93, 118, and 150) ms with TR = 4,100 ms. The other sequence parameters (i.e., acceleration factor, diffusion times, b-value, diffusion directions, number of b0s images, diffusion gradient strength, matrix size, and field of view) were equal to those employed in the previous diffusion-T1-T2 sequence. The number of slices was 46, and the voxel size was 2.5 × 2.5 × 2.5 mm3. The acquisition time per TE was 5 min, and the total scan time was 20 min.

Additionally, a structural T1–weighted (T1w) image was collected for each subject using a 3D MPRAGE sequence with the following parameters: TR = 2,300 ms, TE = 2 ms, TI = 857 ms, voxel size = 1 mm isotropic, and flip angle = 9°, for the purposes of spatial normalization.

The nine diffusion-T1-T2 4D volumes with different TEs and TIs, and the four diffusion-T2 4D volumes with different TEs were preprocessed separately in the following order: (1) noise level estimation and removal using the MP-PCA method (Veraart et al., 2016) by using the matrix centering and patch-based aggregation options (Manjon et al., 2013), as implemented in dipy (Garyfallidis et al., 2014); (2) attenuation of the Rician-noise dependent bias in the signal by implementing the postprocessing correction scheme proposed by Gudbjartsson and Patz (1995) and (3) motion, geometric distortions, and eddy current corrections using the “topup” and “eddy” tools included in FSL (Andersson et al., 2003; Andersson and Sotiropoulos, 2016).

3.3. Estimation of the intra-axonal relaxation times

Diffusion-T1-T2 model: after computing the spherical mean signal for each pair of the preprocessed diffusion-T1-T2 datasets with different TEs and TIs (see Figure 2), the intra-axonal relaxation times were estimated by fitting the diffusion-relaxation model in Eq. (5) using the “L-BFGS-B” method for bound constrained minimization included in the Scipy python library (Virtanen et al., 2020), with the following bounds: 0 ≤ K < ∞, 40≤T2a⁢(m⁢s)≤2000, and 300≤T1a⁢(m⁢s)≤5000. The bounds for the intra-axonal relaxation times were chosen to be higher and lower than those expected for the myelin water and Cerebrospinal fluid (Mackay et al., 1994; Labadie et al., 2014), respectively.

Diffusion-T2 model: the estimation was performed by fitting the diffusion-relaxation model in Eq. (6) to the spherical mean signals estimated from the diffusion-T2 data, using the L-BFGS-B method (Virtanen et al., 2020) with the following bounds: 0 ≤ K < ∞, 40≤T2a⁢(m⁢s)≤2000.

3.4. Histological samples

Two histological datasets were employed. The first one contains two histological samples measured on the same subject. The first sample, which we call “Histology1,” was measured and reported by Caminiti et al. (2009). For completeness, we provide a summary of the histological procedures. Axon radii were measured in four regions of interest (i.e., ROI2, ROI5, ROI8, and ROI10) in the midsagittal CC cross-section of a postmortem human brain (female, 63 years old). These ROIs include axons connecting the prefrontal, motor, parietal and visual cortices, respectively. All analyses were performed with Neurolucida 7 software (MBF Biosciences) and a digital camera-mounted Olympus BX51 microscope. Three sagittal blocks of the CC were removed from the brain. The sample was immersion-fixed in 4% (w/vol) paraformaldehyde in phosphate-buffered saline solution within 27–30 h of death, cryoprotected, cut frozen, and stained for myelin. Axons were sampled within 112 × 87 μm2 frames divided into 25-μm squares. The axonal profiles were chosen for measurement if they presented a dark complete or nearly complete myelin ring with a clear center. Longitudinally cut axons were excluded, and the radius of slightly obliquely cut axons (which appeared as ellipses) was approximated to its smallest radius. Since fixation artifacts were frequent, the sampling was restricted to profiles that could be followed through the thickness of the whole section. Limitations of the optical microscopy prevented measurement of axons radius smaller than ∼0.17 μm. A different number of axons were measured per ROI, ranging from 1,178 (ROI10) to 9,605 (ROI2) axons. No correction for shrinkage effects was applied to the measured radii because accurate shrinkage estimates were unavailable. For more technical details, see Caminiti et al. (2009). The second sample, which we call “Histology2,” was measured by the same team (Prof. Giorgio Innocenti) using the same material and following the same sampling procedure. The main difference was that this time, eleven ROIs (i.e., ROI0-ROI10) encompassing the whole midsagittal CC cross-section were analyzed, and the number of measured axons per ROI was smaller: from 153 (ROI5) to 720 (ROI1) axons. It is important to note that the spatial locations of ROI2, ROI5, ROI8, and ROI10 are the same in both histological samples. However, the sampling procedure employed in the Histology2 sample was repeated without including the axons measured in the Histology1 sample. The anatomical location of the ROIs in both histological samples and the number of measured axons per ROI are displayed in Figure 3.

Figure 3. Anatomical location of the two independent histological samples of the first histological dataset (Histology1 and Histology2) taken from eleven regions of interest (ROIs) in the Corpus Callosum. The number of studied axons per sample and ROI are reported for each case. The second sample (Histology2) consisted of axons not included in the first sample (Histology1).

The second histological dataset, which we call “Histology3,” was reported by Wegiel et al. (2018). This electron microscopic study of the CC included nine control subjects (age-range = 4–52 years old; mean-age = 26.3 ± 15.8 years; postmortem-interval = 15 ± 6.6 h; six males and three females) with well-preserved CC ultrastructure. Each brain was fixed in 10% buffered formalin for at least 3 months, washed for 24 h in water to remove fixer, dehydrated, embedded in celloidin, and cut into 200-μm-thick sections. Samples were oriented to cut axons perpendicularly to the long axon axis and stained with a 2% solution of p-phenylenediamine. Each section was stained with uranyl acetate and photographed at a magnification of 15,000x using a Hitachi H7500 transmission electron microscope with an Advanced Microscopy Technique (AMT) Image Capture Engine (Danvers, MA). Axons from five different segments (i.e., I, II, III, IV, and V) of the midsagittal CC cross-sections of the nine control subjects were measured. The study was limited to myelinated axons, which were better preserved than non-myelinated axons. For each case, 12 electron micrographs were used, and background correction was applied to reduce the risk of distortions during image analysis. Axons were manually delineated, and the Image J software was employed to obtain the axon radius (Feret’s radius, μm) and area (μm2). No correction for shrinkage effects to the measured radii was reported. The total number of axons measured in the nine control subjects was 15,085 (1,676 per subject, and 335 per segment, on average). For additional details, see Wegiel et al. (2018).

We note that the CC segments employed in both histological datasets (i.e., Histology1-Histology2 and Histology3) are related. Segment I (Histology3) approximately corresponds to the union of ROI0, ROI1, and ROI2 (Histology1-Histology2); segment II is located around ROI3 and ROI4; the union of segments III and IV is similar to the union of ROI5, ROI6, and RO7; and segment V covers ROI8, ROI9 and ROI10. These relationships were used to compare the histological estimates from both studies and the MRI-based radius estimates.

3.5. Estimation of the mean histological effective radius

For each ROI of each sample, we estimated the mean histological axon radius. However, as the mean axon radius estimated from MRI generally differs from the mean histological radius (Burcaw et al., 2015; Veraart et al., 2020), we derived an approximate expression for the mean effective radius for our diffusion-relaxation models, finding that reff ≈ ⟨r2⟩/⟨r⟩, which differs from the previous result reported in Burcaw et al. (2015), Veraart et al. (2020) (The complete derivation is reported in the Appendix section). This key result shows that the mean effective radius derived from our model is not heavily weighted by the tail of the axon radius distribution as that in Burcaw et al. (2015), Veraart et al. (2020). Consequently, we used this expression to estimate the mean effective histological axon radius for each CC ROI in both samples (see Figure 3), which was compared with the MRI-predicted mean radius.

In order to estimate the effective radius, knowing both the mean histological axon radius and the mean squared radius is required. For the Histology1 and Hostology2 samples, these values were calculated from the whole radius distribution per ROI. We don’t have access to the radius distributions of the Histology3 sample. Fortunately, in that study, the mean histological radius and the mean axon area were reported (Wegiel et al., 2018). We used the mean axon area to estimate the mean squared radius assuming a circular geometry.

3.6. Spatial registration

The histologist that measured the axons in the Histology2 sample drew the locations of the eleven histological ROIs on the structural T1w image of the subject scanned using the diffusion-T1-T2 sequence, which we used to create a cluster mask. Therefore, we used that T1w image as a reference to spatially register the estimated parameter maps for all the subjects (i.e., intra-axonal relaxation times and K maps). The same affine registration matrix and non-linear deformation field were applied to each subject’s estimated parameter map. These registration parameters were determined by non-linearly registering the estimated K map (whose visual appearance is similar to a T1w image, e.g., see Figure 4 in the results section) to the reference T1w image. The registration was carried out using the state-of-the-art (Klein et al., 2009) Symmetric Normalization (SyN) method (Avants et al., 2008) implemented in the ANTs software (ANTsPy). Before the registration, we corrected the K map and T1w image for spatial intensity variations due to B1-Radiofrequency field inhomogeneities using FSL (Smith et al., 2004). All the registered images were visually inspected to verify the accuracy of the normalization procedure. All the subsequent analyses employed the registered maps. Furthermore, the ROIs were eroded to remove peripheral voxels that do not correspond to the corpus callosum and are affected by partial volume effects with surrounding tissue and CSF.

Figure 4. Axial slices of the T2a, T1a, and K maps estimated from the in vivo diffusion-T1-T2 relaxation MRI data in native space (i.e., before registering the images to the reference T1w image). Note that the intra-axonal relaxation times are only meaningful in the white matter because the assumptions underlying the estimation method are invalid in gray matter or CSF. The values of K (in arbitrary units) are higher in the white matter because this parameter is proportional to the intra-axonal volume. We highlight two regions with different intra-axonal relaxation times: the genu of the Corpus Callosum and the corticospinal tract (CST).

The number of voxels included in each ROI ranges from 170 (ROI0) to 604 (ROI1) in the cluster mask defined in the reference T1w image. The equivalent number of voxels in the native space of the diffusion-T1-T2 MRI data with a lower spatial resolution (obtained after applying the resulting non-linear inverse registration to the cluster mask) ranges from 10 (ROI10) to 20 (ROI1).

3.7. Calibration step to predict axon radii

The first sample of the first histological dataset (Histology1) was employed to estimate the unknown parameters of the surface-based relaxation models in Eqs. (1)-(2). These equations were fitted independently using the mean intra-axonal T2 and T1 times and the mean effective histological radii estimated in the same four ROIs of the Histology1 sample. The fitting allowed us to determine the cytoplasmic T2c and T1c times and the surface relaxivity coefficients ρ2 and ρ1, which best explain the data in these regions. This was done by fitting the linear equation y = mx + n, where y=1⁢/⁢T2c and x = 2/r for values from the four CC ROIs. Note that these parameters can be estimated as ρ2 = m and T2c=1⁢/⁢n. A similar independent linear model was used to fit the T1 data for estimating ρ1 and T1c.

Subsequently, we predicted the mean effective axon radii, using Eq. (3), in the eleven CC ROIs of the second sample of the first dataset (Histology2) and the CC segments defined in the second histological dataset (Histology3). The forecasted and histological axon radii were compared using a linear regression model. The linear relationship among the parameters was quantified and tested by the slope and intercept of the fitted regression line and the Pearson correlation coefficient. It is important to mention that when there are two variables, such as in our study, the p-value of the slope of the regression line and the p-value of the correlation coefficient are identical. Therefore, to avoid redundancy, we have reported only the p-values of the slopes in our findings. In the Results section, we present the raw p-values without applying the correction for multiple comparisons. However, in the Discussion section, we mention the analyses that have survived the Bonferroni correction.

4. Results 4.1. diffusion-T1-T2 and Histology1-Histology2 data

Figure 4 shows the T2a, T1a, and K maps estimated from the in vivo diffusion-T1-T2 MRI data for different brain slices. The estimated relaxation times are within the expected range for white matter. The values in the whole medial part of the CC were distributed in the following ranges: 70<T2a⁢(m⁢s)<130, 650<T1a⁢(m⁢s)<760.

The results of the calibration step are depicted in Figure 5. It shows the regression line fitting the inverse of the mean intra-axonal T2 per ROI to the inverse of the mean histological radius in the four ROIs of the Histology1 sample (for more details, see Figure 1), employing the surface-based relaxation model in Eq. (1), as described in the subsection “Calibration step to predict axon radii.” The correlation coefficient between both variables was 0.97, and the p-value of the slope (i.e., for a hypothesis test whose null hypothesis is that the slope is zero) was p = 0.03. We found the calibrated parameters T2c≈126.97 ms and ρ2 ≈ 1.16nm/ms from the estimated coefficients.

Figure 5. Linear fitting of the inverse of the intra-axonal T2 times (y-axis) estimated from the in vivo diffusion-T1-T2 MRI data to the inverse of the inner axon radius (x-axis) measured from the first histological sample (Histology1) of the first histological dataset. The scatter plot depicts the mean values computed for all the voxels inside four corpus callosum (CC) regions of interest, corresponding to ROI2, ROI5, ROI8, and ROI10 in the Histology1 sample. The number of axons sampled for each CC ROI is displayed in the legend. The intercept and slope of the regression line were 0.0079 ms–1 and 0.00116, respectively. The slope of the regression line was significantly different from zero (p = 0.030).

In Figure 6, we compare the effective histological radii in the eleven ROIs of the Histology2 sample and those predicted using the intra-axonal T2 times estimated from the in vivo diffusion-T1-T2 MRI data [Eq. (3)]. The intercept and slope of the regression line were 0.026 μm and 1.055, respectively; the correlation coefficient was 0.676, and the p-value for the slope and the correlation was p = 0.022. To further investigate the data, we analyzed a subset of seven ROIs, excluding the four ROIs in the same locations as those in the Histology1 sample. We obtained a slightly higher intercept of 0.12 μm and a smaller slope of 0.89 compared to the analysis conducted with eleven ROIs. The resulting correlation coefficient decreased to 0.557, and the p-value for the slope was not significant, p = 0.19.

Figure 6. Linear fitting of the effective histological radius estimated from the second histological sample (Histology2) of the first histological dataset to the predicted radius from the intra-axonal T2 times, calculated from the in vivo diffusion-T1-T2 MRI data. The scatter plot depicts the mean values computed for all the voxels inside the eleven corpus callosum (CC) regions, corresponding to ROI0-ROI10. The number of axons sampled for each CC ROI is displayed in the legend. The slope of the regression line was significantly different from zero (p = 0.022).

Figures 7, 8 show results from similar experiments using the intra-axonal T1. Figure 7 depicts the linear fitting of the inverse of the mean intra-axonal T1 per ROI estimated from in vivo diffusion-T1-T2 MRI data to the inverse of the mean effective radius corresponding to the Histology1 sample [Eq. (2)]. The correlation coefficient between both variables was 0.755 lower than that previously found for the intra-axonal T2 in Figure 5, and the p-value of the slope and the correlation did not reach statistical significance, p = 0.25. From the estimated coefficients, we found the calibrated parameters to be T1c≈870 ms and ρ1 ≈ 0.087 nm/ms.

Figure 7. Linear fitting of the inverse of the intra-axonal T1 times (y-axis) estimated from the in vivo diffusion-T1-T2 MRI data to the inverse of the inner axon radius (x-axis), measured from the first histological sample (Histology1) of the first histological dataset. The scatter plot depicts the mean values computed for all the voxels inside four corpus callosum (CC) regions, corresponding to ROI2, ROI5, ROI8, and ROI10. The number of axons sampled for each CC ROI is displayed in the legend. The intercept and slope of the regression line were 0.0011 and 0.000087. The p-value for the slope was not statistically significant (p = 0.23).

Figure 8. Linear fitting of the effective histological radius determined in the second histological sample (Histology2) of the first histological dataset to the predicted radius from the intra-axonal T1 times estimated from the in vivo diffusion-T1-T2 MRI data. The scatter plot depicts the mean values computed for all the voxels inside the eleven corpus callosum (CC) regions, corresponding to ROI0-ROI10 in the Histology2 sample. The number of axons sampled for each CC ROI is displayed in the legend. The slope of the regression line was significantly different from zero (p = 0.039).

The linear relationship between the effective mean axon radius estimated in the Histology2 sample and the radius predicted by using the intra-axonal T1 times [Eq. (3)] is shown in Figure 8. The intercept and slope of the regression line were 0.064 μm and 1.002, respectively. The correlation coefficient was 0.628, and the slope was significant, p = 0.039. When analyzing the subset of seven ROIs, excluding the four ROIs from the Histology1 sample, we obtained a new slope of 0.962 (p = 0.16), which was not statistically significant. The intercept was 0.065 μm, and the correlation coefficient was 0.598.

Table 1 reports the mean histological effective axon radius per ROI and the predicted values from the intra-axonal T2 and T1 times, respectively. The predicted axon radius from both intra-axonal T2 and T1 times were very similar to each other. A linear fitting between both estimates revealed a slope close to one (0.947) and an intercept close to zero (0.041 μm). The slope was significantly non-zero (p = 4e–5), and the correlation coefficient was 0.927.

Table 1. Mean effective radius (in μm) for each region of interest (ROI) in the corpus callosum.

4.2. Diffusion-T2 and Histology1-Histology2-Histology3 data

We complement the results presented in the previous section by reporting the predicted radii for the subjects scanned with the diffusion-T2 MRI sequence and by including the Histology3 dataset. Notably, the parameters T2c and ρ2 were not recalibrated for these subjects; instead, we used the values estimated in the previous section.

The estimated intra-axonal T2 values in the whole medial part of the CC for the three subjects were distributed in the following ranges 80–130 ms, 90–125 ms, and 85–115 ms, respectively.

In Figure 9, the predicted mean effective radius, derived from the intra-axonal T2 times of the three subjects, is presented for all the CC ROIs. The figure also depicts the mean histological effective radius for the three histological samples (Histology1, Histology2, and Histology3).

Figure 9. Predicted axon radius from intra-axonal T2 times estimated from the in vivo diffusion-T2 MRI data for the eleven ROIs (ROI0-ROI10) of the Histology2 sample. Additionally, as a reference, the mean effective histological radius calculated from the three histological samples (Histology1, Histology2, and Histology3) is also reported. Although the number and location of the ROIs used in the Histology3 sample differ from those employed in the Histology1-Histology2 samples, they can be regrouped to cover similar anatomical areas (see subsection “Histological samples” for more details). The histological and T2-based radii follow the expected “low-high-low” trend in axon radii. The axon radii from the Histology1-Histology2 samples are consistently higher (about 25%) than those in the Histology3 sample.

To assess the validity of the calibrated parameters, which were estimated from the subject scanned with the diffusion-T1-T2 sequence, for the subjects scanned with the diffusion-T2 sequence, we repeated the calibration process using the mean intra-axonal T2 times estimated for the three subjects and the Histology1 sample as a reference, as before. The recalibrated parameters were remarkably similar to those obtained previously: T2c≈127.17 ms and ρ2 ≈ 1.13nm/ms.

We compared the T2-based predicted radii for the subject that underwent two scans, using both diffusion-relaxation MRI sequences, which values are reported in Figure 9 and Table 1 (as Subject 3). The linear fitting between both estimates produced a statistically significant regression line with a slope close to one (0.993, p < 0.001) and an intercept close to zero (−0.029 μm). The correlation coefficient between both estimates was 0.884.

Finally, we employed the calibrated model to predict the axon radius across the whole WM. Axial and sagittal slices of the voxel-wise T2-based inner axon radius estimated for all the subjects are shown in Figure 10. The maps are approximately symmetrical, the spatial variability of the estimated radius is apparent in both slices, and all subjects show a similar pattern of small and big axons in the same anatomical regions.

Figure 10. Axial and sagittal slices of the T2-based inner axon radius for the three scanned subjects. Subject3 underwent two scans, with scan2 (46 slices) and scan1 (10 slices) representing the in vivo diffusion-T2 and diffusion-T1-T2 MRI data, respectively. All maps were normalized to the reference T1w image, where the histological CC ROIs were defined, and the predicted radii were plotted over the reference image. A white matter mask was used to suppress voxels in gray matter or cerebrospinal fluid.

5. Discussion

This proof-of-concept study shows that (1) the intra-axonal T2 and T1 relaxation times are highly modulated by the axon radius (see Figures 5, 7), as measured from histological data (see Figure 3), (2) a simple surface-based relaxation model can explain this dependence (see Figure 1), and (3) the intra-axonal relaxation times