Estimating axon radius using diffusion-relaxation MRI: calibrating a surface-based relaxation model with histology

Axon radius is a potential biomarker for brain diseases and a crucial tissue microstructure parameter that determines the speed of action potentials. Diffusion MRI (dMRI) allows non-invasive estimation of axon radius, but accurately estimating the radius of axons in the human brain is challenging. Most axons in the brain have a radius below one micrometer, which falls below the sensitivity limit of dMRI signals even when using the most advanced human MRI scanners. Therefore, new MRI methods that are sensitive to small axon radii are needed. In this proof-of-concept investigation, we examine whether a surface-based axonal relaxation process could mediate a relationship between intra-axonal T2 and T1 times and inner axon radius, as measured using postmortem histology. A unique in vivo human diffusion-T1-T2 relaxation dataset was acquired on a 3T MRI scanner with ultra-strong diffusion gradients, using a strong diffusion-weighting (i.e., b = 6,000 s/mm2) and multiple inversion and echo times. A second reduced diffusion-T2 dataset was collected at various echo times to evaluate the model further. The intra-axonal relaxation times were estimated by fitting a diffusion-relaxation model to the orientation-averaged spherical mean signals. Our analysis revealed that the proposed surface-based relaxation model effectively explains the relationship between the estimated relaxation times and the histological axon radius measured in various corpus callosum regions. Using these histological values, we developed a novel calibration approach to predict axon radius in other areas of the corpus callosum. Notably, the predicted radii and those determined from histological measurements were in close agreement.


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 1500 mT/m in animal preclinical scanners) are used (Assaf et al., 2004(Assaf et al., , 2008Assaf 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 micrometre 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., 2008Assaf et al., , 2013Alexander et al., 2010Alexander et al., , 2019Dyrby et al., 2013Dyrby et al., , 2018Novikov et al., 2019;Veraart et al., 2020;Fan et al., 2020;Jelescu 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 r 2 ), larger axons contribute more than smaller axons to the measured signal. After considering these two factors together, an approximate expression for the mean 'effective' dMRIbased radius eff r per voxel can be derived, which depends on the higher-order moments of the unknown axon radius distribution. The resulting analytical expression ( ) 1 62 4 eff r r r  (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).

5
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 intraaxonal relaxation times. This might be explained by the fact that approaches to estimating the intraaxonal 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>4000 s/mm 2 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 fibre 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 6 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.

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.
According to the general model provided by (Zimmerman and Brittin, 1957), the inverse of the observed intra-axonal T2 can be modelled 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 where 11 s T  = is the longitudinal surface relaxivity.
Insert Figure 1 (1.5 columns) 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 coloured in green and red for illustrative purposes) in fast exchange (Zimmerman and Brittin, 1957). This model is equivalent to the Brownstein-Tarr model in the fast diffusion limit (Brownstein and Tarr, 1977). 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].

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- where k is a scalar that depends on the MRI machine, pulse sequence, image-reconstruction algorithm, digital converter, etc.; PD is the proton density; a f is the intra-axonal water volume fraction; ( , ) a M 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; 2 a T and 1 a T are the intra-axonal relaxation times.
Following the approach of (Edén, 2003;Lasič et al., 2014;Kaden et al., 2016bKaden et al., , 2016a, Eq. (4) can be simplified by computing the orientation-averaged spherical mean signal M as: where 2 a T and 1 a T 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., 10 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).

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.  (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/mm 2 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 5000ms, and the voxel size was 2.5x2.5x3.5 mm 3 . Ten slices were acquired with matrix size and field of view of 88x88 and 220x220 mm 2 , respectively. The acceleration factor was 2, and the total acquisition time was 42 minutes.
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=4100ms. 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.5x2.5x2.5 mm 3 . The acquisition time per TE was 5 minutes, and the total scan time was 20 minutes.
Additionally, a structural T1-weighted (T1w) image was collected for each subject using a 3D MPRAGE sequence with the following parameters: TR = 2300 ms, TE = 2 ms, TI=857 ms, voxel size = 1 mm isotropic, and flip angle=9°, for the purposes of spatial normalisation.
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 centring and patch-based aggregation options (Manjon et al., 2013), as implemented in dipy (Garyfallidis et al., 2014) (https://dipy.org/); (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).
Insert Figure 2 (1 column) 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).

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 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    ,

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 13 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-30h of death, cryoprotected, cut frozen, and stained for myelin. Axons were sampled within 112×87 μm 2 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 centre. 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 artefacts 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 1178 (ROI10) to 9605 (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 analysed, 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.
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 wellpreserved CC ultrastructure. Each brain was fixed in 10% buffered formalin for at least three months, washed for 24 h in water to remove fixer, dehydrated, embedded in celloidin, and cut into 200-μmthick 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 (μm 2 ). 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 (1676 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.

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 2 eff r r 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.

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 nonlinear deformation field were applied to each subject's estimated parameter map. These registration parameters were determined by nonlinearly 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: https://github.com/ANTsX/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 normalisation procedure. All the subsequent analyses 16 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.
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 nonlinear inverse registration to the cluster mask) ranges from 10 (ROI10) to 20 (ROI1).

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) Figure 4 shows the 2 a T , 1 a T , 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.

diffusion-T1-T2 and Histology1-Histology2 data
The values in the whole medial part of the CC were distributed in the following ranges: 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 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 analysed 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.

19
Insert Figure 6 (1.5 columns) (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 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) Fig 3. The second row lists the radii corresponding to the Histology2 sample. The third and four rows report the predicted axon radii from the intraaxonal T2 and T1 times, respectively, estimated from the in vivo diffusion-T1-T2 MRI data.

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 2 c T 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).
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: 2 127.17 c T  ms and 2 1.13   nm/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 Table 1 and Figure 9 (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.
Insert Figure 9 (1.5 columns) 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.
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.
Insert Figure 10 (1.5-2 columns) 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 normalised 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 grey matter or cerebrospinal fluid.

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 Figs. 5 and 7), as measured from histological data (see Fig. 3), (2) a simple surface-based relaxation model can explain this dependence (see Fig. 1), and (3) the intraaxonal relaxation times may also be sensitive to the smallest axons. Indeed, we did not observe a clear overestimation bias in the estimated axon radius (see Figs. 6,8,and 9) in comparison to the histological values, as reported in previous dMRI studies (Assaf et al., 2008;Alexander et al., 2010;Dyrby et al., 2013;Horowitz et al., 2015) where only the largest radii might have been detected. This result suggests that our new approach may also be sensitive to differences in axon radius below the 'diffusion resolution limit' of ~1-2 μm. One possible explanation for this finding is that the intraaxonal T2 times are not influenced by the strength of the diffusion gradients, as opposed to the intraaxonal radial diffusivities used to estimate axon radii in dMRI. Moreover, we found that the effective mean radius estimated by our approach, i.e., 2 eff r r r  , produces much smaller radii than those 25 from diffusion models heavily weighted by the tail of the axon radius distribution, i.e.,  (Burcaw et al., 2015;Veraart et al., 2020). This important result suggests that, from a modelling point of view, the employed diffusion-relaxation model may be more valuable than previous pure dMRI models for estimating axonal radii. The predicted mean effective radius obtained from the intra-axonal T2 and T1 times fell within a narrow range of 0.52-1.13 µm and 0.51-1.12 µm, respectively, which closely matched the range of histological axon radii (0.57-0.95 µm). The smallest predicted effective radii were observed in ROI1, ROI0, and ROI2, while the largest radii were found in ROI6 and ROI5, followed by ROI4, ROI9 and ROI10 (refer to Table 1). Nevertheless, we cannot rule out the possibility that the calibration step, informed by the histological values, may have reduced any potential overestimation effect.

( )
Inspecting the estimated 2 a T and 1 a T relaxation maps (see Fig. 4), we notice that both relaxation times tend to be smaller in the genu and splenium of the CC than in the corticospinal tract (connecting the motor cortex to the spinal cord). Although these values could be affected by fibre orientation effects with respect to the B0 field (see the subsection "Orientation dependence on relaxation times" in the Appendix), the corticospinal tract is characterised by axons with larger inner radius (Aboitiz et al., 1992;Innocenti et al., 2014;Barakovic et al., 2021a). This observation agrees with multi-echo T2 relaxometry studies showing that the intra-and extra-axonal T2 times (and the myelin content) in the corticospinal tract are larger than in the CC, e.g., see (Yu et al., 2020;Canales-Rodríguez et al., 2021b, 2021c. A consistent trend was observed in the T2-based predicted axon radii for all three subjects, as shown in Figure 10. The voxel-wise maps in Figure 10 and the ROI-based estimates in Figure 9 agree with previous estimates derived from dMRI data acquired using much higher b-values (Veraart et al., 2021).
In agreement with our results, a previous multi-echo T2 relaxometry study found a positive correlation between axon radius and T2 (including both the intra-and extra-axonal compartments) in six samples of an excised and fixed rat spinal cord (Dula et al., 2010). Moreover, two previous experimental studies investigated the microstructural correlates of T1 in white matter (Hofer et al., 2015;Harkins et al., 2016). In line with our findings, a significant correlation between 1/T1 and axon radius was reported by (Harkins et al., 2016) in white matter tracts of a rat spinal cord. Similarly, the analysis performed by (Hofer et al., 2015) found a tendency for the lowest T1 in the genu of the human CC (composed of densely packed smaller axons) and the highest T1 in the somatomotor region (dominated by fibres with large radii). In those studies, however, the estimated relaxation times characterise the relaxation process in the intra-and extra-axonal compartments combined. In contrast, we report a more specific relationship by analysing the intra-axonal relaxation times associated with the inner axon radius.
A multi-gradient-echo MRI model was proposed to estimate axon density based on the susceptibilitydriven non-monotonic time-dependent MRI signal decay (Nunes et al., 2017). They employed a simple (phenomenological) general-linear model to predict the average axonal diameters using four modelled parameters, including the T2 relaxation times of the intra-and extra-axonal compartments.

Impact on previous and future studies
Our study has important implications for previous and future dMRI studies of white matter microstructure. Previous studies, such as those by ( This simplification may affect the estimation of the intra-axonal diffusivities from which the axon radii are derived. Alternatively, this issue could be attenuated by using sufficiently high b-values, as shown in studies by (Veraart et al., 2020;Pizzolato et al., 2022), which helps eliminate the contribution from the extra-axonal signal. However, this may be insufficient in voxels with a broad distribution of intra-axonal T2 times. These multi-compartment models should be extended to include T2 dependence, as discussed in studies by (Veraart et al., 2018;Lampinen et al., 2019;Tax et al., 2021). Recently, (L et al., 2022) demonstrated that more accurate estimates of neurite size could be obtained by investigating the coupling between relaxation rate and diffusivity using multi-TE diffusion-relaxation MRI data. For further discussion on this issue, please refer to the Appendix subsection "Is the intra-axonal relaxation process mono-exponential and time-independent?".

Underlying assumptions and confounding factors
The proposed diffusion-relaxation model specifically applies to WM regions composed of myelinated axons, where the exchange of water molecules and other macromolecules and elements (such as iron/ferritin) between the intra-and extra-axonal spaces is negligible. In the human brain's CC, for example, more than 95% of axons in most regions are myelinated (Aboitiz et al., 1992). While the non-myelinated portions of the axon (i.e., nodes of Ranvier) contain a high density of voltage-gated ion channels that facilitate ion passage across the axonal membrane, including K+ and Na+, which is associated with a concomitant water flux (Badaut et al., 2002), the myelinated portions of the axon (i.e., internodes) are not exposed to the extracellular environment. Although the axonal membrane in the nodes of Ranvier is semipermeable to small diffusing molecules, such as water, the internodes' 27 length is significantly greater (∼100 times the outer axon diameter (Hursh, 1939;Rushton, 1951)) than the nodes of Ranvier (∼1 µm (Arancibia-Cárcamo et al., 2017)). As a result, most multicompartment T2 (Lancaster et al., 2003;Deoni et al., 2013) and 'standard' dMRI models (see (Novikov et al., 2019) for a review) assume that the measured MRI signal is not significantly affected by the inter-compartmental molecular exchange in WM regions composed of myelinated axons.
Therefore, it is important to note that our model is unsuitable for GM or WM regions affected by demyelination processes, such as in Multiple Sclerosis, or any pathological condition with an increase in intra-axonal iron. These conditions can significantly reduce intra-axonal relaxation times and the estimated radii, rendering our model invalid. However, it is worth noting that we use long TEs in our model. If the intra-axonal T2 time of a given axon is significantly reduced (e.g., below 20-40ms) due to external factors, the contribution of this axon to the overall voxel-wise measured signal will be greatly diminished.
However, it is interesting that our calibration approach could also be extended to cases where water exchange between intra-and extra-axonal spaces is non-negligible, provided the exchange is similar across axons with different radii. In such cases, the effect of the exchange on the observed intraaxonal relaxation times can be modelled by a global scaling of the cytoplasmic relaxation time, which is accounted for during calibration.
A more suitable approach for modelling systems that are coupled by means of a relaxation exchange process could be based on the Bloch-McConnell equations (McConnell, 1958), which generalise the relaxation model employed in this study (Eqs. (5) and (6)). However, fitting such a model requires estimating additional parameters, including membrane permeability and extra-axonal relaxation times, which may be prone to numerical degeneracies. Additionally, the MRI acquisition time required for fitting the Bloch-McConnell model (using both high and low b-values) is longer than that required for our proposed model.
A study on human postmortem brains revealed that T2* is more sensitive than T2 to changes in WM iron concentration (Langkammer et al., 2010). While it is established that the macromolecular and iron content is altered in certain pathologies (Stüber et al., 2014), more research is required to understand how these abnormalities affect the intra-axonal space and how they can impact the intraaxonal relaxation times.
We assume that signals measured at very high b-values are primarily attributable to the intra-axonal space, given that the signals from free-water and extra-axonal compartments decay more rapidly with the b-value (Veraart et al., 2020). To further suppress signals from tissue compartments with very short T2s, such as myelin water (Mackay et al., 1994) and other confined water molecules, we also used long TEs. Hence, the resulting signals are expected to come from intra-axonal water molecules.
However, there are other 1D-stick-like structures in the WM, such as the radiating processes of astrocytes, which can have large diameters that might contaminate the resulting signals (Veraart et al., 2020(Veraart et al., , 2021, as well as cell nuclei, vacuoles, and other restricted compartments (Andersson et al., 2020). Therefore, further studies are necessary to understand the potential effects of these compartments on the measured T2 and predicted radii.

Acquisition sequences
The diffusion-T1-T2 sequence was implemented to investigate the impact of axon radius on the intraaxonal T1 and T2 times independently. Our results demonstrate that both relaxation times are sensitive to changes in axon radius, with T2 exhibiting a slightly higher sensitivity. Consequently, we can obtain two separate estimates of axon radius using the relaxation times calculated from this sequence (see Table 1). However, this is not our recommended acquisition protocol due to the long acquisition time required. Alternatively, a more practical approach for estimating axon radius is to use the diffusion-T2 sequence. A faster version of this sequence could be implemented by utilising only two TEs, although the resulting estimates may be more affected by underlying noise. This possibility shall be investigated in future studies.
When implementing these sequences, it is important to identify the optimal b-value to attenuate the extra-axonal signal. Based on in vivo human brain data and numerical experiments using analytical equations, the general rule of thumb is that a b-value in the range of 4000-6000 s/mm 2 is sufficient (Jensen et al., 2016;McKinnon and Jensen, 2019). In our study, we used the highest b-value within this range. However, it is worth noting that determining the optimal b-value involves a trade-off influenced by the SNR, which is affected by other sequence parameters, including the TE and voxel size. Our data were acquired using the Connectom 3T scanner at CUBRIC, which has been previously used to collect data with b-values up to 30000 s/mm 2 (Veraart et al., 2020(Veraart et al., , 2021. Ultra-high b-values with very strong diffusion gradients are necessary for pure dMRI models to improve sensitivity to smaller axon radii (Nilsson et al., 2017). However, our sequences do not require b-values larger than 6000 s/mm 2 because all the necessary information is derived from the relaxation times, which depend on the TEs/TIs.

Main limitations and future studies
29 While our study provides valuable insights into the relationship between axon radii and MRI relaxation times, it is important to acknowledge some limitations. First, the in vivo diffusionrelaxation MRI data and postmortem histological samples were obtained from different subjects of different ages and genders. Although some studies suggest that there are no sex differences in the fibres composition of the corpus callosum (Aboitiz et al., 1992), others have found age-related changes in axon size (Aboitiz et al., 1996), which may affect the comparison between the postmortem and in vivo measurements. Therefore, the estimated relaxation times of the cytoplasmic water and surface relaxivities must be considered as approximated guide values.
Second, the histological analysis of the second and third histological samples (Histology2 and Histology3, covering eleven and five CC sectors, respectively) are based on a reduced number of axons compared to the first (Histology1) sample. This may introduce sampling biases that could affect the accuracy of the histological radius estimates. An extended discussion is provided in the Appendix subsection "Histological tissue shrinkage and sampling issues". As such, a perfect agreement between the effective histological radius and the predicted MRI-based radius was not expected.
Third, the analysis was confined to the mid-sagittal plane of the CC. Therefore, the estimated mean cytoplasmic relaxation times and surface relaxivities are specific to this region. It is possible that different values may be obtained if other white matter tracts were included in the analysis. However, the extension of the analysis to other regions would require modelling the orientation susceptibility effects, which is beyond the scope of this proof-of-concept study. For more details, refer to the Appendix subsection "Orientation dependence on relaxation times".
Fourth, our study had a relatively small number of data points available for computing correlations, with only four ROIs to implement the calibration. This limited sample size restricts the statistical power and precision of the estimated correlations, leading to increased uncertainty and decreased reliability of the findings. While it is generally recommended to account for multiple comparisons to reduce the risk of false-positive findings, we opted not to implement such correction. Given the exploratory nature of our study, we prioritised sensitivity over stringent control of false positives.
Consequently, our findings should be interpreted cautiously, requiring further validation in independent studies. However, for completeness, we report that if we correct our results for multiple comparisons using the Bonferroni method, only two analyses survive the correction: the correlation of the radii estimated using the T2 and T1 relaxation times reported in Table 1 and the T2-based predicted radii for the subject who underwent two scans, using the two diffusion-relaxation MRI sequences employed in this study.

30
Fifth, the proposed model is not suitable for GM or WM regions affected by demyelination processes or any pathological condition increasing intra-axonal iron. These conditions can significantly reduce intra-axonal relaxation times and the estimated radii, rendering the model invalid. A more detailed discussion of these underlying assumptions and confounding factors can be found in the previous subsection, "Underlying assumptions and confounding factors", in the Discussion section.
Sixth, the estimation of intra-axonal T2 from the spherical mean of the strongly diffusion-weighted signal may be subject to bias due to the presence of isotropically-restricted compartments, including cell nuclei and vacuoles (Andersson et al., 2020). However, this issue can be mitigated by utilizing the spherical variance instead (Pizzolato et al., 2022). For more detailed information, please refer to the Appendix subsection "The effect of spherical cells: spherical mean vs spherical variance." Seventh, although the spherical mean signal is not affected by the presence of fibre crossings and orientation dispersion (Lindblom et al., 1977;Callaghan et al., 1979;Kaden et al., 2016a), it is influenced by the orientation susceptibility effect. In other words, the measured signal still depends on the angle between the B0 magnetic field and the fibre orientation. In our study, the regions of interest were located in the medial part of the CC, where the angle between the B0 vector field and the nerve fibres remains relatively constant. More details on this topic can be found in the Appendix subsection "Orientation dependence on relaxation times.
To better assess the generalisability of our approach, further validation studies are necessary. In particular, we plan to test our method using biomimetic phantoms with known ground truth (Hubbard et al., 2015;McHugh et al., 2018;Huang et al., 2021;Zhou et al., 2021) and ex-vivo data from the same brains and multiple white matter regions. Such datasets would allow us to investigate whether the cytoplasmic relaxation times are truly independent of axon radius (see the Appendix subsection "Is the cytoplasmic T2 constant?"). This could be achieved by repeating the calibration process using different subsets of ROIs and comparing the resulting estimates. However, obtaining sufficient histological ROIs and measured axons per ROI will be crucial to minimize sampling bias and get robust results not affected by noise. Additionally, including data from the same brains (e.g., from non-human studies) will enable us to guarantee that we are studying the same axonal bundles.
An interesting future direction would be to utilise bundle-specific intra-axonal T2 values (Barakovic et al., 2021b) to estimate bundle-specific inner axon radius, which could potentially resolve multiple axonal radii per voxel. This approach may potentially predict axon radius beyond the current dMRI resolution limit using clinical scanners. However, one limitation of translating the diffusionrelaxation MRI sequence to clinical scanners is the decreased signal-to-noise ratio resulting from using high b-values and long TEs. One potential solution to mitigate this could be to reduce the bvalue to 4000-5000 s/mm 2 and use numerical simulations to determine the optimal range of TEs, based on the intra-axonal relaxation times reported in this study and the expected noise range.
Despite these limitations, our study provides a promising approach for estimating axon radii and understanding their relationship with MRI relaxation times. Future studies could address these limitations and expand the analysis to other brain regions to further validate the technique.

Code and data availability statement
The datasets used in this study and the Python code can be made available upon request from the corresponding authors, subject to the following terms and conditions. The mean effective histological radii of the Histology1, Histology2, and Histology3 samples are reported in Table 1 and Figure 9, respectively. We can also share any other derived metric from the Histology1-Histology2 samples.
Additional results for the Histology2 and Histology3 samples are available in (Caminiti et al., 2009) and (Wegiel et al., 2018), respectively. The MRI data will be available upon signing a data-sharing agreement with Cardiff University. Finally, we can provide the Python scripts used in this study upon request. the author has applied a CC-BY public copyright licence to any author accepted manuscript version arising from this submission. The presented study is a tribute to Professor Giorgio Innocenti (1946Innocenti ( -2021.

Appendix: effective axon radius
We derive the mean effective radius that can be estimated from the intra-axonal T2 and T1 relaxation times. For simplicity, we will separately analyse the components of the measured signals that exclusively depend on T2 and T1. The signal arising from the intra-axonal compartments is the sum of signals from the spins inside all axons. The measured T2-weighted signal for a given echo time TE is

Axon radius estimated from
where P is the total number of axons, i N is the number of spins inside the i th axon with transverse relaxation time 2 i T , and k is a constant that depends on the sequence/scanner.
Assuming that the proton density (PD) does not depend on the axon radius, then We estimate a single intra-axonal T2 per voxel, which is equivalent to assuming that all the T2s in Eq.
(10) are equal to 2 a T (i.e., and hence that all axons in the voxel have the same radius r ); in that case, Eq. (10) becomes This is the expression that we used to correct the histological radii, which is the mean effective radius estimated from this relaxation model.
Note that we neglected the TR dependence because, in practice, this experimental parameter is much higher than the intra-axonal T1, and its contribution is minor.
As we estimate a single apparent intra-axonal T1 per voxel, our model is equivalent to assuming that all the T1s are equal to 1 a T (i.e., all axons in the voxel have the same radius r ); thus, Eq. (16) After plugging the surface-based relaxation model in Eq.
(2) and removing common terms on both sides of Eq. (18)

Histological tissue shrinkage and sampling issues
The histological datasets were inspected to investigate the trend in axon radii. As expected, the data followed the "low-high-low" pattern in axon radii, as shown in Figure 9. However, the mean effective histological radii differed between the samples. The axon radii from the Histology1-Histology2 samples were about 25% higher than those in the Histology3 sample. These differences could be due to genuine anatomical variations between the postmortem brains or related to the histological procedures and corresponding tissue shrinkage factors. The T2-based predicted radii in all subjects followed a similar "low-high-low" pattern closer to the values measured in the Histology1 sample, as this was the calibration sample.
In this study, the histological samples were not corrected for tissue shrinkage, which can affect the accuracy of the estimated axon radii. Consequently, the in vivo axons may be thicker than the reported histological values (Barazany et al., 2009;Horowitz et al., 2015). The extent of tissue shrinkage can vary widely depending on the used histological preparation techniques, with reported shrinkage factors ranging from 1-65% (Lamantia and Rakic, 1990;Aboitiz et al., 1992;Houzel et al., 1994;Riise and Pakkenberg, 2011). It is also unclear if shrinkage affects all brain axons equally, as previous research has shown varying shrinkage levels in different cellular compartments (Hursh, 1939).
However, there is currently limited knowledge about the effects of shrinkage on CC axons in the human brain (Innocenti et al., 2019). Please refer to (Dyrby et al., 2018) for further information on tissue shrinkage issues.
Sampling biases can impact histological radii measurements. One issue is that only a small amount of tissue is typically sampled, so the microstructure properties of these regions may not accurately represent properties in other regions within the ROIs (Assaf et al., 2008). Another issue is that larger axons can influence the mean effective radius more than the mean radius of the distribution. Since larger axons are less common, accurately detecting their proportions in a sample requires measuring a larger number of axons. We observed this effect in the Histology1 and Histology2 samples, where the effective radii in the four ROIs used in Histology1 (which had denser spatial sampling) were consistently higher than those in the same ROIs measured in the Histology2 sample (see Figure 9).
For these reasons, the presented histological results should not be considered the definitive "ground truth". Future studies should aim to identify optimal histological procedures, such as those suggested by (Sepehrband et al., 2016b), and also explore the use of neural network approaches for automatic measurement of tens of thousands of axons per ROI to reduce sampling biases (Mordhorst et al., 36 2022). It is also worth noting that because the proposed calibration approach uses histological data as a reference, the predicted radii are relative to the specific histological sample employed.

The effect of spherical cells: spherical mean vs spherical variance
A recent study showed that isotropically-restricted compartments might bias the intra-axonal T2 estimated from the spherical mean of the strongly diffusion-weighted signal (Pizzolato et al., 2022).
Thus, our estimates could be partially affected by cell nuclei, vacuoles, and other types of structures in the white matter (Andersson et al., 2020). As a remedy for that problem, it was proposed to use the spherical variance (Pizzolato et al., 2022) as a 'filter' since the spherical variance of an ordered axon bundle would be high, but in an isotropic component would be close to zero. Although the results obtained in that study are promising, the spherical variance is more sensitive to noise than the spherical mean. Moreover, a larger number of diffusion gradient directions than that used in our study is necessary to employ this novel technique (48 vs >96 in (Pizzolato et al., 2022)). In future studies, we plan to acquire dMRI data using a higher angular resolution to compare both techniques' outputs and filter out any contribution from spherical cells.

Is the cytoplasmic T2 constant?
The cytoplasmic T2 may be influenced by the intra-axonal microstructure, such as the number of organelles and the density of cytoskeletal elements such as neurofilaments, microtubules, and actin, as well as the chemical composition, including the type and density of macromolecules and water content.
Numerous morphometric studies have provided evidence of a linear correlation between neurofilament and microtubule numbers and axonal cross-sectional area (Friede and Samorajski, 1970;Hoffman et al., 1987). These studies suggest that myelinated axons contain more neurofilaments than microtubules and that the axon radius adjusts to maintain a constant density of neurofilaments. It was demonstrated that this relationship is regulated by the relative degree of phosphorylation of the mid-sized and heavy neurofilaments (Rao et al., 1998). Furthermore, the myelin-associated glycoprotein is implicated in the signalling cascade controlling neurofilament phosphorylation (Lunn et al., 2002) and axon radius. As neurofilaments are the more abundant cytoskeletal elements and their density is nearly constant in myelinated axons with different radii, we do not anticipate a relationship between cytoplasmic T2 and axon radius mediated by neurofilament density in the axonal cytoskeleton.
However, a previous study using electron probe x-ray microanalysis (LoPachin et al., 1991) measured the concentrations of biologically relevant elements (such as Na, P, S, Ca, CI, K, and Mg, in mmol/kg dry or wet weight) and water content in the axoplasm of rat optic nerve myelinated axons. The study found that dry and wet weight concentrations of Na, P, S and Ca were not dependent on the axonal radius. In contrast, the axoplasmic concentration of K, CI, and Mg was related to axon radius.
Furthermore, the water content in medium and large axons was similar (between 91% and 92%) but slightly reduced in small axons (89%). These findings suggest that the chemical composition of the axoplasm depends on the axon radius (LoPachin et al., 1991). Therefore, until the effect of K, CI, and Mg on intra-axonal relaxation times is clarified, the surface-based relaxation model employed (i.e., Eqs. (1) and (2)) should be regarded as a first-order approximation.
Despite this limitation, our findings (refer to Figures 5 and 7) suggest a linear relationship between the inverse of intra-axonal relaxation times and axon radius, consistent with predictions made by the surface-based relaxation model we employed. Our empirical results demonstrate that the calibration step enables us to estimate the mean axon radius in various regions of the midsagittal CC (refer to Figures 6, 8, and 9). As we did not observe any significant nonlinear relationships between intraaxonal relaxation times and axon radii (within the range of measured radii in the midsagittal CC), we conclude that any nonlinear dependence is weak and can be disregarded. Hence, either the cytoplasmic relaxation times remain constant, as assumed in this study, or they vary linearly with axon radius. In the following, we present some examples where the calibrated model could predict the axon radius accurately, even if the cytoplasmic relaxation times depend on the axon radius.
Let us consider two distinct scenarios: In the first case, the cytoplasmic T2 increases with the radius (similar to the observed trend for intra-axonal T2 time), while in the second case, it decreases. The former corresponds to the model ( ) 38 where 2 2 k  =+ . Althout the new parameters 2 const T and  cannot be interpreted as the cytoplasmic T2 and surface relaxivity, they can be estimated by employing our calibration approach. Therefore, the resulting model would be equally valid for predicting axon radius.
The second case corresponds to a model that predicts a decrease in In our experiments, we observed a net reduction of 2 a T with r . Hence, the surface relaxivity term must dominate the relaxation over k , i.e., It is important to note that the models presented in this section (Eqs. (20)-(23)) are hypothetical and were discussed to illustrate the flexibility and limitations of the calibration approach in cases where the underlying assumptions are not met. Similar results can be obtained by using the intra-axonal T1 time or assuming a surface relaxivity that depends on the radius.

Orientation dependence on relaxation times
By computing the spherical mean of the diffusion signal, the resulting orientation-averaged signal is independent of the fibre orientation distribution and thus is not affected by the presence of fibre crossings and varying levels of fibre orientation dispersion (Lindblom et al., 1977;Callaghan et al., 1979;Kaden et al., 2016a). However, the spherical mean does not eliminate the dependence on the orientation susceptibility effects, i.e., the measured signal still depends on the angle between the B0 magnetic field and the fibre orientation. Some previous studies have reported this orientation dependence for both the T2* and T2 (Oh et al., 2013;Aggarwal et al., 2016;Gil et al., 2016) and T1 (Knight and Kauppinen, 2016;Knight et al., 2017Knight et al., , 2018Schyboll et al., 2018Schyboll et al., , 2020. Notably, while 39 (McKinnon and Jensen, 2019) reported a significant intra-axonal T2 orientation effect, a recent study found that extra-axonal T2 is more affected by this phenomenon than intra-axonal T2 (Tax et al., 2021). Given these inconsistent findings, further research is needed to determine whether the orientation-dependent T2 effect is significant enough to be considered in this model.
In our study, the regions of interest were located in the medial part of the CC, where the angle between the B0 vector field and the nerve fibres remains relatively constant. Therefore, our findings are not likely affected by B0-orientation-related bias. However, the orientation effect should be modelled in brain regions with different fibre orientations, as it may affect the estimation. Despite this potential limitation, in Figure 10, we present T2-based radius images across the entire white matter, showing the spatial variability of estimated radii across different regions, especially in the sagittal slices depicting the midsagittal CC cross-sections. The estimates from all subjects demonstrate a similar concordant pattern, as well as the maps of the same subject (Subject3) obtained from the two diffusion-relaxation MRI sequences, although some differences are noticeable due to the different voxel sizes used in both acquisitions.

Is the intra-axonal relaxation process mono-exponential and time-independent?
A recent theoretical formulation by (Kiselev and Novikov, 2018) demonstrated how the interplay between diffusion and spin dephasing in a heterogeneous environment could produce a non-monoexponential time-dependent transverse relaxation signal. While this effect may be significant for short TEs, our relatively long TEs (i.e., >73 ms) and diffusion times (Δ=22 ms, δ=8 ms) used in this study (compared to the small intra-axonal space where the restricted diffusion process takes place) indicate that a mono-exponential signal relaxation is expected for the spins inside each axon.
In our study, we estimated a single intra-axonal relaxation time per voxel. However, if axon radii are distributed with non-negligible variance, a more complete formulation must consider distributions of relaxation times. Estimating a nonparametric distribution of relaxation times is problematic from a practical point of view because a large number of TEs/TIs would be required. Nevertheless, an approach similar to that introduced in AxCaliber ( Assaf et al., 2008) could be adopted by assuming a parametric form for such distributions, as shown in (Sepehrband et al., 2016a). Future studies should investigate this generalization further.