Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation

Abstract

We review, systematize and discuss models of diffusion in neuronal tissue, by putting them into an overarching physical context of coarse-graining over an increasing diffusion length scale. From this perspective, we view research on quantifying brain microstructure as occurring along three major avenues. The first avenue focusses on transient, or time-dependent, effects in diffusion. These effects signify the gradual coarse-graining of tissue structure, which occurs qualitatively differently in different brain tissue compartments. We show that transient effects contain information about the relevant length scales for neuronal tissue, such as the packing correlation length for neuronal fibers, as well as the degree of structural disorder along the neurites. The second avenue corresponds to the long-time limit, when the observed signal can be approximated as a sum of multiple nonexchanging anisotropic Gaussian components. Here, the challenge lies in parameter estimation and in resolving its hidden degeneracies. The third avenue employs multiple diffusion encoding techniques, able to access information not contained in the conventional diffusion propagator. We conclude with our outlook on future directions that could open exciting possibilities for designing quantitative markers of tissue physiology and pathology, based on methods of studying mesoscopic transport in disordered systems.

Keywords: MRI; brain; coarse-graining; diffusion; mesoscopic transport; microstructure.

FIG. 1.

The mesoscopic scale in brain dMRI, as an intermediate scale between the elementary (molecular) and the macroscopic (resolution).

FIG. 2.. Diffusion as coarse-graining.

An example of a medium where the mesoscopic structure is created by randomly placing black disks of two different radii, r_small = 1 μm and r_large = 20 μm, top left panel. To obtain snapshots of the medium as effectively seen by the diffusing molecules at different time scales, we used a Gaussian filter with width L/2, where $L(t)=\sqrt{2Dt}$, and ignored the time dependence of D(t) in the definition of diffusion length, using a typical value D = 1 μm²/ms for the illustration purposes (cf. Sec. 2.4 below).

FIG. 3.

The parameter space of dMRI is at least twodimensional: By increasing q one accesses the progressively higher-order diffusion cumulants, Sec. 1.8, whereas the dependence along the t-axis reflects their epvolution over an increasing diffusion length scale $L~\sqrt{Dt}$, Eq. (1.12). The b-value alone does not uniquely describe the measurement, unless diffusion in all tissue compartments is Gaussian; contour lines of b = q²t are schematically drawn in beige. Large-q limits: Top-left is high-resolution limit L(t) ≪ l_c, ql_c ≫ 1, Sec. 1.5(i); middle is the $L(t)\gtrsim l_c$, $q{l}_c\gtrsim 1$ limit of probing the pore correlation function, Sec. 1.6. The hierarchy of dMRI models (pictures), cf. Fig. 4, as well as the cumulant representation with different number of terms, cf. Fig. 5, are superimposed. The decrease of the signal from axonal bundles parallel to the increasing gradient is shown by their darkening (top right). In Section 2 we move along the t-axis at low q, and in Section 3 we move along the q-axis at asymptotically long t. Section 4 is devoted to effects beyond this diagram, contained in voxelaveraged products of propagators at different t and q.

FIG. 4.. Models are pictures…

Here they are drawn with coarse-graining occurring, roughly, from left to right. References: Mitra 1993 [32], universal short-t limit; Novikov 2014 [33], universal long-t behavior; Burcaw 2015 [23] and Fieremans 2016 [34], long-t behavior transverse and along WM fibers; Tanner 1978 [49], Powles 1992 [50], Sukstanskii 2004 [35], periodic 1-dimensional lattice; Novikov 2011 [51], random permeable barriers in any dimension, and its application to myofibers (Sigmund 2014 [52] and Fieremans 2016 [53]); Callaghan 1979 [54], first model of diffusion inside random narrow cylinders; Yablonskiy 2002 [55], diffusion in finite-diameter cylinders modeling lung alveoli; Stanisz 1997 [56], first model for WM fiber tracts made of ellipsoids; Assaf 2004 [57], CHARMED; Assaf 2008 [58], AxCaliber; Alexander 2010 [59], ActiveAx; Kroenke 2004 [36], NAA diffusion inside neurites. The widely adopted t → ∞ picture of narrow “sticks” for the neurites, embeded in the extra-neurite space (the Standard Model): Jespersen 2007 [37], Jespersen 2010 [38], Fieremans 2010 [39], Fieremans 2011 (WMTI) [40], Sotiropoulos 2012 (Ball and rackets) [41], Zhang 2012 (NODDI) [42], Reisert 2014 (MesoFT) [43], Jelescu 2016 (NODDIDA) [44], Reisert 2017 [45], Veraart 2017 (TEdDI) [48], Novikov 2018 (LEMONADE [46], RotInv [47]).

FIG. 5.. …while representations are formulas.

References: Le Bihan 1991 [60], first biexponential representation of dMRI signal from brain; Basser 1994 [31], diffusion tensor imaging (DTI); Jensen 2005 [61], diffusion kurtosis imaging (DKI); Kiselev 2011 [18], cumulant expansion; Novikov 2008 and 2010 [27, 62], effective medium theory (transverse relaxation and diffusion, correspondingly); Özarslan 2013 [63], expansion in harmonic oscillator basis; Yablonskiy 2003 [64], inverse Laplace transform (multi-exponential representation). Anisotropic multi-exponential representations: Wang 2011 [65], diffusion basis spectrum imaging (DBSI); White 2013 [66], restriction spectrum imaging (RSI); Scherrer 2016 [67], distribution of anisotropic microstructural environments in diffusion-compartment imaging (DIAMOND).

FIG. 6.

General relations between the basic diffusion metrics:$\mathcal{D}(\omega )$, $\mathcal{D}(t)$, D_inst(t) and D(t), and the signal attenuation up to $\mathcal{O}(q^2)$.

FIG. 7.

OGSE measurements in cortical GM: circles are data from average of 5 rats [123] and squares from 6 neonatal mice at 24 hours after unilateral hypoxic ischemic injury [126]. Red: normal rat brain and contralateral side of mouse brain. Blue: globally ischemic rat and ipsilateral side of hypoxia-ischemia injured mouse brain. PGSE data not shown. Dashed lines are fits from Fig. 4 of ref. [33], dotted lines are ω^1/2 fits (shown as guide to the eye; power-law exponent fit for mouse data was not robust due to narrow frequency range). Note that while the absolute (ω) values differ between rat and mouse, the general features are similar: data is well described with ω^1/2 behavior for normal and ischemic GM (except, possibly, the ischemic mouse, where major structural changes may have occurred in 24h); and the coefficient in front of ω (the slope) increases in ischemia, consistent with short-range structural disorder increase along the neurites (e.g., due to beading).

FIG. 8.. The Standard Model of diffusion in neuronal tissue, Eq. (3.4).

In the t → ∞ regime (iii), elementary fiber segments (fascicles), consisting of intra- and extra-neurite compartments, are described by at least 4 independent parameters: f, D_a, $D_e^{||}$ and $D_e^{\perp }$. CSF can be further added as the third compartment, cf. Eq. (3.5). Within a macroscopic voxel, such segments contribute to the directional dMRI signal according to their ODF $\mathcal{P}(\widehat{\mathbf{n}})$.

FIG. 9.

Comparison of NODDI (Sec. 3.4.2) and WMTI (Sec. 3.4.1) parameter evolution with age in human corpus callosum splenium [232]. A qualitatively (but not quantitatively) similar trend of continued increase in the intra-axonal water fraction f_intra ≡ f was observed for both models, consistent with on-going myelination. WMTI displays a trend of increased fiber alignment (expressed by the orientation dispersion $\langle {\cos }^2\psi \rangle$, derived from the intra-axonal diffusion tensor ), which could be a manifestation of continued pruning in the first year of life, while NODDI does not. The CSF fraction is set, f_iso ≡ f_CSF = 0, in WMTI.

FIG. 10.

(a) Example of a DDE sequence within the framework of a double spin echo. (b) The resulting gradient waveform obtained by multiplying each gradient by ${(-1)}^{n_{\pi }}$, where n_π is the number of π pulses following the given gradient. In the text, we assume narrowpulse approximation, such that δ_i → 0, with the Larmor frequency gradients g_i sufficiently large to yield finite q_i = g_iδ_i (no summation over i).

FIG. 11.

Two-dimensional temporal integration involved in the second-order cumulant, Eq. (4.7) leading to Eq. (4.6) for the DDE measurement. Labels q₁ and q₂ indicate the time interval in which q(t) equals to q₁ and q₂, respectively; for simplicity, the vector indices are not shown. The green shaded area along the diagonal symbolizes ${\mathcal{D}}_{ij}\left(\left|t_1-t_2\right|\right)$, Eq. (4.8), where it significantly deviates from 0, with the width of this region set by the correlation time t_c. The nontrivial cross-term q_1iq_2j in Eq. (4.6) arises from the off-diagonal quadrants. As this contribution is weighted with the velocity autocorrelation function, it tends to zero when the mixing time, τ (indicated by the thin lines along each dimension) becomes larger than the correlation time, τ ≫ t_c. In particular, no non-trivial cross-term is present for Gaussian diffusion, for which t_c → 0.

FIG. 12.

Examples of model systems considered in the text. In (a), a system of identical spherical pores is shown, whereas, in (b), the pores have a distribution of sizes. In (c), an approximately isotropic distribution of ellipsoidal pores is sketched and, in (d), the pores are coherently oriented. Systems (a)–(c) are macroscopically isotropic, system (d) is not. Systems (c) and (d) are microscopically anisotropic. Ensemble heterogeneity is only seen in systems (b) (size) and (c) (orientation). Here, spins contributing to the signal are assumed to only reside within the pores.

FIG. 13.

Analytical structure of a causal (retarded) response function on the complex plane of ω. When calculating the inverse Fourier transform such as Eq. (A5), the original integration contour over the real axis can be closed in the infinite semicircle with Im ω > 0 (light blue dashed line) when t < 0, according to the Jordan’s lemma. Causality then requires that no singularities are present in the upper half of the complex plane, in which case the integration contour can be shrunk to a point. For t > 0, the contour can be closed where Im ω < 0 (light red dashed line). This contour can be shrunk to encircle the singularities of the transformed function (red solid lines).