3560

Tensorial formulation allowing to verify or falsify the microstructural standard model from multidimensional diffusion MRI
Jose M. Pozo1,2, Santiago Coelho1,2, and Alejandro F. Frangi1,2

1Centre for Computational Imaging & Simulation Technologies in Biomedicine (CISTIB), School of Computing, University of Leeds, Leeds, United Kingdom, 2Leeds Institute for Cardiac and Metabolic Medicine (LICAMM), School of Medicine, University of Leeds, Leeds, United Kingdom

Synopsis

The standard model (SM) has been proposed as an appropriate diffusion MRI microstructural model for brain white matter. Using the cumulant expansion up to 2nd order in the b-tensor, it has been recently shown that multidimensional diffusion encoding makes the model parameter estimation problem well-posed. However, the tensorial properties of the expansion and their relationship with the SM assumptions have not been exploited. We reformulate the solution of the SM in an elegant tensorial form and analyse the constraints the SM imposes, showing how they can potentially falsify the model. Some simulations show the test feasibility for high-quality signals ($$$\text{SNR}\geq 50$$$).

Introduction

Biophysical models bear the promise of enabling estimation of specific microstructural information from diffusion MRI (dMRI). The standard model (SM)1 of brain white matter incorporates relevant tissue properties without imposing excessive constraints. The SM signal for a general $$$q$$$-space trajectory (multidimensional dMRI) determining a b-tensor $$$\boldsymbol{B}$$$2 is

$$S(\boldsymbol{B})=S_0\int_{\mathbb{S}^2}{\cal{}K}(\boldsymbol{B};\hat{\boldsymbol{u}})\,{\cal{}P}(\hat{\boldsymbol{u}})\mathop{\rm{}d\!}{}\!{\hat{\boldsymbol{u}}}$$

with a general fibre orientation distribution $$${\cal{}P}(\hat{\boldsymbol{u}})$$$, and fibres composed of two/three impermeable Gaussian compartments with cylindrically symmetric response signal (kernel)

$${\cal{}K}(\boldsymbol{B};\hat{\boldsymbol{u}})=f_{\text{a}}\exp\bigl(-D_{\text{a}}u_iu_jB_{ij}\bigr)+f_{\text{e}}\exp\bigl(-(D_{\text{e}}^\perp\delta_{ij}+\Delta_{\text{e}}u_iu_j)B_{ij}\bigr)+(1-f_{\text{a}}-f_{\text{e}})\exp\bigl(-D_{\text{w}}\delta_{ij}B_{ij}\bigr),$$

representing intra-axonal and extra-axonal environments (and sometimes cerebrospinal fluid), with constant diffusivities ($$$D_{\text{a}},D_{\text{e}}^\perp,\Delta_{\text{e}},D_{\text{w}}$$$) and fractions ($$$f_{\text{a}},f_{\text{e}}$$$) for all orientations $$$\hat{\boldsymbol{u}}$$$.

Multidimensional dMRI accesses information unavailable through single diffusion encoding (SDE), reflected in the cumulant expansion

$$\log{S}=\log{S_0}-D_{ij}B_{ij}+\frac12Z_{ijk\ell}B_{ij}B_{k\ell}+{\mathcal{O}}(\boldsymbol{B}^3).$$

where $$$Z_{ijk\ell}=Z_{(ij)(k\ell)}=Z_{(k\ell)(ij)}$$$ (major and minor symmetries). SDE ($$$B_{ij}=b\,n_in_j$$$) only accesses the fully symmetric part, corresponding to the kurtosis tensor:

$$Z_{ijk\ell}B_{ij}B_{k\ell}=b^2Z_{(ijk\ell)}n_in_jn_kn_\ell=\frac13\bar{D}^2W_{ijk\ell}\,n_in_jn_kn_\ell.$$

Recent work3,4 has shown that additional information in $$$\boldsymbol{Z}$$$ is able to overcome the degeneracy of the SM. However, these works have not fully exploited the tensorial properties of $$$D_{ij}$$$ and $$$Z_{ijk\ell}$$$ and their relationship with the SM assumptions. Herewith, the solution to the SM from $$$D_{ij}$$$ and $$$Z_{ijk\ell}$$$ is reformulated using an elegant tensorial formalism. Then, we analyse the constraints the SM imposes on $$$D_{ij}$$$ and $$$Z_{ijk\ell}$$$ and show how to use these to potentially invalidate the SM.

Methods

Cumulant expansion of the SM Let us consider a generalized SM for a general number of compartments:

$${\cal{}K}(\boldsymbol{B},\hat{\boldsymbol{u}})=\sum_cf_c\exp\bigl(-(D_c^\perp\delta_{ij}+\Delta_c{}u_iu_j)B_{ij}\bigr).$$

Defining $$$\hat{Z}_{ijk\ell}= Z_{ijk\ell}+D_{ij}D_{k\ell}$$$, the SM cumulant expansion can be expressed as

$$D_{ij}=\alpha\,H_{\!\scriptscriptstyle(2)}{}_{ij}+\beta\,\delta_{ij}\qquad\text{and}\qquad\hat{Z}_{ijk\ell}=\gamma\,H_{\!\scriptscriptstyle(4)}{}_{ijk\ell}+\delta\,\bigl(H_{\!\scriptscriptstyle(2)}{}_{ij}\delta_{k\ell}+\delta_{ij}H_{\!\scriptscriptstyle(2)}{}_{k\ell}\bigr)+\epsilon\,\delta_{ij}\delta_{k\ell},$$

where

$$\alpha=\sum_cf_c\Delta_c,\quad\beta=\sum_cf_cD_c^\perp,\quad\gamma=\sum_cf_c\Delta_c{}^2,\quad\delta=\sum_cf_c\Delta_cD_c^\perp,\quad\epsilon=\sum_cf_cD_c^\perp{}^2,$$

and

$$H_{\!\scriptscriptstyle(2)}{}_{ij}=\int_{\mathbb{S}^2}u_iu_j\,\mathcal{P}(\hat{\boldsymbol{u}})\mathop{\rm{}d\!}\hat{\boldsymbol{u}}\qquad\text{and}\qquad{}H_{\!\scriptscriptstyle(4)}{}_{ijk\ell}=\int_{\mathbb{S}^2}u_iu_ju_ku_\ell\,\mathcal{P}(\hat{\boldsymbol{u}})\mathop{\rm{}d\!}{}S_{\hat{\boldsymbol{u}}}.$$

Observe that $$$H_{\!\scriptscriptstyle(4)}{}_{ijkk}=H_{\!\scriptscriptstyle(2)}{}_{ij}$$$ and $$$H_{\!\scriptscriptstyle(2)}{}_{ii}=1$$$.

For the SM with 2 compartments (or 3 with $$$D_{\text{w}}$$$ known), the full set of kernel parameters (diffusivities and fractions) can be determined uniquely3,4 by the 5 variables $$$(\alpha,\beta,\gamma,\delta,\epsilon)$$$. However, even without intending to determine the kernel parameters, the 2 tensors include additional information that has not been properly analysed yet.

Decomposition into irreducible invariant tensor The diffusion tensor is decomposed into the irreducible invariant tensors

$$D_{ij}=D^{\scriptscriptstyle(2)}_{ij}+\frac13D^{\scriptscriptstyle(0)}\delta_{ij},\qquad\text{where}\qquad{}D^{\scriptscriptstyle(2)}_{ii}=0\quad\text{and}\quad{D^{\scriptscriptstyle(0)}}=D_{ii},$$

and the tensor $$$\hat{\boldsymbol{Z}}$$$ in 3D5 into

$$\hat{Z}_{ijk\ell}=S_{ijk\ell}+A_{ijk\ell}=S^{\scriptscriptstyle(4)}_{ijk\ell}+\frac67S^{\scriptscriptstyle(2)}_{(ij}\delta^{\phantom{()}}_{k\ell)}+\frac15\delta_{(ij}\delta_{k\ell)}{S^{\scriptscriptstyle(0)}}+\frac32\Bigl(A^{\scriptscriptstyle(2)}_{ij}\delta^{\phantom{()}}_{k\ell}+A^{\scriptscriptstyle(2)}_{k\ell}\delta^{\phantom{()}}_{ij}-2A^{\scriptscriptstyle(2)}_{(ij}\delta^{\phantom{()}}_{k\ell)}\Bigr)+\frac16(\delta_{ij}\delta_{k\ell}-\delta_{(ij}\delta_{k\ell)}){A^{\scriptscriptstyle(0)}}$$

where $$$S_{ijk\ell}=\hat{Z}_{(ijk\ell)}$$$ is the fully symmetric part and $$$A_{ijk\ell}$$$ is its complement, and each decomposed into the traceless tensors:

$${}S^{\scriptscriptstyle(4)}_{ijkk}=0,\qquad{}S^{\scriptscriptstyle(2)}_{ii}=0,\qquad\text{with}\qquad {}S^{\scriptscriptstyle(2)}_{ij}=S_{ijkk}-\frac13{S^{\scriptscriptstyle(0)}}\delta_{ij},\quad\text{and}\quad{S^{\scriptscriptstyle(0)}}=S_{iikk},$$

and

$${}A^{\scriptscriptstyle(2)}_{ii}=0,\qquad\text{with}\qquad{}A^{\scriptscriptstyle(2)}_{ij}=A_{ijkk}-\frac13{A^{\scriptscriptstyle(0)}}\delta_{ij},\quad\text{and}\quad{A^{\scriptscriptstyle(0)}}=A_{iikk}.$$

In 3D, $$$A^{\scriptscriptstyle(4)}_{ijk\ell}$$$ identically vanishes. This gives 7 irreducible independent tensors: $$$\boldsymbol{{D^{\scriptscriptstyle(2)}}},{D^{\scriptscriptstyle(0)}},\boldsymbol{{S^{\scriptscriptstyle(4)}}},\boldsymbol{{S^{\scriptscriptstyle(2)}}},{S^{\scriptscriptstyle(0)}},\boldsymbol{{A^{\scriptscriptstyle(2)}}},{A^{\scriptscriptstyle(0)}}$$$.

Applying the same decomposition to the fully symmetric tensor $$$H_{\!\scriptscriptstyle(4)}{}_{ijk\ell}$$$ we obtain the traceless tensors $$$H^{\scriptscriptstyle(4)}_{ijk\ell}$$$ and $$$H^{\scriptscriptstyle(2)}_{ij}$$$, satisfying $$$H^{\scriptscriptstyle(4)}_{ijkk}=H^{\scriptscriptstyle(2)}_{ii}=0$$$. They include the same information as spherical harmonics of the corresponding order.

Results

Tensorial constraints from the SM Applying the tensorial decomposition to the cumulant expansion of the SM, we get the equations

$$\hspace{4.5em}{}D^{\scriptscriptstyle(2)}_{ij}=\alpha{}H^{\scriptscriptstyle(2)}_{ij}\hspace{5em}{D^{\scriptscriptstyle(0)}}=\alpha+\beta\\{}S^{\scriptscriptstyle(4)}_{ijk\ell}=\gamma{}H^{\scriptscriptstyle(4)}_{ijk\ell}\hspace{1em}\quad{}S^{\scriptscriptstyle(2)}_{ij}=\bigl(\gamma+\tfrac73\delta\bigr)H^{\scriptscriptstyle(2)}_{ij}\hspace{1em}\quad{S^{\scriptscriptstyle(0)}}=\gamma+5\bigl(\tfrac23\delta+\epsilon\bigr)\\\hspace{6.6em}{}A^{\scriptscriptstyle(2)}_{ij}=\tfrac23\delta{}H^{\scriptscriptstyle(2)}_{ij}\hspace{4.6em}{A^{\scriptscriptstyle(0)}}=4\bigl(\tfrac23\delta+\epsilon\bigr).$$

From this system, we can easily determine the orientation distribution tensors:

$$H^{\scriptscriptstyle(2)}_{ij}=\frac1\gamma\bigl({}S^{\scriptscriptstyle(2)}_{ij}-\tfrac72{}A^{\scriptscriptstyle(2)}_{ij}\bigr),\qquad\text{and}\qquad{}H^{\scriptscriptstyle(4)}_{ijk\ell}=\frac1\gamma{}S^{\scriptscriptstyle(4)}_{ijk\ell},\qquad\text{with}\qquad\gamma={S^{\scriptscriptstyle(0)}}-\tfrac54{A^{\scriptscriptstyle(0)}}.$$

One can sequentially determine the remaining variables $$$\alpha,\beta,\delta,\epsilon$$$, except for $$$H^{\scriptscriptstyle(2)}_{ij}=0$$$, which includes the isotropic case.

For our purpose, the most important observation is that 3 of the tensors, $$$D^{\scriptscriptstyle(2)}_{ij},S^{\scriptscriptstyle(2)}_{ij},A^{\scriptscriptstyle(2)}_{ij}$$$, are proportional between themselves and to $$$H^{\scriptscriptstyle(2)}_{ij}$$$. Thus, the minimum of the residuals

$$R_S=\lVert\boldsymbol{{S^{\scriptscriptstyle(2)}}}-\lambda_S\boldsymbol{D^{\scriptscriptstyle(2)}}\rVert \qquad\text{and}\qquad R_A=\lVert\boldsymbol{{A^{\scriptscriptstyle(2)}}}-\lambda_A\boldsymbol{D^{\scriptscriptstyle(2)}}\rVert,$$

which happens for

$$\lambda_S=D^{\scriptscriptstyle(2)}_{ij}S^{\scriptscriptstyle(2)}_{ij}/\bigl(D^{\scriptscriptstyle(2)}_{k\ell}D^{\scriptscriptstyle(2)}_{k\ell}\bigr)\qquad\text{and}\qquad \lambda_A=D^{\scriptscriptstyle(2)}_{ij}A^{\scriptscriptstyle(2)}_{ij}/\bigl(D^{\scriptscriptstyle(2)}_{k\ell}D^{\scriptscriptstyle(2)}_{k\ell}\bigr),$$

should be zero for noiseless dMRI signals. Hence their value should be comparable to the noise level of the dMRI acquisition.

Simulation test We simulate dMRI from a combination of linear and planar diffusion encoding3 with 3 shells ($$$b=0,1,2\,\text{ms}/\mu\text{m}^2$$$) each with 90 directions, and with a series of signal-to-noise ratios ($$$\text{SNR}=25, 33, 50, 100, 200$$$) at $$$b_0$$$. The signal has been generated for two cases representing fibre crossing:

  • SM: Kernel, $$$f_{\text{a}}=0.3, D_{\text{a}}=2.5, D_{\text{e}}^\perp=1,\Delta_{\text{e}}=1$$$ [$$$\mu\text{m}^2/\text{ms}$$$], and two Watson orientation distributions with equal weight and concentrations, $$$\kappa_1=\kappa_2=4$$$, but perpendicular orientations, $$$\boldsymbol{\mu}_1,\boldsymbol{\mu}_2$$$.
  • Non-SM: Combining two SM with different orientation distribution (same 2 Watson components) but each bundle with a different kernel: $$$f_{\text{a},1}=0.3,f_{\text{a},2}=0.2,D_{\text{a},1}=2.5,D_{\text{a},2}=2.8,D_{\text{e},1}^\perp=1,\Delta_{\text{e},1}=1,D_{\text{e},2}^\perp=1.2,\Delta_{\text{e},2}=1.5$$$ [$$$\mu\text{m}^2/\text{ms}$$$].

The resulting distributions of residuals for each of the cases is presented in Fig. 1. For $$$\text{SNR}\geq 50$$$, the SM and non-SM cases can be reasonably distinguished.

Discussion and Conclusions

From the cumulant expansion of the SM and the irreducible invariant decomposition of the resulting tensors $$$\boldsymbol{D}$$$ and $$$\boldsymbol{Z}$$$ we have presented an elegant tensorial form relating them to the SM parameters. This allows a simple expression of the inverse problem solution and evidencing that the tensors include redundancies that can be used as a test to falsify or verify the SM. We have shown in a simulated experiment the viability of this test, in distinguishing the SM from the non-SM, with achievable, although high, SNRs.

The simulated non-SM case can be interpreted as crossing fibres where each bundle has different microstructural properties, which is a reasonable possibility. We can also go out of the SM by considering non-cylindrically symmetric kernels or non-Gaussian compartments.

Finally, although we have considered multidimensional dMRI with general B-tensor, one of the tests, involving the residual with the symmetric part of $$$\boldsymbol{Z}$$$, can be applied with SDE.

Acknowledgements

This work has been supported by the OCEAN project (EP/M006328/1) and MedIAN Network (EP/N026993/1) both funded by the Engineering and Physical Sciences Research Council (EPSRC) and the European Commission FP7 Project VPH-DARE@IT (FP7-ICT-2011-9-601055).

References

1. Novikov D S, Veraart J, Jelescu I O, Fieremans E. Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion MRI. NeuroImage. 2018; 174:518-538.

2. Westin C-F, Knutsson H, Pasternak O, Szczepankiewicz F, Özarslan E, Van Westen D, Mattisson C, Bogren M, O'Donnell L J, Kubicki M, Topgaard D, Nilsson M. q-space trajectory imaging for multidimensional diffusion MRI of the human brain. NeuroImage. 2016; 135:345-362.

3. Coelho S, Pozo J M, Jespersen S N, Jones D K, Frangi A F. Double diffusion encoding prevents degeneracy in parameter estimation of biophysical models in diffusion MRI. Preprint arXiv. 2018; 1809.05059v1.

4. Reisert M, Kiselev V G, Dhital B. A unique analytical solution of the white matter standard model using linear and planar encodings. Preprint arXiv. 2018; 1808.04389v1.

5. Itin Y and Hehl F W. Irreducible decompositions of the elasticity tensor under the linear and orthogonal groups and their physical consequences. In Journal of Physics: Conference Series, 2015, volume 597.

Figures

Distribution of the residuals $$$R_S$$$ and $$$R_A$$$ for the two simulated cases: SM (red) with two orthogonal Watson distributions with equal weights and kernels, and non-SM (blue) with two orthogonal Watson distributions with equal weights and different kernels. The distributions were estimated from 500 independent noise realizations and these were computed for various SNR levels.

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)
3560