4410

Trueness and precision of statistical descriptors obtained from multidimensional diffusion signal inversion algorithms
Alexis Reymbaut1,2, Paolo Mezzani1,3, João Pedro de Almeida Martins1,2, and Daniel Topgaard1,2
1Physical Chemistry, Lund University, Lund, Sweden, 2Random Walk Imaging AB, Lund, Sweden, 3Physics, Università degli Studi di Milano, Milan, Italy

Synopsis

In first approximation, the diffusion signal writes as the Laplace transform of an intra-voxel diffusion tensor distribution (DTD). Several algorithms have been introduced to estimate the DTD’s statistical descriptors (mean diffusivity, variance of isotropic diffusivities, mean squared diffusion anisotropy, etc.) by inverting data obtained from tensor-valued diffusion encoding schemes. However, the trueness and precision of these estimations have not been systematically assessed and compared across methods. Here, we compare such estimations in silico for a 1D Gamma fit, a generalized two-term cumulant approach, and 2D and 4D Monte-Carlo inversion techniques, using a common and clinically feasible tensor-valued acquisition scheme.

Introduction

In biological tissues, typical MRI voxels comprise multiple microscopic environments, whose local organization can be captured by microscopic diffusion tensors. The measured diffusion MRI signal can, therefore, be written as the multidimensional Laplace transform of an intra-voxel diffusion tensor distribution (DTD).1 While tensor-valued diffusion encoding schemes have been designed to probe specific features of the DTD via the sizes (b-value $$$b$$$), shapes (normalized anisotropy $$$b_\Delta\in[-0.5,1]$$$) and orientations $$$(\Theta,\Phi)$$$ of encoding tensors $$$\mathbf{b}$$$,2-4 several algorithms have been introduced to invert such datasets and estimate statistical descriptors of the DTD. However, the trueness and precision of these estimations have not been systematically assessed and compared across inversion methods.

Methods

We considered the estimations of three statistical descriptors:
  • the mean diffusivity $$$\mathrm{E}[D_\mathrm{iso}]$$$,
  • the mean squared diffusion anisotropy $$$\tilde{\mathrm{E}}[D^2_\mathrm{aniso}]=\mathrm{E}[(D_\mathrm{iso}D_\Delta)^2]\,/\,\mathrm{E}[D_\mathrm{iso}]^2$$$,
  • the variance of isotropic diffusivities $$$\mathrm{V}[D_\mathrm{iso}]$$$,
where $$$D_\mathrm{iso}$$$ and $$$D_\Delta\in[-0.5,1]$$$ are the isotropic diffusivity and the normalized diffusion anisotropy, and $$$\mathrm{E}[\,\cdot\,]$$$ and $$$\mathrm{V}[\,\cdot\,]$$$ denote the voxel-scale average and variance, respectively. We chose an acquisition scheme similar to that of previous in vivo studies of brain tissue microstructure,5 that can be readily implemented in a standard clinical MRI scanner as a 6-7 minutes protocol comprising 30 axial slices:
  • 53 linearly encoded signals ($$$b_\Delta=1$$$) distributed as $$$6,10,16$$$ and $$$21$$$ points at $$$b=100,700,1400$$$ and $$$2000$$$ s/mm2, respectively.
  • 32 spherically encoded signals ($$$b_\Delta=0$$$) sampled as $$$6,6,10$$$ and $$$10$$$ points at $$$b=100,700,1400$$$ and $$$2000$$$ s/mm2, respectively.
  • one $$$b=0$$$ signal.
The linear points were measured for orientation sets derived from platonic solids,6 with a sufficiently high number of directions per $$$b$$$-shell to yield a rotationally invariant powder-averaged (orientationally averaged) signal.5

We compared the performances of four signal inversion techniques:
  • for powder-averaged signals $$$\overline{\mathcal{S}}(b,b_\Delta)$$$, the Gamma fitting (Gamma)7 and the 2D Monte-Carlo inversion (MC-2D).8
  • for non powder-averaged signals $$$\mathcal{S}(\mathbf{b})=\mathcal{S}(b,b_\Delta,\Theta,\Phi)$$$, the covariance tensor approximation (Cov, a tensorial two-term cumulant expansion)6 and the 4D Monte-Carlo inversion (MC-4D).9
While Gamma fits the diffusion signal via the Laplace transform of the 1D Gamma distribution$$\frac{\overline{\mathcal{S}}(b,b_\Delta)}{\mathcal{S}_0}=\left(1+b\,\frac{\mu_2(b_\Delta)}{\overline{\mathit{D}}}\right)^{-\overline{\mathit{D}}^{\,2}/\mu_2(b_\Delta)}\,,$$where $$$\overline{\mathit{D}}$$$ is the average diffusivity and $$$\mu_2(b_\Delta)=\mathrm{V}[D_\mathrm{iso}]+(4b_\Delta^2/5)\times\mathrm{E}[D_\mathrm{iso}]^2\,\tilde{\mathrm{E}}[D^2_\mathrm{aniso}]$$$ is the second central moment of the underlying diffusivity distribution, Cov writes the signal as the two-term cumulant expansion$$\frac{\mathcal{S}(\mathbf{b})}{\mathcal{S}_0}\underset{\mathbf{b}:\langle\mathbf{D}\rangle\to0}{\simeq}\exp\!\left(-\mathbf{b}:\langle \mathbf{D}\rangle+\frac{1}{2}\,\mathbf{b}^{\otimes 2}:\mathbb{C}\right)\,,$$where $$$\langle\mathbf{D}\rangle$$$ is the average diffusion tensor, $$$\mathbb{C}=\langle\mathbf{D}^{\otimes\,\!2}\rangle-\langle\mathbf{D}\rangle^{\otimes\,\!2}$$$ is the fourth-order covariance tensor with the outer tensor product $$$\mathbf{D}^{\otimes\,\!2}=\mathbf{D}\otimes\mathbf{D}$$$, and ":" is the Frobenius inner product. Alternatively, MC-2D considers the general powder-averaged signal3$$\frac{\overline{\mathcal{S}}(b,b_\Delta)}{\mathcal{S}_0}=\sqrt{\frac{\pi}{2}}\int_{-0.5}^1\int_0^{+\infty}\!\frac{\mathrm{erf}(\sqrt{3bb_\Delta\,\!D_\mathrm{iso}D_\Delta})}{\sqrt{3bb_\Delta\,\!D_\mathrm{iso}D_\Delta}}\,\mathrm{exp}(-bD_\mathrm{iso}[1-b_\Delta\,\!D_\Delta])\,\mathcal{P}(D_\mathrm{iso},D_\Delta)\,\mathrm{d}D_\mathrm{iso}\,\mathrm{d}D_\Delta$$and MC-4D considers the general signal$$\frac{\mathcal{S}(b,b_\Delta,\Theta,\Phi)}{\mathcal{S}_0}=\int_0^{2\pi}\int_0^\pi\int_{-0.5}^1\int_0^{+\infty}\!\mathrm{exp}(-bD_\mathrm{iso}[1+2b_\Delta\,\!D_\Delta\,P_2(\cos\beta)])\,\mathcal{P}(D_\mathrm{iso},D_\Delta,\theta,\phi)\,\mathrm{d}D_\mathrm{iso}\,\mathrm{d}D_\Delta\,\mathrm{d}\theta\,\mathrm{d}\phi,$$where $$$P_2(x)=(3x^2-1)/2$$$ is the second Legendre polynomial, $$$\mathcal{P}(\,\cdot\,)$$$ is the associated 2D or 4D DTD, and $$$\beta$$$ is the shortest angle between the diffusion tensor orientation $$$(\theta,\phi)$$$ and the encoding tensor orientation $$$(\Theta,\Phi)$$$. In practice, these last two equations are discretized as $$$\mathbf{S}=\mathbf{K}\cdot\mathbf{w}$$$, where $$$\mathbf{S}$$$ is the column vector containing the signals, $$$\mathbf{K}$$$ is the inversion kernel matrix and $$$\mathbf{w}$$$ is the column vector containing the weights $$$w$$$ of any given configuration $$$\{D_\mathrm{iso},D_\Delta,\theta,\phi\}$$$. The Monte-Carlo approaches randomly sample such configurations and weigh their propensity to fit the acquired signals via non-negative least-squares fitting.

To ensure fair comparison of the different inversion techniques, all readily available as Matlab code on GitHub,10,11 we follow these steps:
  1. A system is simulated by generating a set of ground-truth features $$$\{D_\mathrm{iso},D_\Delta,\theta,\phi,w\}$$$, from which the ground-truth descriptors and the ground-truth signals are computed using the general 4D kernel.
  2. Each signal inversion technique is run on an identical set of signals with added Rician noise at either signal-to-noise ratio (SNR) of 30 or at infinite SNR.
  3. Step 2 is repeated 100 times to build up statistics (medians and interquartile ranges) on estimation of the statistical descriptors.
While trueness is quantified by the difference between the median and the ground truth (bias), precision is quantified by the interquartile range. An "accurate" estimation is defined as presenting both high trueness and high precision.

Results and discussion

Figures 1, 2, 3 and 4 respectively present the effect of the variance of isotropic diffusivities, anisotropy, orientational order, and fraction of isotropic and anisotropic components on the descriptors' estimations for all investigated inversion methods.

We found that these methods all share similar precision (except for a lower precision of MC-2D) but differ in term of trueness. Gamma exhibits infinite-SNR biases when the signal strongly deviates from mono-exponentiality and is unaffected by orientational order. As for Cov, it shows infinite-SNR biases when this deviation originates either from the variance in isotropic diffusivities or from the low orientational order of anisotropic diffusion components, because signals from such distributions contain features that are not captured by a two-term cumulant expansion. MC-2D shows remarkable trueness in all studied systems. However, its trueness in systems that are anisotropic at the voxel scale drastically worsens if the acquisition scheme does not possess enough directions to yield a good rotational invariance of the powder-averaged signal. Finally, MC-4D presents no infinite-SNR bias but significantly suffers from noise in the data, while preserving contrast in finite-anisotropy systems.

Conclusion

Methods that constrain the inversion (Gamma and Cov) exhibit, within the common acquisition scheme we chose, inherent limitations in estimating statistical descriptors, even at infinite SNR. Unconstrained methods (MC-2D and MC-4D) show perfect accuracy at infinite SNR but suffer from data noise, worsening either precision (MC-2D) or trueness (MC-4D) at finite SNR while preserving reasonable contrast. Future denoising efforts may alleviate these problems.

Acknowledgements

This work was financially supported by the Swedish Foundation for Strategic Research (AM13-0090, ITM17-0267) and the Swedish Research Council (2018-03697). Daniel Topgaard owns shares in Random Walk Imaging AB (Lund, Sweden, http://www.rwi.se/), holding patents related to the described methods. We warmly thank Filip Szczepankiewicz and Markus Nilsson for fruitful discussions.

References

1. Jian B, Vemuri BC, Özarslan E et al., A novel tensor distribution model for the diffusion-weighted MR signal. Neuroimage 2007;37(1):164-176.

2. Eriksson S, Lasic S and Topgaard D, Isotropic diffusion weighting in PGSE NMR by magic-angle spinning of the q-vector. Journal of Magnetic Resonance 2013;226:13-18.

3. Eriksson S , Lasic S, Nilsson M et al., NMR diffusion-encoding with axial symmetry and variable anisotropy: Distinguishing between prolate and oblate microscopic diffusion tensors with unknown orientation distribution. The Journal of Chemical Physics 2015;142(10):104201.

4. Topgaard D, Multidimensional diffusion MRI. Journal of Magnetic Resonance 2017; 275:98-113.

5. Szczepankiewicz F, Sjölund J, Ståhlberg F et al. Tensor-valued diffusion encoding for diffusional variance decomposition (DIVIDE): Technical feasibility in clinical MRI systems. PLOS ONE 2019;14(3):1-20.

6. Westin C-F, Knutsson H, Pasternak O et al. Q-space trajectory imaging for multidimensional diffusion MRI of the human brain. NeuroImage 2016;135:345-362.

7. Lasic S, Szczepankiewicz F, Eriksson S et al., Microanisotropy imaging: quantification of microscopic diffusion anisotropy and orientational order parameter by diffusion MRI with magic-angle spinning of the q-vector. Frontiers in Physics 2014;2:11.

8. de Almeida Martins J P and Topgaard D, Two-dimensional correlation of isotropic and directional diffusion using NMR. Physics Review Letters 2016;116:087601.

9. Topgaard D. Diffusion tensor distribution imaging. NMR in Biomedicine 2019;32(5):e4066

10. Nilsson M, Szczepankiewicz F, Lampinen B, et al. An open-source framework for analysis of multidimensional diffusion MRI data implemented in MATLAB. In: International Society for Magnetic Resonance Imaging (ISMRM) 2018.

11. Nilsson M, Szczepankiewicz F, C-F Westin and Topgaard D. MATLAB code for multidimensional diffusion MRI. https://github.com/daniel-topgaard/md-dmri. Accessed: 2019-07-09.

12. Reymbaut A, Gilbert G, Szczepankiewicz F et al., The “Magic DIAMOND" method: probing brain microstructure by combining b-tensor encoding and advanced diffusion compartment imaging. In: International Society for Magnetic Resonance Imaging (ISMRM) 2018.

Figures

Figure 1. Characterization of bimodal distributions of Diso (two Gaussian distributions, identical standard deviation 0.05 μm2/ms) with constant E[Diso]=0.8 μm2/ms and varying V[Diso]. Estimations are shown for 100 realizations of Gamma (orange circles), Cov (green squares), MC-2D (red triangles) and MC-4D (blue stars) with either finite (top) or infinite (bottom) SNR. Symbols/error bars: medians/interquartile ranges of the descriptors across realizations. Black dashed lines: ground truth. Bottom right panels: DTDs of the investigated voxel contents.

Figure 2. Characterization of anisotropic systems consisting of aligned anisotropic components, each of which is described by a unimodal distribution of D$$$_\Delta$$$ (Gaussian distribution) with constant Diso=0.8 μm2/ms and constant standard-deviation-to-the-mean ratio of 0.1. Even though the first and third simulated systems exhibit the same squared anisotropy, the former contains planar anisotropic components (subscript P) and the latter contains linear anisotropic components (subscript L). Symbol/color conventions are identical to those of Figure 1.

Figure 3. Characterization of anisotropic systems as a function of the orientation order parameter OP=E[P2($$$\cos\theta$$$)]. Each anisotropic component is described by a unimodal distribution of D$$$_\Delta$$$ (Gaussian distribution) with constant standard-deviation-to-the-mean ratio of 0.1. Each order parameter is fixed by orientating the anisotropic components following a given Watson distribution. These systems share common E[Diso]=0.8 μm2/ms, $$$\tilde{E}$$$[D2aniso]=0.47 and V[Diso]=0. Symbol/color conventions are identical to those of Figure 1.

Figure 4. Characterization of a composite system with varying signal fraction fiso$$$\in$$$[0,1] of a cerebrospinal-fluid-like structure (Gaussian distribution of Diso with mean 3 μm2/ms and standard deviation 0.1 μm2/ms), starting from a fiber-like structure. This structure is made of Watson-distributed anisotropic components with a reasonable internal distribution of Diso and D$$$_\Delta$$$ inspired by results of the Magic DIAMOND model,12 and an intermediate order parameter mimicking fiber dispersion. Symbol/color conventions are identical to those of Figure 1.

Proc. Intl. Soc. Mag. Reson. Med. 28 (2020)
4410