The dot-compartment revealed? Diffusion MRI with ultra-strong gradients and spherical tensor encoding in the living human brain.

The so-called "dot-compartment" is conjectured in diffusion MRI to represent small spherical spaces, such as cell bodies, in which the diffusion is restricted in all directions. Previous investigations inferred its existence from data acquired with directional diffusion encoding which does not permit a straightforward separation of signals from 'sticks' (axons) and signals from 'dots'. Here we combine isotropic diffusion encoding with ultra-strong diffusion gradients (240 ​mT/m) to achieve high diffusion-weightings with high signal to noise ratio, while suppressing signal arising from anisotropic water compartments with significant mobility along at least one axis (e.g., axons). A dot-compartment, defined to have apparent diffusion coefficient equal to zero and no exchange, would result in a non-decaying signal at very high b-values (b≳7000s/mm2). With this unique experimental setup, a residual yet slowly decaying signal above the noise floor for b-values as high as 15000s/mm2 was seen clearly in the cerebellar grey matter (GM), and in several white matter (WM) regions to some extent. Upper limits of the dot-signal-fraction were estimated to be 1.8% in cerebellar GM and 0.5% in WM. By relaxing the assumption of zero diffusivity, the signal at high b-values in cerebellar GM could be represented more accurately by an isotropic water pool with a low apparent diffusivity of 0.12 μm2/ms and a substantial signal fraction of 9.7%. The T2 of this component was estimated to be around 61ms. This remaining signal at high b-values has potential to serve as a novel and simple marker for isotropically-restricted water compartments in cerebellar GM.

Introduction

Diffusion Magnetic Resonance Imaging (dMRI) (Le Bihan and Breton, 1985) probes structures at much smaller length-scales than the imaging resolution by sensitising the signal to the random molecular motion of water. Biophysical modelling of the contributions to this signal aims to characterise tissue microstructure properties by carefully selecting model compartments (typically multiple non-exchanging water pools) that have a measurable impact on the signal (Stanisz et al., 1997). In healthy white matter (WM), biophysical models typically include anisotropic extra- and intra-axonal compartments (Alexander et al., 2010; Assaf and Basser, 2005; Fieremans et al., 2011; Jespersen et al., 2007; Kroenke et al., 2004; Lampinen et al., 2019; Novikov et al., 2018b; Reisert et al., 2017; Sotiropoulos et al., 2012; Stanisz et al., 1997; Zhang et al., 2012). The inclusion of a so-called “dot-compartment” for WM-modelling is motivated by the observation of an almost constant, non-attenuating signal at very high b-values (e.g., ). This has been hypothesised to arise from the ubiquity of small isotropic spaces (e.g., glial cell-bodies) wherein the diffusion of water molecules is highly restricted in all directions (Alexander et al., 2010; Stanisz et al., 1997), leading to a near-zero apparent diffusivity. A method to measure the signal fraction of such isotropically-restricted components accurately in vivo could thus potentially provide a proxy for the density of cells and enable quantification of cellular pathology in a wide range of neurological and psychiatric disorders (Palombo et al., 2019b).

Previous work investigating compartmental contributions to the dMRI signal from conventional pulsed-gradient encoding – also called Stejskal-Tanner encoding (Stejskal and Tanner, 1965) or linear tensor encoding (LTE (Westin et al., 2016)) – showed that including a dot-compartment provided a more complete description of the WM dMRI signal, both ex vivo (Panagiotaki et al., 2012) and in vivo (Ferizi et al., 2014; Zeng et al., 2018). However, a dot-compartment is not generally included in WM biophysical models, e.g. (Assaf and Basser, 2005; Behrens et al., 2003; Jespersen et al., 2007; Kroenke et al., 2004; Novikov et al., 2018b; Zhang et al., 2012). Moreover, a recent study of the dMRI signal in WM at b-values up to on a clinical MRI system suggested that the WM dot-signal-fraction is negligible (Veraart et al., 2019).

Probing the dot-compartment in anisotropic tissue is challenging with conventional LTE, due to the strong relationship between encoding-direction and orientation-distribution of anisotropic tissue microenvironments. Even when measuring along the dominant axis of a fibre bundle in which there is orientation dispersion, a slow diffusing component can be observed (due to the gradient direction not being perfectly parallel to all of the fibres); it is therefore challenging to disentangle this from the scenario in which a dot-compartment is present (Fig. 1a). Here, we address this problem by the use of spherical tensor encoding (STE, also called isotropic diffusion encoding) to render signals insensitive to orientation and anisotropy (Eriksson et al., 2013; Lasič et al., 2014; Mori and Van Zijl, 1995; Westin et al., 2016; Wong et al., 1995). STE at high b-values can suppress the dMRI signal from water pools that are mobile along at least one axis (Fig. 1a). At sufficiently high b-values only the signal from compartments with very low or zero diffusivity in all directions would remain.

Fig. 1

a) Simulations of LTE and STE data for two different scenarios (schematically represented in the middle): dispersed sticks representing axons surrounded by extra-axonal space (top, blue surround), vs dots ​+ ​dispersed sticks surrounded by extra-axonal space (bottom, green surround). Here, we used a Watson distribution to simulate a stick orientation dispersion (OD) of 0.7 and (blue), and OD ​= ​0.5 and (green). In both simulations, the parallel diffusivity was set to 2.1 for the stick compartment, and the extra-cellular compartment consisted of non-exchanging Watson-distributed ‘zeppelins’ with the same OD as the sticks and parallel and perpendicular diffusivities of 1.9 and 0.8 , respectively. The two scenarios result in very similar signals for LTE across a wide b-value range, and can be disentangled better at high b-values with STE. A linear y-scale is chosen here to not make small differences seem disproportionally large. b) STE simulated as in (a) but with varying , at different SNR levels. The dashed-dotted line represents the rectified noise floor, and the error bars represent the mean and standard deviation over 5000 noise realisations. A logarithmic y-scale is chosen here to improve the visualisation for different . b is given in .

Previous work using STE obtained by a series of pulsed gradients on a clinical system concluded that the dot-signal-fraction is likely lower than 2% in WM, and therefore has a negligible contribution to the dMRI signal (Dhital et al., 2018). However, the gradient amplitude available on clinical MRI scanners (40–80 ​mT/m) limits the maximal b-value per unit signal-to-noise ratio (SNR) – needed for reliable quantification of the dot-signal-fraction – whereas ultra-strong gradients (e.g., 300 ​mT/m) allow much higher b-values per unit SNR (Jones et al., 2018; Setsompop et al., 2013). Furthermore, the previous implementation of STE used waveforms with exceedingly low efficiency (Sjölund et al., 2015). In this work, we leverage the power of ultra-strong gradients and optimised asymmetric STE gradient waveforms to reduce the echo time (TE) significantly, thereby increasing SNR. This allows signal decays to be examined in the living human brain over a much larger range of b-values and TEs typically unachievable using clinical MRI scanners (Tax et al., 2017), and thus provides a more reliable assessment of signal fractions which could result from isotropically-restricted compartments. In addition, we extend the analysis to tissue types beyond cerebral WM, including deep grey matter (GM) and the cerebellum.

Theory

Assuming Gaussian diffusion within a compartment, the signal arising from the ith compartment, represented by diffusion tensor and which contributes a relative signal fraction to the signal, probed by symmetric b-tensor can be described by

The total signal is then the sum of the signals from the individual compartments, with summing to one. The b-tensor is a positive semi-definite tensor which we here design to be axially symmetric; it can then be characterised by its trace – better known as the b-value - and its anisotropy (Eriksson et al., 2013; Topgaard, 2017; Westin et al., 2016), where and are the eigenvalues corresponding to the eigenvectors along and perpendicular to the symmetry axis, respectively. represents the signal at , and denotes the trace of the matrix product between the tensors.

In the case of STE, the b-tensor is isotropic and thus and . For non-exchanging Gaussian compartments, the STE-signal simplifies towhere is the mean apparent diffusivity of each compartment.

An isotropically restricted compartment typically exhibits a very low mean apparent diffusivity. If we index this compartment as and assume , then the only remaining signal when approaching high b-values (beyond a certain b-value, ) is that arising from the isotropic restricted compartment:

For example, for a two-compartment system with and , the signal from the second compartment is reduced to 0.1% for , while the signal from the first compartment is only reduced to 42%. This means that the behaviour of at increasing b-values is increasingly dominated by compartments with lower apparent diffusivity.

In the case of a dot-compartment with zero mean apparent diffusivity, i.e. , Eq. [3] simplifies tosuch that the dot-signal-fraction is equal to the relative signal that remains at high b-values. Fig. 1b shows the simulated signal in the case of non-exchanging compartments of which one is a dot-compartment. Even if the signal does not yet exhibit a plateau, the relative signal at the highest b-value can serve as an upper limit of because . The accuracy of this limit is affected by the presence of the rectified noise floor, with standard deviation of the Gaussian noise added to each of the real and imaginary channels (Jones and Basser, 2004) (Fig. 1b).

Methods

Data

Five healthy adult volunteers were included in the study (3 female), which was approved by the Cardiff University School of Medicine ethics committee. Written informed consent was obtained from all participants.

Participants were scanned on a 3T Connectom MRI system (Siemens Healthcare, Erlangen, Germany) with an ultra-strong 300 ​mT/m gradient set. The acquisition protocol included a structural MPRAGE (Magnetization Prepared RApid Gradient Echo) (de Lange et al., 1991) with voxel size 1 ​× ​1 ​× ​1 and dMRI sequences. The dMRI data were acquired using a prototype spin-echo sequence with an echo-planar imaging (EPI) readout, that enables user-defined gradient waveforms to be used for diffusion encoding (Szczepankiewicz et al., 2019a). For STE we used b ​= ​[250, 1500, 3000, 4500, 6000, 7500, 9000, 10 ​500, 12 ​000, 13 ​500, 15 ​000] , repeated [6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36] times, respectively. The b-values and repetitions were interleaved over volumes to reduce the impact of system drift (Hutter et al., 2018; Vos et al., 2016). For LTE, the b-tensor principal eigenvectors were distributed over the unit sphere for each b-shell. b ​= ​0 ​ (b0) images were acquired every 15th image for monitoring and correction of subject motion. Additional b0 images with reversed phase-encoding were acquired to correct for susceptibility distortions (Chang and Fitzpatrick, 1992). No in-plane acceleration was used, and imaging parameters were: voxel size ​= ​4 ​× ​4 ​× ​4 mm3, matrix ​= ​64 ​× ​64, 34 slices, TR ​= ​4300 ​, partial-Fourier ​= ​6/8, bandwidth ​= ​1594 ​Hz/pixel.

The waveforms used for STE and LTE are shown in Fig. 2, and were optimised numerically (Sjölund et al., 2015) to be Maxwell-compensated (Szczepankiewicz et al., 2019b) and enable a TE as short as . These waveforms render superior encoding efficiency due to their optimised asymmetric trajectory in q-space compared to standard 1-scan-trace imaging (which requires TE ​= ​270 ​ for ).

Fig. 2

Linear tensor encoding (LTE) and spherical tensor encoding (STE) waveforms for , and the corresponding power spectra of the dephasing vector . Timings for the first waveform, temporal gap (180° pulse), and second waveform were [28.6, 6.9, 28.6] for LTE and [35.5, 6.9, 25.6] for STE. The maximum gradient amplitudes along a single axis were 131 and 240 ​mT/m for LTE and STE, respectively.

Two of the volunteers were additionally scanned with the STE sequence using different TE ​= ​[88 ​115 ​140 165] , to obtain estimates of T2. Note that the timings of the waveforms did not change.

Preprocessing

The dMRI data were corrected for Rician noise bias (Koay et al., 2009a; St-Jean et al., 2016) using an estimate of the Gaussian noise standard deviation from PIESNO (Koay et al., 2009b) and an estimate of the true underlying Rician signal from denoising (Veraart et al., 2016), to determine whether or not any plateau arising in the signal decay curve could be attributed to the effects of the noise floor. We proceeded with the debiased (but not denoised) data in further processing. The data were checked for signal intensity errors including slice-wise outliers (Sairanen et al., 2018). The STE data were corrected for subject motion by registering the interleaved b0 images to the first b0 image and applying the corresponding transformations to the diffusion-weighted images (DWIs). The LTE data were corrected for subject motion and eddy-current geometrical distortions using FSL EDDY (Andersson and Sotiropoulos, 2016). Susceptibility geometrical distortions were corrected using TOPUP (Andersson et al., 2003) and for geometrical distortions due to gradient nonlinearities using code kindly provided by colleagues at the Athinoula A. Martinos Center for Biomedical Imaging at Massachusetts General Hospital (Glasser et al., 2013; Jones et al., 2018; Rudrapatna et al., 2018; Setsompop et al., 2013).

The MPRAGE image was segmented into regions using Freesurfer (Fischl et al., 2002) and affinely-registered to the corrected b0 image using FSL FLIRT (Jenkinson et al., 2012). The resulting WM, GM, deep GM (dGM), cerebellar WM (cWM) and cerebellar GM (cGM) segmentations were then used to guide the delineation of regions-of-interest (ROIs) for further analysis. Only voxels in which the tissue probability derived from the Freesurfer segmentations was larger than 90% were considered, and the ROIs were drawn manually to avoid including signal artefacts. For WM, two separate regions were considered: ROIs were drawn on coronal slices in medial WM lateral to the midbody of the corpus callosum (denoted by mWM), and in the occipital regions (denoted by oWM), see Fig. 4.

Fig. 4

STE signal decay for 5 healthy subjects (mean and standard deviation in each ROI), in ROIS in the mWM (red), oWM (orange), cWM (green), dGM (cyan), and cGM (blue); examples of the ROIs are shown for Subject 1. The left column shows the signal before Rician-bias correction plotted with a logarithmic y-scale, to better visualise deviations from mono-exponential behaviour, with a close-up at high b-values. The right column shows the signal after Rician-bias correction plotted with a linear y-scale to be able to visualise negative values. b is given in . The dotted line represents the mean of the estimated noise floor in each ROI (Koay et al., 2009b); the lines of each ROI visually overlap.

Quantitative characterisation of the STE signal at high b-values

Eq. [3] was fitted to the data with  ​= ​10 ​000 ​, (thus including  ​= ​[10 ​500, 12 ​000, 13 ​500, 15 ​000] ) using a nonlinear least-squares trust-region-reflective algorithm implemented in MATLAB (The MathWorks, Natick, USA). The fit was randomly initialised 10 times within bounds [0 max()] and [0 0.3] for and respectively (the fit was constrained within bounds [0 ] and [0 1]), and the solution with the lowest residual norm was selected. In addition, an estimate of was obtained from the two lowest b-value DWIs, to reduce contribution from CSF (Baron and Beaulieu, 2015). From this, the ‘tissue signal fraction’ was estimated (Eq. [3]). Finally, from and an upper limit was derived (Eq. [4]).

An estimate of T2 was obtained by fitting to the high b-value data with a similar fitting procedure as described before, and initialising and constraining between [0 300] . For all fits, gradient-nonlinearities were taken into account by considering the voxel-wise effective B-tensor computed from the spatially-varying coil tensor (Bammer et al., 2003).

Results

ROI delineation

Fig. 3 shows the Freesurfer segmentation results overlaid on individual diffusion-weighted images of one participant. Fig. 4 shows results of the manually delineated ROIs visualised for one of the healthy subjects. The mWM, oWM, cWM, cGM, and dGM ROIs include on average 59, 58, 39, 199, and 52 voxels across participants, respectively. The GM segmentations only include a few voxels that are classified as > 90% GM, which are sparsely distributed. We therefore only consider data in the mWM, oWM, cWM, cGM, and dGM ROIs.

Fig. 3

a) STE signal as a function of b-value (in ), with the Freesurfer tissue segmentations indicated in red ​= ​WM, yellow ​= ​GM, cyan ​= ​deep GM, green ​= ​cerebellar WM and blue ​= ​cerebellar GM. b) Relative STE signal change.

STE signal decay across all b-values

Fig. 3a shows the signal of the image intensity in individual DWIs as a function of b-value in a healthy brain. The signal intensity in most of the cerebral WM has decayed substantially at 2. However, the cerebellar GM retained a remarkably high signal at these high b-values, remaining well above the noise floor even at . Fig. 3b shows the signal in a more quantitative fashion; regions with lower intensity have a higher relative signal change compared to the signal. The cerebellar GM persistently has a high intensity compared to other regions and thus the lowest relative signal change.

Fig. 4 shows the signal decay for each ROI in the five healthy subjects. Both the original and Rician-bias-corrected signal decay curves are shown, accompanied by an estimate of the noise floor. At , the signals from mWM and dGM clearly approach the noise floor. In contrast, oWM, cWM, and cGM exhibit a mean-signal that is above the noise floor for all five subjects. After Rician-noise-bias-correction, the signal is still above zero albeit it can be seen that it continues to decay.

STE signal characterisation at high b-values

Table 1 gives quantitative features related to the STE signal decay at high b-values. For each parameter the median of the voxel-wise fits in each ROI is given, and 10th – 90th percentiles are reported within brackets. The third column presents estimates of the relative rectified noise floor, derived from estimates of the noise standard deviation and b0 signal (i.e. and ). An estimate of the relative noise floor of 0.5% indicates that an SNR of about 250 on the b0 signal could be achieved (SNR estimates are presented in the fourth column of Table 1). The fifth column gives an estimate of the temporal SNR (tSNR), defined by the temporal mean of the images divided by their temporal standard deviation (see also Supplementary Fig. 3). The sixth column presents estimates for the dot-signal-fraction (Eq. [4]) and the last two columns present estimates in the case of an isotropically-restricted compartment with non-zero diffusivity (Eq. [3], i.e. , and ). All estimates are obtained after Rician bias-correction to reduce bias from the least-squares fitting. We will describe characteristics of these features for the different ROIs in the following paragraphs.

Table 1

Parameter estimates (median of the voxel-wise fits in each ROI and 10–90 percentile) for the standard deviation , , and , in different ROIs. has units of . No estimates of and are given for dGM and mWM as these signals are hitting the noise floor.

Region	Subject	σ˜π/2/S˜0⋅100%	SNR	tSNR	f˜_dot⋅100%	f˜_1⋅100%	D˜_1
mWM	1	0.6	220	115	0.0 (−0.2–0.2)
	2	0.4	285	155	0.0 (−0.1–0.2)
	3	0.5	272	140	0.1 (−0.1–0.2)
	4	0.5	246	100	−0.1 (−0.1–0.1)
	5	0.5	232	168	0.0 (−0.2–0.2)
		0.5	251	136	0.2 (-0.1–0.2)
dGM	1	0.9	139	63	−0.2 (−0.5–0.1)
	2	0.8	157	61	−0.2 (−0.5–0.0)
	3	0.8	158	54	−0.2 (−0.5–0.0)
	4	0.0	145	46	−0.3 (−0.5–0.1)
	5	1.1	114	64	−0.3 (−0.6–0.0)
		0.9	143	58	−0.2 (-0.5–0.1)
oWM	1	0.3	404	72	0.2 (0.1–0.4)	2.7 (0.9–7.2)	0.16 (0.09–0.27)
	2	0.3	404	81	0.2 (0.1–0.4)	2.7 (1.1–22.1)	0.18 (0.10–0.37)
	3	0.3	432	89	0.2 (0.1–0.3)	3.2 (1.1–8.4)	0.16 (0.11–0.27)
	4	0.3	417	93	0.2 (0.1–0.3)	2.6 (1.0–9.0)	0.17 (0.11–0.31)
	5	0.3	417	96	0.2 (0.1–0.4)	2.2 (1.1–9.6)	0.14 (0.09–0.31)
		0.3	415	86	0.2 (0.1–0.4)	2.6 (1.1–9.8)	0.16 (0.09–0.31)
cWM	1	0.4	306	96	0.4 (0.1–0.8)	5.4 (0.9–13.5)	0.18 (0.02–0.35)
	2	0.5	236	79	0.1 (−0.2–0.4)	5.6 (0.9–62.6)	0.24 (0.00–0.87)
	3	0.6	216	50	0.1 (−0.3–0.4)	5.6 (0.3–28.4)	0.24 (0.09–0.77)
	4	0.5	236	45	0.2 (−0.1–0.5)	4.8 (0.1–83.7)	0.24 (0.09–1.00)
	5	0.6	216	65	0.1 (−0.1–0.4)	7.4 (0.2–50.7)	0.27 (0.10–0.74)
		0.5	242	67	0.2 (-0.1–0.6)	5.4 (0.2–46.9)	0.23 (0.07–0.68)
cGM	1	0.2	570	78	2.0 (1.5–2.3)	10.9 (8.4–15.7)	0.14 (0.12–0.17)
	2	0.3	448	73	1.6 (1.2–1.8)	9.6 (7.9–12.3)	0.13 (0.11–0.16)
	3	0.3	418	51	1.5 (1.2–1.8)	9.2 (7.2–11.7)	0.12 (0.10–0.15)
	4	0.3	432	41	1.8 (1.3–2.1)	9.8 (8.0–12.4)	0.12 (0.10–0.14)
	5	0.3	392	62	1.7 (1.4–2.0)	9.1 (7.2–11.6)	0.12 (0.10–0.14)
		0.3	452	61	1.8 (1.4–2.2)	9.7 (7.7–12.8)	0.12 (0.10–0.16)

For the mWM and dGM ROI, the mean signal at high b-values converges to the noise floor (Fig. 4). We estimate an upper limit of of 0.5% and 0.9% respectively.

For the oWM and cWM ROI, the estimated upper limits of are 0.3% and 0.5% respectively. The signal at high b-values is still decaying, and can thus be better explained by the presence of a compartment with non-zero apparent diffusivity with estimated signal fractions of 2.6% and 5.4%, and estimated apparent mean diffusivities of 0.16 and 0.23 for oWM and cWM, respectively. Fig. 5a shows scatter plots of these estimates, showing that the spread is large (see also the percentiles in Table 1).

Fig. 5

a) Parameter estimates and , in the oWM (orange), cWM (green), and cGM (blue) ROIs, where each point represents a voxel. b-c) Map of the fits of Eqs. [3] (not assuming zero apparent diffusivity) and [4] (assuming zero apparent diffusivity) in an axial slice of the cGM, respectively; the cerebellar WM is masked out. has units .

For the cGM ROI, we find an upper limit of of 1.8%, with a residual signal that is well above the noise floor. Per Eq. [3], we estimate an average apparent mean diffusivity of 0.12 and an average signal fraction of 9.7%. In some areas, the signal fraction is estimated as high as 15.7%. These estimates are consistent across healthy subjects (Fig. 5a). When visualising the estimates in the cerebellar GM one can observe that they are spatially heterogeneous (Fig. 5b). As a comparison, we show the spatial variability of in Fig. 5c.

Comparison of LTE and STE signals

In Fig. 6 one can readily appreciate the difference between b0-normalised STE and directionally-averaged LTE signals in the different tissue types. These diffusion weightings also give complementary information in GM, where the STE encoding at high b-values has suppressed signal arising from compartments that are mobile along at least one axis (e.g., ‘sticks’ that could represent axons). The overlap of the signal decay curves is high between the healthy controls.

Fig. 6

a) Signal upon LTE and STE () with the same intensity scale. b) LTE (dashed lines) and STE (solid lines) signals, with b in . Colours correspond to Fig. 4.

Characterisation of T2 at high b-values

Fig. 7a shows estimates of T2 and D for the two healthy controls that were scanned with STE at different TE. The T2 estimates are consistent between the two subjects and range between 57 and 69 ​ (10–90 percentile) with a median of 61 ​. Fig. 7b shows the estimated as a function of TE for both subjects (median and 10–90 percentile). A decreasing trend is observable, and the difference between the highest and lowest TE is found to be significant in both subjects using a paired-sample t-test (p ​< ​1.15e-6).

Fig. 7

a) Parameter estimates and , in the cGM ROI. has units and has units . b) Parameter estimates in the cGM ROI (median and 10–90 percentile) as a function of TE. The dotted errorbars represent the median and 10–90 percentile, and the solid (very narrow) errorbars the mean and standard error of the mean. has units .

Discussion

In this study we report a clear depiction in vivo from an isotropically-restricted compartment in dMRI. This compartment is present particularly in the cerebellar GM, but support for its existence can also be found in the WM. Our observations were enabled by ultra-strong gradient hardware (Jones et al., 2018; Setsompop et al., 2013) and recent developments for tensor-valued diffusion encoding (Sjölund et al., 2015; Szczepankiewicz et al., 2019b). STE provides essential complementary information to LTE, but the waveforms generally take up more time than Stejskal-Tanner LTE encoding, leading to long TEs and thereby inferior SNR. With the help of ultra-strong diffusion gradients (240 ​mT/m along a single axis, Fig. 2), a TE as short as 88 ​ could be achieved even for a b-value of 15 ​000 ​. As a result, the SNR was such that we could clearly observe signal amplitudes well above the noise floor.

A plateau of the diffusion-weighted signal (i.e., region of no further signal decay with increasing b-value), even at high b-values, was not observed in any region of interest. This makes the significant contribution of water residing in a dot-compartment with zero apparent diffusivity and no exchange unlikely, in agreement with previous work (Dhital et al., 2018; Veraart et al., 2019). Such a compartment would reflect a spherical compartment with radius that is negligible compared with the diffusion length where is the ‘bulk’ diffusion coefficient and the effective diffusion time. This results in an apparently zero diffusion coefficient, a similar assumption underlying the perpendicular diffusivity of a ‘stick’ compartment (Novikov et al., 2018a; Veraart et al., 2019). Nevertheless, a slowly decaying STE signal was observed in the cerebellar GM and some WM regions.

Signal characterisation and implications

The observation of a slowly decaying STE signal can be supported by two hypotheses: (i) a zero-apparent-diffusivity compartment exists but is not observed as such because it is in exchange with its surroundings; or (ii) the compartment exhibits a low but non-zero apparent diffusivity. The effect of the first hypothesis is illustrated in Supplementary Fig. 1, which shows the noiseless signal decay for different exchange times using a two-compartment Kärger model (Kärger, 1971; Nilsson et al., 2010). At infinite exchange times, the estimated dot-signal fraction approaches its true value. However, at exchange times of e.g. 500 ​ the signal does not exhibit a plateau and the estimated upper limit of is negatively biased. This figure should be interpreted as an illustration of how exchange could affect the signal in a simple scenario with a well-defined effective diffusion time; however free-waveforms give rise to a diffusion-time spectrum and representing the interaction with exchange will be more complex. This is subject to future work.

Regarding the second hypothesis, a slow-diffusing component has not been observed previously in STE data. Previous work has characterised mean apparent diffusivities derived from STE data up to by using a regularised inverse Laplace transform (Avram et al., 2019) or by fitting a finite series of exponentials that could represent different compartments and comparing the fits of the models through the Akaike Information Criterion (AIC) (Dhital et al., 2018). These works showed little deviation from mono-exponential behaviour in WM and single-peak diffusivity distributions in brain parenchyma in the range of b-values used. However, in the logarithmic plots in Fig. 4 one can clearly observe that the signal decay starts deviating from mono-exponential behaviour for in most tissue types, which could explain why this component has not been reported previously. Rather than quantifying the signal across the entire range of b-values and comparing the fit of models with different numbers of compartments, we focus here on quantifying the STE signal at high b-values. Using data from the whole range of b-values, the “cut-off” b-value was visually determined as the lowest b-value beyond which the signal decay approaches again a straight line in a logarithmic plot (Fig. 4). The strategy of working in a regime where the signal of some of the compartments is compressed has been adopted in other studies to focus on the intra-axonal signal with LTE (Kleban et al., 2019; McKinnon and Jensen, 2019; Veraart et al., 2019). The advantage is that the regime of interest can be studied into greater detail, that the number of parameters is reduced, and the fit is not biased by potentially wrongful assumptions about the other compartments. Similarly, we have used STE at high b-values here to suppress signal from compartments with significant mobility in at least one direction.

Using a simple representation based on the often-adopted assumptions of Gaussian diffusion and no exchange, the results provide support for the presence of an isotropic water pool with low diffusivity in the oWM, cWM, and cGM ROIs. In WM, Dhital et al. (2018) found that for a hypothetically small, yet finite, diffusivity of , the relaxed upper limit of the signal fraction was 2.7%. In the present study, we found a similarly low signal fraction, but the diffusivity was estimated to be twice as high () in the oWM and cWM ROIs, albeit with a high variability across voxels (Table 1). In the medial WM the signal converged to the noise floor; this could be caused by the larger distance to the RF receiving coils (and thus lower SNR), or a genuinely lower density of slow-diffusing components compared with the occipital WM, or both.

In the cGM ROI, the signal fraction of the slowly diffusing isotropic water pool was estimated to be as large as 15%, and this component thus makes a significant contribution to the signal. Linking this finding to tissue microstructure derived from histology or realistic numerical simulations of brain cells (Palombo et al., 2019a) is the subject of future work. It has been suggested previously that in cortical GM, the abundance of cell bodies has a significant impact on the LTE signal at high b-values (Palombo et al., 2019b). In that work, the LTE signal at high b-values was considered to be arising from non-exchanging sticks representing neurites, and spheres with a finite radius representing cell bodies. Following this picture, STE at high b-values would nullify the stick-signal and only the signal specific to the cell bodies would remain. Fig. 8a shows simulated signal decays for the STE protocol used in this study, for spherical compartments with different diffusivities and radii . The signal resulting from a sphere acquired with arbitrary waveforms was derived by (Codd and Callaghan, 1999) using the matrix formalism (Callaghan, 1997), and we have used the implementation in the MISST toolbox (Drobnjak et al., 2011, 2010; Ianuş et al., 2013) to generate these results. From the signal patterns it becomes apparent that disentangling and will be challenging. Adopting the same strategy as in (Palombo et al., 2019b), we have explored the fitting with fixed to , cf. Fig. 8b. Specifically, a spherical compartment was fitted to the data with  ​= ​10 ​000 ​, (thus including  ​= ​[10 ​500, 12 ​000, 13 ​500, 15 ​000] ) using the same nonlinear least-squares trust-region-reflective algorithm as was used for other experiments. The fit was constrained within bounds [0 ] and [0 20] for and , respectively. Estimated values for were found to be around 4.5 , which seems to be lower than can be visually derived from (Palombo et al., 2019b). The estimates were consistent across the cerebellar GM and across subjects. Note that the waveforms used here were not optimised for size estimates and for avoiding regimes of exchange. For example, the b-value regime and waveforms used here were such that large spherical compartments with high intra-cellular diffusivity become suppressed. Interestingly, Fig. 8a suggests that the remaining signal could both reflect the existence of an additional compartment with smaller effective radius than reported in Palombo et al., or with a larger radius but more hindered intra-cellular diffusion than previously assumed.

Fig. 8

a) Left: STE signal simulation of a spherical compartment for different values of and . is given in . Right: estimates of in (Eq. [3]) for different values of and . b) Estimates of and in the cerebellar GM of five healthy controls, assuming .

The use of pulsed-gradients allows a more precise definition of the time-scale of diffusion. The STE waveforms in Fig. 2 have broader frequency spectra, affecting the way time-dependent diffusion is encoded (Jespersen et al., 2019; Lundell et al., 2019). Under the assumption of Gaussian (and thus time-independent) diffusion in each compartment (as in Section 2), the net signal becomes non-monoexponential but remains time-independent; as such the signal decay arising from two sets of waveforms with the same B-tensor, but different frequency spectra, would look identical. However, the assumption of compartmental Gaussian diffusion is theoretically only valid for sufficiently short or long diffusion times or low diffusion weightings; beyond these regimes time-dependent diffusion will be encoded differently by waveforms with different frequency spectra. The use of different LTE waveforms with different frequency characteristics and b-values up to 5000 ​ has previously revealed a strong contrast in the cerebellum (Lundell et al., 2017, 2015). In the specific case of STE as studied here, several works (de Swiet and Mitra, 1996; Jespersen et al., 2019; Lundell et al., 2019) have shown that non-Gaussian diffusion within each compartment can lead to anisotropic time-dependence, i.e., probing different time-dependence in different directions. This means that for anisotropic pores, such as cylinders and ellipsoids, the signal decay in STE at lower b-values still depends on the orientation and dispersion of the pores (or the rotation of the waveforms). In the present study, we focused on the high b-value regime to completely suppress the signal from anisotropic compartments that have significant mobility along at least one axis. Therefore, the remaining signal is expected to come from restricted isotropic compartments only, and is as such expected to be rotationally invariant. Non-Gaussian diffusion within these isotropic compartments becomes a contributing factor if one for example tries to estimate the variance of the highly restricted isotropic diffusivities (Jespersen et al., 2019), which is beyond the aim of this study.

In summary, to disentangle the microstructural origins of the slowly decaying STE signal at high b-values (i.e. it being exchange and/or the length scale of the restrictions), additional and independent measurements with different time-characteristics are necessary. Supplementary Fig. 1 (exchange) and Fig. 8 (restrictions) suggest that for the experimental data as presented in this study multiple degenerate solutions could exist, and at this stage the data can be adequately represented by the simple representation used.

The T2 of the low-diffusivity component in cerebellar GM was estimated to be around 61 ​, which is shorter than typical values reported for white matter and grey matter at 3T (Wansapura et al., 1999). The decrease in suggests that the T2 of the remaining compartments are longer; this is illustrated by a simple bi-exponential simulation in Supplementary Fig. 4. The relatively short T2 highlights the importance of achieving a short TE to be able to detect the signal arising from this component.

The cerebellum has an important role in motor coordination, but it is becoming increasingly apparent that it also has an active role in cognition and emotion (O’Halloran et al., 2012; Tedesco et al., 2011). The neurons in the cerebellar cortex are highly organised, consisting of densely-packed granule cells and larger Purkinje cells with a cloud of dendritic spines. The cell bodies of the Purkinje cells are quite large, about 25–40 (Herndon, 1963), whereas the granular cell bodies are much smaller (7–10 (Stuart et al., 2016)) yet larger than the typical diameter of axons. One can speculate that the isotropically-restricted signal comes from within small spaces that may be intra-cellular (e.g. granule cell soma with low intra-cellular diffusivity or dendritic spines) or extra-cellular (e.g. between densely packed granule cells) or both. To elucidate the biological underpinnings of the observed signal, our currents efforts focus on studying the signal in patients with known cerebellar cell loss, such as spinocerebellar ataxia type 2 (SCA2) which affects granule- and Purkinje cells (Tada et al., 2015). Our preliminary results show that parameter estimates of the simple representation used in this study (Eq [3], Fig. 5). are affected beyond the inter- and intra-subject variability in healthy controls, placing some confidence in the use of this representation as a biomarker. Furthermore, future work could correlate changes in the STE signal with cell loss quantified by histology in cerebellar knockout mice. Altogether, studying the STE signal provides exciting avenues for gaining further insight into changes in tissue microstructure in disorders associated with the cerebellum, in addition to the suite of existing dMRI contrasts. dMRI studies have already shown changes in ataxia (Dayan et al., 2016; Salvatore et al., 2014), Parkinson’s disease, and Alzheimer’s disease (Mormina et al., 2017), where metrics such as mean diffusivity and diffusion-tensor (DT)-derived fractional anisotropy (FA) were studied. These studies mostly focused on cerebellar WM (e.g. peduncles). Recently, measures beyond the DT have been derived in cerebellar WM and GM, with the aim of being more specific to different compartments and the underlying neurobiology (Savini et al., 2018). Fig. 5b shows spatial variability in the estimated parameter maps. In future work we aim to look at the variability across and within different lobules, by registration to atlases (Diedrichsen et al., 2009).

SNR and spatial resolution

The spatial resolution used here to achieve the necessary SNR (i.e. voxel size 4 ​× ​4 ​× ​4 ) is relatively coarse, especially if one tries to study highly curved structures such as the cerebellum grey matter. We have aimed to reduce the effect of partial voluming in two ways: 1) registering the T1 segmentation to the diffusion data, and discarding voxels in which the tissue fraction was below 90%, and 2) further reducing partial volume effects with CSF by estimating the S(0) signal from diffusion-weighted images, a strategy proposed by (Baron and Beaulieu, 2015). We note that ‘conventional’ resolutions in diffusion MRI (e.g. 2 ​mm isotropic) are already too large to capture the fine folding in the cerebellum grey matter, and we therefore focus our efforts on achieving high SNR. In addition to coarse voxels, the SNR was further increased (e.g. compared to the study of Dhital et al., which used the same spatial resolution) by shortening TE. Super-resolution and gSlider-SMS (Setsompop et al., 2018) diffusion acquisitions provide exciting future avenues for increasing the spatial resolution while maintaining sufficient SNR. Compared with LTE super-resolution reconstruction, which has been extended recently to incorporate the angular relation between different diffusion measurements (Van Steenkiste et al., 2016), STE super-resolution would theoretically be more straightforward as the need to vary the orientation of the principal eigenvectors of the B-tensor is obviated.

Pre-processing

The low SNR of the STE data at high b-values made the pre-processing of the data challenging. dMRI pre-processing pipelines typically include motion correction and geometric distortion correction. The geometric distortions generally include those resulting from eddy currents and susceptibility differences, and the use of strong gradients requires an additional step to correct for any possible geometric distortions arising from gradient nonlinearities. Subject motion and eddy-current geometric distortions in high b-value data are often corrected for using a prediction-based framework (Andersson et al., 2017; Ben-Amitay et al., 2012); high b-value images are predicted from the corrected low b-value images, and the acquired high b-value images are subsequently registered to the predicted images. Strategies to predict high b-value data with different B-tensors from low b-value data are available (Nilsson et al., 2015), but the deformations allowed at high b-values have to be fairly constrained because only a relatively low signal can be observed in only few regions. When applying tools optimised for LTE images and/or moderate b-value STE images, we observed suspiciously large deformations in the high b-value STE data that could not be verified. In this study, we therefore opted for a conservative strategy where we acquired interleaved b0 images (every 15th image) to correct for subject motion in STE data. This necessarily led to differences in the processing of LTE and STE data; i.e., the STE data were only corrected with a rigid transformation which cannot account for higher order deformations e.g., due to eddy currents. While, theoretically, the eddy current deformations between STE images of the same b-value should be similar, future work should be attributed to optimising the processing of high b-value STE data. Future work will furthermore focus on collecting complementary information by means of real-time motion tracking (Maclaren et al., 2012; Zaitsev et al., 2006) – e.g. optical tracking (Qin et al., 2009)) and dynamic field measurements (De Zanche et al., 2008) – to provide robust correction for subject motion and geometrical distortions in these data.

In this work we corrected the images for geometric distortions arising from gradient nonlinearities, but gradient nonlinearities additionally cause spatiotemporally varying B-tensors. Strategies have been developed to take this into consideration, which were mostly evaluated on data acquired with Stejskal-Tanner encoding (Bammer et al., 2003; Glasser et al., 2013; Jones et al., 2018; Rudrapatna et al., 2018). Here we have taken gradient nonlinearities into account by computing the voxel-wise effective gradients (Fig. 8) and the corresponding effective B-tensor (Fig. 5, Fig. 7, and Table 1) (Bammer et al., 2003). Future work will be attributed to more thoroughly investigating the effect of gradient nonlinearities on the signal arising from free waveforms.

Correcting for the Rician noise bias is of importance here to obtain accurate estimates of the parameters in Eqs. [3-4] when using least-squares optimisation. The data was reconstructed using adaptive combine, which has shown to approach an effective number of coils of 1 (Sakaie and Lowe, 2017). Supplementary Fig. 3b shows an example of the background signal distribution in one subject before and after debiasing, which approximate Rician and Gaussian distributions respectively. The PIESNO method of (Koay et al., 2009b) for noise estimation identified background voxels that were minimally affected by sources of signal instabilities (e.g. ghosting). For debiasing we used the approach of Koay et al. (2009a) which relies only on magnitude data. When phase data are available, this can alternatively be leveraged to obtain Gaussian-distributed data (Eichner et al., 2015; Pizzolato et al., 2016). The process of Rician debiasing can yield signal estimates below the noise floor. To evaluate the accuracy and precision of this approach, we applied the same debiasing step as described in Section 3.2 to the simulated data of Fig. 1 (using the same acquisition protocol as in the in vivo data). Supplementary Fig. 2 shows estimates of for different SNR, before and after Rician debiasing. Indeed, estimation before Rician debiasing results in overestimation of . The expectation value of the error term in the case of nonlinear least squares and Rician-distributed data has been shown to converge to zero relatively slowly as a function of SNR (Veraart et al., 2013), which means that estimates can still be biased even if the SNR is larger than 2. After debiasing, our simulations indicate that signal estimates below the noise floor likely have a negative bias, resulting in a slight underestimation of . This, together with other potential inaccuracies in the signal and noise estimates (e.g. if the noise is non-stationary) could have caused the negative dot fraction estimates in e.g. the dGM (Table 1). The bias from the Rician debiasing step, however, converges to zero faster than the bias if no Rician debiasing would be performed. The noise floor estimates are reported as an upper limit on the estimated dot fraction in Table 1.

In addition to noise, other factors can cause signal variations across DWIs. For example, Nyquist ghosting, incomplete chemical shift suppression, or Gibbs ringing can produce signal errors. While in this study we aimed to correct for Gibbs ringing, other sources of variance may still have affected the signal. In the process of drawing the ROIs, we have carefully avoided regions of ghosting and incomplete chemical shift suppression to minimise their effects on the estimates. The tSNR in Table 1 reflects variability across the b0 images. This variability may be amplified by contributions from e.g. perfusion, subject motion, and physiological pulsations, which are in turn exacerbated in acquisitions with partial Fourier encoding. This may explain the difference between tSNR estimates in the brain and SNR estimates from the background (Table 1). In addition, non-stationarity of the noise may further contribute to this difference; although estimating the noise standard deviation at different spatial neighborhoods in the background did not reveal a strong variation in the estimated noise. The estimation of non-stationary noise is challenging but developments in this field can yield more accurate upper bounds on the dot fraction.

Conclusion

In this work, we combined ultra-strong gradients and efficient spherical tensor encoding to study the isotropic dMRI signal at ultra-high b-values, targeting the dot-compartment. Ultra-strong gradients allowed us to significantly reduce the TE, and therefore increase SNR, when acquiring data at high b-values. We further optimised encoding efficiency and TE by using asymmetric gradient waveforms instead of pulsed-gradients. A dot-compartment with zero diffusivity and no exchange would result in the signal plateauing for sufficiently high b-values; however, we found a signal significantly deviating from zero, yet still decaying across different WM regions and in the cerebellar GM. This observation is not in line with a spherical compartment of negligible size compared to the diffusion length and negligible exchange. We further studied the apparent diffusivity and signal fraction in the cerebellar GM assuming Gaussian diffusion and no exchange, finding these to be remarkably consistent across healthy controls. Future work will investigate the link between this hypothesised compartment and tissue microstructure, and investigate its potential as a biomarker in pathology affecting the cerebellar GM.

Declaration of competing interest

MN declares research support from Random Walk Imaging (formerly Colloidal Resource), and patent applications in Sweden (1250453–6 and 1250452–8), USA (61/642 ​594 and 61/642 ​589), and PCT (SE2013/050 ​492 and SE2013/050 ​493). The remaining authors declare no conflict of interest.