the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Transferring principles of solidstate 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
CarlFredrik Westin
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 millimeterscale 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 fivedimensional relaxation–diffusion distributions by augmenting a conventional diffusionweighted imaging sequence with signal encoding principles from multidimensional solidstate 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.
 Article
(6340 KB)  Fulltext XML

Supplement
(7053 KB)  BibTeX
 EndNote
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 farreaching 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 ^{1}H nuclei of water molecules to produce threedimensional 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 pervoxel 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 diffusionencoding 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 modelbased 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 spinecho (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 modelfree 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 solidstate NMR spectroscopy (SchmidtRohr 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 rotorsynchronized radio frequency pulse sequences (Gan, 1992) in multidimensional solidstate 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 fivedimensional 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 threedimensional 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.
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_{2} (${R}_{\mathrm{2}}=\mathrm{1}/{T}_{\mathrm{2}}$ where T_{2} 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):
where P(R_{2},D) is a joint probability distribution of R_{2} and D, the kernel $K({\mathit{\tau}}_{\mathrm{E}},\mathbf{b},{R}_{\mathrm{2}},\mathbf{D})$ links the analysis space (R_{2},D) to the acquisition space (τ_{E}, b), S_{0} denotes the signal amplitude at (τ_{E}=0, b=0), and Sym${}_{\mathrm{3}}^{+}$ represents the mathematical space containing all 3×3 symmetric positivedefinite matrices. The magnetic field gradient waveforms define an axially symmetric btensor 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 solidstate NMR (Frydman et al., 1992). While diffusion encoding performed by a conventional PGSE sequence is limited to a single btensor “shape” (b_{Δ}=1), we have shown that variation of b_{Δ} enables modelfree 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 btensors 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_{2} > 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 (R_{2}≈70 s^{−1}) – one of the primary focuses of early multiecho 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_{2}, 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_{2}, D_{iso}, D_{Δ}, θ, φ) distributions, whose dimensions directly correspond to those of the 5D acquisition space (τ_{E}, b, b_{Δ}, Θ, Φ):
The relaxation–diffusion encoding kernel is defined as
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 $\mathrm{cos}\mathit{\beta}=\mathrm{cos}\mathrm{cos}\mathit{\theta}+\mathrm{cos}(\mathrm{\Phi}\mathit{\phi})\mathrm{sin}\mathrm{\Theta}\mathrm{sin}\mathit{\theta}$. According to Eq. (3), each (τ_{E}, b, b_{Δ}, Θ, Φ) coordinate establishes correlations across the separate dimensions of the R_{2}–D space. Consequently, sampling various combinations of echo times and btensor parameters facilitates a comprehensive mapping of tissuespecific 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 20channelreceiver 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 spinecho sequence (Lasič et al., 2014) and gradient waveforms shown in Fig. 1. The depicted waveforms give four distinct btensor 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 btensor 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}=\phantom{\rule{0.25em}{0ex}}\mathrm{25}.$1 ms and were separated by 8.0 ms (see Fig. 1b), and a pair of τ= 15.1 ms singlepulsed 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 timedependent diffusion (Woessner, 1963; Callaghan and Stepišnik, 1996).
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^{9} sm^{−2}, 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^{9} 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×10^{9} 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 echoplanar readout with a $\mathrm{220}\times \mathrm{220}\times \mathrm{66}$ mm^{3} field of view, spatial resolution of $\mathrm{2}\times \mathrm{2}\times \mathrm{6}$ mm^{3}, and a partial Fourier factor of 6∕8. Spatial resolution was sacrificed in favor of high signaltonoise ratios (SNRs). The $\mathrm{2}\times \mathrm{2}\times \mathrm{6}$ mm^{3} 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). Motioncorrected and nonmotioncorrected data were then inverted using a quick 12bootstrap 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_{2}, 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_{2}/s^{−1}) < 1.5, −10 < log(D_{}/m^{2}s^{−1}) < −8.5, −10 < log(D_{⊥}/m^{2}s^{−1}) < −8.5, 0 < cosθ < 1, 0 < φ < 2π) space. A discrete P(R_{2},D) distribution is then estimated by solving a discretized version of Eq. (2) via a standard nonnegative 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_{2}, 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_{2}, D_{}, D_{⊥}, θ, φ) configuration, and finally fitting the merged set is repeated a total of 20 times in order to find a P(R_{2}, D_{}, D_{⊥}, θ, φ) distribution yielding a low residual sum of squares. Following 20 rounds, the resulting (R_{2}, 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_{2}, 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_{2}, D_{}, D_{⊥}, θ, φ) configuration corresponding to the lowest sum of squared residuals. A final plausible P(R_{2}, D_{}, D_{⊥}, θ, φ) solution is subsequently estimated at the end of the mutation cycle by selecting the 10 (R_{2}, D_{}, D_{⊥}, θ, φ) points with the highest weights and performing a final NNLS fit.
The procedure described above is performed voxelwise, resulting in an array of spatially resolved P(R_{2}, 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_{2}, D_{}, D_{⊥}, θ, φ) set is fitted to the signal, and the resulting P(R_{2},D) is taken as the solution of the analyzed voxel. Finally, the P(R_{2}, D_{}, D_{⊥}, θ, φ) distribution is mapped onto the (R_{2}, 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 voxelwise ensembles of N distinct P(R_{2}, 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_{2},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/mddmri (last access: 25 February 2020) (Nilsson et al., 2018b). Inversion of the 45 min dataset took ∼72 h on a 12Core Intel Xeon E5 2.7 GHz CPU, with a 64 GB DDR3 memory.
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_{2}–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_{2} 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_{2}D nonparametric distributions are retrieved from the experimental data using the modelfree 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 $\mathrm{log}({D}_{\mathrm{}\mathrm{}}/{D}_{\perp})$ instead of D_{Δ}. The distinctive characters of the raw signal patterns in Fig. 2b result in unique voxelwise distributions that capture the gross microscopic features of the main cerebral components. Namely, CSF is characterized by high D_{iso}, low R_{2}, and D_{} $\sim {D}_{\perp}$; in contrast, GM and WM both exhibit lower D_{iso} and higher R_{2}, 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_{2},D) components covering a significant range of the R_{2}–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 discretecomponent 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 extraaxonal components modeled in previous diffusionrelaxation 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 $\mathrm{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}\phantom{\rule{0.125em}{0ex}}\mathit{>}\phantom{\rule{0.125em}{0ex}}\mathrm{1})$ and oblate (D_{Δ} < 0, ${D}_{\mathrm{}\mathrm{}}/{D}_{\perp}\phantom{\rule{0.125em}{0ex}}$< 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 nonnegligible, 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 voxelaveraged diffusion tensor that in turn gives rise to a positive bias in anisotropy (Pierpaoli and Basser, 1996; Jones and Cercignani, 2010).
3.2 Statistical measures of tissue heterogeneity
The R_{2}–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_{2}–D distributions. A multitude of statistical functionals can be computed from the same distribution, meaning that the pervoxel P(R_{2},D) ensembles generate a comprehensive set of distinct voxelwise parameters. As shown in Fig. 3, the Monte Carlo realizations of P(R_{2},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 intravoxel 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_{2}–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_{2}]〉 and 〈E[D_{iso}]〉 clearly identify CSFrich areas (low 〈E[R_{2}]〉 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 intravoxel 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_{2},D) distributions (displayed in Fig. 2c) comprise two populations at distant (R_{2},D_{iso}) coordinates, and both voxels are thus characterized by high values of Var[R_{2}] 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_{2}–D space. In WM+CSF distributions, for instance, higher values of diffusion anisotropy are correlated with higher R_{2} and lower D_{iso}, which results in positive Cov[R_{2},${{D}_{\mathrm{\Delta}}}^{\mathrm{2}}$] and negative Cov[D_{iso},${{D}_{\mathrm{\Delta}}}^{\mathrm{2}}$]. The elevated 〈Var[R_{2}]〉 and 〈Var[D_{iso}]〉, and negative 〈Cov[R_{2},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_{2},${{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_{2} and D_{iso}, WM+GM voxels result in nearly zero values of Var[R_{2}], Var[D_{iso}], Cov[D_{iso},y], and Cov[R_{2},y]. Instead, WM+GM voxels are signaled by finite values of 〈Var[${{D}_{\mathrm{\Delta}}}^{\mathrm{2}}$]〉, originated by the $\mathrm{log}({D}_{\mathrm{}\mathrm{}}/{D}_{\perp})$ spread observed in the underlying R_{2}–D distribution (see the WM+GM distribution in Fig. 3c).
3.3 Binresolved metrics of tissue heterogeneity
A more detailed picture of intravoxel 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_{2},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 perbin 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_{2} 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.
Global and binresolved averages for all the analyzed voxels of the entire 3D image matrix are compiled in Fig. 5, where pervoxel average means of R_{2}, D_{iso}, and ${{D}_{\mathrm{\Delta}}}^{\mathrm{2}}$ are plotted against their respective uncertainties, σ[E[R_{2}]], σ[E[D_{iso}]], and $\mathit{\sigma}\left[\mathrm{E}\right[{{D}_{\mathrm{\Delta}}}^{\mathrm{2}}\left]\right]$, and average signal amplitudes 〈S_{0}〉. 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_{2}–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_{0}〉, 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_{0}〉.
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_{2} rates can be consistently detected in voxels containing GM+WM mixtures (see Fig. 6a), where conventional 1D R_{2} 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_{2} differences greater than the estimated uncertainties, thus providing evidence that the differentiation between the R_{2} rates of the two bins is indeed a meaningful result.
All binresolved 〈E[R_{2}]〉 plots in Fig. 5 display a secondary cluster at high R_{2} values. Inspection of Fig. 6b reveals that the fast relaxing cluster corresponds to the nonmasked extrameningeal 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 highR_{2} components are observed to be concentrated around R_{2}=30 s^{−1} (see Fig. 6a), the upper R_{2}limit of the Monte Carlo inversion procedure. The “pileup” of fastrelaxing contributions around the maximum allowed R_{2} value is a wellknown artifact of Laplace inversions (Saab et al., 1999).
The 〈E[R_{2}]〉 map of the thick bin features three main R_{2} populations: high R_{2} in the skull region (red voxels), low R_{2} in peripheral brain regions (green voxels), and intermediate R_{2} 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 < $\mathrm{log}({D}_{\mathrm{}\mathrm{}}/{D}_{\perp})$ < 0.6, −10 < log (D_{iso}/m^{2} s^{−1}) < −8.7) subspace in three separate R_{2} regions, and defined the “low” (−0.5 < log(R_{2}/s^{−1}) < 1.2), “medium” (1.2 < log(R_{2}/s^{−1}) < 1.4), and “high” (1.4 < log(R_{2}/s^{−1}) < 2) subbins of Fig. 6c. In T_{2} 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_{2}=30 s^{−1}; the R_{2}=100 s^{−1} boundary is defined simply to render a more aesthetically pleasing plot (see Fig. 6c). The binresolved signal fraction maps were then compared with a highresolution longitudinal relaxationweighted (R_{1}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 binresolved 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_{2}–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_{2} 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 threetissue 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.
3.4 Clinical feasibility of the R_{2}–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 clinicresearch 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_{2}–D distributions and parameter maps are compiled in the Supplement. Figure 7 shows the binresolved 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_{2}–D properties of those same tissue types are accurately captured by the perbin mean parameter maps. The most obvious difference between the two datasets is the lower quality of the R_{2} metrics derived from the abbreviated data. This is evidenced by unreasonably high R_{2} rates in the ventricles (see the 〈E[R_{2}]〉 maps in Fig. 7b), and a higher difficulty in separating between the mean R_{2} rates of the thin and thick bins. Only 65 % of mixed voxels from the abbreviated dataset show a meaningful R_{2} separation, as opposed to the 75 % determined in the previous subsection. The lower resolution along the R_{2} 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_{2}. In future experiments, we plan to address this issue by enforcing a more uniform distribution of data points along the various echo times.
The proposed framework resolves intravoxel heterogeneity on a 5D space of transverse relaxation rates R_{2} and diffusion tensor parameters (D_{iso}, D_{Δ}, θ, φ). Pervoxel brain composition is broken down into a nonpredefined 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_{1}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_{2}–D properties provides complementary information to R_{1}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 fibertracking 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 intravoxel 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 componentspecific relaxation rates. As mentioned before, this information is unattainable with classical multiecho R_{2} 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_{2} differences can be observed despite the limited number and range of echo times sampled in this work; here, the separation between R_{2} components is mostly driven by the excellent resolution in the diffusion dimensions. The measurement of Dresolved transverse relaxation rates may complement previous work on tractspecific R_{1} 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 wellestablished approach of approximating the micrometerscale 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 gradientmodulation 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 monoexponential 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^{9} 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 extraaxonal compartments in brain tissue has been of longstanding interest in the MRI field (Does, 2018). Recently, Veraart et al. (2017) estimated subtle differences in R_{2} and diffusivity parameters for the intra and extraaxonal components of human brain white matter by applying a constrained twocomponent model to data acquired with a conventional relaxation–diffusion correlation protocol relying on the Stejskal–Tanner experiment. The obtained R_{2} 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 extracellular 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 btensor shape dimension. The limited resolution is consistent with the fact that Eq. (2) states an illposed 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_{2}, isotropic diffusivities D_{iso}, normalized diffusion anisotropy D_{Δ}, and diffusion tensor orientation (θ, φ). The correlations allow modelfree estimation of pervoxel relaxation–diffusion distributions P(R_{2},D) that combine the chemical sensitivity of R_{2} with the link between microstructure and the diffusion metrics. The rich information content of P(R_{2},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 ultrahigh fields (7T and above), and ultrastrong field gradients (Setsompop et al., 2013; Jones et al., 2018), can boost SNR and reduce echotimeperunitb 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, multiband 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 echoplanar 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 patternmatching algorithms may considerably boost the data inversion speed.
The software analysis tools discussed in this paper are available for download from a public GitHub repository: https://github.com/JoaoPdAMartins/mddmri (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: https://doi.org/10.5194/mr1272020supplement.
DKJ, CFW, and DT conceived the project. JPdAM, CMWT, and FS designed the acquisition protocol. FS and CFW 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.
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 (AM130090, ITM170267) and the Swedish Research Council (20143910, 201803697). Chantal M. W. Tax is supported by a Rubicon grant (680501527) from the Netherlands Organisation for Scientific Research (NWO). Filip Szczepankiewicz and CarlFredrik 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).
This research has been supported by the Stiftelsen för Strategisk Forskning (grant nos. AM130090 and ITM170267), the Vetenskapsrådet (grant nos. 20143910 and 201803697), the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (grant no. 680501527), the National Institutes of Health (grant no. P41EB015902), and the Wellcome Trust (grant nos. 096646/Z/11/Z and 104943/Z/14/Z).
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, https://doi.org/10.1016/j.neuroimage.2017.12.096, 2018.
Bak, M. and Nielsen, N. C.: REPULSION, a novel approach to efficient powder averaging in solidstate NMR, J. Magn. Reson., 125, 132–139, https://doi.org/10.1006/jmre.1996.1087, 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, https://doi.org/10.1002/mrm.25897, 2016.
Basser, P. J. and Pierpaoli, C.: Microstructural and physiological features of tissues elucidated by quantitativediffusiontensor MRI, J. Magn. Reson. Ser. B, 111, 209–219, https://doi.org/10.1016/j.jmr.2011.09.022, 1996.
Bazin, P.L., Alkemade, A., van der Zwaag, W., Caan, M., Mulder, M., and Forstmann, B. U.: Denoising HighField MultiDimensional MRI With Local Complex PCA, Front. Neurosci.Switz., 13, 1066, https://doi.org/10.3389/fnins.2019.01066, 2019.
Benjamini, D. and Basser, P. J.: Magnetic resonance microdynamic imaging reveals distinct tissue microenvironments, Neuroimage, 163, 183–196, https://doi.org/10.1016/j.neuroimage.2017.09.033, 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, https://doi.org/10.1016/j.mri.2017.08.001, 2017.
de Almeida Martins, J. P. and Topgaard, D.: TwoDimensional Correlation of Isotropic and Directional Diffusion Using NMR, Phys. Rev. Lett., 116, 087601, https://doi.org/10.1103/PhysRevLett.116.087601, 2016.
de Almeida Martins, J. P. and Topgaard, D.: Multidimensional correlation of nuclear relaxation rates and diffusion tensors for modelfree investigations of heterogeneous anisotropic porous materials, Sci. Rep., 8, 2488, https://doi.org/10.1038/s41598018198269, 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 diffusionweighted acquisitions, Magn. Reson. Med., 75, 372–380, https://doi.org/10.1002/mrm.25644, 2016.
Does, M. D.: Inferring brain tissue composition and microstructure via MR relaxometry, Neuroimage, 182, 136–148, https://doi.org/10.1016/j.neuroimage.2017.12.087, 2018.
English, A. E., Whittal, K. P., Joy, M. L. G., and Henkelman, R. M.: Quantitative twodimensional time correlation relaxometry, Magn. Reson. Med., 22, 425–434, https://doi.org/10.1002/mrm.1910220250, 1991.
Eriksson, S., Lasic, S., and Topgaard, D.: Isotropic diffusion weighting in PGSE NMR by magicangle spinning of the qvector, J. Magn. Reson., 226, 13–18, https://doi.org/10.1016/j.jmr.2012.10.015, 2013.
Eriksson, S., Lasic, S., Nilsson, M., Westin, C. F., and Topgaard, D.: NMR diffusionencoding with axial symmetry and variable anisotropy: Distinguishing between prolate and oblate microscopic diffusion tensors with unknown orientation distribution, J. Chem. Phys., 142, 104201, https://doi.org/10.1063/1.4913502, 2015.
Frydman, L., Chingas, G. C., Lee, Y. K., Grandinetti, P. J., Eastman, M. A., Barrall, G. A., and Pines, A.: Variableangle correlation spectroscopy in solidstate nuclear magnetic resonance, J. Chem. Phys., 97, 4800–4808, https://doi.org/10.1063/1.463860, 1992.
Galvosas, P. and Callaghan, P. T.: Multidimensional inverse Laplace spectroscopy in the NMR of porous media, C. R. Physique, 11, 172–180, https://doi.org/10.1016/j.crhy.2010.06.014, 2010.
Gan, Z.: Highresolution chemical shift and chemical shift anisotropy correlation in solids using slow magic angle spinning, J. Am. Chem. Soc., 114, 8307–8309, https://doi.org/10.1021/ja00047a062, 1992.
Haldar, J. P., Liu, Y., Liao, C., Fan, Q., and Setsompop, K.: Fast submillimeter diffusion MRI using gSliderSMS and SNRenhancing joint reconstruction, Magn. Reson. Med., in press, https://doi.org/10.1002/mrm.28172, 2020.
Halle, B.: Molecular theory of fielddependent proton spinlattice relaxation in tissue, Magn. Reson. Med., 56, 60–72, https://doi.org/10.1002/mrm.20919, 2006.
Istratov, A. A. and Vyvenko, O. F.: Exponential analysis in physical phenomena, Rev. Sci. Instrum., 70, 1233–1257, https://doi.org/10.1063/1.1149581, 1999.
Jeurissen, B., Tournier, J. D., Dhollander, T., Connelly, A., and Sijbers, J.: Multitissue constrained spherical deconvolution for improved analysis of multishell diffusion MRI data, Neuroimage, 103, 411–426, https://doi.org/10.1016/j.neuroimage.2014.07.061, 2014.
Jian, B., Vemuri, B. C., Özarslan, E., Carney, P. R., and Mareci, T. H.: A novel tensor distribution model for the diffusionweighted MR signal, Neuroimage, 37, 164–176, https://doi.org/10.1016/j.neuroimage.2007.03.074, 2007.
Jones, D. K.: Diffusion MRI, Oxford University Press, 2010.
Jones, D. K. and Cercignani, M.: Twentyfive pitfalls in the analysis of diffusion MRI data, NMR Biomed., 23, 803–820, https://doi.org/10.1002/nbm.1543, 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, https://doi.org/10.1002/(SICI)15222594(199909)42:3<515::AIDMRM14>3.0.CO;2Q, 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 “superscanner”: 10 key advantages of ultrastrong gradients for diffusion MRI, Neuroimage, 182, 8–38, https://doi.org/10.1016/j.neuroimage.2018.05.047, 2018.
Kärger, J.: Zur Bestimmung der Diffusion in einem Zweibereichsystem mit Hilfe von gepulsten Feldgradienten, Ann. Phys., 479, 1–4, https://doi.org/10.1002/andp.19694790102, 1969.
Klein, S., Staring, M., Murphy, K., Viergever, M. A., and Pluim, J. P.: Elastix: a toolbox for intensitybased 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, https://doi.org/10.1016/j.jpsychires.2005.05.005, 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, https://doi.org/10.1002/hbm.24542, 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 magicangle spinning of the qvector, Front. Phys., 2, 11, https://doi.org/10.3389/fphy.2014.00011, 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, https://doi.org/10.1007/s0041501786096, 2017.
Lawson, C. L. and Hanson, R. J.: Solving least squares problems, PrenticeHall, Englewood Cliffs, NJ, 1974.
Le Bihan, D.: Molecular diffusion, tissue microdynamics and microstructure, NMR Biomed., 8, 375–386, https://doi.org/10.1002/nbm.1940080711, 1995.
Lerch, J. P., van der Kouwe, A. J., Raznahan, A., Paus, T., JohansenBerg, H., Miller, K. L., Smith, S. M., Fischl, B., and Sotiropoulos, S. N.: Studying neuroanatomy using MRI, Nat. Neurosci., 20, 314–326, https://doi.org/10.1038/nn.4501, 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, https://doi.org/10.1038/s41598019452357, 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, https://doi.org/10.1038/nature11971, 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, https://doi.org/10.1002/mrm.1910310614, 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, https://doi.org/10.1016/j.pnmrs.2011.07.002, 2012.
Mitra, P. P.: Multiple wavevector extension of the NMR pulsedfieldgradient spinecho diffusion measurement, Phys. Rev. B, 51, 15074–15078, https://doi.org/10.1103/PhysRevB.51.15074, 1995.
Nilsson, M., Szczepankiewicz, F., van Westen, D., and Hansson, O.: ExtrapolationBased References Improve Motion and EddyCurrent Correction of High BValue DWI Data: Application in Parkinson's Disease Dementia, PLoS One, 10, e0141825, https://doi.org/10.1371/journal.pone.0141825, 2015.
Nilsson, M., Englund, E., Szczepankiewicz, F., van Westen, D., and Sundgren, P. C.: Imaging brain tumour microstructure, Neuroimage, 182, 232–250, https://doi.org/10.1016/j.neuroimage.2018.04.075, 2018a.
Nilsson, M., Szczepankiewicz, F., Lampinen, B., Ahlgren, A., De Almeida Martins, J. P., Lasic, S., Westin, C.F., and Topgaard, D.: An opensource 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, https://doi.org/10.1002/mrm.27101, 2018.
Padhani, A. R., Liu, G., MuKoh, D., Chenevert, T. L., Thoeny, H. C., Takahara, T., DzikJurasz, A., Ross, B. D., Van Cauteren, M., Collins, D., Hammoud, D. A., Rustin, G. J. S., Taouli, B., and Choyke, P. L.: DiffusionWeighted Magnetic Resonance Imaging as a Cancer Biomarker: Consensus and Recommendations, Neoplasia, 11, 102–125, https://doi.org/10.1593/neo.81328, 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, https://doi.org/10.1002/mrm.22055, 2009.
Pierpaoli, C. and Basser, P. J.: Toward a quantitative assessment of diffusion anisotropy, Magn. Res. Med., 36, 893–906, https://doi.org/10.1002/mrm.1910360612, 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, https://doi.org/10.1148/radiology.201.3.8939209, 1996.
Prange, M. and Song, Y. Q.: Quantifying uncertainty in NMR T2 spectra using Monte Carlo inversion, J. Magn. Reson., 196, 54–60, https://doi.org/10.1016/j.jmr.2008.10.008, 2009.
Saab, G., Thompson, R. T., and Marsh, G. D.: Multicomponent T2 relaxation of in vivo skeletal muscle, Magn. Reson. Med., 42, 150–157, https://doi.org/10.1002/(SICI)15222594(199907)42:1<150::AIDMRM20>3.0.CO;25, 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 diffusioncompartment imaging (DIAMOND), Magn. Reson. Med., 76, 963–977, 2016.
SchmidtRohr, K. and Spiess, H. W.: Multidimensional solidstate NMR and polymers, Academic Press, 1994.
Setsompop, K., Kimmlingen, R., Eberlein, E., Witzel, T., CohenAdad, 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, https://doi.org/10.1016/j.neuroimage.2013.05.078, 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, https://doi.org/10.1016/j.jmr.2015.10.012, 2015.
Song, Y. Q.: Magnetic resonance of porous media (MRPM): a perspective, J. Magn. Reson., 229, 12–24, https://doi.org/10.1016/j.jmr.2012.11.010, 2013.
Stejskal, E. O. and Tanner, J. E.: Spin diffusion measurements: Spin echoes in the presence of a timedependent field gradient, J. Chem. Phys., 42, 288–292, https://doi.org/10.1063/1.1695690, 1965.
Szczepankiewicz, F., Westin, C. F., and Nilsson, M.: Maxwellcompensated design of asymmetric gradient waveforms for tensorvalued diffusion encoding, Magn. Reson. Med., 82, 1424–1437, https://doi.org/10.1002/mrm.27828, 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, https://doi.org/10.1016/j.jmr.2016.12.007, 2017.
Topgaard, D.: Diffusion tensor distribution imaging, NMR Biomed., 32, e4066, https://doi.org/10.1002/nbm.4066, 2019.
Topgaard, D. and Söderman, O.: Selfdiffusion in twoand threedimensional powders of anisotropic domains: An NMR study of the diffusion of water in cellulose and starch, J. Phys. Chem. B, 106, 11887–11892, https://doi.org/10.1021/jp020130p, 2002.
Turner, R., Le Bihan, D., and Scott Chesnicks, A.: Echoplanar imaging of diffusion and perfusion, Magn. Reson. Med., 19, 247–253, https://doi.org/10.1002/mrm.1910190210, 1991.
Veraart, J., Novikov, D. S., Christiaens, D., AdesAron, B., Sijbers, J., and Fieremans, E.: Denoising of diffusion MRI using random matrix theory, Neuroimage, 142, 394–406, https://doi.org/10.1016/j.neuroimage.2016.08.016, 2016.
Veraart, J., Novikov, D. S., and Fieremans, E.: TE dependent Diffusion Imaging (TEdDI) distinguishes between compartmental T2 relaxation times, Neuroimage, 182, 360–369, https://doi.org/10.1016/j.neuroimage.2017.09.030, 2017.
Wang, H., Zheng, R., Dai, F., Wang, Q., and Wang, C.: Highfield mr diffusionweighted image denoising using a joint denoising convolutional neural network, J. Magn. Reson. Imaging, 50, 1937–1947, https://doi.org/10.1002/jmri.26761, 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, https://doi.org/10.1093/brain/awr307, 2011.
Whittall, K. P. and MacKay, A. L.: Quantitative interpretation of NMR relaxation data, J. Magn. Reson., 84, 134–152, https://doi.org/10.1016/00222364(89)900115, 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, https://doi.org/10.1002/mrm.1910370107, 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.: Singleshot spiral imaging enabled by an expanded encoding model: Demonstration in diffusion MRI, Magn. Reson. Med., 77, 83–91, https://doi.org/10.1002/mrm.26493, 2017.
Woessner, D. E.: N.M.R. spinecho selfdiffusion measurements on fluids undergoing restricted diffusion, J. Phys. Chem., 67, 1365–1367, https://doi.org/10.1021/j100800a509, 1963.
Zatorre, R. J., Fields, R. D., and JohansenBerg, H.: Plasticity in gray and white: neuroimaging changes in brain structure during learning, Nat. Neurosci., 15, 528–536, https://doi.org/10.1038/nn.3045, 2012.
Zhang, H., Schneider, T., WheelerKingshott, C. A., and Alexander, D. C.: NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain, Neuroimage, 61, 1000–1016, https://doi.org/10.1016/j.neuroimage.2012.03.072, 2012.
Zhang, Y. and Blumich, B.: Spatially resolved DT2 correlation NMR of porous media, J. Magn. Reson., 242, 41–48, https://doi.org/10.1016/j.jmr.2014.01.017, 2014.