4408

General tools for diffusion tensor distributions, and matrix-variate Gamma approximation for multidimensional diffusion MRI data

Alexis Reymbaut^1,2
¹Physical Chemistry, Lund University, Lund, Sweden, ²Random Walk Imaging AB, Lund, Sweden

Synopsis

Either on the voxel scale or within sub-voxel diffusion compartments, tissue microstructure can be described using a diffusion tensor distribution $$$\mathcal{P}(\mathbf{D})$$$. One way to resolve microstructural heterogeneity relies on choosing a plausible parametric functional form to approximate $$$\mathcal{P}(\mathbf{D})$$$. However, such a high-dimensional mathematical object is usually intractable. Here, we define matrix moments enabling the computation of diffusion metrics for any arbitrary functional choice approximating $$$\mathcal{P}(\mathbf{D})$$$. Applying these general tools to the matrix-variate Gamma distribution on the voxel scale, we obtain a new signal representation, the matrix-variate Gamma approximation, that we validate in vivo and in silico.

Introduction

The millimeter-scale resolution of diffusion MRI encompasses multiple microstructural tissue types, which implies that the diffusion signal probes a collection of microscopic diffusion profiles. Mathematically, this collection is described by a diffusion tensor distribution (DTD) $$$\mathcal{P}(\mathbf{D})$$$,¹ yielding the diffusion signal$$\frac{\mathcal{S}(\mathbf{b})}{\mathcal{S}_0}=\int_{\mathrm{Sym}^+(3)}\!\mathcal{P}(\mathbf{D})\,\exp(-\mathbf{b}:\mathbf{D})\,\mathrm{d}\mathbf{D}=\left\langle\exp(-\mathbf{b}:\mathbf{D})\right\rangle\,,$$where $$$\mathbf{b}$$$ is the encoding tensor,^2-4 $$$\mathrm{Sym}^+(3)$$$ denotes the space of 3x3 symmetric positive-definite matrices, ":" is the Frobenius inner product, and $$$\langle\cdot\rangle$$$ is the voxel-scale average. While nonparametric estimation of the DTD's features may be achievable,^5-7 an alternative approach consists in approximating $$$\mathcal{P}(\mathbf{D})$$$ with a plausible parametric functional form whose parameters can be fitted against the acquired signals. This second approach, first implemented with distributions of scalar diffusivities,^8-10 has recently been extended to DTDs in two ways:

the covariance tensor approximation,¹¹ a two-term cumulant approach equivalent to considering a voxel-scale matrix-variate Gaussian (mv-Gaussian) DTD.¹²
the DIAMOND model,^13-15 a compartmental model wherein each anisotropic sub-voxel compartment is described by a matrix-variate Gamma (mv-Gamma) DTD.¹²

The former approach considers a distribution that, unlike the mv-Gamma distribution, allows for unphysical negative-definite diffusion tensors. However, the latter approach relies on assumptions regarding the voxel content.

In this work, we first derive general tools facilitating the implementation of any functional choice for $$$\mathcal{P}(\mathbf{D})$$$. We then apply these tools to the mv-Gamma distribution and develop a new signal representation that describes the voxel content with a single mv-Gamma DTD: the "matrix-variate Gamma approximation". Finally, we validate this approximation in vivo and in silico.

Theory - Matrix moments of diffusion tensor distributions

We denote by "$$$\cdot$$$" and "$$$\otimes$$$" the matrix product and outer tensor product, respectively, with $$$\mathbf{D}^{\otimes\,\!2}\equiv\mathbf{D}\otimes\mathbf{D}$$$. Considering an arbitrary $$$\mathcal{P}(\mathbf{D})$$$, its moment-generating function writes¹²$$\mathcal{M}(\mathbf{Z})=\int\!\mathcal{P}(\mathbf{D})\,\exp\!\left[\mathrm{Tr}(\mathbf{Z}\cdot\mathbf{D})\right]\,\mathrm{d}\mathbf{D}=\left\langle\exp\!\left[\mathrm{Tr}(\mathbf{Z}\cdot\mathbf{D})\right]\right\rangle\,,$$for $$$\mathbf{Z}$$$ such that $$$\mathcal{M}(\mathbf{Z})$$$ is well-defined. $$$\mathcal{M}(\mathbf{Z})$$$ relates to the diffusion signal as¹⁵$$\mathcal{S}(\mathbf{b})=\mathcal{S}_0\,\mathcal{M}(-\mathbf{b})\,.$$By analogy with the moment-generating function of any scalar distribution, we showed that the first and second matrix derivatives of $$$\mathcal{M}(\mathbf{Z})$$$ yield the mean diffusion tensor $$$\langle\mathbf{D}\rangle$$$ and the fourth-order covariance tensor $$$\mathbb{C}=\langle\mathbf{D}^{\otimes\,\!2}\rangle-\langle\mathbf{D}\rangle^{\otimes\,\!2}$$$ associated to $$$\mathcal{P}(\mathbf{D})$$$. Indeed, we computed these derivatives using previous works^16-18 as $$$\partial\mathcal{M}/\partial\mathbf{Z}=\left\langle\mathrm{exp}\!\left[\mathrm{Tr}(\mathbf{Z}\cdot\mathbf{D})\right]\,\mathbf{D}\right\rangle$$$ and $$$\partial^2\mathcal{M}/\partial\mathbf{Z}^2=\left\langle\mathrm{exp}\!\left[\mathrm{Tr}(\mathbf{Z}\cdot\mathbf{D})\right]\,\mathbf{D}\otimes\mathbf{D}\right\rangle$$$, deducing$$\langle\mathbf{D}\rangle=\left.\frac{\partial\mathcal{M}}{\partial\mathbf{Z}}\right|_{\mathbf{Z}=\mathbf{0}}\qquad\mathrm{and}\qquad\mathbb{C}=\left.\frac{\partial^2\mathcal{M}}{\partial\mathbf{Z}^2}\right|_{\mathbf{Z}=\mathbf{0}}-\left.\frac{\partial\mathcal{M}}{\partial\mathbf{Z}}\right|_{\mathbf{Z}=\mathbf{0}}^{\otimes\!\,2}\,.$$Additional efforts would yield the sixth-order skewness tensor¹¹ $$$\mathbb{S}=\langle\mathbf{D}^{\otimes\,\!3}\rangle-3\langle\mathbf{D}\rangle\otimes\langle\mathbf{D}^{\otimes\!\,2}\rangle+2\langle\mathbf{D}\rangle^{\otimes\!\,3}$$$ via $$$\partial^3\mathcal{M}/\partial\mathbf{Z}^3$$$. While measures of voxel-scale diffusivity and anisotropy are extracted from $$$\langle\mathbf{D}\rangle$$$, measures of sub-voxel anisotropy and sub-voxel variance in isotropic diffusivities are obtained from $$$\mathbb{C}$$$ according to previous work.¹¹

Validation: We used the above formulae to retrieve the known mean tensor $$$\mathbf{M}$$$ and covariance tensor $$$\mathbf{\Sigma}\otimes\mathbf{\Psi}$$$ of the mv-Gaussian distribution¹²$$\mathcal{P}_\mathrm{Gauss}(\mathbf{X})=\frac{1}{(2\pi)^{9/2}\,\mathrm{Det}(\mathbf{\Sigma})^{3/2}\,\mathrm{Det}(\mathbf{\Psi})^{3/2}}\,\exp\left[\mathrm{Tr}\left(-\frac{1}{2}\,\mathbf{\Sigma}^{-1}\cdot(\mathbf{X}-\mathbf{M})\cdot\mathbf{\Psi}^{-1}\cdot(\mathbf{X}-\mathbf{M})^\mathrm{T}\right)\right]\,,$$for $$$\mathbf{X},\mathbf{M}$$$ arbitrary 3x3 real matrices, and $$$\mathbf{\Sigma},\mathbf{\Psi}\in\mathrm{Sym}^+(3)$$$, given its moment-generating function¹²$$\mathcal{M}_\mathrm{Gauss}(\mathbf{Z})=\exp\left[\mathrm{Tr}\left(\mathbf{Z}\cdot\mathbf{M}+\frac{1}{2}\,\mathbf{Z}^\mathrm{T}\cdot\mathbf{\Sigma}\cdot\mathbf{Z}\cdot\mathbf{\Psi}\right)\right]\,.$$

Theory - Matrix-variate Gamma approximation

We considered the non-central mv-Gamma distribution¹² used in DIAMOND:^14,15$$\mathcal{P}_\Gamma(\mathbf{D})=\frac{\mathrm{Det}(\mathbf{D})^{\kappa-2}}{\mathrm{Det}(\mathbf{\Psi})^\kappa\,\Gamma_3(\kappa)}\,\exp\!\left[-\mathrm{Tr}(\mathbf{\Theta}+\mathbf{\Psi}^{-1}\cdot\mathbf{D})\right]\,\mathcal{F}_{0,1}(\kappa,\mathbf{\Theta}\cdot\mathbf{\Psi}^{-1}\cdot\mathbf{D})\,,$$where $$$\mathbf{D}\in\mathrm{Sym}^+(3)$$$, $$$\kappa\in\;]1,+\infty[$$$, $$$\mathbf{\Psi}\in\mathrm{Sym}^+(3)$$$, and $$$\mathbf{\Theta}\in \mathrm{Sym}(3)$$$, with $$$\mathrm{Sym}(3)$$$ denoting the space of 3x3 symmetric matrices. $$$\Gamma_3$$$ is the multivariate Gamma function and $$$\mathcal{F}_{0,1}$$$ is the hypergeometric function of matrix argument of order $$$(0,1)$$$. The moment-generating function $$$\mathcal{M}_\Gamma(\mathbf{Z})$$$ of $$$\mathcal{P}_\Gamma$$$ is defined as^12,15$$\mathcal{M}_\Gamma(\mathbf{Z})=\left[\mathrm{Det}(\mathbf{I}_3-\mathbf{Z}\cdot\mathbf{\Psi})\right]^{-\kappa}\exp\!\left[\mathrm{Tr}\!\left(\left[(\mathbf{I}_3-\mathbf{Z}\cdot\mathbf{\Psi})^{-1}-\mathbf{I}_3\right]\cdot\mathbf{\Theta}\right)\right]$$for $$$\mathbf{Z}\in\mathrm{Sym}(3)$$$ such that $$$(\mathbf{I}_3-\mathbf{Z}\cdot\mathbf{\Psi})\in\mathrm{Sym}^+(3)$$$, with the 3x3 identity matrix $$$\mathbf{I}_3$$$. The diffusion signal associated to $$$\mathcal{P}_\Gamma$$$ is given by¹⁵$$\mathcal{S}_\Gamma(\mathbf{b})=\mathcal{S}_0\,\mathcal{M}_\Gamma(-\mathbf{b})=\mathcal{S}_0\,[\mathrm{Det}(\mathbf{I}_3+\mathbf{\Psi}\cdot\mathbf{b})]^{-\kappa}\exp\!\left[-\mathbf{b}:[(\mathbf{I}_3+\mathbf{\Psi}\cdot\mathbf{b})^{-1}\cdot\mathbf{\Psi}\cdot\mathbf{\Theta}]\right]\,.$$Using our expressions for the matrix moments, we retrieved the mean diffusion tensor of $$$\mathcal{P}_\Gamma$$$ (employed in DIAMOND),^14,15 and computed for the first time its covariance tensor, yielding$$\langle\mathbf{D}\rangle=\mathbf{\Psi}\cdot[\kappa\mathbf{I}_3+\mathbf{\Theta}]\qquad\mathrm{and}\qquad\mathbb{C}=\langle\mathbf{D}\rangle\otimes\mathbf{\Psi}+\mathbf{\Psi}\otimes\langle\mathbf{D}\rangle-\kappa\,\mathbf{\Psi}\otimes\mathbf{\Psi}\,.$$ The expression of $$$\langle\mathbf{D}\rangle$$$ and the symmetry of the different tensors imply that $$$\mathbf{\Psi}$$$ and $$$\mathbf{\Theta}$$$ share the same eigenvectors. Therefore, fitting $$$\mathcal{S}_\Gamma$$$ requires 11 parameters: $$$\mathcal{S}_0$$$, $$$\kappa$$$, the eigenvalues for $$$\mathbf{\Psi}$$$ and $$$\mathbf{\Theta}$$$, and the Euler angles orienting $$$\mathbf{\Psi}$$$'s eigenvectors. If numerically non positive-semidefinite, $$$\mathbb{C}$$$ is substituted with the nearest symmetric positive-semidefinite tensor.¹⁹

Validation of the matrix-variate Gamma approximation

Statistical descriptors of interest: The mean diffusivity $$$\mathrm{E}[D_\mathrm{iso}]$$$, variance of isotropic diffusivities $$$\mathrm{V}[D_\mathrm{iso}]$$$, and normalized mean anisotropy $$$\tilde{\mathrm{E}}[D_\mathrm{aniso}^2]=\mathrm{E}[(D_\mathrm{iso}D_\Delta)^2]/\mathrm{E}[D_\mathrm{iso}]^2$$$ (similar to the microscopic $$$\mathrm{FA}$$$),¹⁰ with the isotropic diffusivity $$$D_\mathrm{iso}$$$, normalized diffusion anisotropy $$$D_\Delta\in[-0.5,1]$$$, and the voxel-scale expectation and variance $$$\mathrm{E}[\cdot]$$$ and $$$\mathrm{V}[\cdot]$$$, respectively.

In vivo human brain: After implementing an iterative fitting of $$$\mathcal{S}_\Gamma$$$ in Matlab,²⁰ we applied it to a dataset available online,²¹ which signal-to-noise ratio (SNR) was around 30 in the corona radiata.²² Fig.1 presents the acquisition scheme and shows that mv-Gamma reasonably fits the signals arising from diverse voxel contents. Fig.2 features consistent estimated parameter maps of $$$\mathcal{S}_0$$$, $$$\mathrm{E}[D_\mathrm{iso}]$$$, $$$\mathrm{V}[D_\mathrm{iso}]$$$ and $$$\mathrm{FA}$$$. However, $$$\tilde{\mathrm{E}}[D_\mathrm{aniso}^2]$$$ vanishes in areas of crossing fibers.

In silico: Various voxel contents were simulated in terms of distributions of $$$D_\mathrm{iso}$$$, $$$D_\Delta$$$ and orientations, with identical in vivo acquisition scheme and SNR. Medians and interquartile ranges on estimating $$$\mathrm{E}[D_\mathrm{iso}]$$$, $$$\tilde{\mathrm{E}}[D_\mathrm{aniso}^2]$$$ and $$$\mathrm{V}[D_\mathrm{iso}]$$$ were assessed across 100 Rician-noise realizations of the covariance tensor approximation (Cov)¹¹ and the matrix-variate Gamma approximation (MV-Gamma), and compared to the simulated ground truth in isotropic (Fig.3) and anisotropic (Fig.4) systems. MV-Gamma adequately captures voxel-scale anisotropy and shows superior performance to Cov in accurately capturing $$$\mathrm{V}[D_\mathrm{iso}]$$$, but is inferior in estimating $$$\mathrm{E}[D_\mathrm{iso}]$$$. Besides, it cannot capture sub-voxel anisotropy ($$$\tilde{\mathrm{E}}[D_\mathrm{aniso}^2]$$$) in anisotropic systems with low orientational order, thereby confirming the previous in vivo observations.

Conclusion

We defined matrix moments enabling the computation of diffusion metrics for any functional choice approximating a diffusion tensor distribution $$$\mathcal{P}(\mathbf{D})$$$. Applying these tools to the mv-Gamma distribution up to its second matrix moment yields a new signal representation, the matrix-variate Gamma approximation, that performs adequately in orientationally ordered anisotropic systems (further justifying its use within DIAMOND),^14,15 but misses sub-voxel anisotropy in systems with low orientational order. This missing information is likely contained in the sixth-order skewness tensor, which can be obtained via a third matrix moment.

Acknowledgements

This work was financially supported by the Swedish Foundation for Strategic Research (AM13-0090, ITM17-0267) and the Swedish Research Council (2018-03697).

References

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

2. Eriksson S et al., 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 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. de Almeida Martins J P and Topgaard D, Two-dimensional correlation of isotropic and directional diffusion using NMR. Physics Review Letters 2016;116:087601.

6. 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. Scientific Reports 2018;8,2488.

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

8. Yablonskiy D A et al., Statistical model for diffusion attenuated MR signal, Magnetic Resonance in Medicine 2003;50:664.

9. Jbabdi S et al., Model-based analysis of multishell diffusion MR data for tractography: How to get over fitting problems. Magnetic Resonance in Medicine 2012;68:1846.

10. Lasic 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.

11. Westin C-F et al. Q-space trajectory imaging for multidimensional diffusion MRI of the human brain. NeuroImage 2016;135:345-362.

12. Guptar A K and Nagar D K, Matrix Variate Distributions. Edited by Chapman and Hall/CRC, Boca Raton, Florida, 2000.

13. Scherrer B et al., Characterizing brain tissue by assessment of the distribution of anisotropic microstructural environments in diffusion-compartment imaging (DIAMOND). Magnetic Resonance in Medicine 2016;76:963.

14. Scherrer B et al., Decoupling axial and radial tissue heterogeneity in diffusion compartment imaging. In: Information Processing in Medical Imaging (IPMI) 2017;440-452.

15. Reymbaut A 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.

16. Brewer J, Kronecker products and matrix calculus in system theory. IEEE Transactions on Circuits and Systems 1978;25:772.

17. Laue S et al., Computing higher order derivatives of matrix and tensor expressions. In: Advances in Neural Information Processing Systems (NIPS) 2018.

18. Laue S et al., Matrix calculus, www.matrixcalculus.org. Accessed: 2019-07-30.

19. Higham N J, Computing a nearest symmetric positive semidefinite matrix. Linear algebra and its applications 1988;103:103-118.

20. Nilsson M 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.

21. Szczepankiewicz F et al., 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.

22. Szczepankiewicz F et al., Tensor-valued diffusion encoding for diffusional variance decomposition (DIVIDE): Technical feasibility in clinical MRI systems. PLoS ONE 2019;14(3):e0214238.

Figures

Figure 1. Left: Acquisition scheme ($$$\mathbf{b}$$$'s trace $$$b$$$, normalized anisotropy $$$b_\Delta$$$ and orientation $$$(\Theta,\Phi)$$$) with sorted acquisition number $$$n_\mathrm{acq}$$$, and measured (black) and fitted (color) normalized signals $$$\tilde{S}$$$ in selected voxels (cerebrospinal fluid in blue, grey matter in green, white matter in red, and crossing area in orange). Right: $$$\mathcal{S}_0$$$ map indicating the location of the selected voxels.

Figure 2. Left: Estimated axial parameter maps of $$$\mathcal{S}_0$$$ (compared to the measured $$$\mathcal{S}_{0,\mathrm{measured}}$$$), $$$\mathrm{E}[D_\mathrm{iso}]$$$ and $$$\mathrm{V}[D_\mathrm{iso}]$$$. Grey scale codes for magnitude. Right: Estimated axial parameter maps of $$$\mathrm{FA}$$$ and $$$\tilde{E}[D_\mathrm{aniso}^2]$$$. Color codes for orientations and opacity codes for magnitude.

Figure 3. Characterization of isotropic systems with varying $$$\mathrm{V}[D_\mathrm{iso}]$$$ at low $$$\mathrm{E}[D_\mathrm{iso}]$$$ (top), and at high $$$\mathrm{E}[D_\mathrm{iso}]$$$ (bottom). Each system consists of two identical Gaussian distributions of $$$D_\mathrm{iso}$$$. While the left panels feature both SNR=30 and infinite-SNR results, the right panels illustrate the voxel contents for both cases.

Figure 4. Characterization of anisotropic components with constant orientation and varying $$$\tilde{E}[D_\mathrm{aniso}^2]$$$ (top), and with varying orientational order and constant $$$\tilde{E}[D_\mathrm{aniso}^2]$$$ (bottom). The order parameter $$$\mathrm{OP}=\mathrm{E}[P_2(\cos\theta)]$$$ (with second Legendre polynomial $$$P_2$$$) is set by a Watson distribution of these identical components (Gaussian distributions of $$$D_\Delta$$$ with constant $$$\mathrm{E}[D_\mathrm{iso}]$$$). Other conventions are identical to those of Figure 3.

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

4408