Axon Diameters and Myelin Content Modulate Microscopic Fractional Anisotropy at Short Diffusion Times in Fixed Rat Spinal Cord

Introduction

Diffusion Magnetic Resonance Imaging (MRI) has become a mainstay of contemporary microstructural imaging in biomedical applications. Diffusion MRI can provide rich information on the sample's microstructure by interrogating micron-scale dimensions within millimeter-scale voxels [1]. In the hierarchical scaling of dimensions in biological systems, the micron-scale is fortuitously a characteristic length scale of many (sub)cellular structures of interest, such as axons, dendrites or cell bodies, which cannot be accessed using routine spatial resolutions in MRI. Most diffusion MRI methods utilize variants of Stejskal and Tanner's [2] Single Diffusion Encoding (SDE) technique [3], which probes diffusion using a single diffusion epoch spanned by diffusion-sensitizing gradient waveforms. The flexibility of SDE in terms of parameter space led to numerous variants [4], as well as diffusion models [5, 6], that have been devised to probe different aspects of the microstructure. For example, Diffusion Tensor Imaging (DTI) models diffusion using a single tensor [7, 8] under the assumption of (time-dependent) Gaussian diffusion, and the tensor's rotationally invariant properties can then report on diffusion anisotropy and parallel/perpendicular diffusivities. Other methods, such as q-space imaging [9, 10] or diffusion spectrum imaging [11] utilize Fourier relationships between the diffusion propagator and signal decay with the q-value (where $\text{q}=\frac{1}{2\pi }\gamma \delta \text{G}$ is the wavevector, γ is the gyromagnetic ratio, δ represents the gradient duration, and |G| is the gradient amplitude) to extract information on pore size or orientation distributions, respectively. Diffusion time- and/or frequency-dependence can also provide much insight into the restricting geometry by probing the way in which the diffusion path is modulated with time and/or the diffusion spectrum, respectively [12–19]. Furthermore, more advanced biophysical modeling has been recently put forth to characterize specific microstructural components such as neurite density [20, 21], or water fractions tentatively associated with axons in white matter [19] from specific acquisition schemes. Such SDE methods have been widely useful in neuroscience [22] and biomedical applications, typically targeting longitudinal processes such as stroke, learning, or chronic disease progression [1].

One interesting metric that can be probed by diffusion is the microscopic diffusion anisotropy (μA) [23–25], from which its normalized counterpart – the microscopic fractional anisotropy (μFA) – can be derived. μFA defines a single compartment's anisotropy in its own eigenframe [26], e.g., for a sphere μFA = 0 while for an infinite cylinder μFA can approach 1. However, in practice, the MRI signal will always originate from an ensemble, thereby making it necessary to account for orientation dispersion within the ensemble [27]. In systems comprising coherently-aligned anisotropic objects where orientation dispersion is ideally zero, μFA would be equivalent to the fractional anisotropy (FA) derived from DTI. However, in conventional SDE methods, when orientation dispersion is significant, estimated FA values typically do not represent the true anisotropy, or μFA, as they are conflated with orientation dispersion [28, 29]. For example, in ideal randomly oriented infinite cylinders, the averaging of anisotropic compartments results in FA = 0, which—without a-priori knowledge or extensive modeling—would suggest that the microscopic geometry is spherical.

In recent years, the Double Diffusion Encoding (DDE) methodology (Figure 1) has been gaining increasing attention for its potential to refine and identify microstructural aspects not so easily probed by SDE [25, 30]. Unlike SDE, DDE probes diffusion correlations using—as its name suggests—two diffusion encoding periods, spanned by two independent gradient wavevectors, which are separated by a mixing time (τ_m). Comparing q-space-like signal decays using parallel and perpendicular relative gradient orientations, Cheng and Cory have been able to measure the sizes of randomly oriented elongated (anisotropic) yeast cells, and distinguish them from spherical cells [24]. Similarly, Callaghan and Komlosh have shown that diffusivities extracted from parallel vs. perpendicular DDE experiments could provide insight into μFA in randomly oriented liquid crystals characterized by Gaussian diffusion [23]. Such measurements provided the first clues that μFA (termed using many divergent terms [3]) could be recovered from DDE irrespective of orientation dispersion.

Mitra [25], and later Özarlsan [31] derived exact solutions for DDE signals, and have identified the importance of the mixing time in decoupling μA from other effects. In the short mixing time regime, interesting diffusion-diffraction phenomena can be produced [32–36], and angular dependencies can provide insight into pore sizes as shown experimentally first by Koch and Finsterbusch [37, 38] and then by others [39–41]; however, by analyzing the displacement correlation tensor [42], the short τ_m angular DDE experiment aiming to measure compartment sizes was found by Jespersen to be equivalent to a time-dependent SDE experiment [43]. By contrast, in the long mixing time regime, the second order term in the displacement correlation tensor, from which sizes are measured, is decoupled from μA, making its measurement much less complicated [25, 31]. The ability to measure accurate μA values was validated in Shemesh et al. [34] and its importance was shown in biological systems such as ex-vivo neural tissues [44], yeast cells [45], and preclinical in-vivo experiments [46], where the orientational variance of the measurements was highlighted. Lawrenz et al. have proposed rotationally invariant schemes for mapping an index of μA [47, 48], and Jespersen et al. subsequently generalized rotationally invariant DDE measurements up to 5th order (in q-values) via a measurement scheme termed DDE 5-design [26]. Numerous promising studies have also been performed on human scanners [37, 38, 48–51], suggesting quite promising potential for disentangling μFA from the underlying orientation dispersion. Additional recent experiments have even extended the DDE methodology toward MR spectroscopy, aiming to impart specificity toward specific cell populations via cellular-specific metabolites [52, 53].

As alluded to above, the diffusion process in biological tissues is highly time-dependent, and thus the filter with which the diffusion experiment is performed can be important. Oscillating Diffusion Encoding (ODE) experiments [14, 54, 55] have been widely used in SDE to enhance contrast in neural tissue, likely since they access shorter diffusion time than could be reached using pulsed-gradient-spin-echo methods [56]. Additionally, ODE has been shown to be highly beneficial for mapping axonal sizes in rat spinal cord [57, 58] as well as for contrasting malignancy in tissues [59, 60]. More recently, the DDE framework was extended toward accommodation of oscillating gradients, termed Double Oscillating Diffusion Encoding (DODE, Figure 1A), first in theory [61], and more recently, in experiment [62]. Importantly, DODE enables the time/frequency-dependence of μFA to be studied. Furthermore, DODE sequences reach the long mixing time regimes much more easily than their DDE counterparts, thereby making the experiments less mixing-time dependent [61], and, as a result, offering the benefit of reduced echo times. This property is likely due to the mixing beginning already from the first gradient pair, and accumulating over the entire gradient waveform. Such DODE experiments were recently reported for the first time in the ex-vivo mouse brain, and μFA maps derived from DODE indeed showed richer contrast than those of their DDE-derived counterparts [62].

Many studies have investigated the underlying microstructural correlates of FA, mainly in white matter (for a classical review, the reader is referred to Beaulieu [63]). It is clear that although myelin strongly modulates FA, it is not necessary for detection of anisotropy in biological systems. Axonal membranes, for example, can impede the diffusion processes with orientational preference and thus can contribute to FA. However, in most studies attempting to investigate the origins of restriction in tissues, orientation dispersion was conflated with SDE-driven metrics; an interesting question is therefore whether μFA, which should not suffer from orientation dispersion effects, could be associated with microstructural features to different extents than FA. The goal of this study was therefore to investigate how μFA and FA correlate underlying microstructural features such as myelin water fraction (MWF) or axonal diameters. As well, we aimed to investigate whether these parameters are differently correlated, and to qualitatively assess the importance of orientation dispersion, especially in the white matter. The final goal of this study was to determine whether μFA is modulated when different length scales are probed via DODE and DDE sequences. A well-characterized system, namely, fixed spinal cord—which has been extensively used in the past to study diffusion [15, 58, 64–66] or relaxation [67–70] microstructural correlates—was used for these investigations. Our findings demonstrate interesting differences in correlations between μFA and FA and MWF, as well as with the a-priori known axonal sizes in white matter, when measured using DODE or DDE. Interesting findings in gray matter tissues are also presented. Implications for D(O)DE contrasts and future routes for investigations of the origin of μFA in neural tissue, are discussed.

Theory

Most DDE-MRI studies up to date have used only a single b-value to extract μFA. However, very recently, Ianus et al. showed that for most plausible microstructural scenarios, μFA obtained in such a way can be highly biased due to neglecting the higher-order terms in the signal decay [62]. Ianus et al. proposed to more accurately estimate μFA in both DDE and DODE methodologies by performing D(O)DE experiments at multiple b-values, and fitting both μA (from which μFA is then calculated) and the higher-order term via polynomial fits. That is, the D(O)DE signal decay at long mixing times can be expanded with b-value as:

where $\mu A^2=\frac{3}{5}\text{var}({\sigma }_i),{\sigma }_{i=1,2,3}$ are the diffusion tensor eigenvalues, S_‖ and S_⊥ represent the D(O)DE signals acquired using parallel and perpendicular gradients, respectively, and P₃ contains the higher-order terms up to third order (even higher-order terms are neglected). Ianus et al. showed that polynomial fitting can be used to estimate μA² and P₃ from Equation 1. When the mean diffusivity (MD) is additionally measured at lower b-values (e.g., from fitting a tensor to the 12 parallel orientations in the 5-design), μFA can be directly calculated from Eq. 2:

Although polynomial fitting probably yields more accurate estimates of μA², it should be noted that ideally, many b-value shells would be required for robust fitting. An alternative approach would be to acquire a much more minimalistic dataset and still be able to quantify μA² and P₃. Setting $\frac{1}{12}\sum log(S_{\Vert }(b))-\frac{1}{60}\sum log(S_{\perp }(b))\equiv \widetilde{ϵ}(b)$, Eq. 1 can be rewritten for two different b-values b₁ and b₂:

It is then straightforward to show that from two measurements at different b-values, μA² can be directly obtained from

which can then be plugged into Equation 2 to obtain μFA directly. Note that we use the tilde to distinguish the extracted $\widetilde{\mu A^2}$ from the real μA². This approach for accurate μFA extraction thus requires, in principle, only two measurements, one at low b-value, from which MD and $\widetilde{ϵ}(b_1)$ would be obtained, and another at higher b-value, where $\widetilde{ϵ}(b_2)$ would be obtained. However, since at low b-values required for accurate estimation of MD, $\widetilde{ϵ}(b_1)$ may be very small and comparable to noise levels, it is more appropriate to acquire $\widetilde{ϵ}(b_1)$ and $\widetilde{ϵ}(b_2)$ at somewhat higher b-values (where the b² terms are more dominant) and perform a separate, third acquisition for extracting MD at lower b-values. This 3-shell approach was thus preferred in this study.

Materials and Methods

This study was carried out in accordance with the recommendations of the directive 2010/63/EU of the European Parliament of the Council, authorized by the Champalimaud Centre for the Unknown's Animal Welfare Body, and approved by the national competent authority (Direcção Geral de Alimentação e Veterinária, DGAV).

Specimen Preparation

Spinal cord specimens were obtained from adult male Wistar rats (N = 2) weighing ~300 gr. The rats underwent standard transcardial perfusion under deep pentobarbital anesthesia. Cervical spinal cords were extracted, washed in PBS, and kept in 4% paraformaldehyde (PFA) for 24 h at 4°C. The samples were then placed in freshly prepared phosphate buffer saline (PBS) for at least 48 h prior to MRI experiments. The samples were cut to ~1 cm segments and placed in a 5 mm NMR tube filled with fluorinert (Sigma Aldrich, Lisbon, Pt).

MRI Experiments

All MRI experiments were performed on a vertical 16.4T (700 MHz ¹H frequency) Aeon Ascend scanner (Bruker, Karlsruhe, Germany) interfaced with a Bruker AVANCE IIIHD console. A Micro5 probe equipped with a 5 mm birdcage coil for transmit and receive functions and a gradient system capable of producing amplitudes of up to 3T/m isotropically was used. The sample was kept at a constant temperature of 23°C throughout the experiments by means of air flow, and the samples were allowed to equilibrate with the surrounding temperature for at least 4 h before acquiring any diffusion or relaxation experiments.

All diffusion sequences were written in-house and were based on an Echo Planar Imaging (EPI) readout. For both DODE and DDE, the same acquisition parameters were used, namely, two-shot and double-sampled EPI with a readout bandwidth of 555.555 kHz, Field of View (FOV) of 6 × 4 mm² and in-plane matrix size of 120 × 80, leading to an isotropic in-plane resolution of 50 × 50 μm². The slice thickness was 500 μm, and TR/TE = 2,500/52 ms. For both DODE and DDE acquisitions, Jespersen's 5-design sampling scheme [26] was used for the diffusion weighted images, and, additionally, eight images with zero b-value were acquired, such that the total number of images acquired in a given scan was 80. For both DODE and DDE, three separate acquisitions were performed with different b-values, namely, 2b = 1.2, 2.4 and 3.0 ms/μm² (where the factor of 2 reflects the accumulated diffusion weighting along the two diffusion epochs). The specific b-values were chosen based on signal-to-noise and contrast considerations: on the one hand, they have to be sufficiently low such that even higher-order terms do not contribute, but on the other hand, they have to be high enough for μFA contrast to be detectable. The lowest b-value scans were acquired with 12 averages, while the other two b-value shells were acquired with 32 averages each. The DODE diffusion parameters were: T_DODE = 13 ms, N = 5, τ_s = 2 ms. The DDE diffusion parameters were Δ/δ = 12/1 ms, τ_m = 12 ms, see Figure 1 for definitions of the parameters.

Additional experiments were performed for mapping myelin water fraction. Those consisted of a Carr-Purcell-Meiboom-Gill (CPMG)-based acquisition performed using a modified pulse multi-slice-multi-echo (MSME) sequence. The same slice was acquired as in the diffusion images with identical in-plane resolution and FOV. The acquisition bandwidth for the pulse sequence was 100 kHz, and the pulses used for slice-selective excitation and refocusing had durations of 1.16 ms (Shinnar-Le-Roux design) and 50 μs (Gaussian shape), respectively. The respective bandwidths of the excitation and refocusing pulses were 3625 and 32,100 Hz, respectively, such that the refocusing pulse provided complete refocusing on the entire slice. The ΔTE that could be achieved using these parameters was 2.85 ms, and 96 echoes were acquired from 2.85 to 273.6 ms. The repetition time was 2500 ms and two averages were acquired.

Diffusion Data Preprocessing

All preprocessing and analyses were performed using MatLab^® (The MathWorks, Inc., Natick, Massachusetts, United States). Raw images were registered using an implementation of Guizar-Sicairos et al. [71] found in https://goo.gl/3bGU8b. The images were then denoised using Veraart's method based on Marchenko-Pastur distributions in Principal Component Analysis of redundant data [72]. Gibbs unringing was performed using Kellner's method [73] implemented in Matlab. Finally, the denoised and unrung images were very slightly smoothed using a [2 2] median filter.

Relaxation Data Preprocessing

The preprocessing steps for the relaxation data were identical to the diffusion data preprocessing steps, except for an additional step in the very beginning of the pipeline whereby the magnitude data was converted to real data using Eichner's method [74]. All steps listed above including denoising, unringing and median filter smoothing were then executed in sequence.

Diffusion Data Analysis

The first analysis step for D(O)DE data was to fit the diffusion tensor. Diffusivities were computed using a simple linear fitting of S_‖ data acquired at the lowest b-value experiments followed by diagonalization and extraction of the diffusion tensor eignevalues and eigenvectors. The mean diffusivity and fractional anisotropy were then calculated from the tensor eigenvalues as $MD=\frac{1}{3}({\lambda }_1+{\lambda }_2+{\lambda }_3)$ and $FA=\sqrt{\frac{3}{2}\frac{{\left({\lambda }_1-MD\right)}^2+{({\lambda }_2-MD)}^2+{({\lambda }_3-MD)}^2}{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}}$, where λ_i represent the tensor eigenvalues.

The second step in the analysis was to use the data from the two higher b-values to extract μFA. First, $\widetilde{\mu A^2}$ was extracted directly from Eq. 4; the mean diffusivity estimate was then used along with the extracted $\widetilde{\mu A^2}$ to obtain μFA via Equation 2.

Relaxation Data Analysis

Following the preprocessing steps listed above, the filtered relaxation data were subject to a voxelwise inverse Laplace Transform (iLT) using 150 T₂ components log-spaced between 2.1 and 328.3 ms. The T₂ spectra were smoothed by minimum-curvature constraint as in Dula et al. [75] and extended phase graph analysis was performed to account for any ${\text{B}}_1^{+}$ inhomogeneity and ensuing stimulated echoes [76]. The myelin water fraction (MWF) was computed from each spectrum as the fraction of signal originating from components with peak T₂ smaller than 17 ms. ROIs were drawn manually on the raw data closely following Dula et al. [75], and the ROI data underwent the same analysis using the mean signal decay in each ROI.

Statistical Analysis

Gray matter and white matter masks were created by thresholding MWF maps with MWF < 0.22 for gray matter and MWF > 0.25 for white matter. The histograms in Figure 4 were then generated for each metric/method using Matlab's histogram function which automatically selects the bin width to represent the underlying distribution in the most accurate way. Parameter means and standard deviations are reported in the text and Tables.

Correlation analyses between different diffusion metrics were performed using automatic outlier rejection (Grubbs test for outliers) followed by calculation of Spearman's ρ (μFA and FA data from all methods were not normally distributed). An analysis of variance (ANOVA) was performed to compare μFA and FA arising from DODE and DDE methods, with post-hoc Bonferroni tests corrected for multiple comparisons.

To correlate MWF with μFA or FA extracted from the different methods, the diffusion maps were registered to the MWF using Matlab's imregister function using a multimodal configuration, initial radius of 1e-5, maximum number of iterations = 1,000, and allowing for affine transformations due to the small differences in image geometry arising from EPI-based (diffusion) and line-by-line (relaxation) acquisitions.

When linear fits are presented (Figure 8), Matlab's robustfit function was used to extract the coefficients.

Results

Diffusion data quality can be appraised in Figure 2, which plots representative raw data from one of the spinal cords, obtained from experiments with zero b-value (Figure 2A), parallel (Figure 2B), and perpendicular (Figure 2C) diffusion orientations at the highest b-value used in this study. Before denoising, the worst-case signal to noise ratio (SNR)—measured at the highest b-value and with significant diffusion weighting gradients in the direction parallel to the spinal cord's principal axis—was ~20 in white matter. The middle column in Figure 2 shows the corresponding preprocessed data and the ensuing enhancement of image quality from denoising and Gibbs unringing (Figures 2D–F). Figures 2G–I show the result of subtracting raw and denoised images. The lack of structure in the subtracted images suggest that indeed only noise was removed and that no significant signal components were lost during denoising [72]. The SNR of the preprocessed images was enhanced by a factor of ~2.

To assess the different maps obtained in this study, representative μFA and FA maps derived from DODE as well as DDE experiments (hereafter referred to as μFA_DODE and μFA_DDE or FA_DODE and FA_DDE, respectively) are shown in Figure 3. Several interesting qualitative features can be highlighted from these images: (1) both μFA_DODE and μFA_DDE maps (Figures 3A,C) have higher values than their FA_DODE and FA_DDE counterparts (Figures 3B,D) in white matter, as well as in gray matter; (2) μFA_DDE is higher and less tract-specific when compared with μFA_DODE (for approximate definitions of tract locations and spinal cord anatomy, the reader is referred to Figure 3E); (3) μFA_DDE appears quite homogeneous in the WM while μFA_DODE shows more variation within WM; (4) similarly, FA_DDE is more homogeneous in white matter compared with FA_DODE, which shows a greater variance in different tracts. To provide a more quantitative view on these features, Figure 4 plots histograms of μFA and FA in white matter and gray matter (c.f. Figures 4A,B for the ROI masks). In white matter, μFA_DODE is higher than its FA_DODE counterpart (Figure 4C), while in gray matter, μFA_DODE is distributed at much higher values compared to FA_DODE (Figure 4D). Similar trends were observed for DDE but with μFA or FA shifted toward somewhat higher values (Figures 4E,F).

It is also interesting to compare differences between methods within the same tissue type (e.g., comparing same-color distributions down the columns of Figure 4). μFA_DODE is clearly lower and more widely distributed compared with μFA_DDE in white matter. In gray matter, μFA_DDE is high, while μFA_DODE is somewhat smaller. Another interesting finding in gray matter, is that FA_DODE and FA_DDE values are only slightly different. The means and standard deviations of μFA and FA for each method are tabulated in Table 1.

A statistical analysis of these data is given in Figure 5, which presents box plots of the data. A one-way ANOVA revealed that in each tissue type (e.g., white matter or gray matter), all four metrics are highly statistically significantly different from each other (corrected p < 1e-12, post-hoc Bonferroni test). However, it should be noted that although the metrics are different, they are not completely uncorrelated. Table 1 reports Spearman's ρ and its significance levels when comparing μFA and FA (extracted by the same method) in each ROI. While μFA_DODE and FA_DODE are correlated in white matter (Spearman's ρ = ~0.41), μFA_DDE and FA_DDE metrics are only weakly correlated (Spearman's ρ = ~0.19). In gray matter, the correlations between μFA and FA were weak for both methods and (Spearman's ρ = 0.22 and −0.10 for DODE and DDE, respectively). Note that although outlier rejection was used, in all cases < ~1% of the data were identified as outliers and rejected.

To establish whether and how myelin modulates the anisotropy metrics, Carr-Purcell-Meiboom-Gill (CPMG) MRI experiments were performed on the same slice with the same resolution as the diffusion experiments. To assess the quality of the data, Figures 6A,B show the preprocessed data at short and very short TE of 2.9 ms and very long TE of 142.5 ms, respectively, in a representative spinal cord. Even at the very long TE, the SNR remains very high, especially after denoising. Denoising and unringing procedures were validated and found to have no negative impact on the quality of T₂ fitting procedure (data not shown), while improving the fits significantly. Figure 6C shows ROIs drawn in the major tracts of the spinal cord, while Figures 6D,E show the T₂ decays (with the ordinate drawn in log scale) and the resultant T₂ spectra (with the abscissa drawn in linear scale), respectively. The decays in white matter are clearly non-linear, and the myelin water can be seen as an early peak in the T₂ spectrum with its peak T₂ around ~10 ms.

A representative myelin water fraction (MWF) map arising from pixel-by-pixel quantification of the spectra is shown in Figure 7A. Note the sharp contrast between the different tracts in MWF: for example, the dCST shows the lowest MWF (MWF~0.30) while VST and FC exhibit the highest MWF (MWF~0.45). Scatter plots between MWF and μFA or FA in white matter are shown in Figure 7 for DODE (Figure 7B) and DDE (Figure 7C), respectively. Table 2 summarizes the correlation coefficients and associated statistics. A moderate anticorrelation between MWF and μFA_DODE is observed in the white matter (Spearman's ρ = ~ −0.36), while FA_DODE did not correlate with MWF in a statistically significant manner. The DDE counterparts μFA_DDE and FA_DDE exhibited weak anti-correlation and correlation, respectively. Figures 7D,E show similar plots as described above, but for gray matter. Notably, correlations between MWF and FA_DODE, as well as FA_DDE were very weak and their statistical significance not very high; on the contrary, μFA_DODE was found to correlate somewhat with MWF, while μFA_DDE correlated moderately with MWF, with very high statistical significance (c.f. Table 2).

Finally, the correlation of the mean μFA in the different tracts with literature regional averaged axon diameter was assessed. Figures 8A,B plot mean μFA and FA against the axon diameters reported in Dula et al. [75] for the different spinal cord tracts. These data, along with the values tabulated in Table 3, demonstrate that μFA_DODE exhibits very strong anticorrelation with axon diameters (Spearman's ρ = −0.96, p = 0.0028). All other metrics are not significantly correlated with axon diameter.

Discussion

μFA has been recently gaining increasing attention as a potentially useful source of contrast in microstructural MRI due to its ability to disentangle anisotropy from orientation dispersion. Methods other than D(O)DE, targeting μFA such as tailoring b-tensor shapes are emerging, with many potential applications [77–80]. However, such methods may be confounded by time-dependent diffusion effects [27, 81–83], whereas D(O)DE at long mixing times naturally avoids these confounds [43]. It is therefore imperative to investigate how μFA may be correlated with underlying microstructural features such as axon dimensions and myelin, much like the early studies aiming to understand the sources for FA [63, 67, 84, 85]. In general, perhaps the most significant findings of prior studies on FA (conducted nearly invariably with SDE) were that (1) anisotropy in white matter depends on axonal membranes; and (2) the presence of myelin can further modulate FA metrics [63]. The application of oscillating gradients has also been shown to generate more contrast and more accurate estimations of small dimensions as compared to long diffusion time experiments, presumably due to the more efficient probing of smaller dimensions via the shorter diffusion times [58, 86, 87].

The present study aimed to investigate how μFA differs from FA in terms of correlations with myelin water and axonal diameters, and to compare those metrics when measured with DDE or DODE sequences. We first focus attention to our results arising from white matter tissue. Notably, μFA was always larger than FA (Figures 3–5 and Table 2), in agreement with previous DDE experiments in fixed tissues [26] and in-vivo [88]. Since the μFA and FA metrics were extracted from the same acquisition, it is unlikely that other effects such as exchange or relaxation contributed to μFA > FA. Thus, our finding supports the notion that that orientation dispersion is significant even in highly structured tissues, such as spinal cord white matter. This is in excellent agreement with a recent study of SDE-derived diffusion tensor and kurtosis time-dependencies which also pointed to the same conclusion in pig spinal cord [15], as well as with histological studies attempting to measure the dispersion directly in white matter [89]. It is difficult to draw conclusions on whether the orientation dispersion arises within intra- or extra-axonal spaces (or both), or, whether undulations [90] or passing collateral fibers [91] can contribute to these observations. Performing similar spectroscopic measurements utilizing cell-specific markers such as NAA or mI [52, 53], or performing much more extensive time/frequency/b-value-dependent measurements on water [19, 59, 92], or on metabolites [93, 94] may further assist in addressing this question in the future.

Another interesting aspect when comparing μFA with FA in white matter, is that the two metrics are only moderately correlated when measured with DODE, and very weakly correlated when measured with DDE (c.f. Table 1). This finding suggests that when diffusion is encoded using oscillating gradients, spins experience less orientation dispersion than when they are probed using block gradients, since μFA would be perfectly correlated (and identical) to FA for perfectly aligned fibers. Hence, our findings point to specific length scales for orientation dispersion that are probed differently using the different sequences.

Next, we consider the relationships between myelin and μFA. Akin to its FA counterpart–μFA is ambiguous in that a compartment with length “L” and radius “R” can give rise to the same μFA as a compartment with length 2L and radius 2R. The axial path length could be restricted due to nodes of Ranvier, non-ideal cylindrical structure, varicosities, etc. However, if the path length parallel to the (assumingly) ellipsoids is constant, then one could predict that when larger amounts of myelin surround an axon, the μFA will be smaller as the restriction will increase in the perpendicular direction. However, in our study, a moderate negative correlation was observed between MWF and μFA_DODE in white matter (Figure 7 and Table 2). This can be explained by considering the dependence of MWF and axon diameter via the g-ratio [95]: the larger the axon, the thicker the myelin around it in (healthy) mammalian white matter [96]. Hence, the negative correlation between μFA_DODE and MWF would reflect indirectly the approximately constant g-ratio in healthy tissue, rather than enhanced restriction. Interestingly, μFA_DDE showed a much weaker, yet still negative correlation with MWF. Since the microstructure has not changed between measurements, this likely reflects that DODE and DDE probe different path lengths parallel to the spinal cord's major axis: the larger the diffusion time, the longer path will be probed in the unrestricted dimension, and thence the μFA will be larger and less reflective of axon diameter or, by proxy, its myelin thickness. FA_DDE showed a small positive correlation with MWF, which perhaps reflects the ambiguity of probing restriction and orientation distribution at the same time. Extracellular space contributions again cannot be neglected here, but for coherently aligned systems the arguments are similar as one could potentially treat the space between densely packed axons as potentially even more restricted than the intra-axonal space itself [45]. It is also worth mentioning that MWF extracted from multiexponential T₂ measurements, as performed in this study, have been shown in the past to reflect microstructural metrics such as axon size and myelin thickness very faithfully in white matter [67, 75, 97].

Our most striking findings in this study, perhaps, is that μFA_DODE showed an extremely high, and statistically significant, negative correlation with axon diameters reported by Dula et al. [75] and Harkins et al. [98] for the different tracts (Table 3). This observation lends further credence to the explanation above: the finite parallel length scale probed by DODE makes the measurement strongly dependent on the perpendicular restriction, which in this case is reflected through axon sizes. Although the axon diameters were obtained from literature, it is worth stressing that axon diameter dependence in healthy spinal cords is highly reproducible and that the tracts analyzed were obtained from very similar cervical slices as in Dula et al. [75]. Such a strong correlation is also highly unlikely to be obtained randomly. It is very interesting to also note that all other metrics did not correlate in a statistically significant fashion with axon diameters: μFA_DDE likely due to its probing of longer parallel lengths, and the FA from both methods due to its inherent conflation or restriction with orientation dispersion.

In the spinal cord gray matter, very low FA_DODE and FA_DDE values were measured, suggesting a much lower degree of restriction compared to white matter diffusion. However, the μFA_DODE and μFA_DDE metrics in gray matter were still very high in the gray matter. In fact, the values reported in Table 1 also reflect $\frac{\mu F{A}_{DODE}^{GM}}{\mu F{A}_{DODE}^{WM}}=~0.92$ and $\frac{\mu F{A}_{DDE}^{GM}}{\mu F{A}_{DDE}^{WM}}=~0.89$. Combined with the low FA values in the gray matter, our findings suggest that a significant component of gray matter tissue experiences restricted diffusion but with a large degree of orientation dispersion. This finding is also in good agreement with previous literature demonstrating significant angular DDE modulations in ex-vivo gray matter [44, 45]. Further studies are needed to establish which underlying biological components give rise to such high μFA in gray matter, but dendrites, astrocyte branches, and nonmyelinated or myelinated axons traversing gray matter could be suspected [93, 99]. Time-dependent or spectroscopic experiments on metabolites could provide insight into such questions in the future.

Several limitations can be identified in this study. First, we have introduced a new way of measuring $\widetilde{\mu A^2}$ harnessing the 5-design acquisition at two b-values to reduce the recently-reported bias in μA² estimation due to higher order terms. Our new method is likely inferior to a sampling of a large range of b-values and the ensuing polynomial fitting as done in Ianus et al. [62]. However, the advantage of the current approach is that it manages to avoid a prohibitively long experiment duration. Future studies will identify the accuracy and precision of the method proposed above vis-à-vis the ground-truth, and attempt to find optimal b-values for measuring $\widetilde{\mu A^2}$ as accurately and with as little bias as possible. Second, to compute μFA, we executed a third measurement at lower b-value to extract MD, which is then input into Eq. 2 along with $\widetilde{\mu A^2}$. However, MD itself may be conflated with higher-order terms, as pointed out recently by Chuhutin et al. [100]; in this study, this issue was not accounted for, and may induce minor biases in the measurements of μFA. Better estimation of MD could probably be performed by sampling one or more low b-values and fitting kurtosis and MD at the same time from spherically averaged data. In addition, we have not explored the impact of specific b-value selection. At too low b-values, the difference in the log signals is very small, while at higher b-values, even higher-order terms may come into play. Third, the sample size was quite small (N = 2 spinal cords, only a single slice per cord), such that the variability across animals was not very well sampled. However, it is worth noting that the results were actually very consistent between both spinal cords: the mean μFA and FA, for both DODE and DDE, varied <10% between the cords (both in gray and white matter tissues), and the MWF varied <6% between the tissues. Although this consistency is promising for the robustness of the approach, the small number of samples renders this study perhaps more exploratory. Fourth, the experiments were performed at a relatively long TE of 52 ms. Given that the MWF was associated with T₂ < 20 ms and that the other water T₂s were distributed between ~20 and 60 ms, the diffusion experiments can be considered completely filtered for (directly contributing) myelin water, as $e^{\frac{-TE}{T_{2myelin}}}~$ 0.005. Exchange between myelin water and intra/extra-axonal water is likely to occur, which may also confound the measurements, although it should be noted that at least for conventional DODE MRI, the relatively long TE is nearly unavoidable due to the necessity of non-negligible diffusion gradient waveform durations. Double-stimulated-echo approaches [101, 102] would thus be nearly impossible to execute for DODE, even before considering the significant SNR reduction associated with such sequences, (1/2)^N, where N is the number of stimulated echoes. Finally, a histological study was not here performed, and the study relies on literature reports of correlations between MRI-derived MWF and myelin thickness and the values for axon diameters. Future studies can expand the findings here and perform more direct correlations with histology, although it should be pointed out that big differences in these parameters are unlikely to be observed for healthy tissues. In addition, it would be fruitful to modulate the microstructure actively and to observe how μFA varies, e.g., using genetic mutations that alter myelin content. All these highly interesting avenues will be pursued in the future, but the present study provides the first steps in this direction.

Conclusions

This study investigated the microstructural correlates of μFA and FA using high resolution D(O)DE experiments in fixed spinal cords at 16.4 T. Our results indicate very strong anticorrelations of μFA_DODE with axon size, and moderate anticorrelations of μFA_DODE with MWF, whereas μFA_DDE, FA_DODE and FA_DDE correlate to a much lesser or no extent with those microstructural features. These findings shed light on the mechanisms of restriction in spinal cord white matter when investigate without conflation by orientation dispersion. The correlations of μFA_DODE with axon diameters and myelin water fraction are thus promising for future investigations of longitudinal variations in these properties, e.g., in disease or with learning.

Author Contributions

NS designed the study, collected and analyzed data, and wrote the paper.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This study was funded in part by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 679058-DIRECT-fMRI). The author would like to thank Prof. Sune N. Jespersen and Mr. Jonas Lynge Olesen (Aarhus University) for providing the codes for denoising and Gibbs unringing, as well as for many insightful discussions. The author also thanks Dr. Daniel Nunes for extracting the tissues used in this study, Dr. Andrada Ianus and Ms. Teresa Serradas Duarte for providing parts of code used in the analyses performed here, and Prof. Mark D. Does from Vanderbilt University for the REMMI pulse sequence and its analysis tools, that were supported through grant number NIH EB019980.

References