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

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 ¹H 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 1 Acquisition 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 R₂ by varying the value of τ_E. (b) Numerically optimized gradient waveforms (Sjölund et al., 2015) yielding four distinct b-tensor shapes ($b_{\mathrm{\Delta }}=-\mathrm{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 R₂ ($R_{\mathrm{2}}=\mathrm{1}/T_{\mathrm{2}}$ where T₂ 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):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{S\left({\mathit{\tau }}_{\mathrm{E}},\mathbf{b}\right)}{S_{\mathrm{0}}}}\\ & =\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\int }_{\mathbf{D}\in {\mathrm{Sym}}_{\mathrm{3}}^{+}}P(R_{\mathrm{2}},\mathbf{D})K({\mathit{\tau }}_{\mathrm{E}},\mathbf{b},R_{\mathrm{2}},\mathbf{D})\mathrm{d}\mathbf{D}\mathrm{d}{R}_{\mathrm{2}},\end{array}\end{array}$$

where P(R₂,D) is a joint probability distribution of R₂ and D, the kernel $K({\mathit{\tau }}_{\mathrm{E}},\mathbf{b},R_{\mathrm{2}},\mathbf{D})$ links the analysis space (R₂,D) to the acquisition space (τ_E, b), S₀ denotes the signal amplitude at (τ_E=0, b=0), and Sym${}_{\mathrm{3}}^{+}$ 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 R₂ > 60 s⁻¹ 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 (R₂≈70 s⁻¹) – 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 (R₂, 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, $D_{\mathrm{iso}}=(D_{\mathrm{|}\mathrm{|}}+\mathrm{2}{D}_{\perp })/\mathrm{3}$, and normalized diffusion anisotropy, $D_{\mathrm{\Delta }}=(D_{\mathrm{|}\mathrm{|}}-D_{\perp })/\mathrm{3}{D}_{\mathrm{iso}}$ (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(R₂, D_iso, D_Δ, θ, φ) distributions, whose dimensions directly correspond to those of the 5D acquisition space (τ_E, b, b_Δ, Θ, Φ):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\displaystyle \frac{S\left({\mathit{\tau }}_{\mathrm{E}},b,b_{\mathrm{\Delta }},\mathrm{\Theta },\mathrm{\Phi }\right)}{S_{\mathrm{0}}}}\\ & =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{1}/\mathrm{2}}{\overset{\mathrm{1}}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}K({\mathit{\tau }}_{\mathrm{E}},b,b_{\mathrm{\Delta }},\mathrm{\Theta },\mathrm{\Phi },R_{\mathrm{2}},D_{\mathrm{iso}},D_{\mathrm{\Delta }},\mathit{\theta },\mathit{\varphi })\\ & \times P\left(R_{\mathrm{2}},D_{\mathrm{iso}},D_{\mathrm{\Delta }},\mathit{\theta },\mathit{\varphi }\right)\mathrm{d}\mathit{\varphi }\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}{D}_{\mathrm{\Delta }}\mathrm{d}{D}_{\mathrm{iso}}\mathrm{d}{R}_{\mathrm{2}}.\end{array}\end{array}$$

The relaxation–diffusion encoding kernel is defined as

$$\begin{array}{}\text{(3)}& K\left(\mathrm{\dots }\right)=\exp \left(-{\mathit{\tau }}_{\mathrm{E}}{R}_{\mathrm{2}}\right)\exp \left(-b{D}_{\mathrm{iso}}\left[\mathrm{1}+\mathrm{2}{b}_{\mathrm{\Delta }}{D}_{\mathrm{\Delta }}{P}_{\mathrm{2}}\left(\cos \mathit{\beta }\right)\right]\right),\end{array}$$

where $P_{\mathrm{2}}\left(x\right)=(\mathrm{3}{x}^{\mathrm{2}}-\mathrm{1})/\mathrm{2}$ denotes the second Legendre polynomial, and β is the arc angle between the major symmetry axes of b and D, given by $\cos \mathit{\beta }=\cos \cos \mathit{\theta }+\cos (\mathrm{\Phi }-\mathit{\phi })\sin \mathrm{\Theta }\sin \mathit{\theta }$. According to Eq. (3), each (τ_E, b, b_Δ, Θ, Φ) coordinate establishes correlations across the separate dimensions of the R₂–D 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⁻¹ at the maximum slew rate of 200 T m⁻¹ s⁻¹. 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_{\mathrm{\Delta }}=\mathit{\{}-\mathrm{0.5},\mathrm{0.0},\mathrm{0.5},\mathrm{1.0}\mathit{\}}$), which were probed at varying combinations of echo times, b values, and b-tensor orientations. The waveforms giving $b_{\mathrm{\Delta }}=-\mathrm{0.5}$, 0.0, and 0.5 (see Fig. 1b) were calculated with a numerical optimization package (Sjölund et al., 2015) (https://github.com/jsjol/NOW, 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 $\mathit{\tau }=\mathrm{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 2 Representative 5D relaxation–diffusion encoded signals S(τ_E,b) and distributions P(R₂,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 R₂-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 R₂, isotropic diffusivity D_iso, and axial–radial diffusivity ratio $D_{\mathrm{|}\mathrm{|}}/D_{\perp }$. An auxiliary relaxation time T₂ scale was included along the log (R₂) axis to aid the inspection of the P(R₂,D) plots. The diffusion tensor orientation (θ, φ) is color-coded as [R,G,B] $=[\cos \mathit{\phi }\sin \mathit{\theta },\sin \mathit{\phi }\sin \mathit{\theta },\cos \mathit{\theta }]\cdot |D_{\mathrm{|}\mathrm{|}}-D_{\perp }|/max(D_{\mathrm{|}\mathrm{|}},D_{\perp })$, 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(R₂,D) distribution onto the respective 2D planes. Panels (b) and (c) display the signals S(τ_E,b) and corresponding P(R₂,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×10⁹ sm⁻², respectively), both $b_{\mathrm{\Delta }}=-\mathrm{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×10⁹ sm⁻²), 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×10⁹ sm⁻²). 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 $\mathrm{220}\times \mathrm{220}\times \mathrm{66}$ mm³ field of view, spatial resolution of $\mathrm{2}\times \mathrm{2}\times \mathrm{6}$ mm³, and a partial Fourier factor of 6∕8. Spatial resolution was sacrificed in favor of high signal-to-noise ratios (SNRs). The $\mathrm{2}\times \mathrm{2}\times \mathrm{6}$ mm³ 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(R₂, D_iso, 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(R₂/s⁻¹) < 1.5, −10 < log(D_||/m²s⁻¹) < −8.5, −10 < log(D_⊥/m²s⁻¹) < −8.5, 0 < cosθ < 1, 0 < φ < 2π) space. A discrete P(R₂,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 (R₂, 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 (R₂, D_||, D_⊥, θ, φ) configuration, and finally fitting the merged set is repeated a total of 20 times in order to find a P(R₂, D_||, D_⊥, θ, φ) distribution yielding a low residual sum of squares. Following 20 rounds, the resulting (R₂, D_||, D_⊥, θ, φ) configuration is selected, split, and subjected to a small random mutation. The original and mutated configurations are merged and a new P(R₂, 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 (R₂, D_||, D_⊥, θ, φ) configuration corresponding to the lowest sum of squared residuals. A final plausible P(R₂, D_||, D_⊥, θ, φ) solution is subsequently estimated at the end of the mutation cycle by selecting the 10 (R₂, 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(R₂, 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 (R₂, D_||, D_⊥, θ, φ) set is fitted to the signal, and the resulting P(R₂,D) is taken as the solution of the analyzed voxel. Finally, the P(R₂, D_||, D_⊥, θ, φ) distribution is mapped onto the (R₂, D_iso, 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(R₂, 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(R₂,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: https://github.com/JoaoPdAMartins/md-dmri (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 R₂–D 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 R₂ 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 R₂-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_{\mathrm{|}\mathrm{|}}/D_{\perp })$ 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 D_iso, low R₂, and D_|| $\sim D_{\perp }$; in contrast, GM and WM both exhibit lower D_iso and higher R₂, with WM being differentiated by its high $D_{\mathrm{|}\mathrm{|}}/D_{\perp }$. 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(R₂,D) components covering a significant range of the R₂–D 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_{\mathrm{|}\mathrm{|}}/D_{\perp })=\mathrm{0}$ plane. The bimodality of the GM distributions is an artifact attributed to the fact that prolate (D_Δ > 0, $D_{\mathrm{|}\mathrm{|}}/D_{\perp }\mathit{>}\mathrm{1})$ and oblate (D_Δ < 0, $D_{\mathrm{|}\mathrm{|}}/D_{\perp }$< 1) diffusion tensors with similar D_iso yield signal patterns that are only clearly discerned when D_Δ > 0.5 or, equivalently, $D_{\mathrm{|}\mathrm{|}}/D_{\perp }$ > 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 3 Statistical measures derived from the relaxation–diffusion distributions. The ensemble of 96 distinct P(R₂,D) solutions was used to calculate means E[x], variances Var[x], and covariances Cov[x,y] of all combinations of transverse relaxation rate R₂, isotropic diffusivity D_iso, and squared anisotropy ${D_{\mathrm{\Delta }}}^{\mathrm{2}}$. The statistical measures were all derived from the entire R₂-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 R₂–D 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 R₂–D distributions. A multitude of statistical functionals can be computed from the same distribution, meaning that the per-voxel P(R₂,D) ensembles generate a comprehensive set of distinct voxel-wise parameters. As shown in Fig. 3, the Monte Carlo realizations of P(R₂,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 R₂–D space (see Fig. 3). All of the statistical measures derived in this work parameterize diffusion tensor anisotropy with ${D_{\mathrm{\Delta }}}^{\mathrm{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[R₂]〉 and 〈E[D_iso]〉 clearly identify CSF-rich areas (low 〈E[R₂]〉 and high 〈E[D_iso]〉), while high $\langle \mathrm{E}\left[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\right]\rangle$ 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(R₂,D) distributions (displayed in Fig. 2c) comprise two populations at distant (R₂,D_iso) coordinates, and both voxels are thus characterized by high values of Var[R₂] and Var[D_iso] (see histograms of Fig. 3). As CSF and GM are both characterized by a low anisotropy, GM+CSF exhibits low values of Var[${D_{\mathrm{\Delta }}}^{\mathrm{2}}$]; in contrast, WM+CSF displays a significant dispersion along ${D_{\mathrm{\Delta }}}^{\mathrm{2}}$, which results in high Var[${D_{\mathrm{\Delta }}}^{\mathrm{2}}$] values. Covariance measures provide information about the correlations across the various dimensions of the R₂–D space. In WM+CSF distributions, for instance, higher values of diffusion anisotropy are correlated with higher R₂ and lower D_iso, which results in positive Cov[R₂,${D_{\mathrm{\Delta }}}^{\mathrm{2}}$] and negative Cov[D_iso,${D_{\mathrm{\Delta }}}^{\mathrm{2}}$]. The elevated 〈Var[R₂]〉 and 〈Var[D_iso]〉, and negative 〈Cov[R₂,D_iso]〉 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_{\mathrm{\Delta }}}^{\mathrm{2}}$]〉, positive 〈Cov[R₂,${D_{\mathrm{\Delta }}}^{\mathrm{2}}$]〉, and negative 〈Cov[D_iso,${D_{\mathrm{\Delta }}}^{\mathrm{2}}$]〉 locate WM+CSF voxels in those same regions, while low values of 〈Var[${D_{\mathrm{\Delta }}}^{\mathrm{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 R₂ and D_iso, WM+GM voxels result in nearly zero values of Var[R₂], Var[D_iso], Cov[D_iso,y], and Cov[R₂,y]. Instead, WM+GM voxels are signaled by finite values of 〈Var[${D_{\mathrm{\Delta }}}^{\mathrm{2}}$]〉, originated by the $\log (D_{\mathrm{|}\mathrm{|}}/D_{\perp })$ spread observed in the underlying R₂–D distribution (see the WM+GM distribution in Fig. 3c).

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

Figure 4 Parameter maps with bin-resolved means of the relaxation–diffusion distributions. (a) Division of the R₂-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 R₂-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[R₂]〉, isotropic diffusivity 〈E[D_iso]〉, squared anisotropy $\langle \mathrm{E}\left[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\right]\rangle$, and diffusion tensor orientation 〈E[Orientation]〉. The orientation maps (column 4) are color-coded as [R,G,B] $=[D_{xx},D_{yy},D_{zz}]/max(D_{xx},D_{yy},D_{zz})$, where D_ii 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(R₂,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 R₂ 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 5 Uncertainty estimation of the statistical measures derived from the relaxation–diffusion distributions. 3D density (color) scatter plots show the relationship between average initial signal intensity 〈S₀〉, the average of mean values derived from the R₂-D distributions 〈E[x]〉, and their corresponding uncertainties σ[E[x]]. For display purposes, signal intensity values were normalized to the maximum recorded 〈S₀〉, max(〈S₀〉). 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[R₂]〉 (row 1), isotropic diffusivity 〈E[D_iso]〉 (row 2), and squared anisotropy $\langle \mathrm{E}\left[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\right]\rangle$ (row 3) were computed from all voxels whose 〈S₀〉 was greater than 5 % of max(〈S₀〉). The resulting dataset comprises 55 327 voxels spread throughout all slices of the acquired 3D volume. The uncertainties of 〈E[R₂]〉, 〈E[D_iso]〉, and $\langle \mathrm{E}\left[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\right]\rangle$ correspond to the median absolute deviation between measures extracted from 96 independent solutions of Equation (2): σ[E[R₂]], σ[E[D_iso]), and $\mathit{\sigma }\left[\mathrm{E}\right[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\left]\right]$, respectively. All displayed data were derived from both the entire R₂-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 R₂, D_iso, and ${D_{\mathrm{\Delta }}}^{\mathrm{2}}$ are plotted against their respective uncertainties, σ[E[R₂]], σ[E[D_iso]], and $\mathit{\sigma }\left[\mathrm{E}\right[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\left]\right]$, and average signal amplitudes 〈S₀〉. Although the displayed statistical analysis is restricted to mean values, similar calculations can be done using any other scalar measure derived from the 5D R₂–D 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 $\langle S_{\mathrm{0}}\rangle /max(\langle S_{\mathrm{0}}\rangle )$ < 0.1 (see, for example, the D_iso scatterplots for the thin and thick populations). Despite the lower precision at low 〈S₀〉, the various average mean values are observed to be nearly constant throughout the $\langle S_{\mathrm{0}}\rangle /max(\langle S_{\mathrm{0}}\rangle )$ > 0.1 region; the only exception is $\langle \mathrm{E}\left[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\right]\rangle$ for the thin fraction, which shows a higher susceptibility to noise as evidenced by its positive correlation with 〈S₀〉.

Figure 6 Per-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[R₂]〉 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[R₂]〉 values exhibit either significant (second column) or nonsignificant (third column) differences. (b) Bin-resolved signal fractions (brightness) and average per-bin means (color) of R₂, and squared anisotropy ${D_{\mathrm{\Delta }}}^{\mathrm{2}}$. Regions 1 and 2 identify microstructural properties singled-out in the Results section. (c) Subdivision of the thick bin into three different R₂ subspaces. The contributions from different sub-bins are compared with a high-resolution R₁-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 R₂: 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 R₁-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 R₂ rates can be consistently detected in voxels containing GM+WM mixtures (see Fig. 6a), where conventional 1D R₂ 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 R₂ differences greater than the estimated uncertainties, thus providing evidence that the differentiation between the R₂ rates of the two bins is indeed a meaningful result.

All bin-resolved 〈E[R₂]〉 plots in Fig. 5 display a secondary cluster at high R₂ 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-R₂ components are observed to be concentrated around R₂=30 s⁻¹ (see Fig. 6a), the upper R₂-limit of the Monte Carlo inversion procedure. The “pile-up” of fast-relaxing contributions around the maximum allowed R₂ value is a well-known artifact of Laplace inversions (Saab et al., 1999).

The 〈E[R₂]〉 map of the thick bin features three main R₂ populations: high R₂ in the skull region (red voxels), low R₂ in peripheral brain regions (green voxels), and intermediate R₂ 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_{\mathrm{|}\mathrm{|}}/D_{\perp })$ < 0.6, −10 < log (D_iso/m² s⁻¹) < −8.7) subspace in three separate R₂ regions, and defined the “low” (−0.5 < log(R₂/s⁻¹) < 1.2), “medium” (1.2 < log(R₂/s⁻¹) < 1.4), and “high” (1.4 < log(R₂/s⁻¹) < 2) sub-bins of Fig. 6c. In T₂ 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 R₂=30 s⁻¹; the R₂=100 s⁻¹ 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 (R₁-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 R₂–D 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 R₂ 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 7 15 min protocol – bin-resolved signal contributions and mean parameter maps. (a) Map of average initial signal intensity 〈S₀〉 (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[R₂]〉, isotropic diffusivity 〈E[D_iso]〉, squared anisotropy $\langle \mathrm{E}\left[{D_{\mathrm{\Delta }}}^{\mathrm{2}}\right]\rangle$, 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 R₂–D 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 R₂–D 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 R₂–D 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 R₂ metrics derived from the abbreviated data. This is evidenced by unreasonably high R₂ rates in the ventricles (see the 〈E[R₂]〉 maps in Fig. 7b), and a higher difficulty in separating between the mean R₂ rates of the thin and thick bins. Only 65 % of mixed voxels from the abbreviated dataset show a meaningful R₂ separation, as opposed to the 75 % determined in the previous subsection. The lower resolution along the R₂ 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 R₂. 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 R₂ and diffusion tensor parameters (D_iso, 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 R₁-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 R₂–D properties provides complementary information to R₁-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 R₂ distribution protocols (Whittall et al., 1997), and its extraction is facilitated by the vast correlations across the full (D_iso, D_Δ, θ, φ) space (de Almeida Martins and Topgaard, 2018). We would like to reinforce that small R₂ differences can be observed despite the limited number and range of echo times sampled in this work; here, the separation between R₂ 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 R₁ 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 $\left|D_{\mathrm{\Delta }}\right|$ < 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_{\mathrm{\Delta }}}^{\mathrm{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×10⁹ sm⁻² (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 R₂ 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 R₂ values differ by less than a factor of 2 while the D_iso 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 R₂, isotropic diffusivities D_iso, normalized diffusion anisotropy D_Δ, and diffusion tensor orientation (θ, φ). The correlations allow model-free estimation of per-voxel relaxation–diffusion distributions P(R₂,D) that combine the chemical sensitivity of R₂ with the link between microstructure and the diffusion metrics. The rich information content of P(R₂,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⁻¹) results in an improvement in SNR of approximately 50 % compared to that achievable on the system used in this study (80 mT m⁻¹ 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.