Strong diffusion gradients allow the separation of intra- and extra-axonal gradient-echo signals in the human brain.

The quantification of brain white matter properties is a key area of application of Magnetic Resonance Imaging (MRI), with much effort focused on using MR techniques to quantify tissue microstructure. While diffusion MRI probes white matter (WM) microstructure by characterising the sensitivity of Brownian motion of water molecules to anisotropic structures, susceptibility-based techniques probe the tissue microstructure by observing the effect of interaction between the tissue and the magnetic field. Here, we unify these two complementary approaches by combining ultra-strong (300mT/m) gradients with a novel Diffusion-Filtered Asymmetric Spin Echo (D-FASE) technique. Using D-FASE we can separately assess the evolution of the intra- and extra-axonal signals under the action of susceptibility effects, revealing differences in the behaviour in different fibre tracts. We observed that the effective relaxation rate of the ASE signal in the corpus callosum decreases with increasing b-value in all subjects (from 17.1±0.7s-1 at b=0s/mm2 to 14.6±0.7s-1 at b=4800s/mm2), while this dependence on b in the corticospinal tract is less pronounced (from 12.0±1.1s-1 at b=0s/mm2 to 10.7±0.5s-1 at b=4800s/mm2). Voxelwise analysis of the signal evolution with respect to b-factor and acquisition delay using a microscopic model demonstrated differences in gradient echo signal evolution between the intra- and extra-axonal pools.

Introduction

Methods that are sensitive to inter-individual differences in human brain microstructure in vivo hold great promise for understanding the function of the healthy brain and the pathophysiology of brain development or degeneration. While in the past characterisation of tissue properties has relied on postmortem, histological examination (Aboitiz et al., 1992), the advent of Magnetic Resonance Imaging (MRI) has had a transformative impact on research by allowing non-invasive assessment of tissue in vivo at different stages of brain development and degeneration. Although MRI research is dominated by examination of the bulk properties of tissue, in recent years the field of microstructural MRI has blossomed giving rise to “virtual in vivo histology”.

The characterisation of brain microstructure using MRI remains challenging, despite progress recently made in this field. At the core of this challenge lies the fact that the typical size of an MRI image voxel is about a thousand times larger than the diameter of an average neuron. This means that the bulk signal from each imaging voxel is a superposition of the signals from the inside and outside of many cells. While the signal decay of each individual compartment can often be approximated by a mono-exponential, this rarely holds for the sum of the signals from the different compartments. Indeed, a large number of studies have demonstrated that the MR signal variation from brain white matter (WM) as a function of a single experimental variable (e.g. echo time, inversion time, or the strength of the diffusion-weighting) cannot be accurately described by a single exponential (MacKay et al., 1994; Whittall et al., 1997; Assaf and Cohen, 2000; van Gelderen et al., 2012; Wharton and Bowtell, 2012; Sati et al., 2013; Oh et al., 2013; Labadie et al., 2014; Kim et al., 2015). For example, diffusion MRI (dMRI) reveals two diffusing components: the intra-axonal signal (restricted space) experiences a much weaker signal attenuation compared to the extra-axonal signal, where diffusion is assumed to be largely hindered and the apparent diffusivity is comparatively larger than in the intra-axonal space (Assaf and Cohen, 2000; Assaf et al., 2004). On the other hand, various relaxometry studies have demonstrated that the myelin water shows a much faster transverse (R2) (MacKay et al., 1994; Whittall et al., 1997), effective transverse (van Gelderen et al., 2012; Wharton and Bow-tell, 2012; Sati et al., 2013) and longitudinal (R1) (Oh et al., 2013; Labadie et al., 2014; Kim et al., 2015) relaxation rate than the intra- and extra-axonal signals. For example Oh et al. (2013) used a double inversion recovery in conjunction with a gradient echo acquisition to suppress long T1-signals so as to probe the evolution of the myelin water signal with echo time. Finally, tissue placed in a strong magnetic field will interact with it, producing small changes in this field at a microscopic level (Lee et al., 2010). This interaction is characterised by the magnetic susceptibility. Recent work has shown that the magnetic susceptibility properties of the myelin sheath produce fibre-orientation-dependent differences in the nuclear magnetic resonance (NMR) frequencies and relaxation rates of the intra- and extra-axonal signals in white matter (Wharton and Bowtell, 2012; Sati et al., 2013; van Gelderen et al., 2012). Experimental evidence for this behaviour is based on analysis of subtle effects on the evolution of multi-echo gradient echo signals at late echo times, when the signal from the myelin water has decayed away (Wharton and Bowtell, 2012; Sati et al., 2013; Nam et al., 2015; Thapaliya et al., 2018; Tendler and Bowtell, 2019; Alonso-Ortiz et al., 2018).

The advantage of varying more than one experimental parameter (multi-dimensional MRI) to improve the estimation of NMR responses has long been recognised (Peemoeller and Pintar, 1969; Bernin and Topgaard, 2013). Variation of the spin-echo time in combination with diffusion weighting has proven to be effective in separating the R2 relaxation rates of the intra- and extra-axonal signals in vivo (Tax et al., 2017; Veraart et al., 2018; McKinnon and Jensen, 2019; Lampinen et al., 2019).

In this study we report a technique which can be used to disentangle the gradient echo signals from the intra- and extra-axonal compartments. This involved the development and implementation of a diffusion-filtered, asymmetric spin echo (D-FASE) EPI sequence, which we subsequently used to measure the evolution of the signal from white matter under the simultaneous action of diffusion filtering and dephasing due to local frequency variations. By applying the diffusion weighting perpendicular to the WM fibres we can strongly suppress the extra-axonal signal while maintaining intra-axonal signal: in other words we can filter out the extra-axonal signal by means of a diffusion filter. By exploiting the very latest in ultra-strong magnetic field gradient technology using a Connectom scanner (Setsompop et al., 2013; Jones et al., 2018), it was possible to produce the level of filtering needed to differentially attenuate the signal from the two compartments, while maintaining reasonable spin echo times and thus a sufficient signal-to-noise ratio (SNR).

Methods and materials

Data acquisition

The D-FASE sequence was created by adding diffusion gradients to an asymmetric spin echo sequence (Stables et al., 1998; Streicher et al., 2014; Hutter et al., 2018), yielding images that are sensitive to dephasing due to both diffusion and microscopic susceptibility effects. Setting the acquisition time delay to zero results in a standard Stejskal-Tanner diffusion-weighted spin echo sequence (Stejskal and Tanner, 1965), while setting the gradient amplitude to zero and the acquisition time delay to non-zero, yields a conventional asymmetric spin-echo sequence (Stables et al., 1998).

Data were collected from four healthy subjects (one female, age: 29 ± 2 years) on a 3 T MR system equipped with Connectom gradients with a maximum gradient strength of 300 mT/m (Connectom, Siemens Healthcare, Erlangen, Germany) and 32 channel receive array head coil (Nova Medical, Wilmington, United States). The study was performed with ethics approval from the Cardiff University School of Psychology ethics review board and written informed consent was obtained from all participants. Diffusion gradients were applied along the superiorinferior (SI), anterior-posterior (AP) and left-right (LR) axes with b-values = [1200, 2400, 4800] s/mm2 at each of 6 acquisition delay times Δt = [0, 10, 20, 30, 40, 50] ms. Diffusion scans acquired for each of the three diffusion-weighting orientations were repeated 20 times and interleaved with b0 scans. 19 sagittal slices with 2 mm isotropic resolution were acquired using the above D-FASE sequence parameters. Other acquisition parameters included spin echo time TESE = 57 ms, repetition time TR = 2.9 s, bandwidth 1976 Hz/Px, and parallel imaging approach GRAPPA (Griswold et al., 2002) with acceleration factor of 2 in the phase encoding direction. Fig. 1 shows the D-FASE sequence and a summary of the acquisition parameters. The protocol additionally included a single-shell scan with 64 diffusion-weighting directions (b = 3000 s/mm2, TESE = 57 ms, TR = 2.9 s, bandwidth 1976 Hz/Px) for estimation of fibre-orientation to the magnetic field and ROI calculation.

Fig. 1

Data acquision: A. Diffusion-filtered asymmetric spin echo (D-FASE) sequence. In a diffusion-weighted spin echo sequence, acquisition is delayed by Δt relative to the spin echo to explore dephasing due to both diffusion and susceptibility effects. B. Scanning parameters. For each Δt, diffusion gradients were applied along each of the main radiological axes (SI, AP, and LR). At non-zero b-values the measurements were performed 20 times to assess the spread of data. Deployment of the ultrastrong diffusion gradients allowed us to keep the spin echo time short (TESE = 57 ms) and therefore to suppress the extra-axonal signal while maintaining a high signal-to-noise ratio.

Data processing

The k-space data for all datasets were reconstructed using the GRAPPA-kernel (Griswold et al., 2002) estimated from the b = 0 s/mm2 data without acquisition delay (Δt = 0 ms) to limit any inconsistencies between the reconstructed images due to differing SNR (Ding et al., 2015). The complex data from each coil element k were combined using the coil elements’ sensitivities m and phases ϕ as follows (Roemer et al., 1002; Bydder et al., 2002): where m and ϕ from the kth coil element were approximated from the data acquired at b = 0 s/mm2 Δt = 0 ms:

The reconstructed D-FASE data were corrected for misalignment due to subject motion and eddy currents (Andersson and Sotiropoulos, 2016), Gibbs ringing (Kellner et al., 2016) and outlier detection (Sairanen et al., 2018). An example of reconstructed and corrected magnitude data is shown in Fig. 2.

Fig. 2

The diffusion-filtered asymmetric spinecho (D-FASE) signal was acquired by simultaneously varying the strength of the diffusion filtering, the b-value (images from left to right), and by delaying the acquisition by the time Δt (images from top to bottom). While at time Δt = 0 ms (the first row at the top), the signal is equivalent to a standard diffusion-weighted spinecho signal, the signal at b = 0 (the first column on the left) is equivalent to an asymmetric spinecho signal. The D-FASE data shown here are from a mid-sagittal slice and were acquired with diffusion filtering applied along the anterior-posterior (AP) direction, hence the signal intensity is highest in the corpus callosum in the images acquired at high b-values, where white matter fibres are predominantly oriented left-right. The logarithm of the data is displayed.

Multi-directional diffusion-weighted data (64 diffusion-gradient directions, b = 3000 s/mm2) were processed with MRtrix3 (Tournier et al., 2019) to identify peaks in the fibre orientation distribution (Dhollander et al., 2016, 2018; Tournier et al., 2004, 2007). A single fibre mask was estimated by selecting those voxels in which the amplitude of the second peak was less than 5% of that of the first peak (Tax et al., 2014). Three white matter regions-of-interest (ROIs) with the dominant fibre orientation within 15° of the main radiological axes were selected: along the right-left (RL) axis (mostly corpus callosum, ), along the anterior-posterior (AP) axis (mostly cingulum, ), and along the superior-inferior (SI) axis (mostly corticospinal tract, ). Finally, a more conservative threshold of 0.5 was applied to the amplitude of the first peak to reduce partial volume effects from cerebrospinal fluid in the ROIs. Fig. 3 shows ROIs superimposed over the peaks data in a subset of slices for each subjects.

Fig. 3

The ROIs in the corpus callosum (red, LR), the corticospinal tract (blue, SI) and cingulum (green, AP) are highlighted on a single-fibre-orientation colormap for a subset of slices. The peaks in fibre orientation distribution and the regions of interest were calculated from the multi-directional diffusion-weighted data (64 diffusion-gradient directions, b = 3000 s/mm2).

For each b-value a mono-exponential function was fitted to the data from the corpus callosum, the cingulum and the corticospinal tract when the diffusion weighting was applied perpendicular to the main fibre direction. This yielded the initial characterisation of the signal evolution with acquisition delay Δt and the average relaxation rate as a function of b. Microscopic signal analysis was then performed on data from the corpus callosum.

Microscopic signal analysis

We modelled the D-FASE signal evolution as a superposition of intra- and extra-axonal signals and as a function of two parameters: the acquisition delay Δt and the strength of the diffusion filtering b. The difference in the susceptibility effects experienced by the intra- and extra-axonal compartments is manifested in differences in the average resonance frequency and effective transverse relaxation rate of the intra- and extra-axonal signals (van Gelderen et al., 2012; Wharton and Bowtell, 2012; Sati et al., 2013). The non-zero difference in frequency is an essential indicator of the anisotropic susceptibility of the myelin sheath (Wharton and Bowtell, 2012). Moreover, under the action of diffusion filtering perpendicular to WM fibres, the extra-axonal signal decays at a faster rate (Dfast) with b than the intra-axonal signal. The apparent diffusivity of the intra-axonal compartment (Dslow) is assumed to be non-zero, predominantly due to orientational dispersion and the diffusion gradients not being perfectly perpendicular to the main fibre direction. Therefore, the magnitude of the D-FASE signal can be modelled according to Equation (3):

Here, is the effective transverse relaxation rate of the intra-axonal signal, is the additional effective transverse relaxation rate of the extra-axonal signal, and ω is the difference in angular frequency of the two signals. Three cases were considered: i) variable ω and , ii) variable and ω = 0 Hz, iii) variable ω and . For each case, the Bayesian information criterion (BIC) (Schwarz, 1978) was calculated.

Nonlinear least-squares fitting was performed using the Matlab lsqnonlin function; fitting parameters and parameters bounds are listed in Table 1.

Table 1

The absolute value of the signal from the two-compartment model in Equation (3) was fitted to the magnitude data from the corpus callosum for all Δt and b-values for AP and SI diffusion directions, simultaneously. Fitting parameters and parameter bounds for fitting the 2-pool, 2-parameter model (Equation (3)) to D-FASE data. A1 and A2 are signal amplitudes and Dslow and Dfast are the apparent diffusivities for the intra- and extra-axonal pools, respectively, is the effective transverse relaxation rate of the intra-axonal signal, is the additional effective transverse relaxation rate of the extra-axonal signal, and ω is the difference in angular frequency of the two signals. Intra-axonal signal fraction f was estimated from the intra- and extra-axonal signal amplitudes: f = A1/(A1 + A2).

Fitting parameters	A 1	A 2	R2*	ΔR2*	ω	D slow	D fast
	[a.u]	[a.u]	[s−1]	[s−1]	[Hz]	[10−3mm2/s]	[10−3mm2/s]
Lower bound	0	0	0	0	0	0	0
Upper bound	max(S) × 10	max(S) × 10	50	30	10	0.5	5
Estimated value	f = A1/(A1 + A2)

Results

Macroscopic effects

First, we investigated how the macroscopic signal evolution as a function of the acquisition delay time Δt changed with different levels of diffusion filtering. To address this question we estimated the macroscopic (apparent) effective transverse relaxation rate, , of the signal decay at each b-value, under the assumption that the signal evolution is mono-exponential. The data from the corpus callosum show a reduction of with increasing b-value when the diffusion filtering is applied either along the AP or SI directions (Fig. 4): from 17.4 ± 3.4 s−1 without diffusion filtering to 14.6 ± 2.9 s−1 at the maximum diffusion-filtering strength (μ ± σ from the voxels in the corpus callosum in all subjects). As the b-value increases, so does the relative fraction of the intra-axonal signal. Therefore, the reduction in the macroscopic relaxation rate with increasing diffusion filtering strength reveals that there are differences in signal evolution within the intra- and extra-axonal compartments.

Fig. 4

At each diffusion filtering (DF) level set by the b-values, the macroscopic effective transverse relaxation, , rate was estimated from the D-FASE signal by fitting a mono-exponential function of the acquisition delay Δt. The reduction of the macroscopic relaxation rate with increasing b is evident and hints at there being a slower decaying intra-axonal signal. The maps were created from the data with diffusion gradients applied along the anterior-posterior (AP) axis, which is approximately perpendicular to the main fibre direction in the corpus callosum and the corticospinal tract. The data from the mid-sagittal slice are shown.

Another important question is whether the change of the effective transverse relaxation rate with diffusion filtering is different between white matter fibre tracts. For all subjects, we compared the relaxation rates from ROIs within the corticospinal tract (parallel to ), and in the corpus callosum, and cingulum (perpendicular to ) (Fig. 5). The same figure shows the apparent intra-axonal signal fractions estimated at Δt = 0 ms assuming a bi-exponential diffusivity with negligible diffusivity in the intra-axonal compartment (Fig. 5F). The macroscopic transverse relaxation rate of the signal decay at b = 0 s/mm2 is highest in the cingulum (19.9 ± 0.8 s−1) followed by the value in the corpus callosum (17.1 ± 0.7 s−1) and then that in the corticospinal tract (12.0± 1.1 s in all subjects. Moreover, the macroscopic transverse relaxation rates in the corpus callosum and the cingulum decrease faster with increasing b in all subjects than in the corticospinal tract, resulting in -values of 17.0 ± 0.9s−1 in the cingulum, 14.6 ± 0.7s−1 in the corpus callosum, and 10.7 ± 0.5 s−1 in the corticospinal tract at the highest b-value of 4800 s/mm2.

Fig. 5

At each diffusion-filtering (DF) level set by the b-values, the macroscopic effective transverse relaxation rate was estimated from the D-FASE signal by fitting a mono-exponential function of the acquisition delay Δt within the ROIs in corpus callosum (red), cingulum (green) and corticospinal tract (blue). Results from monoexponential fitting to the data from a single subject are shown in A-D. The distributions of the mean macroscopic effective transverse relaxation rates among all subjects as a function of b for each ROI are shown in E and the corresponding apparent intra-axonal signal fractions (estimated with a bi-exponential in diffusivity at Δt = 0) are shown in F. The variation of the macroscopic -values with diffusion gradient strength and their differences between the three ROIs are evident.

Extracting microscopic effects

The advantage of diffusion filtering for the ASE signal becomes obvious when the two-compartment two-parameter model in Equation (3) is fitted to the 2D data from the corpus callosum (the residual myelin-water signal fraction was expected to be less than 5% of the total signal at echo-time of 57 ms under the assumption of myelin water fraction of 0.15 and MW T2 ≈ 20 ms and thus was assumed negligible). While the transverse relaxation rates of the intra- and extra-axonal signals overlap (the range, i.e. lower-upper adjacent, of the intra-axonal is 9– 20 s−1 with a mean value of 14.3 s−1, and of the extra-axonal is 4– 30 s−1 with a mean value of 17.5 s−1), the diffusive properties of these signals are clearly separable (the range of Dslow is 0 – 0.2 · 10−3mm2/s with a mean value of 0.07 · 10−3mm2/s, the range of Dfast is 0.1 – 2.2 · 10−3mm2/s with a mean value of 1.03 · 10−3mm2/s). Thus, with the added diffusion filtering, the ASE signals can be assigned to the intra- and extra-axonal compartments unambiguously (Fig. 6). The fitting results show combined data from all subjects, estimated using the case of Equation (3) leading to the smallest BIC, which was predominantly either ω = 0 (Case ii) or (Case iii).

Fig. 6

Combined results from all subjects when fitting the two-parameter two-pool model to the D-FASE data from the corpus callosum are shown. Three cases of Equation (3) were considered: i) ω and were allowed to vary, ii) ω = 0, iii) . For each voxel the fitting results were selected from the case which delivered the lowest BIC and the best fitting results were merged together for all subjects and voxels within the ROI in the corpus callosum. Dslow and Dfast values are plotted against and , respectively, the colormap represents the absolute value of the frequency offset between the intra- and the extra-axonal compartments. The smaller dots in the scatterplot correspond to the intra-axonal compartment, the size of the dots is proportional to the mean diffusivities in the intra-/extra-axonal compartments, respectively. Additionally, the distribution of diffusivities, relaxation rates, frequency offsets and intra-axonal signal fractions are displayed in the histograms. Blue and red correspond to the values from the intra- and extra-axonal compartments, respectively. Axis limits correspond to the fitting boundaries of the two-pool model. Case ii) can be distinguished from cases i) and iii) in the scatterplot and in the histogram of frequency offset, as it only represents ω = 0 Hz.

The differences between the ASE signals from the two compartments are reflected in differences in the transverse relaxation rates and in the frequency offsets ω/2π. The values within the corpus callosum range between 0 and 20 s−1 in the four subjects (when led to the smallest BIC). The magnitude of the frequency offset values range from 0 to 7 Hz in all subjects (when S(ω ≠ 0) led to the smallest BIC) with the median value being 3.6 Hz [2.8 Hz, 4.5 Hz] (the values in square brackets are the corresponding lower and upper quartiles; mean: 3.7 Hz, std: 1.8 Hz) which would correspond to 8.4 Hz [6.5 Hz, 10.5 Hz] at a magnetic field strength of 7 T and is in good agreement with previous studies, e.g. in (Sati et al., 2013; Tendler and Bowtell, 2019). Intra- and extra-axonal -values estimated in this work are notably in correspondence with previously reported T2-values estimated in diffusion-T2 studies: intra-axonal T2,a ranging 70 – 100 ms (≙ 10-14 s−1) (Veraart et al., 2018; McKinnon and Jensen, 2019) and extra-axonal T2,e ranging 40 – 50 ms (≙ 20-25 s−1) (Veraart et al., 2018). The distributions of values within the corpus callosum for the mid-sagittal slice are shown for each subject in Fig. 7. The colormaps in the bottom row of Fig. 7 indicate, which variation of Equation (3) gave the smallest BIC.

Fig. 7

A subset (mid-sagittal slice) of fitted values from Fig. 6 shown across the corpus callosum for each subject (columns), rows from top to bottom: apparent intra-axonal signal fraction f, intra-axonal relaxation rate , extra-axonal relaxation rate , absolute frequency offset value |ω|/2π [Hz], fast diffusivity Dfast [10−3mm2/s], and slow diffusivity Dslow [10−3mm2/s]. The last row shows which case of Equation (3) had the lowest BIC.

Discussion

In this work we have introduced an approach by which the NMR signal can be independently sensitised to both diffusion and relaxation effects in the intra- and extra-axonal compartments.

Diffusion filtering

We used strong diffusion gradients applied (nearly) perpendicular to the main fibre axis as a diffusion (Dslow-pass) filter, since this allowed us to attenuate the extra-axonal signal at a faster rate, Dfast than that of the intra-axonal signal, Dslow. The largest differential signal attenuation occurs when the diffusion-encoding gradient is applied perpendicular to the long axis of the axon. This can be difficult to achieve for a population of axons when there is orientational dispersion (Sotiropoulos et al., 2012; Ronen et al., 2014). However, we considered the signal as being diffusion-filtered, as long as the distributions of the (slow) intra- and (fast) extra-axonal apparent diffusivities within a ROI are well separated. In our experiments, diffusion gradients applied along the AP/SI, AP/RL, and RL/SI directions, served as good filters for the corpus callosum (fibre orientation predominantly RL), the corticospinal tract (fibre orientation predominantly SI), and the cingulum (fibre orientation predominantly AP), respectively.[1] The diffusion-filter strength has to be high enough to suppress the extra-axonal signal, but this must be achieved with a spin-echo time that is short enough to provide the SNR needed for the subsequent analysis. Such filter strengths can be obtained by deploying an MR system equipped with ultra-strong gradients (maximum gradient strength of 300 mT/m was deployed in this work): at the spin-echo times of 57 ms we could achieve a b-value of 4800 s/mm2 along the main radiological axes, which allowed us to suppress the extra-axonal signal by more than 95% under the assumption of exponential decay of the extra-axonal signal and a radial diffusivity of 0.75 · 10−3mm2/s.

Deployment of the diffusion filtering allowed us to separate the intra- and extra-axonal ASE signals unambiguously. Previous studies have evaluated the complex multi-echo gradient-echo signal evolution by fitting a multi-exponential function of the gradient-echo time. The signal fraction with a negative frequency offset (Sati et al., 2013; Thapaliya et al., 2018; Tendler and Bowtell, 2019; Alonso-Ortiz et al., 2018), or the long- signal fraction with the smaller amplitude (van Gelderen et al., 2012) was assumed to be intra-axonal. Diffusion filtering allowed us to assign the intra- and extra-axonal ASE signals from the magnitude data unambiguously based on the differences in the diffusive properties of the two compartments rather than on modelling of the subtle differences in the evolution of the gradient-echo signal at long echo times (i.e., when the myelin water signal has decayed away).

Origin of macroscopic anisotropy

We observed differences in the macroscopic -values between the three white matter fibre tracts in all subjects. In particular, is lower and shows less variation with b in the corticospinal tract than in the corpus callosum and the cingulum. Fibres in the corpus callosum and in the cingulum are predominantly perpendicular to , while those in the corticospinal tract are predominantly parallel to . The higher at b = 0 s/mm2 in the corpus callosum and the cingulum than in the corticospinal tract is consistent with the anisotropy that was identified in previous gradient echo experiments (Bender and Klose, 2010; Denk et al., 2011; Rudko et al., 2014; Gil et al., 2016). While some of the observed differences may arise from microscopic differences between the fibre populations (such as differences in intra-axonal signal fraction, axon diameter, orientation dispersion, or myelination), results from (Gil et al., 2016) showed that although the orientation-independent contribution to was higher in the corticospinal tract than in the cingulum, the size of the orientation-dependent contribution meant that was greater in the cingulum when oriented perpendicular to than in the CST when oriented parallel to .

In our work, the diffusion weighting perpendicular to the main fibre axis relatively attenuates the signal from the extra-axonal compartment, thus the macroscopic can also be seen as a function of the (b-value-dependent) intra-axonal signal fraction f. In a first-order approximation (when and ωΔt ≪ 1) can be represented as a weighted sum of the intra-/extra-axonal relaxation rates with the resonance frequency difference ω(θ), as a function of fibre orientation to , taken into account:

For instance high intra-axonal signal fraction f at b = 0s/mm2 would mean that is predominantly intra-axonal and would change little with b. In all subjects intra-axonal signal fractions of the corticospinal tract as estimated at Δt = 0ms (Fig. 5F) are similar or smaller than those in the corpus callosum and the cingulum and therefore do not explain the differences in between the tracts. Previous work has shown that nerve fibre orientation with respect to has a strong effect on the macroscopic (Bender and Klose, 2010; Denk et al., 2011; Rudko et al., 2014; Gil et al., 2016), and this is consistent with our findings. However, we cannot rule out the possibility that the microscopic differences between the fibre tracts or partial volume effects due to the presence of cerebrospinal fluid could contribute to a portion of the observed effects. Therefore, it would be beneficial to include measurements with different orientations of the head with respect to in future experiments in order to control for effects arising from differences between fibre populations. It would also be beneficial to add in low-b-value scans to provide measurements in which the signals from intra- and extra-axonal signals are relatively unaffected, but the CSF signal is suppressed.

Anisotropic magnetic susceptibility

The observed variation of macroscopic with diffusion filtering shows that there are differences in and/or frequency offset between the intra- and extra-axonal compartments. The difference in behaviour with fibre orientation, with the corticospinal tract showing a lower in the absence of diffusion weighting than the cingulum or corpus callosum, and a lesser change with the diffusion weighting, is consistent with the susceptibility effects predicted by the hollow cylinder model of the myelin sheath (Wharton and Bowtell, 2012). However, the differences in the macroscopic measured in the corpus callosum and cingulum indicate that factors other than fibre orientation also effect , and further work, ideally involving measurements made with the head held at different orientations to are needed to unambiguously demonstrate that the effects of the myelin sheath are the dominant cause of the inter-compartmental differences in gradient echo signal evolution.

When modelling the signal evolution, we found that the inclusion of a frequency offset between the intra- and extra-axonal compartments provided a better fit to the experimental data (judged by lower Bayesian Information Criterion (Schwarz, 1978)) in 50% of the corpus callosum voxels in all subjects. Such a frequency offset is explained by the hollow cylinder model in which the myelin sheath has anisotropic susceptibility with a radially oriented principal axis at each location (Wharton and Bowtell, 2012). The frequency offsets estimated in this study are consistent with previously reported values from work performed at 7 T (Wharton and Bowtell, 2012; Sati et al., 2013; Thapaliya et al., 2018; Tendler and Bowtell, 2019). Moreover, the intra-axonal frequency offset has been reported to be a function of the anisotropic susceptibility, fibre orientation to , and the ratio of the intra- and extra-axonal radii (Wharton and Bowtell, 2012). Through this technique we have access to a quantitative index of myelin properties that is inaccessible through dMRI alone.

Although our method provides reasonable estimates of the frequency offsets between the intra- and extra-axonal compartments, the frequency offset and the relaxation rate difference will have similar effects on the magnitude of the signal evolution. It has been already shown that frequency offsets are much better represented in the evolution of the signal’s phase rather than its magnitude, (Wharton and Bowtell, 2012), therefore to achieve a reliable quantification of myelin properties using the D-FASE sequence, phase-data-processing methods that can deal with the challenges arising from low SNR and the use of strong diffusion gradients need to be developed. In order to see the effect of frequency offset between the intra-/extra-axonal signals we need to access the local phase change due to suppression of the extra-axonal signal. However, phase generally also varies with b at large-length scale due to eddy current effects and head motion. Therefore, in order to assess the local phase change it is not sufficient to compare the unprocessed phases of the signals acquired at b = 0 s/mm2 and b = 4800 s/mm2. The existing approaches used to process gradient-echo phase data usually deal with a single dataset and are not able to take into account the large potential phase errors due to the strong gradients applied in this study. Calculation of the phase-change-difference between the corpus callosum and a reference region, such as CSF or grey matter could not be applied due to low SNR at high b-values in those regions. The alternative comparison between phase-changes in the corpus callosum and the corticospinal tract would also result in a partially non-local phase evolution due to spatially non-uniform phase distribution induced by eddy-current-related field perturbations. Bipolar diffusion gradients (Alexander et al., 1997) or complementary field measurements could potentially be used to estimate and correct for the diffusion-gradient-induced phase effects (Chan et al., 2014; Wilm et al., 2015).

Pulsation artifacts and partial volume effects from CSF also need to be taken into account (Enzmann and Pelc, 1991; Poncelet et al., 1992; Skare and Andersson, 2001; Jones and Pierpaoli, 2005; Bopp et al., 2018) and can be eliminated by introducing cardiac triggering and higher spatial resolution at a cost of longer acquisition and echo times.

Potential advantage for dMRI

We demonstrated the benefits of using diffusion filtering to separate the gradient echo signals arising from the intra- and extra-axonal pools. It would be interesting to consider whether the modulation of the EPI readout delay has an effect on estimated diffusion measures. It has been recently demonstrated that a two-compartment white matter model fitted to dMRI data acquired with a standard pulsed-gradient spin echo sequence does not provide a unique solution (degeneracy occurs) (Jelescu et al., 2016). It was also empirically demonstrated that adding extra dimensions to the measurement space (e.g. by varying the echo time of the measurement) the degeneracy can be resolved (Veraart et al., 2018). Based on the observations made in (Veraart et al., 2018) and due to the similar nature of R2- and -decays, our method has the potential to combat the degeneracy of dMRI parameter estimation. Moreover, our implementation allows us to fix the diffusion time whose variability is a potential confounding factor in recently reported dMRI-R2 studies (Veraart et al., 2018; McKinnon and Jensen, 2019). Still, the obtained diffusivity properties have to be interpreted with caution, since we do not explicitly model orientation dispersion of WM fibres and the intra-axonal signal fraction is echo-time-dependent due to differences in R2 between the intra- and extra-axonal compartments.

Challenges and future directions

We have touched upon some limitations of this study in the sections above. In this section we summarise the main challenges and discuss future directions for D-FASE work.

Challenge of separating from Δω using magnitude data. Although fitting results favoured the model with {, Δω = 0} over {, Δω = 0} in about 50% of the voxels in the corpus callosum ROIs for each subject, our data also show that using magnitude data does not allow us to separate the effects of a frequency offset and relaxation rate difference unambiguously. Therefore, it is important to develop a robust method for background phase removal in corresponding phase data so that these data can be included in the signal analysis. A background phase removal method will need to fulfill the following requirements: i) it should be robust against SNR reduction with increasing diffusion weighting; ii) it should account for temporally non-linear phase effects due to eddy currents from the strong diffusion gradients; iii) it needs to account for temporal variations of due to physiological effects, such as respiratory motion.

Coil combination method. Our coil combination method uses uncombined phase and magnitude data from one of three volume data sets acquired at b = 0 s/mm2, Δt = 0 ms to approximate the coil elements’ sensitivities and phases. When this method is applied to that volume itself, the coil combination method will be effectively sum-of-squares, which could potentially have different noise characteristics, than other volumes. We expect that the effect of the Rician noise characteristic is negligible at the high SNR levels of the image data acquired with b = 0 s/mm2. Nevertheless, an alternative method to estimate coil elements’ sensitivities and phases, such as acquiring gradient-echo data with the receiver array coil and repeating the same measurement with the body coil, could improve data combination from the coil array.

Susceptibility or microstructural features? We observed consistent difference in the macroscopic relaxation rates depending on the applied b-values and in the evolution of signals arising from different tracts in all subjects. Previous studies suggest that some of the observed effects could arise from the differences in fibre orientation θ to in the selected ROIs (corticospinal tract vs corpus callosum and cingulum). However, we cannot rule out that the differences between the fibre tracts partially arise from other microscopic and fibre specific measures, such as axon diameter, fibre density and orientational dispersion effects. For example, we observed differences in macroscopic between the corpus callosum and the cingulum (which are both mainly perpendicular to ). These differences could also reflect differences in intra-axonal R2-values between the fibre tracts, which potentially arise from the axon-diameter-dependent surface relaxation (Kaden and Alexander, 2013) or fibre dispersion. On the other hand, these differences might also reflect g-ratio-dependent frequency offsets between the intra-/and extra-axonal signals, as proposed in (Wharton and Bowtell, 2012): where γ is the gyromagnetic ratio, G is the ratio between the intra- and extra-axonal radii, χ is the anisotropic susceptibility and θ is fibre orientation to .

In order to thoroughly investigate and separate the different sources of contrast between different fibre tracts, further experiments are required. In order to separate out orientation-dependent relaxation effects one will need to include experiments involving re-orientation of the participant’s head in the scanner to modulate fibre orientation to for the same fibre tract. Additional modulation of the spin-echo time at fixed diffusion weighting may provide supplementary information on transverse relaxation rate and its link to axonal diameter. However, some investigation is required to show that this information is not redundant.

Conclusion

In conclusion, by exploiting the ultra-strong gradients of a Connectom scanner, we have been able to integrate diffusion-encoding gradients into an asymmetric spin-echo (ASE) sequence to produce selective suppression of the extra-axonal signal, while maintaining sufficiently short echo time to measure the ASE signal amplitude reliably from the intra-axonal space. This approach, for the first time, enables direct measurement of compartmental-specific gradient echo signals and provides a powerful complement to approaches in diffusion-only, and diffusion-relaxometry studies, where parallel efforts have been made to obtain compartmental (e.g. intra-vs extra-axonal) specific microstructural measurements.