Articles | Volume 1, issue 1
Research article
28 Feb 2020
Research article |  | 28 Feb 2020

Transferring principles of solid-state and Laplace NMR to the field of in vivo brain MRI

João P. de Almeida Martins, Chantal M. W. Tax, Filip Szczepankiewicz, Derek K. Jones, Carl-Fredrik Westin, and Daniel Topgaard

Magnetic resonance imaging (MRI) is the primary method for noninvasive investigations of the human brain in health, disease, and development but yields data that are difficult to interpret whenever the millimeter-scale voxels contain multiple microscopic tissue environments with different chemical and structural properties. We propose a novel MRI framework to quantify the microscopic heterogeneity of the living human brain as spatially resolved five-dimensional relaxation–diffusion distributions by augmenting a conventional diffusion-weighted imaging sequence with signal encoding principles from multidimensional solid-state nuclear magnetic resonance (NMR) spectroscopy, relaxation–diffusion correlation methods from Laplace NMR of porous media, and Monte Carlo data inversion. The high dimensionality of the distribution space allows resolution of multiple microscopic environments within each heterogeneous voxel as well as their individual characterization with novel statistical measures that combine the chemical sensitivity of the relaxation rates with the link between microstructure and the anisotropic diffusivity of tissue water. The proposed framework is demonstrated on a healthy volunteer using both exhaustive and clinically viable acquisition protocols.

1 Introduction

The structure of the brain is affected by both disease and normal development over a wide range of length scales. To measure and map the cellular architecture and molecular composition of the living human brain is a challenging experimental endeavor that promises far-reaching implications for both clinical diagnosis and our understanding of normal brain function. Over recent decades, magnetic resonance imaging (MRI) methods have been crucial for the progress of neuroanatomical studies (Lerch et al., 2017). Most clinical MRI applications rely on detecting 1H nuclei of water molecules to produce three-dimensional images with a spatial resolution on the millimeter scale. Even though the attainable resolution is clearly insufficient for direct observation of individual cells, chemical and microstructural features can be investigated by probing their effect on magnetic resonance observables such as nuclear relaxation rates (Halle, 2006) and the translational diffusivity (Le Bihan, 1995) of water. Relaxation and diffusion parameters can thus indirectly report on various microscopic properties, including cell density (Padhani et al., 2009), orientation of nerve fibers (Basser and Pierpaoli, 1996), and the presence of nutrients (Daoust et al., 2017). Current quantitative relaxation (Tofts, 2003) and diffusion (Jones, 2010) MRI observables are exquisitely sensitive to the cellular processes associated with knowledge acquisition (Zatorre et al., 2012), neuropsychiatric disorders (Kubicki et al., 2007), and different tumor types (Nilsson et al., 2018a), but suffer from poor specificity, and the same experimental data may support several distinct biological scenarios (Zatorre et al., 2012).

More detailed information can be obtained by taking into account that each MRI voxel comprises hundreds of thousands of cells with potentially different properties, implying that the per-voxel signal may include contributions from multiple microenvironments with distinct values of the MRI observables. To resolve the various microenvironments within a single voxel remains a highly challenging problem of vital importance for the progression of quantitative MRI studies. The signals from heterogeneous materials are often approximated as integral transformations of nonparametric distributions of relaxation rates or diffusivities (Istratov and Vyvenko, 1999), which may be estimated by Laplace inversion of data acquired as a function of the relevant experimental variable (Whittall and MacKay, 1989). Within the context of human brain MRI, the components of the distributions have been assigned to water populations residing in specific tissue microenvironments such as myelin (Mackay et al., 1994) and tumors (Laule et al., 2017). The power to resolve and individually characterize the different components can be boosted by combining multiple relaxation- and diffusion-encoding blocks and analyzing the data as joint probability distributions of the relevant observables (English et al., 1991). These ideas follow the principles of multidimensional nuclear magnetic resonance (NMR) spectroscopy and form the basis for multidimensional Laplace NMR which has become routine in the field of porous media (Galvosas and Callaghan, 2010; Song, 2013) and is now being combined with MRI (Zhang and Blumich, 2014; Benjamini and Basser, 2017). Recently, similar relaxation–diffusion correlation protocols have been translated to in vivo studies using model-based rather than nonparametric data inversion (De Santis et al., 2016; Veraart et al., 2017). So far, relaxation–diffusion correlation studies have relied on the Stejskal–Tanner experiment (Stejskal and Tanner, 1965), a pulsed gradient spin-echo (PGSE) sequence that has been in use for more than 50 years and where the signal is encoded for diffusion along a single axis using a pair of collinear gradient pulses. The limitations of the conventional experimental design become apparent when considering a white matter voxel comprising anisotropic domains with multiple orientations. When projected onto the measurement axis defined by the magnetic field gradients, the combination of diffusion anisotropy and orientation dispersion gives rise to a broad distribution of effective diffusivities (Topgaard and Söderman, 2002) that is challenging to retrieve with nonparametric Laplace inversion and, most importantly, impossible to differentiate from a spread of isotropic diffusivities (Mitra, 1995). Consequently, despite the fact that the relaxation–diffusion correlation yields more detailed information than conventional quantitative MRI, the inherent limitations of the Stejskal–Tanner experiment prevent unambiguous discrimination between isotropic and anisotropic contributions to the diffusivity distributions as well as model-free resolution of tissue microenvironments for heterogeneous anisotropic materials such as brain tissue.

We have recently shown that data acquisition and processing schemes for correlating isotropic and anisotropic nuclear interactions in multidimensional solid-state NMR spectroscopy (Schmidt-Rohr and Spiess, 1994) can be translated to diffusion NMR (de Almeida Martins and Topgaard, 2016), relaxation–diffusion correlation NMR (de Almeida Martins and Topgaard, 2018), and diffusion MRI (Topgaard, 2019), yielding nonparametric diffusion tensor distributions (Jian et al., 2007) with resolution of multiple isotropic and anisotropic diffusion components. These “multidimensional diffusion MRI” methods (Topgaard, 2017) rely on varying both the amplitude and orientation of the magnetic field gradients within a single encoding block in order to mimic the effects of sample reorientation (Frydman et al., 1992) and rotor-synchronized radio frequency pulse sequences (Gan, 1992) in multidimensional solid-state NMR to target specific aspects of the tensorial property being investigated. Here, we incorporate these ideas into a clinically feasible relaxation–diffusion correlation MRI protocol to quantify the microscopic heterogeneity of the living human brain. The suggested acquisition and analysis protocols resolve tissue heterogeneity on a five-dimensional space of transverse relaxation rates and axisymmetric diffusion tensors that report on the underlying chemical composition and microscopic geometry. Nonparametric relaxation–diffusion distributions are obtained for each voxel in the three-dimensional image using Monte Carlo data inversion to deal with the nonuniqueness of the Laplace inversion and estimate the uncertainty of quantitative parameters derived from the distributions (Prange and Song, 2009). Subvoxel tissue environments are resolved without limiting assumptions on the number or properties of the individual components and are characterized with statistical measures that have intuitive relations with the local microstructure.

Figure 1Acquisition protocol for 5D relaxation–diffusion MRI. (a) Pulse sequence for acquiring images encoded for relaxation and diffusion in a 5D space defined by the echo time τE, and b-tensor trace b, normalized anisotropy bΔ, and orientation (Θ, Φ). An EPI image readout block acquires the spin echo produced by slice-selective 90 and 180 radio-frequency pulses. The 180 pulse is encased by a pair of gradient waveforms allowing for diffusion encoding according to principles from multidimensional solid-state NMR (Topgaard, 2017) (red, green, and blue lines). The signal is encoded for the transverse relaxation rate R2 by varying the value of τE. (b) Numerically optimized gradient waveforms (Sjölund et al., 2015) yielding four distinct b-tensor shapes (bΔ=-0.5, 0.0, 0.5, and 1) (Eriksson et al., 2015).


2 Methods

2.1 Multidimensional relaxation–diffusion encoding

Figure 1a displays a pulse sequence wherein the signal S(τE,b) from a given voxel is encoded for information about the transverse relaxation rate R2 (R2=1/T2 where T2 is the transverse relaxation time) and diffusion tensor D by the experimental variables echo time τE and diffusion encoding tensor b according to de Almeida Martins and Topgaard (2018):

(1) S τ E , b S 0 = 0 + D Sym 3 + P ( R 2 , D ) K ( τ E , b , R 2 , D ) d D d R 2 ,

where P(R2,D) is a joint probability distribution of R2 and D, the kernel K(τE,b,R2,D) links the analysis space (R2,D) to the acquisition space (τE, b), S0 denotes the signal amplitude at (τE=0, b=0), and Sym3+ represents the mathematical space containing all 3×3 symmetric positive-definite matrices. The magnetic field gradient waveforms define an axially symmetric b-tensor that is parameterized by its trace (b), orientation (Θ, Φ), and normalized anisotropy (bΔ) (Eriksson et al., 2015), the latter controlling the influence of diffusion anisotropy on the detected signal in a manner corresponding to the effect of the angle between the main magnetic field and the rotor spinning axis in solid-state NMR (Frydman et al., 1992). While diffusion encoding performed by a conventional PGSE sequence is limited to a single b-tensor “shape” (bΔ=1), we have shown that variation of bΔ enables model-free separation and quantification of the isotropic and anisotropic contributions to the diffusion tensors (de Almeida Martins and Topgaard, 2016). In this work, we used the numerically optimized gradient waveforms displayed in Fig. 1b (Sjölund et al., 2015) to generate b-tensors at four distinct values of bΔ. In common with conventional diffusion MRI, our method requires a minimum echo time of ∼50 ms to accommodate diffusion encoding, causing the signal contributions from components with R2 > 60 s−1 to be reduced to less than 5 % of their initial amplitude. This means that the proposed protocol would require substantial signal averaging in order to quantify the fractions of fast relaxing components, thus precluding a mapping of myelin water (R2≈70 s−1) – one of the primary focuses of early multi-echo MRI methods (Mackay et al., 1994) – within a time compatible with either clinical practice or research.

Throughout the signal encoding process, the relaxation and diffusion of water are both affected by molecular exchange between chemically different sites and interactions with cell membranes. Averaging all these complex effects into sets of effective relaxation rates and apparent diffusion tensors, subvoxel composition can be reported as a collection of independent tissue microenvironments, each of which is characterized by a set of (R2, D) coordinates (de Almeida Martins and Topgaard, 2018). Assuming axial symmetry, the various microscopic diffusion tensors are parameterized by four independent dimensions: two eigenvalues corresponding to the axial and radial diffusivities, D|| and D, and the polar and azimuthal angles, θ and φ, describing the orientation of D relative to the laboratory frame of reference. The D|| and D diffusivities can be combined to define measures of isotropic diffusivity, Diso=(D||+2D)/3, and normalized diffusion anisotropy, DΔ=(D||-D)/3Diso (Eriksson et al., 2015), which report on the “size” and “shape” of the corresponding microscopic diffusion patterns (Topgaard, 2017). Tissue microscopic heterogeneity is therefore characterized with P(R2, Diso, DΔ, θ, φ) distributions, whose dimensions directly correspond to those of the 5D acquisition space (τE, b, bΔ, Θ, Φ):

(2) S τ E , b , b Δ , Θ , Φ S 0 = 0 0 - 1 / 2 1 0 π 0 2 π K ( τ E , b , b Δ , Θ , Φ , R 2 , D iso , D Δ , θ , ϕ ) × P R 2 , D iso , D Δ , θ , ϕ d ϕ sin θ d θ d D Δ d D iso d R 2 .

The relaxation–diffusion encoding kernel is defined as

(3) K ( ) = exp - τ E R 2 exp - b D iso 1 + 2 b Δ D Δ P 2 cos β ,

where P2(x)=(3x2-1)/2 denotes the second Legendre polynomial, and β is the arc angle between the major symmetry axes of b and D, given by cosβ=coscosθ+cos(Φ-φ)sinΘsinθ. According to Eq. (3), each (τE, b, bΔ, Θ, Φ) coordinate establishes correlations across the separate dimensions of the R2D space. Consequently, sampling various combinations of echo times and b-tensor parameters facilitates a comprehensive mapping of tissue-specific relaxation and diffusion properties.

2.2 MRI measurements

A healthy volunteer (female, 31 years) was scanned on a Siemens Magnetom Prisma 3T system equipped with a 20-channel-receiver head coil, and capable of delivering gradients of 80 mT m−1 at the maximum slew rate of 200 T m−1 s−1. The measurements were approved by a local Institutional Review Board (Partners Healthcare System), and the research subject provided written informed consent prior to participation.

Experimental data were acquired using the prototype spin-echo sequence (Lasič et al., 2014) and gradient waveforms shown in Fig. 1. The depicted waveforms give four distinct b-tensor anisotropies (bΔ={-0.5,0.0,0.5,1.0}), which were probed at varying combinations of echo times, b values, and b-tensor orientations. The waveforms giving bΔ=-0.5, 0.0, and 0.5 (see Fig. 1b) were calculated with a numerical optimization package (Sjölund et al., 2015) (, last access: 1 November 2019), including compensation for the effects of concomitant gradients (Szczepankiewicz et al., 2019). This procedure yielded a pair of asymmetric gradient waveforms lasting 30.8 and 25.0 ms, separated by approximately 8.0 ms. Linear encoding (bΔ=1) was implemented with two separate gradient waveforms: a symmetric bipolar gradient waveform whose encoding blocks lasted τ=25.1 ms and were separated by 8.0 ms (see Fig. 1b), and a pair of τ= 15.1 ms single-pulsed gradients bracketing a time period of 13.7 ms. The spectral profile of the bipolar gradient waveform was tuned to that of the asymmetric gradient waveforms in order to reduce the influence of time-dependent diffusion (Woessner, 1963; Callaghan and Stepišnik, 1996).

Figure 2Representative 5D relaxation–diffusion encoded signals S(τE,b) and distributions P(R2,D) for selected voxels in a living human brain. (a) Acquisition scheme showing τE, b, bΔ, Θ, and Φ as a function of acquisition point. (b) Experimental (gray circles) and fitted (black points) S(τE,b) signals from three representative voxels containing white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF). The presented signal data were acquired according to the scheme shown in panel (a) and is drawn with the same horizontal axis. (c) Nonparametric R2-D distributions obtained for both pure (WM, GM, CSF) and mixed (WM+GM, WM+CSF, GM+CSF) voxels. The discrete distributions are reported as scatter plots in a 3D space of the logarithms of the transverse relaxation rate R2, isotropic diffusivity Diso, and axial–radial diffusivity ratio D||/D. An auxiliary relaxation time T2 scale was included along the log (R2) axis to aid the inspection of the P(R2,D) plots. The diffusion tensor orientation (θ, φ) is color-coded as [R,G,B] =[cosφsinθ,sinφsinθ,cosθ]|D||-D|/max(D||,D), and the circle area is proportional to the statistical weight of the corresponding component. The contour lines on the sides of the plots represent projections of the 5D P(R2,D) distribution onto the respective 2D planes. Panels (b) and (c) display the signals S(τE,b) and corresponding P(R2,D), respectively, for the same WM, GM, and CSF voxels.


A total of 852 individual images were recorded at different combinations of (τE, b, bΔ, Θ, Φ) throughout the entire scan time of 45 min. The acquisition protocol is summarized in Fig. 2a. Briefly, bΔ=1 was acquired over 72 directions distributed over four b values (6, 10, 16, and 40 directions at b=0.1, 0.7, 1.4, and 2×109 sm−2, respectively), both bΔ=-0.5, and 0.5 were collected across 64 directions spread out over four b values (6, 10, 16, and 32 directions at, respectively, b=0.1, 0.7, 1.4, and 2×109 sm−2), and bΔ=0 was acquired for a single gradient waveform orientation, repeated 6 times over six b values (b=0.1, 0.3, 0.7, 1, 1.4, and 2×109 sm−2). For each (b, bΔ) coordinate, the set of directions was optimized using an electrostatic repulsion scheme (Bak and Nielsen, 1997; Jones et al., 1999). The various (b,bΔ, Θ, Φ) sets were then repeatedly acquired at three different echo times (τE=80, 110, and 150 ms) using the spectrally tuned waveforms. The nontuned Stejskal–Tanner waveform was used to acquire bΔ=1 data at τE=60 and 80 ms. Comparison between data acquired with the bipolar and the Stejskal–Tanner gradient waveforms at τE=80 ms allowed us to assess the validity of the Gaussian diffusion approximation (Callaghan and Stepišnik, 1996).

All images were recorded using a repetition time of 3 s, and an echo-planar readout with a 220×220×66 mm3 field of view, spatial resolution of 2×2×6 mm3, and a partial Fourier factor of 6∕8. Spatial resolution was sacrificed in favor of high signal-to-noise ratios (SNRs). The 2×2×6 mm3 anisotropic voxel configuration enables a large coverage with a minimal number of slices and yields axial maps with a high spatial resolution wherein anatomical features of interest can be easily identified. The acquired images were corrected for subject motion in ElastiX (Klein et al., 2009), using the extrapolated reference method detailed in Nilsson et al. (2015). Motion-corrected and non-motion-corrected data were then inverted using a quick 12-bootstrap procedure (see the following subsection for more details on the inversion), and the resulting parameter maps were subsequently compared. As no substantial differences were found between the results from the corrected and noncorrected datasets, we opted to not use motion correction in our final analysis. No denoising approaches were used prior to data inversion.

2.3 Nonparametric Monte Carlo inversion

Algorithms designed to solve Eq. (2) have been reviewed in both general (Istratov and Vyvenko, 1999) and magnetic resonance (Mitchell et al., 2012) literature. While classical inversion methods can be successfully used to estimate the 5D P(R2, Diso, DΔ, θ, φ) distribution, they become costly in terms of memory at the high dimensionality of our protocol. To circumvent this difficulty, we introduced an inversion approach wherein our correlation space is explored through a directed iterative algorithm, as explained in de Almeida Martins and Topgaard (2018). The algorithm starts by randomly selecting 200 points from the (0 < log(R2/s−1) < 1.5, −10 < log(D||/m2s−1) < −8.5, −10 < log(D/m2s−1) < −8.5, 0 < cosθ < 1, 0 < φ < 2π) space. A discrete P(R2,D) distribution is then estimated by solving a discretized version of Eq. (2) via a standard non-negative least squares (NNLS) algorithm (Lawson and Hanson, 1974). Points with nonzero weights are stored and merged with a new randomly generated set of 200 (R2, D||, D, θ, φ) points, and the weights of the merged set of points are found through a NNLS fit (Lawson and Hanson, 1974). The process of selecting points with nonzero weights, subsequently merging them with a random (R2, D||, D, θ, φ) configuration, and finally fitting the merged set is repeated a total of 20 times in order to find a P(R2, D||, D, θ, φ) distribution yielding a low residual sum of squares. Following 20 rounds, the resulting (R2, D||, D, θ, φ) configuration is selected, split, and subjected to a small random mutation. The original and mutated configurations are merged and a new P(R2, D||, D, θ, φ) distribution is determined by fitting the merged set to the data using the NNLS algorithm (Lawson and Hanson, 1974). The mutation and fitting procedure is repeated 20 times to find the local (R2, D||, D, θ, φ) configuration corresponding to the lowest sum of squared residuals. A final plausible P(R2, D||, D, θ, φ) solution is subsequently estimated at the end of the mutation cycle by selecting the 10 (R2, D||, D, θ, φ) points with the highest weights and performing a final NNLS fit.

The procedure described above is performed voxel-wise, resulting in an array of spatially resolved P(R2, D||, D, θ, φ) discrete distributions. Owing to the stochastic nature of the inversion protocol, we may fail at retrieving a nontrivial solution, which produces a small number of randomly located black voxels in the parameter maps. To correct for this, we combine the points from each voxel with the ones from its six nearest neighbors, subsequently fitting the set of 7×10 points to the underlying signal in order to find the 10 most likely points. The new (R2, D||, D, θ, φ) set is fitted to the signal, and the resulting P(R2,D) is taken as the solution of the analyzed voxel. Finally, the P(R2, D||, D, θ, φ) distribution is mapped onto the (R2, Diso, DΔ, θ, φ) space.

Following the works of Prange and Song (Prange and Song, 2009), we replace traditional regularization constraints (Whittall and MacKay, 1989) with an unconstrained Monte Carlo approach that estimates voxel-wise ensembles of N distinct P(R2, D) solutions consistent with the primary data (de Almeida Martins and Topgaard, 2018). In this study, we estimated ensembles of N=96 solutions per voxel. The level of dispersion within a given solution set characterizes the uncertainty of the inversion procedure and can thus be used to estimate the uncertainty of any quantities derived from P(R2,D) (Prange and Song, 2009; de Almeida Martins and Topgaard, 2018).

The nonparametric Monte Carlo inversion procedure was implemented in MATLAB and is publicly available in our GitHub repository: (last access: 25 February 2020) (Nilsson et al., 2018b). Inversion of the 45 min dataset took ∼72 h on a 12-Core Intel Xeon E5 2.7 GHz CPU, with a 64 GB DDR3 memory.

3 Results

3.1 Spatially resolved 5D relaxation–diffusion distributions

The proposed acquisition protocol translates into distinctive signal decay curves for each of the main components of the human brain. Indeed, voxels encompassing either white matter (WM), gray matter (GM), or cerebrospinal fluid (CSF) are all characterized by clearly distinct signal patterns (see Fig. 2b). The observed differences can be used to infer the gross R2D properties of the various cerebral constituents: WM signals are highly sensitive to both bΔ and (Θ, Φ), indicative of anisotropic diffusion along coherently aligned microscopic domains; GM signal patterns are rather insensitive to bΔ and (Θ, Φ), consistent with isotropic diffusion; and CSF data decays quickly with increasing b while remaining mostly unaffected by the other acquisition variables, features that suggest an isotropic medium characterized by relatively low R2 values. Voxels comprising mixtures of WM, GM, and/or CSF generate patterns that can be interpreted as a superposition of the signal data from the pure components.

Spatially resolved 5D R2-D nonparametric distributions are retrieved from the experimental data using the model-free inversion approach described in the Methods section. Figure 2c displays the solution ensembles for voxels containing WM, GM, and CSF, as well as combinations of those components: WM+GM, WM+CSF, and GM+CSF. Brain tissue possesses various microscopic components, whose relaxation and diffusion properties differ over various orders of magnitude. Therefore, tissue heterogeneity is more suitably described with logarithmic distributions, where pore anisotropy is parameterized with log(D||/D) instead of DΔ. The distinctive characters of the raw signal patterns in Fig. 2b result in unique voxel-wise distributions that capture the gross microscopic features of the main cerebral components. Namely, CSF is characterized by high Diso, low R2, and D|| D; in contrast, GM and WM both exhibit lower Diso and higher R2, with WM being differentiated by its high D||/D. As expected, voxels comprising mixtures of WM, GM, and CSF yield a linear combination of the distributions from the individual components.

Voxels containing pure GM or WM are characterized by clusters of P(R2,D) components covering a significant range of the R2D space. Because both tissue types comprise a plethora of cells with varying geometries or chemical compositions (e.g. axons with various amounts of myelin, dendrites, or glial cells), the observed spread may be interpreted as a direct consequence of the underlying cellular heterogeneity. However, similar broad distributions were also observed in spectroscopic multidimensional diffusion correlation measurements of discrete-component phantoms (de Almeida Martins and Topgaard, 2016, 2018), hinting that the solution spread additionally reflects the measurement and inversion uncertainty. This intrinsic uncertainty masks the effects of finer cellular details like the intra- and extra-axonal components modeled in previous diffusion-relaxation correlation MRI methods (Veraart et al., 2017).

As evidenced by Fig. 2c, pure GM voxels yield bimodal distributions that feature a nearly symmetric spread of components around the log(D||/D)=0 plane. The bimodality of the GM distributions is an artifact attributed to the fact that prolate (DΔ > 0, D||/D>1) and oblate (DΔ < 0, D||/D< 1) diffusion tensors with similar Diso yield signal patterns that are only clearly discerned when DΔ > 0.5 or, equivalently, D||/D > 4 (Eriksson et al., 2015). Diffusion tensor imaging (DTI) studies of the human cortex have revealed a low, yet non-negligible, diffusion anisotropy in cortical GM tissue (Assaf, 2018). The observation of both oblate and prolate components in the pure GM voxel is consistent with those findings, with the intrinsically low anisotropy preventing an unambiguous distinction between DΔ > 0 or DΔ < 0 solutions. The artifactual spread of anisotropic components is expected to worsen with the increase in experimental noise. Random signal fluctuations create small differences between data acquired at different bΔ values and consequently introduce a preference for anisotropic components with arbitrary DΔ sign. This effect is similar to the “eigenvalue repulsion” artifact in conventional DTI, where noise introduces a discrepancy in the eigenvalues of the voxel-averaged diffusion tensor that in turn gives rise to a positive bias in anisotropy (Pierpaoli and Basser, 1996; Jones and Cercignani, 2010).

Figure 3Statistical measures derived from the relaxation–diffusion distributions. The ensemble of 96 distinct P(R2,D) solutions was used to calculate means E[x], variances Var[x], and covariances Cov[x,y] of all combinations of transverse relaxation rate R2, isotropic diffusivity Diso, and squared anisotropy DΔ2. The statistical measures were all derived from the entire R2-D distribution space on a voxel-by-voxel basis. Histograms are used to represent the parameter sets calculated for three voxels containing binary mixtures of white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF). Each histogram comprises 96 estimates of a single statistical measure. The averages of statistical measures, 〈E[x]〉, 〈Var[x]〉 and 〈Cov[x,y]〉, are displayed as parameter maps whose color scales are given by the bars along the abscissas of the histograms. The crosses and arrows identify the heterogeneous voxels analyzed in the histograms; notice that the signaled points correspond to the average (as measured by the median) of the ensembles of plausible solutions shown in the histograms.


3.2 Statistical measures of tissue heterogeneity

The R2D distribution ensembles provide a wealth of information that is challenging to visualize in spatially resolved datasets with large image matrices. Drawing inspiration from the field of porous media, where ensembles of distributions have been converted into ensembles of scalar parameters such as total porosity or a fraction of bound fluid (Prange and Song, 2009), we extract statistical measures from the R2D distributions. A multitude of statistical functionals can be computed from the same distribution, meaning that the per-voxel P(R2,D) ensembles generate a comprehensive set of distinct voxel-wise parameters. As shown in Fig. 3, the Monte Carlo realizations of P(R2,D) are translated into ensembles of statistical measures, with 96 individual estimates being extracted for each measure. For compactness, the ensembles of statistical parameters are reduced to an average 〈⋅〉 and a dispersion measure σ[⋅] that is interpreted as the uncertainty of the estimated functional (Prange and Song, 2009). To render the results more robust to outliers, we report 〈⋅〉 as the ensemble median and estimate σ[⋅] as a median absolute deviation. The calculation of averages (as measured by the median) reduces the underlying ensemble of solutions into a single scalar and allows us to convey intra-voxel composition with parameter maps of average mean values 〈E[x]〉, average variances 〈Var[x]〉, and average covariances 〈Cov[x,y]〉 of all the relevant dimensions of the 5D R2D space (see Fig. 3). All of the statistical measures derived in this work parameterize diffusion tensor anisotropy with DΔ2 rather than DΔ; this is motivated by the intrinsic difficulty of distinguishing between prolate and oblate tensors (Eriksson et al., 2015).

The three maps in the first column of Fig. 3 provide a rough spatial overview of the principal tissue types: 〈E[R2]〉 and 〈E[Diso]〉 clearly identify CSF-rich areas (low 〈E[R2]〉 and high 〈E[Diso]〉), while high E[DΔ2] values separate WM from the two other main cerebral tissues. However, mean parameter maps alone cannot identify or characterize intra-voxel heterogeneity, and their use should be complemented with dispersion measures including, but not limited to, the (co)variance elements displayed in columns 2 and 3 of Fig. 3. For example, voxels surrounding the ventricles do not show a truly distinctive feature in maps of mean values but are characterized by nonzero covariance matrix elements. To understand the origin of the nonzero values, let us focus on the WM+CSF and GM+CSF voxels indicated in Fig. 3. The corresponding P(R2,D) distributions (displayed in Fig. 2c) comprise two populations at distant (R2,Diso) coordinates, and both voxels are thus characterized by high values of Var[R2] and Var[Diso] (see histograms of Fig. 3). As CSF and GM are both characterized by a low anisotropy, GM+CSF exhibits low values of Var[DΔ2]; in contrast, WM+CSF displays a significant dispersion along DΔ2, which results in high Var[DΔ2] values. Covariance measures provide information about the correlations across the various dimensions of the R2D space. In WM+CSF distributions, for instance, higher values of diffusion anisotropy are correlated with higher R2 and lower Diso, which results in positive Cov[R2,DΔ2] and negative Cov[Diso,DΔ2]. The elevated Var[R2] and Var[Diso], and negative Cov[R2,Diso] values found in the ventricular regions are thus interpreted as a product of subvoxel combinations of CSF with other components. A combination of high Var[DΔ2], positive Cov[R2,DΔ2], and negative Cov[Diso,DΔ2] locate WM+CSF voxels in those same regions, while low values of Var[DΔ2] indicate the existence of deep gray matter in the vicinity of the ventricles.

The maps displayed in Fig. 3 can also be used to identify voxels containing WM+GM mixtures. Because WM and GM distributions are characterized by similar values of R2 and Diso, WM+GM voxels result in nearly zero values of Var[R2], Var[Diso], Cov[Diso,y], and Cov[R2,y]. Instead, WM+GM voxels are signaled by finite values of Var[DΔ2], originated by the log(D||/D) spread observed in the underlying R2D distribution (see the WM+GM distribution in Fig. 3c).

Table 1R2-D limits of the “big”, “thin”, and “thick” bins.

Download Print Version | Download XLSX

Figure 4Parameter maps with bin-resolved means of the relaxation–diffusion distributions. (a) Division of the R2-D distribution space into different bins. The distribution space was separated into three bins (gray volumes) named “big”, “thin”, and “thick” that loosely capture the diffusion features of cerebrospinal fluid CSF, white matter WM, and gray matter GM, respectively. The 3D scatter plots display the nonparametric R2-D distributions corresponding to the CSF (top), WM (middle), and GM (bottom) voxels selected in Fig. 2. Superquadric tensor glyphs are used to illustrate the representative D captured by each bin. (b) Parameter maps of average per-bin means (color) of transverse relaxation rate 〈E[R2]〉, isotropic diffusivity 〈E[Diso]〉, squared anisotropy E[DΔ2], and diffusion tensor orientation 〈E[Orientation]. The orientation maps (column 4) are color-coded as [R,G,B] =[Dxx,Dyy,Dzz]/max(Dxx,Dyy,Dzz), where Dii are the diagonal elements of laboratory-framed average diffusion tensors estimated from the various distribution bins. Brightness indicates the signal fractions corresponding to the big (row 1), thin (row 2), and thick (row 3) bins. The white arrows identify deep gray-matter structures.


3.3 Bin-resolved metrics of tissue heterogeneity

A more detailed picture of intra-voxel heterogeneity is obtained by dividing the distribution space into smaller subspaces (“bins”). In line with early diffusion MRI works (Pierpaoli et al., 1996), we define three bins that loosely capture the diffusion properties of the P(R2,D) distributions from the main brain components (see Table 1 and Fig. 4a). The big bin contains CSF contributions, whereas the “thin” and “thick” bins capture the signal fractions from WM and GM, respectively. The names big, thin, and thick are inspired by the geometric properties of the microscopic diffusion tensors that are captured by each individual bin. Visual inspection of Fig. 4b reveals that the spatial distributions of the three bins are consistent with the expected distributions of the corresponding tissues, providing more evidence that the coarsely defined bins allow a separation of the main cerebral constituents. Parameter maps of the per-bin means of the relaxation and diffusion properties are more straightforwardly interpreted than the heterogeneity measures derived from the entire distribution space: for example, the deep gray matter inferred in the previous paragraph is easily identifiable at the center (white arrows) of the thick maps of Fig. 4b. Further, the correlations across the various dimensions of the diffusion space allow the resolution of subtle differences in relaxation rates. Focusing on the first column of Fig. 4b, we notice that the thick fraction exhibits a slightly lower R2 rate than that of the thin fraction. This behavior is in accordance with the previous literature (Tofts, 2003) and is consistently observed across the entire slice.

Figure 5Uncertainty estimation of the statistical measures derived from the relaxation–diffusion distributions. 3D density (color) scatter plots show the relationship between average initial signal intensity S0, the average of mean values derived from the R2-D distributions 〈E[x]〉, and their corresponding uncertainties σ[E[x]]. For display purposes, signal intensity values were normalized to the maximum recorded S0, max(〈S0〉). The contour lines on the side planes show 2D projections of the point density function defining the distribution of data points. The average mean values of transverse relaxation rate 〈E[R2]〉 (row 1), isotropic diffusivity 〈E[Diso]〉 (row 2), and squared anisotropy E[DΔ2] (row 3) were computed from all voxels whose S0 was greater than 5 % of max(〈S0〉). The resulting dataset comprises 55 327 voxels spread throughout all slices of the acquired 3D volume. The uncertainties of 〈E[R2]〉, 〈E[Diso]〉, and E[DΔ2] correspond to the median absolute deviation between measures extracted from 96 independent solutions of Equation (2): σ[E[R2]], σ[E[Diso]), and σ[E[DΔ2]], respectively. All displayed data were derived from both the entire R2-D space (column 1), and the “big” (column 2), “thin” (column 3), and “thick” (column 4) bins defined in Fig. 4a.


Global and bin-resolved averages for all the analyzed voxels of the entire 3D image matrix are compiled in Fig. 5, where per-voxel average means of R2, Diso, and DΔ2 are plotted against their respective uncertainties, σ[E[R2]], σ[E[Diso]], and σ[E[DΔ2]], and average signal amplitudes S0. Although the displayed statistical analysis is restricted to mean values, similar calculations can be done using any other scalar measure derived from the 5D R2D distributions. Examination of the scatter plots in Fig. 5 shows that microscopic populations with low signal fractions generate statistical measures with significantly higher uncertainties. While no immediate correlation is discerned between the estimated mean values and their corresponding uncertainty, the negative correlation between uncertainty and signal fractions introduces a significant dispersion of 〈E[x]〉 at S0/max(S0) < 0.1 (see, for example, the Diso scatterplots for the thin and thick populations). Despite the lower precision at low S0, the various average mean values are observed to be nearly constant throughout the S0/max(S0) > 0.1 region; the only exception is E[DΔ2] for the thin fraction, which shows a higher susceptibility to noise as evidenced by its positive correlation with S0.

Figure 6Per-bin relaxation properties and tissue composition. (a) Transverse relaxation properties specific to each of the “thin” (red) and “thick” (green) bins defined in Fig. 4a. The color-coded composite images (top) and histograms (bottom) display the fractional populations and average mean transverse relaxation values 〈E[R2]〉 of the two bins. The first column displays all of the thin and thick voxels, while the two other columns focus on thin+thick mixtures wherein the bin-specific 〈E[R2]〉 values exhibit either significant (second column) or nonsignificant (third column) differences. (b) Bin-resolved signal fractions (brightness) and average per-bin means (color) of R2, and squared anisotropy DΔ2. Regions 1 and 2 identify microstructural properties singled-out in the Results section. (c) Subdivision of the thick bin into three different R2 subspaces. The contributions from different sub-bins are compared with a high-resolution R1-weighted image segmented into four different tissues: white matter WM, cortical gray matter (GM), deep GM, and cerebrospinal fluid (CSF). Additive color maps display the spatial distribution of sub-bin fractions (from low to high R2: green, red, blue), and of cortical (green) and deep (red) GM. (d) Color-coded composite images showing the contributions of different bins (red = thin, green = thick, blue = big) and conventional R1-based segmentation labels (red = WM, green = cortical + deep GM, blue = CSF).


The minor differences between the relaxation rates of the thin and thick components are also observed in the scatter plots of Fig. 5. A more detailed analysis shows that distinct R2 rates can be consistently detected in voxels containing GM+WM mixtures (see Fig. 6a), where conventional 1D R2 distributions fail to resolve the subtle differences between components (Whittall et al., 1997). The second and third columns of Fig. 6a display mixed voxels, where the thin and thick populations each account for at least 30 % of the total measured signal. Approximately 75 % of the mixed voxels exhibit R2 differences greater than the estimated uncertainties, thus providing evidence that the differentiation between the R2 rates of the two bins is indeed a meaningful result.

All bin-resolved 〈E[R2]〉 plots in Fig. 5 display a secondary cluster at high R2 values. Inspection of Fig. 6b reveals that the fast relaxing cluster corresponds to the nonmasked extra-meningeal tissues and, for the thin fraction, to the pallidum (region 1 in Fig. 6b), a major component of the basal ganglia structures located deep in the brain. The contributions from the high-R2 components are observed to be concentrated around R2=30 s−1 (see Fig. 6a), the upper R2-limit of the Monte Carlo inversion procedure. The “pile-up” of fast-relaxing contributions around the maximum allowed R2 value is a well-known artifact of Laplace inversions (Saab et al., 1999).

The 〈E[R2]〉 map of the thick bin features three main R2 populations: high R2 in the skull region (red voxels), low R2 in peripheral brain regions (green voxels), and intermediate R2 values in the inner brain regions (yellow voxels). To more easily inspect the spatial distribution of the various populations within the thick bin we delimited the (−3.5 < log(D||/D) < 0.6, −10 < log (Diso/m2 s−1) < −8.7) subspace in three separate R2 regions, and defined the “low” (−0.5 < log(R2/s−1) < 1.2), “medium” (1.2 < log(R2/s−1) < 1.4), and “high” (1.4 < log(R2/s−1) < 2) sub-bins of Fig. 6c. In T2 units, the low, medium, and high bins correspond to 63 ms to 3.16 s, 40 to 63 ms, and 10 to 40 ms. Note that the true upper boundary of the high bin is set by the limits of the Monte Carlo inversion and is equal to R2=30 s−1; the R2=100 s−1 boundary is defined simply to render a more aesthetically pleasing plot (see Fig. 6c). The bin-resolved signal fraction maps were then compared with a high-resolution longitudinal relaxation-weighted (R1-weighted) image segmented in four tissue classes: WM, cortical GM, deep GM, and CSF. Figure 6c shows that the spatial distributions of the low and medium subfractions roughly correspond to the expected distributions of cortical GM and deep GM structures, respectively. Despite the similarities between bin-resolved and segmentation maps, the former possesses a grainier appearance and seems to miss a significant portion of deep GM tissue at the center of the slice. While the grainier aspect is caused by the higher noise of the R2D correlation dataset, the absence of central GM is explained by the presence of anisotropic tissues in structures such as the pallidum (region 1 in Fig. 6b) and the thalamus (region 2 in Fig. 6b). Those two deep GM structures are then contained within the thin bin, and not within the thick bin from which we defined the R2 subspaces. Joining the contributions of cortical and deep GM within a single tissue class offers further insight into the link between microscopic tissue composition and binning (see Fig. 6d). Comparing the three-tissue segmentation with maps of the big, thin, and thick fractions confirms that the pallidum and part of the thalamus are captured by the thin bin.

Figure 715 min protocol – bin-resolved signal contributions and mean parameter maps. (a) Map of average initial signal intensity S0 (top); subdivision of the diffusion space into the “big”, “thin”, and “thick” bins (middle); color-coded composite map of per-bin signal contributions (bottom). The colors in the bottom identify the fractions from different bins: [R,G,B] = [thin,thick,big]. (b) Parameter maps of average per-bin means (color) of transverse relaxation rate 〈E[R2]〉, isotropic diffusivity 〈E[Diso]〉, squared anisotropy E[DΔ2], and diffusion tensor orientation 〈E[Orientation]. The color and brightness of the various maps follows the same convention as Fig. 4b.


3.4 Clinical feasibility of the R2D correlation approach

The acquisition protocol discussed thus far can be inserted without further alterations in research studies of brain disease, where subjects are recruited for long scan sessions. However, the associated 45 min scan time impedes its use outside of a clinic-research setting. To assess the potential for clinical translation of the proposed framework, we compare the performance of the exhaustive 45 min protocol with that of an abbreviated protocol, compatible with the time frame of most clinical applications. To this end, we included two different 5D relaxation–diffusion MRI protocols in a single imaging session: the 45 min protocol described in the Methods section, and an abbreviated 15 min protocol whose details are contained in the Supplement. The two acquisition protocols were consecutively used without repositioning the volunteer.

The abbreviated dataset was inverted with the Monte Carlo algorithm described above. The resulting 5D R2D distributions and parameter maps are compiled in the Supplement. Figure 7 shows the bin-resolved parameter maps obtained with the 15 min acquisition protocol. Overall, the parameter maps derived from the abbreviated data resemble slightly noisier reproductions of the maps computed from the exhaustive protocol and provide the same conclusions. Namely, the big, thin, and thick bins demarcate the signal contributions from CSF, WM, and GM, respectively, and the main R2D properties of those same tissue types are accurately captured by the per-bin mean parameter maps. The most obvious difference between the two datasets is the lower quality of the R2 metrics derived from the abbreviated data. This is evidenced by unreasonably high R2 rates in the ventricles (see the 〈E[R2]〉 maps in Fig. 7b), and a higher difficulty in separating between the mean R2 rates of the thin and thick bins. Only 65 % of mixed voxels from the abbreviated dataset show a meaningful R2 separation, as opposed to the 75 % determined in the previous subsection. The lower resolution along the R2 dimension is most likely explained by the fact that the abbreviated protocol concentrates 85 % of its measurements within two unique values of τE, an acquisition scheme that is quite unspecific to dispersion along R2. In future experiments, we plan to address this issue by enforcing a more uniform distribution of data points along the various echo times.

4 Discussion and conclusions

The proposed framework resolves intra-voxel heterogeneity on a 5D space of transverse relaxation rates R2 and diffusion tensor parameters (Diso, DΔ, θ, φ). Per-voxel brain composition is broken down into a non-predefined number of microscopic environments with clearly distinct relaxation and diffusion properties. The heterogeneity within a voxel is thus resolved as linear combinations of independent microscopic components that can be assigned to local tissue environments; on a global scale, the subvoxel environments can be grouped into more general tissue classes. For healthy brain tissue, the detected microenvironments were classified into three broad bins whose diffusion properties respectively match those of the main constituents of the brain: WM, GM, and CSF. The separation between contributions from the three bins was observed to provide a clean 3D mapping of WM, GM, and CSF that agrees well with a conventional R1-based tissue segmentation. This demonstrates that the proposed protocol can indeed separate subvoxel tissue environments with different relaxation and diffusion properties; in the healthy human brain, the resolved environments can be coarsely assigned to contributions from CSF, WM, and GM (see Fig. 6d). The distinction between microscopic tissue environments with different R2D properties provides complementary information to R1-weighted segmentation and enables the resolution of tissue heterogeneity within a single anatomical structure, e.g. resolving anisotropic and isotropic regions within the thalamus.

The protocol presented in this work shows promise for neuroanatomy studies dealing with the resolution of specific microscopic features such as nerve fiber-tracking through heterogeneous voxels (Jeurissen et al., 2014) or free water mapping (Pasternak et al., 2009). Within a clinical setting, disentangling different tissue signals is expected to be useful for pathological conditions associated with intra-voxel tissue heterogeneity, e.g. tumor infiltration in surrounding brain tissue, inflammation of cerebral tissue, or replacement of myelin with free water. In the latter example, the proposed echo times lead to an almost complete decay of the signal contributions from myelin domains, meaning that the effects of axonal demyelination would have to be probed indirectly by tracking a reduction of the signal fraction from anisotropic subvoxel components.

Several approaches have been introduced in the diffusion MRI literature where subvoxel composition is investigated by devising signal models with increasingly complex priors and constraints (Wang et al., 2011; Zhang et al., 2012; Scherrer et al., 2016). While such models can be used to investigate the conditions mentioned in the above paragraph, the attained conclusions will be heavily dependent on the assumptions used to construct the model (Novikov et al., 2018). Hence, erroneous conclusions may be derived whenever the presupposed MR properties differ from the underlying microstructure (Lampinen et al., 2019). This limitation is alleviated in the present framework where subvoxel heterogeneity is quantified with nonparametric distributions that are retrieved from the data with minimal assumptions on the underlying tissue properties. Moreover, the vast majority of diffusion MRI models has been so far implemented with conventional Stejskal–Tanner sequences, which are known to convolve the signal contributions from DΔ and D orientation. Acquiring data at various bΔ has been shown to disentangle the effects of anisotropy and dispersion in D orientations (Eriksson et al., 2013, 2015), meaning that our 5D (τE, b) acquisition space is expected to provide a more clear component resolution whenever orientation dispersion is present.

Besides resolving the various microscopic domains within a voxel, we were also capable of observing subtle differences in component-specific relaxation rates. As mentioned before, this information is unattainable with classical multi-echo R2 distribution protocols (Whittall et al., 1997), and its extraction is facilitated by the vast correlations across the full (Diso, DΔ, θ, φ) space (de Almeida Martins and Topgaard, 2018). We would like to reinforce that small R2 differences can be observed despite the limited number and range of echo times sampled in this work; here, the separation between R2 components is mostly driven by the excellent resolution in the diffusion dimensions. The measurement of D-resolved transverse relaxation rates may complement previous work on tract-specific R1 rates (De Santis et al., 2016).

At the cellular level, the translational motion of water inside the human brain is influenced by interactions with macromolecules and partially permeable membranes forming compartments with barrier spacings ranging from nanometers for synaptic vesicles and myelin sheaths to micrometers for the plasma membranes of the axons. The diffusion of water during the 0.1 s timescale of MRI signal encoding is thus affected by a myriad of complex phenomena that are not explicitly accounted for in Eq. (2). Instead, we use the well-established approach of approximating the micrometer-scale water displacements as a distribution of anisotropic Gaussian contributions (Jian et al., 2007). The measured diffusivities may depend on the exact choice of experimental variables if the timing parameters of the gradient waveforms match the characteristic timescales of displacements between cellular barriers (Woessner, 1963) or molecular exchange between tissue environments with distinctly different diffusion properties (Kärger, 1969). By augmenting our acquisition protocol with an experimental dimension in which the spectral profiles of the gradient waveforms are comprehensively varied (Callaghan and Stepišnik, 1996; Lundell et al., 2019), microscopic barrier spacings could in principle be estimated by explicitly including the effects of restricted diffusion in the kernel of Eq. (2). Here we chose to minimize the influence of time dependence by designing waveforms with similar gradient-modulation spectra.

In the previous section, we mentioned that prolate (DΔ > 0) and oblate (DΔ < 0) diffusion tensors with |DΔ| < 0.5 result in similar signal decays (Eriksson et al., 2015). In the absence of orientational order, diffusion tensor anisotropy is detected as a deviation from a mono-exponential signal decay, which, to first order, is proportional to DΔ2 (Eriksson et al., 2015). Consequently, the magnitude of DΔ can be easily determined at moderate b values while the sign may require data acquired with b values up to 4×109 sm−2 (Eriksson et al., 2015) and echo times comparable to the ones registered in this work; currently, such acquisition parameters can only be achieved with a specialized scanner (Setsompop et al., 2013; Jones et al., 2018).

Resolving and separately characterizing intra- and extra-axonal compartments in brain tissue has been of long-standing interest in the MRI field (Does, 2018). Recently, Veraart et al. (2017) estimated subtle differences in R2 and diffusivity parameters for the intra- and extra-axonal components of human brain white matter by applying a constrained two-component model to data acquired with a conventional relaxation–diffusion correlation protocol relying on the Stejskal–Tanner experiment. The obtained R2 values differ by less than a factor of 2 while the Diso values are nearly identical and the DΔ values are 1 (by constraint) and approximately 0.5 for the intra- and extra-cellular compartments, respectively. Comparing with the nonparametric distributions in Fig. 2, we note that components with such similar properties would be virtually impossible to resolve in our minimally constrained approach despite the additional information added by the b-tensor shape dimension. The limited resolution is consistent with the fact that Eq. (2) states an ill-posed inverse problem accommodating multiple nonunique solutions – probably also including the one with two thin components as assumed by Veraart et al. We suggest that the unconstrained inversion could be used as a first analysis tool to define the boundaries of a more ambitious model incorporating additional information, e.g. from microanatomy studies that is not directly observable in the MRI data.

This work introduces and demonstrates a novel MRI framework, in which the microscopic heterogeneity of the living human brain is characterized via 5D correlations between the transverse relaxation rate R2, isotropic diffusivities Diso, normalized diffusion anisotropy DΔ, and diffusion tensor orientation (θ, φ). The correlations allow model-free estimation of per-voxel relaxation–diffusion distributions P(R2,D) that combine the chemical sensitivity of R2 with the link between microstructure and the diffusion metrics. The rich information content of P(R2,D) is reported through a set of 21 unique maps obtained by binning and parameter calculation in the 5D distribution space. Being specific to different tissue types while relying on few assumptions, the presented protocol shows promise for explorative neuroscience and clinical studies in which microscopic tissue composition cannot be presumed a priori. While the spatial resolution of the data acquired in this work was relatively limited, sacrificing resolution for SNR, there are several avenues to explore in the future, in hardware, acquisition, and analysis that will boost the SNR per unit time, thereby increasing the potential for improved resolution. From the hardware perspective, the use of ultra-high fields (7T and above), and ultra-strong field gradients (Setsompop et al., 2013; Jones et al., 2018), can boost SNR and reduce echo-time-per-unit-b value, respectively. For example, as noted in Jones et al. (2018), for bΔ=0 encoding, the shorter τE afforded by stronger gradients such as those available on a Connectom scanner (300 mT m−1) results in an improvement in SNR of approximately 50 % compared to that achievable on the system used in this study (80 mT m−1 gradients). From the acquisition perspective, multi-band acquisition schemes (Barth et al., 2016) can speed up overall acquisition times and facilitate a wide brain coverage with smaller voxel sizes. Moreover, replacing the rectilinear echo-planar readout (Turner et al., 1991) with a spiral readout (Wilm et al., 2017) can help to further reduce the echo time, boosting SNR which could be traded for higher spatial resolution. From the analysis side, as noted in the Methods section, no denoising approaches were applied here. Recent advances in denoising and/or joint reconstruction (Veraart et al., 2016; Bazin et al., 2019; Wang et al., 2019; Haldar et al., 2020) could further enhance the SNR, allowing resolution to be pushed higher. Finally, the presented framework can be merged with MRI fingerprinting methodology (Ma et al., 2013), whose pattern-matching algorithms may considerably boost the data inversion speed.

Code and data availability

The software analysis tools discussed in this paper are available for download from a public GitHub repository: (last access: 25 February 2020) (Nilsson et al., 2018b). The presented in vivo data may be directly requested from the authors.


The supplement related to this article is available online at:

Author contributions

DKJ, C-FW, and DT conceived the project. JPdAM, CMWT, and FS designed the acquisition protocol. FS and C-FW acquired the data. The nonparametric Monte Carlo algorithm was designed by JPdAM and DT, and the data analysis was performed by JPdAM, CMWT, and DT. JPdAM and DT wrote the manuscript, and all authors read and reviewed the manuscript.

Competing interests

Daniel Topgaard owns shares in and João P. de Almeida Martins is partially employed by the private company Random Walk Imaging AB (Lund, Sweden), which holds patents related to the described method. All other authors declare no competing interests.


The authors thank Scott Hoge for his assistance with the MRI measurements. João P. de Almeida Martins and Daniel Topgaard were financially supported by the Swedish Foundation for Strategic Research (AM13-0090, ITM17-0267) and the Swedish Research Council (2014-3910, 2018-03697). Chantal M. W. Tax is supported by a Rubicon grant (680-50-1527) from the Netherlands Organisation for Scientific Research (NWO). Filip Szczepankiewicz and Carl-Fredrik Westin are both supported by a National Institutes of Health grant (P41EB015902). Derek K. Jones is supported by a Wellcome Investigator Award (096646/Z/11/Z) and a Wellcome Strategic Award (104943/Z/14/Z).

Financial support

This research has been supported by the Stiftelsen för Strategisk Forskning (grant nos. AM13-0090 and ITM17-0267), the Vetenskapsrådet (grant nos. 2014-3910 and 2018-03697), the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant no. 680-50-1527), the National Institutes of Health (grant no. P41EB015902), and the Wellcome Trust (grant nos. 096646/Z/11/Z and 104943/Z/14/Z).

Review statement

This paper was edited by Markus Barth and reviewed by two anonymous referees.


Assaf, Y.: Imaging laminar structures in the gray matter with diffusion MRI, Neuroimage, 197, 677–688,, 2018. 

Bak, M. and Nielsen, N. C.: REPULSION, a novel approach to efficient powder averaging in solid-state NMR, J. Magn. Reson., 125, 132–139,, 1997. 

Barth, M., Breuer, F., Koopmans, P. J., Norris, D. G., and Poser, B. A.: Simultaneous multislice (SMS) imaging techniques, Magn. Reson. Med., 75, 63–81,, 2016. 

Basser, P. J. and Pierpaoli, C.: Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI, J. Magn. Reson. Ser. B, 111, 209–219,, 1996. 

Bazin, P.-L., Alkemade, A., van der Zwaag, W., Caan, M., Mulder, M., and Forstmann, B. U.: Denoising High-Field Multi-Dimensional MRI With Local Complex PCA, Front. Neurosci.-Switz., 13, 1066,, 2019. 

Benjamini, D. and Basser, P. J.: Magnetic resonance microdynamic imaging reveals distinct tissue microenvironments, Neuroimage, 163, 183–196,, 2017. 

Callaghan, P. T. and Stepišnik, J.: Generalized analysis of motion using magnetic field gradients, in: Advances in magnetic and optical resonance, Elsevier, 325–388, 1996. 

Daoust, A., Dodd, S., Nair, G., Bouraoud, N., Jacobson, S., Walbridge, S., Reich, D. S., and Koretsky, A.: Transverse relaxation of cerebrospinal fluid depends on glucose concentration, Magn. Reson. Imaging, 44, 72–81,, 2017. 

de Almeida Martins, J. P. and Topgaard, D.: Two-Dimensional Correlation of Isotropic and Directional Diffusion Using NMR, Phys. Rev. Lett., 116, 087601,, 2016. 

de Almeida Martins, J. P. and Topgaard, D.: Multidimensional correlation of nuclear relaxation rates and diffusion tensors for model-free investigations of heterogeneous anisotropic porous materials, Sci. Rep., 8, 2488,, 2018. 

De Santis, S., Barazany, D., Jones, D. K., and Assaf, Y.: Resolving relaxometry and diffusion properties within the same voxel in the presence of crossing fibres by combining inversion recovery and diffusion-weighted acquisitions, Magn. Reson. Med., 75, 372–380,, 2016. 

Does, M. D.: Inferring brain tissue composition and microstructure via MR relaxometry, Neuroimage, 182, 136–148,, 2018. 

English, A. E., Whittal, K. P., Joy, M. L. G., and Henkelman, R. M.: Quantitative two-dimensional time correlation relaxometry, Magn. Reson. Med., 22, 425–434,, 1991. 

Eriksson, S., Lasic, S., and Topgaard, D.: Isotropic diffusion weighting in PGSE NMR by magic-angle spinning of the q-vector, J. Magn. Reson., 226, 13–18,, 2013. 

Eriksson, S., Lasic, S., Nilsson, M., Westin, C. F., and Topgaard, D.: NMR diffusion-encoding with axial symmetry and variable anisotropy: Distinguishing between prolate and oblate microscopic diffusion tensors with unknown orientation distribution, J. Chem. Phys., 142, 104201,, 2015. 

Frydman, L., Chingas, G. C., Lee, Y. K., Grandinetti, P. J., Eastman, M. A., Barrall, G. A., and Pines, A.: Variable-angle correlation spectroscopy in solid-state nuclear magnetic resonance, J. Chem. Phys., 97, 4800–4808,, 1992. 

Galvosas, P. and Callaghan, P. T.: Multi-dimensional inverse Laplace spectroscopy in the NMR of porous media, C. R. Physique, 11, 172–180,, 2010. 

Gan, Z.: High-resolution chemical shift and chemical shift anisotropy correlation in solids using slow magic angle spinning, J. Am. Chem. Soc., 114, 8307–8309,, 1992. 

Haldar, J. P., Liu, Y., Liao, C., Fan, Q., and Setsompop, K.: Fast submillimeter diffusion MRI using gSlider-SMS and SNR-enhancing joint reconstruction, Magn. Reson. Med., in press,, 2020. 

Halle, B.: Molecular theory of field-dependent proton spin-lattice relaxation in tissue, Magn. Reson. Med., 56, 60–72,, 2006. 

Istratov, A. A. and Vyvenko, O. F.: Exponential analysis in physical phenomena, Rev. Sci. Instrum., 70, 1233–1257,, 1999. 

Jeurissen, B., Tournier, J. D., Dhollander, T., Connelly, A., and Sijbers, J.: Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data, Neuroimage, 103, 411–426,, 2014. 

Jian, B., Vemuri, B. C., Özarslan, E., Carney, P. R., and Mareci, T. H.: A novel tensor distribution model for the diffusion-weighted MR signal, Neuroimage, 37, 164–176,, 2007. 

Jones, D. K.: Diffusion MRI, Oxford University Press, 2010. 

Jones, D. K. and Cercignani, M.: Twenty-five pitfalls in the analysis of diffusion MRI data, NMR Biomed., 23, 803–820,, 2010. 

Jones, D. K., Horsfield, M. A., and Simmons, A.: Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging, Magn. Reson. Med., 42, 515–525,<515::AID-MRM14>3.0.CO;2-Q, 1999. 

Jones, D. K., Alexander, D. C., Bowtell, R., Cercignani, M., Dell'Acqua, F., McHugh, D. J., Miller, K. L., Palombo, M., Parker, G. J. M., Rudrapatna, U. S., and Tax, C. M. W.: Microstructural imaging of the human brain with a “super-scanner”: 10 key advantages of ultra-strong gradients for diffusion MRI, Neuroimage, 182, 8–38,, 2018. 

Kärger, J.: Zur Bestimmung der Diffusion in einem Zweibereichsystem mit Hilfe von gepulsten Feldgradienten, Ann. Phys., 479, 1–4,, 1969. 

Klein, S., Staring, M., Murphy, K., Viergever, M. A., and Pluim, J. P.: Elastix: a toolbox for intensity-based medical image registration, IEE Trans. Med. Imaging, 29, 196–205, 2009. 

Kubicki, M., McCarley, R., Westin, C.-F., Park, H.-J., Maier, S., Kikinis, R., Jolesz, F. A., and Shenton, M. E.: A review of diffusion tensor imaging studies in schizophrenia, J. Psychiatr. Res., 41, 15–30,, 2007. 

Lampinen, B., Szczepankiewicz, F., Noven, M., van Westen, D., Hansson, O., Englund, E., Martensson, J., Westin, C. F., and Nilsson, M.: Searching for the neurite density with diffusion MRI: Challenges for biophysical modeling, Hum. Brain Mapp., 40, 2529–2545,, 2019. 

Lasič, S., Szczepankiewicz, F., Eriksson, S., Nilsson, M., and Topgaard, D.: Microanisotropy imaging: quantification of microscopic diffusion anisotropy and orientational order parameter by diffusion MRI with magic-angle spinning of the q-vector, Front. Phys., 2, 11,, 2014. 

Laule, C., Bjarnason, T. A., Vavasour, I. M., Traboulsee, A. L., Moore, G. W., Li, D. K., and MacKay, A. L.: Characterization of brain tumours with spin–spin relaxation: pilot case study reveals unique T2 distribution profiles of glioblastoma, oligodendroglioma and meningioma, J. Neurol., 264, 2205–2214,, 2017. 

Lawson, C. L. and Hanson, R. J.: Solving least squares problems, Prentice-Hall, Englewood Cliffs, NJ, 1974. 

Le Bihan, D.: Molecular diffusion, tissue microdynamics and microstructure, NMR Biomed., 8, 375–386,, 1995. 

Lerch, J. P., van der Kouwe, A. J., Raznahan, A., Paus, T., Johansen-Berg, H., Miller, K. L., Smith, S. M., Fischl, B., and Sotiropoulos, S. N.: Studying neuroanatomy using MRI, Nat. Neurosci., 20, 314–326,, 2017. 

Lundell, H., Nilsson, M., Dyrby, T. B., Parker, G. J. M., Cristinacce, P. L. H., Zhou, F. L., Topgaard, D., and Lasič, S.: Multidimensional diffusion MRI with spectrally modulated gradients reveals unprecedented microstructural detail, Sci. Rep., 9, 9026,, 2019. 

Ma, D., Gulani, V., Seiberlich, N., Liu, K., Sunshine, J. L., Duerk, J. L., and Griswold, M. A.: Magnetic resonance fingerprinting, Nature, 495, 187–192,, 2013. 

Mackay, A., Whittall, K., Adler, J., Li, D., Paty, D., and Graeb, D.: In vivo visualization of myelin water in brain by magnetic resonance, Magn. Reson. Med., 31, 673–677,, 1994. 

Mitchell, J., Chandrasekera, T. C., and Gladden, L. F.: Numerical estimation of relaxation and diffusion distributions in two dimensions, Prog. Nucl. Magn. Reson. Spectrosc., 62, 34–50,, 2012. 

Mitra, P. P.: Multiple wave-vector extension of the NMR pulsed-field-gradient spin-echo diffusion measurement, Phys. Rev. B, 51, 15074–15078,, 1995. 

Nilsson, M., Szczepankiewicz, F., van Westen, D., and Hansson, O.: Extrapolation-Based References Improve Motion and Eddy-Current Correction of High B-Value DWI Data: Application in Parkinson's Disease Dementia, PLoS One, 10, e0141825,, 2015. 

Nilsson, M., Englund, E., Szczepankiewicz, F., van Westen, D., and Sundgren, P. C.: Imaging brain tumour microstructure, Neuroimage, 182, 232–250,, 2018a. 

Nilsson, M., Szczepankiewicz, F., Lampinen, B., Ahlgren, A., De Almeida Martins, J. P., Lasic, S., Westin, C.-F., and Topgaard, D.: An open-source framework for analysis of multidimensional diffusion MRI data implemented in MATLAB, in: Proceedings of the 26th Annual Meeting of ISMRM, Paris, France, 16–21 June 2018, 5355, 2018b. 

Novikov, D. S., Kiselev, V. G., and Jespersen, S. N.: On modeling, Magn. Reson. Med., 79, 3172–3193,, 2018. 

Padhani, A. R., Liu, G., Mu-Koh, D., Chenevert, T. L., Thoeny, H. C., Takahara, T., Dzik-Jurasz, A., Ross, B. D., Van Cauteren, M., Collins, D., Hammoud, D. A., Rustin, G. J. S., Taouli, B., and Choyke, P. L.: Diffusion-Weighted Magnetic Resonance Imaging as a Cancer Biomarker: Consensus and Recommendations, Neoplasia, 11, 102–125,, 2009. 

Pasternak, O., Sochen, N., Gur, Y., Intrator, N., and Assaf, Y.: Free Water Elimination and Mapping from Diffusion MRI, Magn. Reson. Med., 62, 717–730,, 2009. 

Pierpaoli, C. and Basser, P. J.: Toward a quantitative assessment of diffusion anisotropy, Magn. Res. Med., 36, 893–906,, 1996. 

Pierpaoli, C., Jezzard, P., Basser, P. J., Barnett, A., and Di Chiro, G.: Diffusion tensor MR imaging of the human brain, Radiology, 201, 637–648,, 1996. 

Prange, M. and Song, Y. Q.: Quantifying uncertainty in NMR T2 spectra using Monte Carlo inversion, J. Magn. Reson., 196, 54–60,, 2009. 

Saab, G., Thompson, R. T., and Marsh, G. D.: Multicomponent T2 relaxation of in vivo skeletal muscle, Magn. Reson. Med., 42, 150–157,<150::AID-MRM20>3.0.CO;2-5, 1999. 

Scherrer, B., Schwartzman, A., Taquet, M., Sahin, M., Prabhu, S. P., and Warfield, S. K.: Characterizing brain tissue by assessment of the distribution of anisotropic microstructural environments in diffusion-compartment imaging (DIAMOND), Magn. Reson. Med., 76, 963–977, 2016. 

Schmidt-Rohr, K. and Spiess, H. W.: Multidimensional solid-state NMR and polymers, Academic Press, 1994. 

Setsompop, K., Kimmlingen, R., Eberlein, E., Witzel, T., Cohen-Adad, J., McNab, J. A., Keil, B., Tisdall, M. D., Hoecht, P., Dietz, P., Cauley, S. F., Tountcheva, V., Matschl, V., Lenz, V. H., Heberlein, K., Potthast, A., Thein, H., Van Horn, J., Toga, A., Schmitt, F., Lehne, D., Rosen, B. R., Wedeen, V., and Wald, L. L.: Pushing the limits of in vivo diffusion MRI for the Human Connectome Project, Neuroimage, 80, 220–233,, 2013. 

Sjölund, J., Szczepankiewicz, F., Nilsson, M., Topgaard, D., Westin, C.-F., and Knutsson, H.: Constrained optimization of gradient waveforms for generalized diffusion encoding, J. Magn. Reson., 261, 157–168,, 2015. 

Song, Y. Q.: Magnetic resonance of porous media (MRPM): a perspective, J. Magn. Reson., 229, 12–24,, 2013. 

Stejskal, E. O. and Tanner, J. E.: Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient, J. Chem. Phys., 42, 288–292,, 1965. 

Szczepankiewicz, F., Westin, C. F., and Nilsson, M.: Maxwell-compensated design of asymmetric gradient waveforms for tensor-valued diffusion encoding, Magn. Reson. Med., 82, 1424–1437,, 2019. 

Tofts, P.: Quantitative MRI of the Brain: Measuring Changes Caused by Disease, John Wiley & Sons, 2003. 

Topgaard, D.: Multidimensional diffusion MRI, J. Magn. Reson., 275, 98–113,, 2017. 

Topgaard, D.: Diffusion tensor distribution imaging, NMR Biomed., 32, e4066,, 2019. 

Topgaard, D. and Söderman, O.: Self-diffusion in two-and three-dimensional powders of anisotropic domains: An NMR study of the diffusion of water in cellulose and starch, J. Phys. Chem. B, 106, 11887–11892,, 2002. 

Turner, R., Le Bihan, D., and Scott Chesnicks, A.: Echo-planar imaging of diffusion and perfusion, Magn. Reson. Med., 19, 247–253,, 1991. 

Veraart, J., Novikov, D. S., Christiaens, D., Ades-Aron, B., Sijbers, J., and Fieremans, E.: Denoising of diffusion MRI using random matrix theory, Neuroimage, 142, 394–406,, 2016. 

Veraart, J., Novikov, D. S., and Fieremans, E.: TE dependent Diffusion Imaging (TEdDI) distinguishes between compartmental T2 relaxation times, Neuroimage, 182, 360–369,, 2017. 

Wang, H., Zheng, R., Dai, F., Wang, Q., and Wang, C.: High-field mr diffusion-weighted image denoising using a joint denoising convolutional neural network, J. Magn. Reson. Imaging, 50, 1937–1947,, 2019. 

Wang, Y., Wang, Q., Haldar, J. P., Yeh, F.-C., Xie, M., Sun, P., Tu, T.-W., Trinkaus, K., Klein, R. S., Cross, A. H., and Song, S.-K.: Quantification of increased cellularity during inflammatory demyelination, Brain, 134, 3590–3601,, 2011.  

Whittall, K. P. and MacKay, A. L.: Quantitative interpretation of NMR relaxation data, J. Magn. Reson., 84, 134–152,, 1989. 

Whittall, K. P., Mackay, A. L., Graeb, D. A., Nugent, R. A., Li, D. K., and Paty, D. W.: In vivo measurement of T2 distributions and water contents in normal human brain, Magn. Reson. Med., 37, 34–43,, 1997. 

Wilm, B. J., Barmet, C., Gross, S., Kasper, L., Vannesjo, S. J., Haeberlin, M., Dietrich, B. E., Brunner, D. O., Schmid, T., and Pruessmann, K. P.: Single-shot spiral imaging enabled by an expanded encoding model: Demonstration in diffusion MRI, Magn. Reson. Med., 77, 83–91,, 2017. 

Woessner, D. E.: N.M.R. spin-echo self-diffusion measurements on fluids undergoing restricted diffusion, J. Phys. Chem., 67, 1365–1367,, 1963. 

Zatorre, R. J., Fields, R. D., and Johansen-Berg, H.: Plasticity in gray and white: neuroimaging changes in brain structure during learning, Nat. Neurosci., 15, 528–536,, 2012. 

Zhang, H., Schneider, T., Wheeler-Kingshott, C. A., and Alexander, D. C.: NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain, Neuroimage, 61, 1000–1016,, 2012. 

Zhang, Y. and Blumich, B.: Spatially resolved D-T2 correlation NMR of porous media, J. Magn. Reson., 242, 41–48,, 2014. 

Short summary
Detailed interpretation of brain MRI data is hampered by the fact that each imaging voxel comprises several types of cells and tissues. To address this, we adapt signal encoding and data inversion strategies from solid-state and low-field NMR to quantify the sub-voxel heterogeneity of the human brain with 5D relaxation–diffusion distributions wherein distinct tissue components are resolved, individually characterized, and subsequently mapped throughout the volume of the brain.