Multidimensional Diffusion MRI Methods With Confined Subdomains

Diffusion Magnetic Resonance Imaging (dMRI) is an imaging technique with exquisite sensitivity to the microstructural properties of heterogeneous media. The conventionally adopted acquisition schemes involving single pulsed field gradients encode the random motion of water molecules into the NMR signal, however typically conflating the effects of different sources contributing to the water motion. Time-varying magnetic field gradients have recently been considered for disentangling such effects during the data encoding phase, opening to the possibility of adding specificity to the recovered information about the medium’s microstructure. Such data is typically represented via a diffusion tensor distribution (DTD) model, thus assuming the existence of several non-exchanging compartments in each of which diffusion is unrestricted. In this work, we consider a model that takes confinement into account and possesses a diffusion time-dependence closer to that of restricted diffusion, to replace the free diffusion assumption in multidimensional diffusion MRI methods. We first demonstrate how the confinement tensor model captures the relevant signal modulations impressed by water diffusing in both free and closed spaces, for data simulated with a clinically feasible protocol involving time-varying magnetic field gradients. Then, we provide the basis for incorporating this model into two multidimensional dMRI methods, and attempt to recover a confinement tensor distribution (CTD) on a human brain dataset.

1 Introduction

Diffusion Magnetic Resonance Imaging (dMRI) is a method used for investigating the microstructural organization of various heterogeneous media. This is achieved by sensitizing the MR signal to the random motion of water molecules inside the scanned substrate. To interpret and extract relevant information from the water motion, several models and signal representations have been developed.

At spatial resolutions achievable with current MR scanners, the scanned sample comprises several compartments within, outside, and possibly in between which diffusion is taking place. A general strategy for capturing this complexity without assuming any specific combination of compartments (see for example [1, 2] for reviews of multi-compartment models for brain white matter), considers modeling the medium as a collection of isolated pores, each represented by a diffusion tensor [3, 4]. This approach leads to a diffusion tensor distribution (DTD), which could also be represented parametrically via normal [5] and Wishart [4] distributions as well as other related distributions [6–8]. Advances in diffusion encoding schemes [9–13] provided ways of disentangling confounding signal contributions, thus possibly enabling the extraction of relevant information about the medium’s structure and composition via such modeling [14]. However, it is rather paradoxical to have free diffusion within isolated compartments, as the cellular membranes have a strong effect on diffusion, making them the primary determinant of diffusion anisotropy [15]. If this picture involving multiple isolated compartments is to be employed, it would be natural to represent the individual subdomains by accounting for confined diffusion within them [16].

A viable alternative to the diffusion tensor representation of individual subdomains utilizes confinement tensors [17] instead. In this case, the molecules are envisioned to be diffusing under the influence of an Hookean restoring force, i.e., according to the Ornstein-Uhlenbeck process [18]. Just like in restricted diffusion, the particle trajectories have limited extent, which has made the Ornstein-Uhlenbeck process a simple toy problem in earlier theoretical works on characterizing the influence of restricted diffusion on the NMR signal [19–21].

Following a series of developments [22, 23], the confinement tensor model has been noted to provide an alternative representation of diffusion anisotropy, very-well suited for studying heterogeneous media [17, 24, 25]. Furthermore, for NMR experiments involving long diffusion encoding pulses, the harmonic confinement becomes the effective model of restricted diffusion, giving rise to an approximately linear dependence of the effective stochastic force on the center-of-mass position of the particles during the application of the gradient pulses [26].

Similarly to the better-known diffusion tensor model, the model proposed by Yolcu et al. [17] captures the pore’s geometry/anisotropy with a tensorial object, which can be visualized as an ellipsoid. However, the confinement tensor model offers an extra parameter to encode diffusivity. This parameter can either be a scalar bulk diffusivity, or another tensorial quantity. In either case, this represents the diffusivity when there is no impediment to the particles’ motion, i.e., when the confinement value approaches 0. Therefore, the confinement tensor model can accommodate both restricted and unrestricted diffusion. In a recent study, the orientationally-averaged signal was studied for confined diffusion measured via single- and double diffusion encoding measurements demonstrating that certain features of the NMR signal [27, 28] that cannot be predicted by diffusion tensors are reproduced by the confinement tensor model [29].

These findings suggest the confinement tensor model as a plausible alternative for representing non-exchanging microscopic domains in multicompartment specimen models. In this work, we therefore propose to incorporate this model into the so-called multi-dimensional MRI methods [14]. In particular, we replace the diffusion with the confinement tensor in Diffusion Tensor Distribution Imaging (DTDI) [4, 30], and illustrate that the moments of the DTD estimated using Q-space trajectory Imaging (QTI) [31] would have a different interpretation for confined diffusion. We start by assessing the capabilities of the confinement model in representing single pores on data simulated using a typical protocol involving general time-varying diffusion gradient fields [32], and then proceed with first attempts at recovering distributions of confinement tensors in a human brain dataset.

2 Background and Theory

2.1 Diffusion Under a Hookean Restoring Force

In a diffusion NMR experiment, diffusing molecules acquire a phase shift depending on their trajectory x(t) and on the time-varying magnetic field gradient G(t). The signal from all molecules can be expressed as

$E=\langle e^{-i\gamma \int \mathrm{d}t\mathbf{x}(t)\cdot \mathbf{G}(t)}\rangle ,(1)$

where γ is the gyromagnetic ratio, and the averaging is performed over all particle trajectories.

For the case of diffusion under a Hookean restoring force, we shall denote by C the confinement tensor, which, upon multiplication by the Boltzmann constant k_B and absolute temperature T, gives the tensorial force constant f = k_BTC defining the Hookean potential V through $V(\mathbf{r})=\frac{1}{2}{\mathbf{r}}^{⊺}\mathbf{f}\mathbf{r}$. Furthermore, we denote the possibly anisotropic free diffusion tensor by D, and assume that D and C commute, i.e., they share the same eigendirections. Finally, we introduce Ω = DC for brevity.

Statistical quantities, such as the signal, can be calculated using the path weight

$\mathrm{Pr}[\mathbf{x}()]\propto \exp \left(-\frac{1}{4}\int \mathrm{d}t{\left(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}+\mathbf{\Omega }\mathbf{x}(t)\right)}^{⊺}{\mathbf{D}}^{-1}\left(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}+\mathbf{\Omega }\mathbf{x}(t)\right)\right),(2)$

which represents the differential probability for a particle to follow the trajectory x(). The NMR signal in (Eq. 1) can thus be evaluated, up to a constant, through the path integral

$E\propto \int \mathcal{D}\mathbf{x}()\exp (-\int \mathrm{d}t(\frac{1}{4}{\left(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}+\mathbf{\Omega }\mathbf{x}(t)\right)}^{⊺}{\mathbf{D}}^{-1}\left(\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}+\mathbf{\Omega }\mathbf{x}(t)\right)+i\gamma \mathbf{x}(t)\cdot \mathbf{G}(t))).(3)$

Thanks to stationarity, the time integration can be taken from −∞ to ∞, in which case employing the substitutions

$\begin{array}{ll}\hfill \mathbf{x}(t)& =\int \frac{\mathrm{d}\omega }{2\pi }{e}^{i\omega t}\widehat{\mathbf{x}}(\omega )\hfill \end{array}(\mathrm{4a})$

$\begin{array}{ll}\hfill \mathbf{G}(t)& =\int \frac{\mathrm{d}\omega }{2\pi }{e}^{i\omega t}\widehat{\mathbf{G}}(\omega )\hfill \end{array}(\mathrm{4b})$

yields

$E=\exp \left(-\int \frac{\mathrm{d}\omega }{2\pi }{\widehat{\mathbf{G}}}^{†}(\omega )\mathbf{K}(\omega )\widehat{\mathbf{G}}(\omega )\right)(5)$

with

$\mathbf{K}(\omega )=2{\gamma }^2\mathbf{D}{({\omega }^2\mathbf{I}+{\mathbf{\Omega }}^2)}^{-1}.(6)$

Converting Eq. 5 to the time domain yields

$E=\exp \left(-\frac{{\gamma }^2}{2}\int \mathrm{d}t\int \mathrm{d}{t}^{\prime }{\mathbf{G}}^{⊺}(t^{\prime })\mathbf{D}{\mathbf{\Omega }}^{-1}{e}^{-\mathbf{\Omega }|t-t^{\prime }|}\mathbf{G}(t)\right).(7)$

2.2 The Confinement Tensor Model

Figure 1A shows the parameters of the confinement tensor model [17] and how they represent a generic pore. The shape of the subdomain is captured by an effective confinement tensor C with units of inverse squared length, like in the case of diffusion under a Hookean force as described above. On the other hand, the rate of water diffusivity is captured by a scalar effective isotropic diffusivity D_eff.

FIGURE 1

FIGURE 1. (A) Confinement tensor model parameters. For a subdomain of generic shape, the geometry is captured by the effective confinement tensor, while the water diffusivity within is represented by an effective isotropic diffusivity. The tensor representing the shape is a second order tensor, while the effective diffusivity is a scalar. (B) Signal for different confinement values for linear tensor encoding experiments at different b-values. For C → 0 the signal converges to the free diffusion regime, where the value is determined as exp ( − b D₀), with b being the b-value [33]. This is represented by the colored circles in the plot. For C → ∞, the signal converges to 1 indicating complete water immobility.

While the expression given in Eq. 7 is the natural way for defining the signal implied by the confinement tensor model, that is not optimal for the actual computation of the signal. To avoid potential numerical issues with the inversion of the Ω tensor within the integral, we find more advantageous to use the equivalent expression given by Yolcu et al. [17] for a gradient waveform applied between time points 0 and t_f,

$E=\exp \left(-D_{\text{eff}}{\int }_0^{t_f}\mathrm{d}t|\mathbf{Q}(t){|}^2\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{D_{\text{eff}}}{2}{\mathbf{Q}}^{⊺}(0){\mathbf{\Omega }}^{-1}\mathbf{Q}(0)\right)(8)$

with

$\mathbf{Q}(t)=\gamma {\int }_t^{t_f}\mathrm{d}{t}^{\prime }{e}^{-\mathbf{\Omega }(t^{\prime }-t)}\mathbf{G}(t^{\prime }).(9)$

Note that for C → 0, the signal in Eq. 8 reduces to the NMR signal expression for isotropic free diffusion (the proof is provided in [17]), while for C → ∞, the signal converges to 1, indicating particles’ immobility. Both these scenarios are shown in Figure 1B, where the signal for confinement values in the range [0, ∞) is shown.

2.3 md—dMRI With Confinement

The expressions derived thus far concern the MR signal for a single confinement tensor. Here, we instead consider the case where a distribution of such tensors is collectively giving rise to the signal. In particular, we provide the signal expression for a confinement tensor distribution (CTD), which could be used for performing Confinement Tensor Distribution Imaging (CTDI), and discuss employing QTI for locally confined diffusion.

2.3.1 Confinement Tensor Distribution

The NMR signal expression for a distribution $\mathcal{P}(\mathbf{D})$ of diffusion tensors is given by [4]

$S(\mathbf{b})=S_0\int \mathcal{P}(\mathbf{D})e^{-\mathbf{b}:\mathbf{D}}\mathrm{d}\mathbf{D},(10)$

where b is the measurement tensor [34], D is the diffusion tensor, and “ :^″ indicates the generalized scalar product between tensors. A similar expression can be introduced to include the confinement tensor model. Due to the extra parameter D_eff, the considered distribution becomes the joint distribution of effective confinement tensors and effective diffusivities $\mathcal{P}(\mathbf{C},D_{\text{eff}})$. The signal expression for an experiment determined by a general time varying magnetic field gradient G(t) is

$S(\mathbf{G}(t))=S_0\int \mathcal{P}(\mathbf{C},D_{\text{eff}})E(\mathbf{G}(t),\mathbf{C},D_{\text{eff}})\mathrm{d}\mathbf{C}\mathrm{d}{D}_{\text{eff}},(11)$

where E (G(t), C, D_eff) is as defined in Eq. 7 or Eq. 8. Eqs 10, 11 can be considered to be generalizations of the Laplace transforms of $\mathcal{P}(\mathbf{D})$ and $\mathcal{P}(\mathbf{C},D_{\text{eff}})$, respectively. Recovering the $\mathcal{P}(\mathbf{D})$ (or $\mathcal{P}(\mathbf{C},D_{\text{eff}})$) from a series of measurements, i.e., numerically inverting the Laplace transform, is well known to be an ill-posed problem [30, 35, 36].

2.3.2 QTI for Locally Confined Diffusion

The QTI technique exploits the sensitivity of the detected signal to the statistical moments of the structural parameters describing the specimen [37]. When the DTD model is employed for the latter, the high signal (low diffusion sensitivity) regime reveals the first few moments of the diffusivities [31, 38]. For a DTD characterized by the distribution P(D) the signal decay in (10) can be expressed by

$E(\mathbf{b})=⟨e^{-b_{k\ell }{D}_{k\ell }}⟩,(12)$

where we employed the Einstein summation convention. At low diffusion sensitivity, the natural logarithm of the signal decay is approximated by the Maclaurin series of the above expression around b ≈ 0, yielding

$\ln E(\mathbf{b})\approx -b_{k\ell }⟨D_{k\ell }⟩+\frac{1}{2}{b}_{k\ell }{b}_{mn}{⟨D_{k\ell }{D}_{mn}⟩}_c,(13)$

where the last quantity is the second cumulant, i.e., ${⟨D_{k\ell }{D}_{mn}⟩}_c=⟨D_{k\ell }{D}_{mn}⟩-⟨D_{k\ell }⟩⟨D_{mn}⟩$. We remind that the components of the measurement tensor are given through [34, 39].

$b_{k\ell }={\gamma }^2{\int }_0^{t_f}\mathrm{d}t{\int }_0^t\mathrm{d}{t}^{\prime }{\int }_0^t\mathrm{d}{t}^{″}{G}_k(t^{\prime })G_{\ell }(t^{″}).(14)$

In the case of a CTD, the averaging takes the form of an integration over D_eff as well as C; see (Eq. 11). Applying the same procedure for the signal in frequency domain, (Eq. 5), yields

$\begin{array}{ll}\hfill \ln E(\widehat{\mathbf{G}}(\omega ))\approx & -\int \frac{\mathrm{d}\omega }{2\pi }{\widehat{G}}_k(\omega ){\widehat{G}}_{\ell }(\omega )\langle K_{k\ell }(\omega )\rangle \hfill \\ \hfill & +\frac{1}{2}\int \frac{\mathrm{d}\omega }{2\pi }\int \frac{\mathrm{d}{\omega }^{\prime }}{2\pi }{\widehat{G}}_k(\omega ){\widehat{G}}_{\ell }(\omega ){\widehat{G}}_m({\omega }^{\prime }){\widehat{G}}_n({\omega }^{\prime })\langle K_{k\ell }(\omega )K_{mn}{({\omega }^{\prime })\rangle }_c,\hfill \end{array}(15)$

while the same is given in the time-domain by

$\begin{array}{ll}\hfill \ln E(\mathbf{G}(t))\approx & -\int \mathrm{d}t\int \mathrm{d}{t}^{\prime }{G}_k(t)G_{\ell }(t^{\prime })\langle H_{k\ell }(t,t^{\prime })\rangle \hfill \\ \hfill & +\frac{1}{2}\int \mathrm{d}t\int \mathrm{d}{t}^{\prime }\int \mathrm{d}{t}^{″}\int \mathrm{d}t‴G_k(t)G_{\ell }(t^{\prime })G_m(t^{″})G_n(t‴)\langle H_{k\ell }(t,t^{\prime })H_{mn}{(t^{″},t‴)\rangle }_c,\hfill \end{array}(16)$

where

$H_{k\ell }(t,t^{\prime })=\frac{{\gamma }^2}{2}{\mathbf{C}}^{-1}{e}^{-\mathbf{D}\mathbf{C}|t^{\prime }-t|}(17)$

since C and D commute.

Note that, the shape of the waveform G(t) is inextricably linked to the signal in the CTD case. Furthermore, the interpretations of the signal decay rate are substantially different for the CTD and DTD assumptions. Thus, when QTI is performed, one can quantify only apparent moments of a DTD while the same analysis employing the CTD model would provide a more meaningful description of the low diffusion sensitivity regime.

3 Implementation

The confinement tensor model was implemented in Matlab (The Mathworks Inc., Natick, Massachussets) according to Eqs 8, 9. Numerical integration was performed using the trapezoidal rule. The signal computation for a given confinement tensor and effective isotropic diffusivity was carried out in a reference frame in which Ω is diagonal. This is achieved by rotating the measurement waveforms G(t) with the rotation matrix determined by the eigenvectors of Ω. This allows for computations to be carried out separately for each of the eigenvalues of Ω thanks to the separability of the model [17]. This approach mitigates numerical issues that arise for small confinement values, in which case inverting Ω may become problematic. Possible issues can be alleviated by considering a Taylor expansion for the second exponential factor in Eq. 8 to remove the dependency of that part of the signal on Ω. The derivation of the expression for computing the approximation of the signal using the Taylor expansion is provided in Appendix.

To fit the confinement model to the data, we used the Matlab function lsqnonlin with Levenberg-Marquardt as algorithm. The unknown estimated quantities consist of the signal without diffusion weighting S₀, the effective diffusivity D_eff, and the six unique components of the confinement tensor C. During the fitting, the tensor Ω is replaced by its Cholesky factorization to ensure the positive semidefiniteness of the estimated confinement tensor [40].

To estimate the joint distribution of confinement tensors and effective isotropic diffusivities, we adapted the existing technology implementing a Monte Carlo inversion of Eq. 10, as detailed in [36] and retrieved from https://github.com/markus-nilsson/md-dmri. As for the original implementation, we limit ourselves to the case of axisymmetric tensors. These can be represented using 4 parameters: the parallel and perpendicular confinement (C_para and C_perp) capture the pore’s geometry, while the other two define the pore orientation through the azimuthal (ϕ) and polar (θ) angles. Altogether, each pore is represented by 5 parameters (C_para, C_perp, ϕ, θ, and D_eff). While performing the inversion, these parameters are searched within the limits 8 ≤ log₁₀(C_para/m⁻²), log₁₀(C_perp/m⁻²) ≤ 12, 0.1 ≤ (D_eff/μm²ms⁻¹) ≤ 3.2, 0 ≤ ϕ ≤ 2π, 0 ≤   cos(θ) ≤ 1. For each voxel, the recovered $\mathcal{P}(\mathbf{C},D_{\text{eff}})$ can be visualized in 3D plots where C_para and C_perp vary along the x and y axes, while D_eff varies along the z axis, respectively. The pore direction is encoded using the RGB color scale. We adopted the convention of displaying the color according to the main diffusion direction, not according to the direction of maximum confinement.

4 Results

4.1 Signal for Single Compartments

In this section we investigate the capabilities of the confinement tensor model in capturing features of both free and restricted diffusion in data where the diffusion sensitization is achieved with general time varying magnetic field gradients. We employ a protocol featuring 217 measurements comprising Linear, Planar, and Spherical Tensor Encoding (LTE, PTE, and STE respectively) [32]. We refer to this protocol as tensor encoding. Signals for diffusion taking place in closed and open geometries were computed using this protocol. The pore shapes and the respective defining parameters were as follows:

• Free isotropic diffusion, D₀ = 3 μm²/ms

• Stick compartment (uni-directional free diffusion), D₀ = 2.5 μm²/ms

• Infinite cylinder, r = 5 μm, D₀ = 3 μm²/ms

• Capped cylinder 1, l = 12 μm, r = 2 μm, D₀ = 2 μm²/ms

• Capped cylinder 2, l = 10 μm, r = 1.5 μm, D₀ = 2.5 μm²/ms,

• Sphere, r = 5 μm, D₀ = 2 μm²/ms

where D₀ is the bulk diffusivity, r is the radius, and l is the length. The diffusion tensor model was used to generate the signals for the free diffusion and the stick compartments, while the method described in [12] was used to generate the signals for the cylinders and the sphere.

Figure 2A shows the results obtained by fitting the confinement tensor model and the diffusion tensor to the simulated signals. Only a subset of the measurements is shown for easier visual inspection. For the considered protocol, the confinement tensor model seems to capture well the features of both free and restricted diffusion, suggesting that the model given in Eq. 8 has sufficient degrees of freedom and there is no need for employing a tensorial diffusivity. Note that as illustrated for the 1D problem of restricted diffusion between two parallel plates [26], in the ideal scenario involving only very long pulses and simple geometries, one would expect the measurement to be sensitive only to the product of squared confinement and diffusivity, in which case there is no need to employ an effective diffusivity in the model. However, in practice there is such sensitivity, and the relationship between the parameters of the model (D_eff, and C) and those of the geometry requires further investigations. What is remarkable however is that having only one additional parameter (D_eff) offers sufficient complexity to capture the information in the signal for the considered acquisition scenario. This is also evident in Figure 2B, which better illustrates how the assumption of free diffusion fails [3, 41] while the confinement tensor model fully captures the signal modulations due to restricted diffusion probed by STE measurements.

FIGURE 2

FIGURE 2. Confinement and diffusion tensor fit to signals arising from particles diffusing in free (Free Diffusion, Stick) and restricted (Infinite Cylinder, Capped Cylinder 1, Capped Cylinder 2, Sphere) geometries. In both panels (A,B) all the 217 data points were used to fit the model. For ease of interpretation, in panel (A) only 52 random points out of the 217 are shown. In panel b) all the data points produced with Spherical Tensor Encoding are shown. Note how the confinement tensor model captures features of diffusion in both free and restricted scenarios. In particular, observe in panel (B) how the signature of restricted diffusion imprinted on the signals produced with STE can be captured by the confinement tensor model while are inevitably missed by the diffusion tensor.

The recovered values of the effective diffusivity D_eff coincided with the bulk diffusivity D₀ for the stick and free water compartments. For the compartments in which diffusion is fully restricted, the estimated values were lower than the nominal D₀. Respectively, 1.65 μm²/ms for Capped Cylinder 1, 2.0 μm²/ms for Capped Cylinder 2, and 1.7 μm²/ms for Sphere.

4.2 Signal for a Distribution of Compartments

Illustrative results obtainable with the Monte Carlo inversion method described in Section 3 were produced on a publicly available brain dataset [32] collected via the tensor encoding protocol used in the simulations. Figure 3 shows Monte Carlo inversion results on a few selected voxels on the bottom, and a bird’s-eye view on the plane displaying the geometric information about the tensors in the distribution on the top. In this last, the colored areas indicate what shape each tensor would have for different values of C_para and C_perp. Pores with isotropic geometries are found along the diagonal where C_para = C_perp. Free isotropic diffusion is found for low C_para = C_perp, while confined isotropic diffusion is found for high C_para = C_perp. Stick-like pores are located at the C_para ≪C_perp corner, while pancake-like pores are found at the C_para ≫C_perp corner.

FIGURE 3

FIGURE 3. Example results of inverting Eq. 11 on real data collected with LTE, PTE, and STE waveforms [32]. The panel on top shows a bird’s eye view on the geometry plane, and how to associate the location of points to a shape according to the C_para and C_perp coordinates. The contours in the 3D plot show the projections of the pores’ clusters onto the various 2D planes.

In the same spirit of what was shown in [42] for the DTD model, the 3D plots in Figure 3 show what the $\mathcal{P}(\mathbf{C},D_{\text{eff}})$ for voxels containing either single or multiple types of brain tissues, as obtained from the data, could be. For example, voxels containing pure CSF would have a $\mathcal{P}(\mathbf{C},D_{\text{eff}})$ of only free diffusion geometry with D_eff ≈ 3.1 μm²/ms. Pure white matter voxels would only contain collections of stick-like geometries (see the voxel from the Corpus Callosum), while, for the considered dataset, gray matter could contain isotropic free water at lower D_eff compared to free water. Voxels with mixed tissue types could build their $\mathcal{P}(\mathbf{C},D_{\text{eff}})$ based on those from single tissues. Note how having separate components encoding for the pore geometry (C) and the water diffusivity (D_eff) allows for clearly identifying scenarios where pores could have the same shape, but different water mobility. See for example the voxel “WM, 2 fiber bundles?“, where the distribution seems to suggest the presence of two differently oriented fiber bundles, which can be teased apart also by looking at their water diffusivity. The same specificity could not be achieved by only considering a distribution of diffusion tensors, where the information about the pore geometry is inextricably entangled to that of water diffusivity.

Note however that all what is presented and discussed in this section are simply initial conjectures, which may very well be the results of falling into the temptation of over-interpreting the outcomes of the DTD or CTD estimation. As shown by [36], inverting Eq. 10, and by extension Eq. 11, is already challenging even at infinite SNR. The situation worsens in real data where the validity of the solutions proves to be very sensitive to the presence of noise. Moreover, as we present and discuss later, the results, and their interpretation, strongly depend upon the adopted acquisition scheme. For example, we could expect to find sphere-like compartments in gray matter in data encoded with different diffusion times and higher diffusion sensitivity, possibly indicating that a relevant fraction of the signal is due to cell bodies. Moreover, while the considered protocol (and data) should encode sufficient information for accurately recovering the pores’ geometry, other waveforms could prove beneficial to study the time-dependence of the diffusion process, augmenting the reliability of the D_eff dimension.

5 Discussion

The results in Figure 2 illustrate how well the confinement model captures the features of both free and restricted diffusion, for data simulated with a clinically-feasible protocol including typical time-varying magnetic field gradients. The signal’s modulation due to restrictions is, under the considered experimental set-up, fully described by studying the problem of diffusion occurring in a potential landscape. This shows that the considered approach retains the right number of degrees of freedom to characterize diffusion processes within individual compartments. This finding is consistent with what was reported by Özarslan et al. [26] for experiments involving long duration pulses, and supports the idea of adopting the confinement tensor for representing isolated pores in multicompartment models.

Having a single model covering both restricted and unrestricted diffusion in different geometries could be advantageous when defining multi-compartmental models based on such shapes as building blocks. Biophysical models, such as the composite hindered and restricted model of diffusion (CHARMED) [43] and neurite orientation dispersion and density imaging (NODDI) [44] strive for modeling specifically the neural tissue, therefore are not suitable for different tissues and other heterogeneous media. The confinement tensor representation of each compartment could be integrated into such models and could provide a convenient means for accounting for restricted diffusion. On the other hand, the confinement tensor distribution model is far more general than such specific models as one would not need to make a priori assumptions on the specimen composition, apart from limiting its representation to numerous non-exchanging and possibly confined domains. The results in Figure 3 exemplified the specificity achievable by modeling a specimen with a joint distribution of confinement tensors-isotropic effective diffusivities. Other information about the water pools, such as T₁ and T₂ relaxations [42, 45], could be added to increase specificity to the tissue heterogeneity. Similarly to what was presented there, the confinement tensor model could also be considered for diffusion—relaxation studies [42, 46–49].

Note that on the specific dataset used in this work we did not observe striking modulation in the signal for isotropic measurements at constant b-value. This could be explained by the experiments not being sensitive to finite-sized anisotropic restrictions, i.e., axons could effectively be pictured as sticks. Under these conditions, the fit to signal for both the DTD and CTD would yield very close results. Having two fundamentally different models exhibiting good fits to the data suggests that the data is possibly not descriptive enough. Another factor contributing to equal performance could be found in both DTD and CTD being overly-parameterized, thus effectively having the capabilities to fit the data equally successfully. This should not however be interpreted as both models being acceptable and providing informative results. In addition, based on the results in Figure 2, we expect the CTD to provide meaningful information on data where restrictions have imprinted a clearer signature.

We would also care to iterate once more on the limitation of performing CTDI (or DTDI) using the technology implemented in https://github.com/markus-nilsson/md-dmri, due to the mathematically ill-posed problem that is being attempted. Different $\mathcal{P}(\mathbf{C},D_{\text{eff}})$ in Eq. 11 will represent the signal equally well, thus possibly leading to wrong interpretations of the microstructural characteristics of the scanned specimen. A similar issue is referred to as the “degeneracy problem” [50] in recovering the brain microstructure. Multi-compartment models present flat fitting landscapes with multiple local minima located in different parts of the parameter space, each of which providing equally sound biological explanation for the signal. One approach to alleviate the problem involves including additional measurements, e. g, diffusion measurements having different temporal profiles, with the goal of disentangling the contribution of different parameters to the model interpretation [51].

When attempting at recovering the joint distribution of confinement tensors and isotropic effective diffusivities, we found from simulations that the pores’ geometry can be obtained relatively faithfully using the tensor encoding protocol only. However, determining the pore diffusivity relatively to the restriction size from data encoded exclusively in such manner seems to be more challenging. We provide examples of this in Figure 4. In Figure 4A we show a modified version of the tensor encoding protocol, where half of the waveforms are replaced with Trapezoidal-Cosine Oscillating Gradient Spin-Echo (TC-OGSE). We refer to this protocol as mixed. The goal is to achieve higher sensitivity to molecules’ diffusivity within restrictions by using waveforms with well defined encoding frequency, and by matching the frequency of such waveform to that of the diffusion process [52–54]. Retaining part of the original protocol should ensure accurate pore geometry estimation. In Figures 4B,C we show the results obtained for a simple scenario where the specimen consists of two pools of water in which molecules are diffusing at different rates. When the data are simulated with only LTE, PTE, and STE, it is possible to accurately recover the expected pore shapes but not the water diffusivity. Conversely, by introducing oscillating gradients, the diffusivity estimates, although still uncertain, converge to the correct values.

FIGURE 4

FIGURE 4. (A) Layout of the mixed protocol. 102 measurements consist of LTE, PTE, and STE waveforms as defined in the tensor encoding protocol, but with their directions redistributed over the sphere to achieve more uniform directional sampling. Another 102 measurements consist of Trapezoidal Cosine Oscillating Gradient Spin Echo waveforms oscillating at 4 different frequencies. The last 13 measurements had null diffusion gradient strength. (B) Results of inverting Eq. 11 for a simple system consisting of two equally weighted water pools with identical geometries but different water diffusivities. Left, data encoded with the tensor encoding protocol. Right, data encoded with the mixed protocol. The red filled dots depict the ground truth values in the 3D plot, while the red crosses show the ground truth values for the projections onto the various planes. (C) The distributions of D_eff obtained in (B). (D) Results inverting Eq. 11 for a substrate comprising a free water compartment, a sphere compartment, a stick compartment, and an intra-axonal compartment, with the following weights in the distribution: 0.25, 0.25, 0.4, 0.1. The substrate ground truth is visualized in the small 3D plot on top. On the left, data encoded with the tensor encoding protocol. On the right, data encoded with the mixed protocol.

In Figure 4D we show the results on a more complex substrate consisting of a stick compartment, an extra-axonal compartment, a sphere compartment and a compartment with free diffusion. As for the simple scenario described above, the estimation of the pores’ diffusivity, in particular one of the sphere compartments, improves when the mixed protocol is used. This corroborates the idea of including measurements not only exploring different shapes of the encoding tensors, but also probing different frequencies [53]. The results in Figure 4D also exemplify the uncertainty around inverting the Laplace transform, even for infinite SNR.

6 Conclusion

In this work we incorporated the confinement tensor model for individual subdomains of heterogeneous media into multidimensional diffusion MRI frameworks. We demonstrated how considering Brownian motion as taking place under the influence of a Hookean potential provides sufficient degrees of freedom to capture the signal modulations arising from water diffusing in restricted geometries. We argued that the confinement tensor distribution (CTD) model is a viable alternative to the diffusion tensor distribution model as CTD relies on the effective model of restricted diffusion, which makes it more consistent with the multicompartmental organization of complex tissue when examined via commonly performed diffusion MRI measurements. Despite its challenges, incorporating this model into multidimensional diffusion MRI methods could provide new insights regarding the structural composition of complex media.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/filip-szczepankiewicz/Szczepankiewicz_DIB_2019.

Author Contributions

DB: Methodology, Software, Validation, Formal analysis, Investigation, Writing—Original Draft, Visualization. CY: Methodology, Theory, Writing—Review and Editing. EÖ: Conceptualization, Methodology, Writing—Review and Editing, Supervision, Project administration, Funding acquisition.

Funding

This project was financially supported by Linköping University Center for Industrial Information Technology (CENIIT), Analytic Imaging Diagnostics Arena (AIDA), and Sweden’s Innovation Agency (VINNOVA) ASSIST.

Conflict of Interest

EÖ is a shareholder of Spin Nord AB.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank Magnus Herberthson for stimulating discussions.

References

1. Panagiotaki E, Schneider T, Siow B, Hall MG, Lythgoe MF, Alexander DC. Compartment Models of the Diffusion MR Signal in Brain white Matter: A Taxonomy and Comparison. NeuroImage (2012) 59:2241–54. doi:10.1016/j.neuroimage.2011.09.081

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Jelescu IO, Budde MD. Design and Validation of Diffusion MRI Models of White Matter. Front Phys (2017) 5:61. doi:10.3389/fphy.2017.00061

PubMed Abstract | CrossRef Full Text | Google Scholar

3. de Swiet TM, Mitra PP. Possible Systematic Errors in Single-Shot Measurements of the Trace of the Diffusion Tensor. J Magn Reson Ser B (1996) 111:15–22. doi:10.1006/jmrb.1996.0055

CrossRef Full Text | Google Scholar

4. Jian B, Vemuri BC, Özarslan E, Carney PR, Mareci TH. A Novel Tensor Distribution Model for the Diffusion-Weighted MR Signal. NeuroImage (2007) 37:164–76. doi:10.1016/j.neuroimage.2007.03.074

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Basser PJ, Pajevic S. A normal Distribution for Tensor-Valued Random Variables: Applications to Diffusion Tensor MRI. IEEE Trans Med Imaging (2003) 22:785–94. doi:10.1109/TMI.2003.815059

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Magdoom KN, Pajevic S, Dario G, Basser PJ. A New Framework for MR Diffusion Tensor Distribution. Sci Rep (2021) 11:2766. doi:10.1038/s41598-021-81264-x

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Shakya S, Batool N, Özarslan E, Knutsson H. Multi-fiber Reconstruction Using Probabilistic Mixture Models for Diffusion MRI Examinations of the Brain. In: T Schultz, E Özarslan, and I Hotz, editors. Modeling, Analysis, and Visualization of Anisotropy. Cham: Springer International Publishing (2017). 283–308. doi:10.1007/978-3-319-61358-1_12

CrossRef Full Text | Google Scholar

8. Herberthson M, Yolcu C, Knutsson H, Westin C-F, Özarslan E. Orientationally-averaged Diffusion-Attenuated Magnetic Resonance Signal for Locally-Anisotropic Diffusion. Sci Rep (2019) 9:4899. doi:10.1038/s41598-019-41317-8

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Cory DG, Garroway AN, Miller JB. Applications of Spin Transport as a Probe of Local Geometry. Abstr Pap Am Chem Soc (1990) 199:105.

Google Scholar

10. Wong EC, Cox RW, Song AW. Optimized Isotropic Diffusion Weighting. Magn Reson Med (1995) 34:139–43. doi:10.1002/mrm.1910340202

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Caprihan A, Wang LZ, Fukushima E. A Multiple-Narrow-Pulse Approximation for Restricted Diffusion in a Time-Varying Field Gradient. J Magn Reson Ser A (1996) 118:94–102. doi:10.1006/jmra.1996.0013

CrossRef Full Text | Google Scholar

12. Özarslan E, Shemesh N, Basser PJ. A General Framework to Quantify the Effect of Restricted Diffusion on the NMR Signal with Applications to Double Pulsed Field Gradient NMR Experiments. J Chem Phys (2009) 130:104702. doi:10.1063/1.3082078

CrossRef Full Text | Google Scholar

13. Westin C-F, Szczepankiewicz F, Pasternak O, Özarslan E, Topgaard D, Knutsson H, et al. Measurement Tensors in Diffusion MRI: Generalizing the Concept of Diffusion Encoding. In: P Golland, N Hata, C Barillot, J Hornegger, and R Howe, editors. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2014. Cham: Springer International Publishing (2014). p. 209–16. doi:10.1007/978-3-319-10443-0_27

CrossRef Full Text | Google Scholar

14. Topgaard D. Multidimensional Diffusion MRI. J Magn Reson (2017) 275:98–113. doi:10.1016/j.jmr.2016.12.007

CrossRef Full Text | Google Scholar

15. Beaulieu C, Allen PS. Determinants of Anisotropic Water Diffusion in Nerves. Magn Reson Med (1994) 31:394–400. doi:10.1002/mrm.1910310408

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Woessner DE. N.M.R. Spin-echo Self-Diffusion Measurements on Fluids Undergoing Restricted Diffusion. J Phys Chem (1963) 67:1365–7. doi:10.1021/j100800a509

CrossRef Full Text | Google Scholar

17. Yolcu C, Memiç M, Şimşek K, Westin CF, Özarslan E. NMR Signal for Particles Diffusing under Potentials: From Path Integrals and Numerical Methods to a Model of Diffusion Anisotropy. Phys Rev E (2016) 93:052602. doi:10.1103/PhysRevE.93.052602

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Uhlenbeck GE, Ornstein LS. On the Theory of the Brownian Motion. Phys Rev (1930) 36:823–41. doi:10.1103/physrev.36.823

CrossRef Full Text | Google Scholar

19. Stejskal EO. Use of Spin Echoes in a Pulsed Magnetic-Field Gradient to Study Anisotropic, Restricted Diffusion and Flow. J Chem Phys (1965) 43:3597–603. doi:10.1063/1.1696526

CrossRef Full Text | Google Scholar

20. Le Doussal P, Sen PN. Decay of Nuclear Magnetization by Diffusion in a Parabolic Magnetic Field: An Exactly Solvable Model. Phys Rev B (1992) 46:3465–85. doi:10.1103/PhysRevB.46.3465

CrossRef Full Text | Google Scholar

21. Mitra P, Halperin B. Effects of Finite Gradient-Pulse Widths in Pulsed-Field-Gradient Diffusion Measurements. J Magn Reson Ser A (1995) 113:94–101. doi:10.1006/jmra.1995.1060

CrossRef Full Text | Google Scholar

22. Afzali M, Yolcu C, Özarslan E. Characterizing Diffusion Anisotropy for Molecules under the Influence of a Parabolic Potential: A Plausible Alternative to DTI. Proc Intl Soc Mag Reson Med (2015) 23:2795.

Google Scholar

23. Zucchelli M, Afzali M, Yolcu C, Westin CF, Menegaz G, Özarslan E. The Confinement Tensor Model Improves Characterization of Diffusion-Weighted Magnetic Resonance Data with Varied Timing Parameters. In: 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI), Prague, Czech Republic, April 13–16, 2016. IEEE (2016). 1093–6. doi:10.1109/isbi.2016.7493456

CrossRef Full Text | Google Scholar

24. Liu C, Özarslan E. Multimodal Integration of Diffusion MRI for Better Characterization of Tissue Biology. NMR Biomed (2019) 32:e3939. doi:10.1002/nbm.3939

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Afzali M, Pieciak T, Newman S, Garyfallidis E, Özarslan E, Cheng H, et al. The Sensitivity of Diffusion MRI to Microstructural Properties and Experimental Factors. J Neurosci Methods (2021) 347:108951. doi:10.1016/j.jneumeth.2020.108951

CrossRef Full Text | Google Scholar

26. Özarslan E, Yolcu C, Herberthson M, Westin CF, Knutsson H. Effective Potential for Magnetic Resonance Measurements of Restricted Diffusion. Front Phys (2017) 5:68. doi:10.3389/fphy.2017.00068

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Mitra PP. Multiple Wave-Vector Extensions of the NMR Pulsed-Field-Gradient Spin-echo Diffusion Measurement. Phys Rev B (1995) 51:15074–8. doi:10.1103/physrevb.51.15074

CrossRef Full Text | Google Scholar

28. Özarslan E, Basser PJ. Microscopic Anisotropy Revealed by NMR Double Pulsed Field Gradient Experiments with Arbitrary Timing Parameters. J Chem Phys (2008) 128:154511. doi:10.1063/1.2905765

CrossRef Full Text | Google Scholar

29. Yolcu C, Herberthson M, Westin CF, Özarslan E. Magnetic Resonance Assessment of Effective Confinement Anisotropy with Orientationally-Averaged Single and Double Diffusion Encoding. In: E Özarslan, T Schultz, E Zhang, and A Fuster, editors. Anisotropy across Fields and Scales. Cham: Springer Nature (2021). 203–23. Mathematics and Visualization. doi:10.1007/978-3-030-56215-1_10

CrossRef Full Text | Google Scholar

30. Topgaard D. Diffusion Tensor Distribution Imaging. NMR Biomed (2019) 32:e4066. doi:10.1002/nbm.4066

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Westin CF, Knutsson H, Pasternak O, Szczepankiewicz F, Özarslan E, van Westen D, et al. Q-space Trajectory Imaging for Multidimensional Diffusion MRI of the Human Brain. NeuroImage (2016) 135:345–62. doi:10.1016/j.neuroimage.2016.02.039

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Szczepankiewicz F, Hoge S, Westin CF. Linear, Planar and Spherical Tensor-Valued Diffusion MRI Data by Free Waveform Encoding in Healthy Brain, Water, Oil and Liquid Crystals. Data in Brief (2019) 25:104208. doi:10.1016/j.dib.2019.104208

PubMed Abstract | CrossRef Full Text | Google Scholar

33. LeBihan D, Breton E. Imagerie de diffusion In Vivo par résonance magnétique nucléaire. C R Acad Sc (1985) 15:1109–12.

Google Scholar

34. Mattiello J, Basser PJ, LeBihan D. Analytical Expressions for the B-Matrix in NMR Diffusion Imaging and Spectroscopy. J Magn Reson A (1994) 108:131–41. doi:10.1006/jmra.1994.1103

CrossRef Full Text | Google Scholar

35. Galvosas P, Callaghan PT. Multi-dimensional Inverse Laplace Spectroscopy in the NMR of Porous media. Comptes Rendus Physique (2010) 11:172–80. doi:10.1016/j.crhy.2010.06.014

CrossRef Full Text | Google Scholar

36. Reymbaut A, Mezzani P, Almeida Martins JP, Topgaard D. Accuracy and Precision of Statistical Descriptors Obtained from Multidimensional Diffusion Signal Inversion Algorithms. NMR Biomed (2020) 33, e4267. doi:10.1002/nbm.4267

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Özarslan E, Shemesh N, Koay CG, Cohen Y, Basser PJ. Nuclear Magnetic Resonance Characterization of General Compartment Size Distributions. New J Phys (2011) 13:15010. doi:10.1088/1367-2630/13/1/015010

CrossRef Full Text | Google Scholar

38. Lasič S, Szczepankiewicz F, Eriksson S, Nilsson M, 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 (2014) 2:11. doi:10.3389/fphy.2014.00011

CrossRef Full Text | Google Scholar

39. Özarslan E, Westin CF, Mareci TH. Characterizing Magnetic Resonance Signal Decay Due to Gaussian Diffusion: the Path Integral Approach and a Convenient Computational Method. Concepts Magn Reson A (2015) 44:203–13. doi:10.1002/cmr.a.21354

CrossRef Full Text | Google Scholar

40. Koay CG, Chang LC, Carew JD, Pierpaoli C, Basser PJ. A Unifying Theoretical and Algorithmic Framework for Least Squares Methods of Estimation in Diffusion Tensor Imaging. J Magn Reson (2006) 182:115–25. doi:10.1016/j.jmr.2006.06.020

CrossRef Full Text | Google Scholar

41. Jespersen SN, Olesen JL, Ianuş A, Shemesh N. Effects of Nongaussian Diffusion on “Isotropic Diffusion” Measurements: An Ex-Vivo Microimaging and Simulation Study. J Magn Reson (2019) 300:84–94. doi:10.1016/j.jmr.2019.01.007

CrossRef Full Text | Google Scholar

42. Topgaard D. Multiple Dimensions for Random Walks. J Magn Reson (2019) 306:150–4. doi:10.1016/j.jmr.2019.07.024

CrossRef Full Text | Google Scholar

43. Assaf Y, Basser PJ. Composite Hindered and Restricted Model of Diffusion (Charmed) Mr Imaging of the Human Brain. NeuroImage (2005) 27:48–58. doi:10.1016/j.neuroimage.2005.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. Noddi: Practical In Vivo Neurite Orientation Dispersion and Density Imaging of the Human Brain. NeuroImage (2012) 61:1000–16. doi:10.1016/j.neuroimage.2012.03.072

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Callaghan PT, Arns CH, Galvosas P, Hunter MW, Qiao Y, Washburn KE. Recent Fourier and Laplace Perspectives for Multidimensional NMR in Porous media. Magn Reson Imaging (2007) 25:441–4. doi:10.1016/j.mri.2007.01.114

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Benjamini D, Basser PJ. Multidimensional Correlation MRI. NMR Biomed (2020) 33. doi:10.1002/nbm.4226

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Slator PJ, Palombo M, Miller KL, Westin CF, Laun F, Kim D, et al. Combined Diffusion-Relaxometry Microstructure Imaging: Current Status and Future Prospects. Magn Reson Med (2021) 86:2987–3011. doi:10.1002/mrm.28963

PubMed Abstract | CrossRef Full Text | Google Scholar

48. de Almeida Martins JP, Tax CMW, Reymbaut A, Szczepankiewicz F, Chamberland M, Jones DK, et al. Computing and Visualising Intra-voxel Orientation-specific Relaxation-Diffusion Features in the Human Brain. Hum Brain Mapp (2021) 42:310–28. doi:10.1002/hbm.25224

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Martin J, Reymbaut A, Schmidt M, Doerfler A, Uder M, Laun FB, et al. Nonparametric D-R1-R2 Distribution MRI of the Living Human Brain. Neuroimage (2021) 245:118753. doi:10.1016/j.neuroimage.2021.118753

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Jelescu IO, Veraart J, Fieremans E, Novikov DS. Degeneracy in Model Parameter Estimation for Multi-Compartmental Diffusion in Neuronal Tissue: Degeneracy in Model Parameter Estimation of Diffusion in Neural Tissue. NMR Biomed (2016) 29:33–47. doi:10.1002/nbm.3450

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Coelho S, Pozo JM, Jespersen SN, Jones DK, Frangi AF. Resolving Degeneracy in Diffusion MRI Biophysical Model Parameter Estimation Using Double Diffusion Encoding. Magn Reson Med (2019) 82:395–410. doi:10.1002/mrm.27714

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Stepišnik J. Time-dependent Self-Diffusion by NMR Spin-echo. Physica B (1993) 183:343–50.

Google Scholar

53. Lundell H, Nilsson M, Dyrby TB, Parker GJM, Cristinacce PLH, Zhou FL, et al. Multidimensional Diffusion MRI with Spectrally Modulated Gradients Reveals Unprecedented Microstructural Detail. Scientific Rep (2019) 9:9026. doi:10.1038/s41598-019-45235-7

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Drobnjak I, Cruz G, Alexander DC. Optimising Oscillating Waveform-Shape for Pore Size Sensitivity in Diffusion-Weighted MR. Microporous Mesoporous Mater (2013) 178:11–4. doi:10.1016/j.micromeso.2013.03.004

CrossRef Full Text | Google Scholar

Appendix: Numerical Approximation of the Signal in C → 0 Regime

The matrix inversion in the second exponential in equation (Eq. 8) can become numerically unstable for C → 0. When working in the coordinate system determined by the eigenvectors of C, the matrix Ω is diagonal

$\mathbf{\Omega }=\left[\begin{array}{ccc}\hfill {\lambda }_1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\lambda }_2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\lambda }_3\hfill \end{array}\right].(18)$

We take ${\vec{\mathbf{v}}}_1$, ${\vec{\mathbf{v}}}_2$, and ${\vec{\mathbf{v}}}_3$ to be the eigenvectors of C, defining the new coordinate system for the experiment. The matrix R having ${\vec{\mathbf{v}}}_1$, ${\vec{\mathbf{v}}}_2$, and ${\vec{\mathbf{v}}}_3$ as columns, can be used to determine the gradient waveforms used to collect the data in the new coordinate system through

${\mathbf{G}}^{\prime }(t)={[G_1^{\prime }(t),G_2^{\prime }(t),G_3^{\prime }(t)]}^{⊺}={\mathbf{R}}^{⊺}\mathbf{G}(t).(19)$

Then, the signal contribution from the second exponential in Eq. 8 can be written as

$\begin{array}{l}\hfill \text{exp}\left(-\frac{D_{\text{eff}}}{2}\left[\begin{array}{ccc}\hfill Q_1(0)\hfill & \hfill Q_2(0)\hfill & \hfill Q_3(0)\hfill \end{array}\right]\left[\begin{array}{ccc}\hfill {\lambda }_1^{-1}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\lambda }_2^{-1}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\lambda }_3^{-1}\hfill \end{array}\right]\left[\begin{array}{c}\hfill Q_1(0)\hfill \\ \hfill Q_2(0)\hfill \\ \hfill Q_3(0)\hfill \end{array}\right]\right)\\ \hfill =\text{exp}\left(-\frac{D_{\text{eff}}}{2}{Q}_1^2(0){\lambda }_1^{-1}\right)\text{exp}\left(-\frac{D_{\text{eff}}}{2}{Q}_2^2(0){\lambda }_2^{-1}\right)\text{exp}\left(-\frac{D_{\text{eff}}}{2}{Q}_3^2(0){\lambda }_3^{-1}\right)\end{array}(20)$

with

$Q_i(0)=\gamma {\int }_0^{t_f}\mathrm{d}{t}^{\prime }{e}^{-{\lambda }_i{t}^{\prime }}{G}_i^{\prime }(t^{\prime }),i=\mathrm{1,2,3}.(21)$

In the case when the jth eigenvalue λ_j is small, so is Q_j (0) due to the gradient echo condition, and the evaluation of the corresponding exponent on the right hand side of (Eq. 20) is numerically difficult. In this case, one can make use of the Taylor expansion of the exponential in (Eq. 21) yielding

$\text{exp}\left(-\frac{D_{\text{eff}}}{2}{Q}_j^2(0){\lambda }_j^{-1}\right)\approx \text{exp}\left(-\frac{D_{\text{eff}}}{2}\left[{\lambda }_j{\alpha }_j^2-{\lambda }_j^2{\alpha }_j{\beta }_j+{\lambda }_j^3\left(\frac{1}{4}{\beta }_j^2+\frac{1}{3}{\alpha }_j{\delta }_j\right)\right]\right),(22)$

where

$\begin{array}{ll}\hfill {\alpha }_j& =\gamma {\int }_0^{t_f}\mathrm{d}\tau G_j^{\prime }(\tau )\tau \hfill \end{array}(\mathrm{23a})$

$\begin{array}{ll}\hfill {\beta }_j& =\gamma {\int }_0^{t_f}\mathrm{d}\tau G_j^{\prime }(\tau ){\tau }^2.\hfill \end{array}(\mathrm{23b})$

$\begin{array}{ll}\hfill {\delta }_j& =\gamma {\int }_0^{t_f}\mathrm{d}\tau G_j^{\prime }(\tau ){\tau }^3.\hfill \end{array}(\mathrm{23c})$