Matrix formalism

The Matrix Formalism Paul Callaghan (1997), Alexandre V Barzykin (1999) representation of the solution to the Bloch-Torrey PDE is based on the generalized Laplace eigenfunctions of the given geometrical configuration. These are the eigenfunctions of the generalized Laplace operator $-\nabla \cdot \mat{D} \nabla$ subject to the same boundary and interface conditions as the BTPDE. It which does not depend on the gradient sequence. The Matrix Formalism truncates the number of Laplace eigenmodes at an acceptable level, by setting a minimum length scale for the eigenfunctions. The gradient sequence dependent Bloch-Torrey operator can then be expressed in the truncated Laplace eigenfunction basis. This allows for a simple analytical expression which approximates the solution to the BTPDE down to an arbitrary length scale.

Eigenvalue decomposition of the generalized Laplace operator

Let $\{(\phi, \lambda)\}$ be the $L^2$-normalized eigenfunctions and eigenvalues of the generalized Laplace operator $\nabla \cdot \mathbf{D} \nabla$, defined almost everywhere on the domain $\Omega = \bigcup_{i = 1}^{N_\text{cmpt}} \Omega_i$. Denoting $\phi^i$ the restriction of $\phi$ to compartment $\Omega_i$, the eigenfunctions respect the following set of equations, for $(i, j) \in \{1, \dots, N_\text{cmpt}\}^2$:

\[-\nabla \cdot \mathbf{D}_i \nabla \phi^i(\vec{x}) = \lambda \phi^i(\vec{x}), \quad \vec{x} \in \Omega_i, \quad i \in \{1, \dots, N_\text{cmpt}\}.\]

where the boundary conditions are the same as for the BTPDE:

\[\begin{alignedat}{3} \mathbf{D}_i \nabla \phi^i(\vec{x}) \cdot \vec{n}_i(\vec{x}) & = -\mathbf{D}_j \nabla \phi^j(\vec{x}) \cdot \vec{n}_j(\vec{x}), \quad & & \vec{x} \in \Gamma_{i j}, \quad (i, j) \in \{1, \dots, N_\text{cmpt}\}^2, & \\ \mathbf{D}_i \nabla \phi^i(\vec{x}) \cdot \vec{n}_i(\vec{x}) & = \kappa_{i j} \left(c_{i j} \phi^j(\vec{x}) - c_{j i}\phi^i(\vec{x})\right), \quad & & \vec{x} \in \Gamma_{i j}, \quad (i, j) \in \{1, \dots, N_\text{cmpt}\}^2, &\\ \mathbf{D}_i \nabla \phi^i(\vec{x}) \cdot \vec{n}_i(\vec{x}) & = -\kappa_i \phi^i(\vec{x}), \quad & & \vec{x} \in \Gamma_i, \quad i \in \{1, \dots, N_\text{cmpt}\}. & \end{alignedat}\]

Let the solutions $(\phi, \lambda)$ to the above equations be denoted by $\left\{(\phi_n, \lambda_n)\right\}_{n \in \mathbb{N}^*}$. We assume the non-negative real-valued eigenvalues are ordered in non-decreasing order:

\[0 = \lambda_1 \leq \lambda_2 \leq \lambda_3 \leq \dots\]

If the domain $\Omega$ consists of only one contiguous group of compartments connected through a chain of permeable membranes, only the first eigenvalue will be zero, and the corresponding eigenfunction will be the only constant function. If there are $N_\text{group} \geq 2$ groups of connected compartments completely separated by interior hard wall membranes, the first $N_\text{group}$ eigenvalues will be zero:

\[0 = \lambda_1 = \dots = \lambda_{N_\text{group}} < \lambda_{N_\text{group} + 1} \leq \dots,\]

and there will be $N_\text{group}$ corresponding groupwise constant eigenfunctions. In the latter case, the equations may be rewritten separately for each connected subdomain to obtain a set of eigenvalues with a multiplicity of one, but the formulation is also valid in the global form with multiple zero eigenvalues. These two formulations will lead to an identical eigenfunction basis (up to a linear combination for the eigenvalues with multiplicity higher than one), where the basis of each subdomain is a subset of the global eigenfunction basis.

Let $\mathbf{L}$ be the diagonal matrix containing the first $N_\text{eig}$ Laplace eigenvalues:

\[\mathbf{L} = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_{N_\text{eig}})\in \R^{N_\text{eig} \times N_\text{eig}}.\]

Then the matrix $\mathbf{L}$ represents the generalized Laplace operator $-\nabla \cdot \mathbf{D} \nabla$ in the truncated Laplace eigenfunction basis.

Bloch-Torrey operator in the Laplace eigenfunction basis

Let $\mathbf{A}(\vec{g})$ be the $N_\text{eig} \times N_\text{eig}$ matrix defined by:

\[\mathbf{A}(\vec{g}) = g_x \mathbf{A}^x + g_y \mathbf{A}^y + g_z \mathbf{A}^z,\]

where $\vec{g} = (g_x, g_y, g_z)^\mathsf{T}$ is the gradient vector and $\mathbf{A}^x$, $\mathbf{A}^y$, and $\mathbf{A}^z$ are three symmetric $N_\text{eig} \times N_\text{eig}$ matrices whose entries are the first order moments in the coordinate directions of the product of pairs of eigenfunctions:

\[\begin{alignedat}{3} A_{mn}^x & = \int_\Omega x \, \phi_m(\vec{x}) \phi_n(\vec{x}) \, \mathrm{d} \Omega(\vec{x}), \quad & & (m, n) \in \{1, \dots, N_\text{eig}\}^2, \\ A_{mn}^y & = \int_\Omega y \, \phi_m(\vec{x}) \phi_n(\vec{x}) \, \mathrm{d} \Omega(\vec{x}), \quad & & (m, n) \in \{1, \dots, N_\text{eig}\}^2, \\ A_{mn}^z & = \int_\Omega z \, \phi_m(\vec{x}) \phi_n(\vec{x}) \, \mathrm{d} \Omega(\vec{x}), \quad & & (m, n) \in \{1, \dots, N_\text{eig}\}^2. \end{alignedat}\]

Similarly, let $\mathbf{T}$ be the $N_\text{eig} \times N_\text{eig}$ Laplace relaxation matrix defined by

\[T_{mn} = \int_\Omega \frac{1}{T_2(\vec{x})} \phi_m(\vec{x}) \phi_n(\vec{x}) \, \mathrm{d} \Omega(\vec{x}), \quad (m,n) \in \{1, \dots, N_\text{eig}\}^2.\]

Then the general time dependent Bloch-Torrey operator

\[-\nabla \cdot \mathbf{D} \nabla + \frac{1}{T_2} + \underline{\mathrm{i}} \gamma \vec{g}(t) \cdot \vec{x}, \quad t \in [0, T_\text{echo}]\]

in the truncated Laplace eigenfunction basis $(\phi_j)_{j = 1, \dots, N_\text{eig}}$ is given by the complex-valued matrix

\[\mathbf{K}(\vec{g}(t)) = \mathbf{L} + \mathbf{T} + \underline{\mathrm{i}} \gamma \mathbf{A}(\vec{g}(t)), \quad t \in [0, T_\text{echo}].\]

Matrix Formalism approximation

Denoting $\vec{\phi} = (\phi_1, \dots, \phi_{N_\text{eig}})^\mathsf{T}$ the vector of the first $N_\text{eig}$ Laplace eigenfunctions, we may project the solution to to the Bloch-Torrey PDE onto the truncated Laplace eigenfunction basis:

\[M^\text{MF}(\vec{x}, t) = \vec{\phi}^\mathsf{T}\!(\vec{x}) \vec{\nu}(t),\]

where the vector of coefficients $\vec{\nu} = (\nu_1, \dots, \nu_{N_\text{eig}})^\mathsf{T} : [0, T_\text{echo}] \to \mathbb{C}^{N_\text{eig}}$ is the solution to

\[\frac{\mathrm{d} \vec{\nu}}{\mathrm{d} t} = -\mathbf{K}(\vec{g}(t)) \vec{\nu}(t), \quad t \in [0, T_\text{echo}].\]

The initial coefficients are given by

\[\vec{\nu}(0) = \int_\Omega \rho(\vec{x}) \vec{\phi}(\vec{x}) \, \mathrm{d} \Omega(\vec{x}).\]

By using a piece-wise constant approximation of the gradient $\vec{g}$ Alexandre V Barzykin (1999), we obtain

\[\vec{\nu}(T_\text{echo}) \approx \left(\prod_{i = 1}^{N_\text{int}} \mathrm{e}^{-\delta_i \mathbf{K}_i}\right) \vec{\nu}(0) = \mathrm{e}^{-\delta_{N_\text{int}} \mathbf{K}_{N_\text{int}}} \dots \mathrm{e}^{-\delta_2 \mathbf{K}_2} \mathrm{e}^{-\delta_1 \mathbf{K}_1} \vec{\nu}(0),\]

where $\{I_i\}_{i = 1, \dots, N_\text{int}}$ are intervals such that $[0, T_\text{echo}] = \bigcup_{i = 1}^{N_\text{int}} I_i$, $\vec{g}(t) = \vec{g}_i$ for $t \in I_i$, $\delta_i = |I_i|$, and $\mathbf{K}_i = \mathbf{K}(\vec{g}_i)$. The constants may be computed through quadrature:

\[\vec{g}_i = \frac{1}{\delta_i}\int_{I_i} \vec{g}(t) \, \mathrm{d} t \approx \frac{1}{2} \left(\vec{g}(\min I_i) + \vec{g}(\max I_i)\right).\]

For the PGSE and double-PGSE sequences $\vec{g}(t) = f(t) g \vec{d}$, where $f(t) \in \{-1, 0, 1\}$ for all $t$, the coefficients of the final magnetization are given by

\[\vec{\nu}(T_\text{echo}) = \mathrm{e}^{-\delta\mathbf{K}(g \vec{d})^*} \mathrm{e}^{-(\Delta - \delta)(\mathbf{L} + \mathbf{T})} \mathrm{e}^{-\delta \mathbf{K}(g \vec{d}} \vec{\nu}(0)\]

and

\[\vec{\nu}(T_\text{echo}) = \mathrm{e}^{-\delta \mathbf{K}(g \vec{d})^*} \mathrm{e}^{-(\Delta - \delta)(\mathbf{L} + \mathbf{T})} \mathrm{e}^{-\delta \mathbf{K}(g \vec{d})} \mathrm{e}^{-t_\text{pause}(\mathbf{L} + \mathbf{T})} \mathrm{e}^{-\delta \mathbf{K}(g \vec{d})^*} \mathrm{e}^{-(\Delta - + \delta)(\mathbf{L} + \mathbf{T})} \mathrm{e}^{-\delta \mathbf{K}(g \vec{d})} \vec{\nu}(0)\]

respectively. These are the exact solutions, although there may still be truncation errors from $N_\text{eig}$.

The signal is given by

\[S^\text{MF}(f, \vec{g}) = \vec{\nu}^\mathsf{T}\!(T_\text{echo}) \int_\Omega \vec{\phi}(\vec{x}) \, \mathrm{d} \Omega(\vec{x}).\]

Apparent diffusion coefficient

From the Matrix Formalism signal, the analytical expression of its ADC for a gradient direction given by a unit vector $\vec{d}$ and a sequence $f$ is the following:

\[D^\text{MF}(\vec{d}, f) = \frac{1}{|\Omega|} \sum_{n = 1}^{N_\text{eig}} (\vec{d} \cdot \vec{a}_n)^2 \lambda_n \frac{\int_0^{T_\text{echo}} F(t) \int_0^t \mathrm{e}^{-\lambda_n(t - s)} f(s) \, \mathrm{d} s \, \mathrm{d} t}{ \int_0^{T_\text{echo}} F^2(t) \, \mathrm{d} t},\]

where the coefficients of the three-row matrix $\mathbf{a} = (\vec{a}_1^, \dots, \vec{a}_{N_\text{eig}}) \in \mathbb{R}^{3 \times N_\text{eig}}$ are given by

\[\begin{split} a_n^x & = \int_{\Omega} x \phi_n(\vec{x}) \, \mathrm{d} \Omega(\vec{x}), \\ a_n^y & = \int_{\Omega} y \phi_n(\vec{x}) \, \mathrm{d} \Omega(\vec{x}), \\ a_n^z & = \int_{\Omega} z \phi_n(\vec{x}) \, \mathrm{d} \Omega(\vec{x}), \end{split} \quad n \in \{1, \dots, N_\text{eig}\}.\]

To clarify the relationship between the ADC and the diffusion encoding direction $\vec{d}$, we rewrite the Matrix Formalism ADC as

\[D^\text{MF}(\vec{d}, f) = \vec{d}^\mathsf{T} \mathbf{D}^{\text{MF}}(f) \vec{d},\]

where the Matrix Formalism effective diffusion tensor is given by

\[\mathbf{D}^{\text{MF}}(f) = \frac{1}{|\Omega|} \sum_{n = 1}^{N_\text{eig}} j_n(f) \vec{a}_n \vec{a}_n^\mathsf{T} = \frac{1}{|\Omega|} \mathbf{a} \mathbf{j}(f) \mathbf{a}^\mathsf{T},\]

with $\mathbf{j}(f) = \operatorname{diag}\left(j_1(f), \dots, j_{N_\text{eig}}(f)\right) \in \mathbb{R}^{N_\text{eig} \times N_\text{eig}}$ depending on $f$:

\[j_n(f) = \lambda_n \frac{\int_0^{T_\text{echo}} F(t) \int_0^t e^{-\lambda_n(t - s)} f(s) \, \mathrm{d} s \, \mathrm{d} t}{\int_0^{T_\text{echo}} F^2(t) \, \mathrm{d} t}, \quad n \in \{1, \dots, N_\text{eig}\}.\]

In the case of a constant initial spin density $\rho$, we also allow for computing the Matrix Formalism Gaussian Approximation (MFGA) signal, given as

\[S^{\text{MFGA}}(g, f, \vec{d}) = \rho |\Omega| \exp\left(-\vec{d}^\mathsf{T} \mathbf{D}^{\text{MF}}(f) \vec{d} \; b(g, f) \right) }.\]

Eigenfunction length scale and orientation

On a line segment of length $L$ and diffusivity $D_0$, the eigenvalues $(\lambda_1, \lambda_2, \dots)$ of the Laplace operator with Neumann boundary conditions are

\[\lambda_n = \left(\frac{\pi (n - 1)}{L}\right)^2 D_0, \quad n = 1, 2, \dots\]

To make the link between the computed eigenvalue and the spatial scale of the eigenmode, we will convert the computed $\lambda_n$ into a length scale:

\[\ell(\lambda) = \begin{cases} + \infty, \quad & \lambda = 0, \\ \pi \sqrt{\frac{\bar{\sigma}}{\lambda}}, \quad & \lambda > 0, \end{cases}\]

and characterize the computed eigenmode by $\ell(\lambda_n)$ instead of $\lambda_n$. The reference diffusivity $\bar{\sigma}$ is taken as a volume weighted mean of the trace average (spherical part) of $(\mathbf{D}_i)_{1 \leq i \leq N_\text{cmpt}}$:

\[\bar{\sigma} = \frac{1}{|\Omega|} \int_\Omega \frac{1}{3}\operatorname{tr} \mathbf{D}(\vec{x}) \, \mathrm{d} \Omega(\vec{x}).\]

To characterize the directional contribution of the eigenmode we use the fact that its contribution to the ADC in the direction $\vec{d}$ is $j_n(f) (\vec{d} \cdot \vec{a}_n)^2$. We thus call $\vec{a}_n = (a_n^x, a_n^y, a_n^z)^\mathsf{T}$ the "diffusion direction" of the $n$th eigenmode. We remind that the three components of $\vec{a}_n$ are the first moments in the 3 canonical basis directions of the associated eigenfunction.