the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### João P. Almeida Martins

### Chantal M. W. Tax

### Filip Szczepankiewicz

### Derek K. Jones

### Carl-Fredrik 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 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.

- Article
(6340 KB) -
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 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 ^{1}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.

## 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
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*_{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
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*_{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 *b*-tensor parameters facilitates a
comprehensive mapping of tissue-specific relaxation and diffusion
properties.

## 2.2 MRI measurements

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

Experimental data were acquired using the prototype spin-echo sequence
(Lasič et al., 2014) and gradient waveforms shown in
Fig. 1. The depicted waveforms give four
distinct *b*-tensor anisotropies (${b}_{\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}=\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 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).

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 echo-planar
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
signal-to-noise 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). 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*_{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 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*_{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 voxel-wise, 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 voxel-wise 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/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.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 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 $\mathrm{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*_{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
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 $\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 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).

## 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 per-voxel *P*(*R*_{2},**D**) ensembles generate a comprehensive
set of distinct voxel-wise 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 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*_{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
CSF-rich 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 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*_{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 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*_{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 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*_{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 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*_{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 bin-resolved 〈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 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*_{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 “pile-up” of
fast-relaxing contributions around the maximum allowed *R*_{2} value is a
well-known 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) sub-bins 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 bin-resolved signal
fraction maps were then compared with a high-resolution longitudinal
relaxation-weighted (*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
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*_{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 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.

## 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 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*_{2}–**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*_{2}–**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*_{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 intra-voxel heterogeneity on a 5D space of
transverse relaxation rates *R*_{2} 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*_{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 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*_{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 **D**-resolved transverse
relaxation rates may complement previous work on tract-specific *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 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^{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 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*_{2} 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*_{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 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*_{2}, 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*_{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 ultra-high
fields (7T and above), and ultra-strong field gradients (Setsompop et al.,
2013; Jones et al., 2018), can boost SNR and reduce
echo-time-per-unit-*b* value, respectively. For example, as noted in
Jones et al. (2018), for *b*_{Δ}=0
encoding, the shorter *τ*_{E} afforded by stronger gradients such as
those available on a Connectom scanner (300 mT m^{−1}) results in an improvement
in SNR of approximately 50 % compared to that achievable on the system
used in this study (80 mT m^{−1} gradients). From the acquisition perspective,
multi-band acquisition schemes (Barth et al., 2016) can
speed up overall acquisition times and facilitate a wide brain coverage with
smaller voxel sizes. Moreover, replacing the rectilinear echo-planar readout
(Turner et al., 1991) with a spiral readout (Wilm
et al., 2017) can help to further reduce the echo time, boosting SNR which
could be traded for higher spatial resolution. From the analysis side, as
noted in the Methods section, no denoising approaches were applied here. Recent
advances in denoising and/or joint reconstruction (Veraart et al.,
2016; Bazin et al., 2019; Wang et al., 2019; Haldar et al., 2020) could further
enhance the SNR, allowing resolution to be pushed higher. Finally, the
presented framework can be merged with MRI fingerprinting methodology
(Ma et al., 2013), whose pattern-matching algorithms
may considerably boost the data inversion speed.

The software analysis tools discussed in this paper are available for download from a public GitHub repository: https://github.com/JoaoPdAMartins/md-dmri (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/mr-1-27-2020-supplement.

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

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

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

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

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 solid-state 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 quantitative-diffusion-tensor 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 High-Field Multi-Dimensional 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.: Two-Dimensional 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 model-free investigations of heterogeneous anisotropic porous materials, Sci. Rep., 8, 2488, https://doi.org/10.1038/s41598-018-19826-9, 2018.

De Santis, S., Barazany, D., Jones, D. K., and Assaf, Y.: Resolving relaxometry and diffusion properties within the same voxel in the presence of crossing fibres by combining inversion recovery and diffusion-weighted acquisitions, Magn. Reson. Med., 75, 372–380, 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 two-dimensional 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 magic-angle spinning of the q-vector, 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 diffusion-encoding with axial symmetry and variable anisotropy: Distinguishing between prolate and oblate microscopic diffusion tensors with unknown orientation distribution, J. Chem. Phys., 142, 104201, 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.: Variable-angle correlation spectroscopy in solid-state nuclear magnetic resonance, J. Chem. Phys., 97, 4800–4808, https://doi.org/10.1063/1.463860, 1992.

Galvosas, P. and Callaghan, P. T.: Multi-dimensional 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.: High-resolution 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 gSlider-SMS and SNR-enhancing joint reconstruction, Magn. Reson. Med., in press, https://doi.org/10.1002/mrm.28172, 2020.

Halle, B.: Molecular theory of field-dependent proton spin-lattice 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.: Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell 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 diffusion-weighted 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.: Twenty-five 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)1522-2594(199909)42:3<515::AID-MRM14>3.0.CO;2-Q, 1999.

Jones, D. K., Alexander, D. C., Bowtell, R., Cercignani, M., Dell'Acqua, F., McHugh, D. J., Miller, K. L., Palombo, M., Parker, G. J. M., Rudrapatna, U. S., and Tax, C. M. W.: Microstructural imaging of the human brain with a “super-scanner”: 10 key advantages of ultra-strong gradients for diffusion MRI, Neuroimage, 182, 8–38, 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 intensity-based medical image registration, IEE Trans. Med. Imaging, 29, 196–205, 2009.

Kubicki, M., McCarley, R., Westin, C.-F., Park, H.-J., Maier, S., Kikinis, R., Jolesz, F. A., and Shenton, M. E.: A review of diffusion tensor imaging studies in schizophrenia, J. Psychiatr. Res., 41, 15–30, 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 magic-angle spinning of the q-vector, 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/s00415-017-8609-6, 2017.

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

Le Bihan, D.: Molecular diffusion, tissue microdynamics and microstructure, NMR Biomed., 8, 375–386, https://doi.org/10.1002/nbm.1940080711, 1995.

Lerch, J. P., van der Kouwe, A. J., Raznahan, A., Paus, T., Johansen-Berg, H., Miller, K. L., Smith, S. M., Fischl, B., and Sotiropoulos, S. N.: Studying neuroanatomy using MRI, Nat. Neurosci., 20, 314–326, 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/s41598-019-45235-7, 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 wave-vector extension of the NMR pulsed-field-gradient spin-echo 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.: Extrapolation-Based References Improve Motion and Eddy-Current Correction of High B-Value DWI Data: Application in Parkinson's Disease Dementia, PLoS One, 10, e0141825, 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 open-source framework for analysis of multidimensional diffusion MRI data implemented in MATLAB, in: Proceedings of the 26th Annual Meeting of ISMRM, Paris, France, 16–21 June 2018, 5355, 2018b.

Novikov, D. S., Kiselev, V. G., and Jespersen, S. N.: On modeling, Magn. Reson. Med., 79, 3172–3193, https://doi.org/10.1002/mrm.27101, 2018.

Padhani, A. R., Liu, G., Mu-Koh, D., Chenevert, T. L., Thoeny, H. C., Takahara, T., Dzik-Jurasz, A., Ross, B. D., Van Cauteren, M., Collins, D., Hammoud, D. A., Rustin, G. J. S., Taouli, B., and Choyke, P. L.: Diffusion-Weighted Magnetic Resonance Imaging as a Cancer Biomarker: Consensus and Recommendations, Neoplasia, 11, 102–125, 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)1522-2594(199907)42:1<150::AID-MRM20>3.0.CO;2-5, 1999.

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

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

Setsompop, K., Kimmlingen, R., Eberlein, E., Witzel, T., Cohen-Adad, J., McNab, J. A., Keil, B., Tisdall, M. D., Hoecht, P., Dietz, P., Cauley, S. F., Tountcheva, V., Matschl, V., Lenz, V. H., Heberlein, K., Potthast, A., Thein, H., Van Horn, J., Toga, A., Schmitt, F., Lehne, D., Rosen, B. R., Wedeen, V., and Wald, L. L.: Pushing the limits of in vivo diffusion MRI for the Human Connectome Project, Neuroimage, 80, 220–233, 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 time-dependent field gradient, J. Chem. Phys., 42, 288–292, https://doi.org/10.1063/1.1695690, 1965.

Szczepankiewicz, F., Westin, C. F., and Nilsson, M.: Maxwell-compensated design of asymmetric gradient waveforms for tensor-valued diffusion encoding, Magn. Reson. Med., 82, 1424–1437, 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.: Self-diffusion in two-and three-dimensional powders of anisotropic domains: An NMR study of the diffusion of water in cellulose and starch, J. Phys. Chem. B, 106, 11887–11892, https://doi.org/10.1021/jp020130p, 2002.

Turner, R., Le Bihan, D., and Scott Chesnicks, A.: Echo-planar 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., Ades-Aron, 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.: High-field mr diffusion-weighted 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/0022-2364(89)90011-5, 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.: Single-shot 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. spin-echo self-diffusion 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 Johansen-Berg, 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., Wheeler-Kingshott, C. A., and Alexander, D. C.: NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain, Neuroimage, 61, 1000–1016, https://doi.org/10.1016/j.neuroimage.2012.03.072, 2012.

Zhang, Y. and Blumich, B.: Spatially resolved D-T2 correlation NMR of porous media, J. Magn. Reson., 242, 41–48, https://doi.org/10.1016/j.jmr.2014.01.017, 2014.