Introduction

Diffusion magnetic resonance imaging is an imaging modality that can probe the tissue microstructure^(1-3). However, the relationship between the diffusion MRI signals and cellular structures is still intricate. Numerical simulations can help investigate the relationship, and the Bloch-Torrey equation is the widely accepted theory that models the diffusion MRI signal^(4,5). Nevertheless, existing simulation frameworks such as finite element methods cannot provide a Fourier representation of the signal. This study proposes the Fourier Potential Method (FPM) as a new spectral method to solve the Bloch-Torrey equation.
The main challenge of spectral methods involves the singularity of the heat kernel. In theory, infinite Fourier modes are required due to the singularity. This limitation may lead to the Gibbs phenomenon that could degrade the approximation accuracy⁽⁶⁾.
We formulate the solution of the Bloch-Torrey equation as a convolution with the heat kernel, which is called the single-layer potential. The convolution allows us to split the heat kernel into a singular part and a smooth part to overcome the singularity. We derived an asymptotic formula to approximate the singular part. Moreover, the recursive method proposed by Greengard et al.^(7,8) can efficiently evaluate the smooth part. Finally, we numerically validated the convergence of our method and showed the error behavior in several simulation parameters.

Method

We restrict ourselves to the 2D diffusion setting with impermeable interfaces and PGSE diffusion-encoding sequence⁽⁹⁾. The main steps of our method are 1) transforming the Bloch-Torrey equation to the diffusion equation using the narrow pulse assumption; 2) formulating the solution of the diffusion equation using the single-layer potential, and splitting it into two parts; 3) approximating the singular part of the single-layer potential using an asymptotic expansion, and storing the smooth part using the Fourier coefficients; 4) computing the diffusion MRI signal using the above representation. We explain each step below. Some mathematical details are omitted due to the length constraint of this abstract.
The narrow pulse assumption allows us to transform the Bloch-Torrey equation to a diffusion equation:
\begin{align}\frac{\partial}{\partial t}{w(\textbf{x},t)}&=\nabla\cdot(\mathcal{D}\nabla w(\textbf{x},t)),&\textbf{x}\in\Omega,t\in[0,\Delta-\delta],\\\mathcal{D}\nabla w(\textbf{x},t)\cdot\textbf{n}&=\mathcal{D}N(\textbf{x},t,\textbf{q}),&\textbf{x}\in\Gamma,t\in[0,\Delta-\delta],\\w(\textbf{x},0)&=0,&\textbf{x}\in\Omega,\end{align}where $$$\Omega$$$ represents the diffusion domain whose boundary is $$$\Gamma$$$ and $$$w$$$ describes the magnetization distribution. The vector $$$\textbf{q}$$$ is $$$\frac{\gamma\delta\textbf{g}}{2\pi}$$$ where $$$\textbf{g}$$$ is the magnetic gradient. The Neumann term $$$N$$$ is complex-valued and decays exponentially in time:\begin{equation}N(\textbf{x},t,\textbf{q})= 2\pi\jmath\rho\,\textbf{q}\cdot\textbf{n}e^{-2\pi\jmath\textbf{q}\cdot\textbf{x}}e^{-4\pi^2\mathcal{D}\|\textbf{q}\|^2t},\textbf{x}\in\Gamma,\end{equation}where $$$\rho$$$ is the initial density of the magnetization. According to Haddar et al.⁽¹⁰⁾, the solution of the diffusion equation is:\begin{equation}w(\textbf{x},t)=\int_{0}^{t}\int_{\partial\Omega}\mathcal{D}G(\textbf{x}-\textbf{y},t-\tau)u(\textbf{y},\tau)ds_\textbf{y}d\tau,\textbf{x}\in\Omega,\end{equation}where $$$u$$$ is an unknown function. The heat kernel G is singular when $$$t$$$ is close to 0. To overcome the singularity, we split $$$w$$$ into a singular part and a smooth part. We derived an asymptotic formula for the singular part on the domain boundary $$$\Gamma$$$ and the smooth part has a Fourier decomposition. In short, the two parts are:\begin{equation}w_{singular}(\textbf{x},t)=\int_{t-\eta}^{t}\int_{\partial \Omega}\mathcal{D}G(\textbf{x}-\textbf{y},t-\tau)u(\textbf{y},\tau)ds_\textbf{y}d\tau\simeq\sqrt{\dfrac{\mathcal{D}\eta}{\pi}}u(\textbf{x}, t)+O(\eta^{3/2}),\end{equation}\begin{equation}w_{smooth}(\textbf{x},t)=\int_{0}^{t-\eta}\int_{\partial\Omega}\mathcal{D}G(\textbf{x}-\textbf{y},t-\tau)u(\textbf{y},\tau)ds_\textbf{y}d\tau=\mathcal{D}\sum_{\nu}f(\nu,t)e^{2\pi\jmath\nu\cdot\textbf{x}}.\end{equation}The parameter $$$\eta$$$ is a controlling factor to separate the singular part from the heat kernel $$$G$$$. Finally, the diffusion MRI signal is\begin{equation}s=\mathcal{D}\rho^{-1}\int_{\partial \Omega}\int_{0}^{\Delta-\delta}\overline{N}(\textbf{x},\Delta-\delta-\tau,\textbf{q}) w(\textbf{x},\tau)d\tau ds_{\textbf{x}}+|\Omega|\rho e^{-4\pi^2\mathcal{D}\|\textbf{q}\|^2(\Delta-\delta)},\end{equation}and the magnetization $$$w$$$ is\begin{equation}w(\textbf{x},t)\simeq\mathcal{D}\sum_{\nu=-\nu_{max}}^{\nu_{max}}f(\nu, t)e^{2\pi\jmath\nu\cdot\textbf{x}}\Delta\nu^2+\sqrt{\dfrac{\mathcal{D}\eta}{\pi}}u(\textbf{x},t),\textbf{x}\in\Gamma.\end{equation}The Fourier-type decomposition of the smooth part of diffusion MRI signal appears in the equation above. We use the recursive method proposed by Greengard et al.⁽⁷⁾ to efficiently evaluate the unknown function $$$u$$$ and the Fourier coefficients $$$f$$$. To numerically implement this method, we discretize the time domain $$$\left[0,\Delta-\delta\right]$$$ with a time step $$$\Delta t$$$ and the domain boundary $$$\Gamma$$$ is discretized with a spatial step $$$\Delta x$$$. In addition, we truncate the Fourier domain at the frequency $$$\nu_{max}$$$. The spectrum is then discretized with a frequency step $$$\Delta\nu$$$.

Results and Discussion

We validate the convergence of our method for some main discretization parameters. The Matrix Formalism method⁽¹¹⁾ can compute analytical signals for circles, so we use the Matrix Formalism signals as reference signals $$$s_{ref}$$$. We note the signal simulated by our method as $$$s$$$. We study the relative error, defined by $$$\varepsilon = |s-s_{ref}| /s_{ref}$$$. The geometry on which we conduct the convergence study is a circle of radius $$$r$$$. The default values for the physical parameters are below:

• $$$r=\{1,2,4\}\mu m$$$, $$$\mathcal{D}_0=0.002\mu m^{2}/\mu s$$$, $$$u_{\mathbf{g}}=[1,0]^{T}$$$
• $$$\delta=10^{-3}\mu s$$$, $$$\Delta=5,000\mu s$$$, $$$b=\{1,4\}ms/\mu m^{2}$$$

We first study the convergence curves regarding the controlling factor $$$\eta$$$, which quantifies the duration of $$$w_{singular}$$$. In Figure 1, the curves show a clear convergence order of $$$3/2$$$ in $$$\eta$$$.
Our method can give a spectral representation for the diffusion MRI signal. To numerically compute the spectrum, we truncate it at $$$\nu_{max}$$$. We present the convergence curves regarding $$$\nu_{max}$$$ in Figure 2. We note that the error decreases exponentially due to the fast decay of Fourier coefficients.
Since our method contains several boundary integrations, we discretize the boundaries by regular polygons. We call the discretized segment length of the boundary $$$\Delta x$$$. Figure 3 illustrates the convergence curves in $$$\Delta x$$$.

Conclusion

This study uses potential theory to derive a Fourier-type representation of the diffusion MRI signal. The splitting of the single-layer potential into two parts avoids the singularity of the heat kernel. It also allows the decomposition of the diffusion MRI signal into the Fourier basis. We numerically validated the convergence of our method and showed the error behavior in several simulation parameters. With proper discretization parameters, the relative error can be less than $$$10^{-4}$$$.

Acknowledgements

We want to thank the French National Institute for Research in Computer Science and Automation (Inria) for providing the PhD scholarship of Chengran Fang and all necessary support.