# Fourier representation of the diffusion MRI signal using layer potentials

Chengran Fang, Demian Wassermann, Jing-Rebecca Li

► To cite this version:

Chengran Fang, Demian Wassermann, Jing-Rebecca Li. Fourier representation of the diffusion MRI signal using layer potentials. SIAM Journal on Applied Mathematics, 2023, 83 (1), pp.99-121. 10.1137/21M1439572. hal-03940100

HAL Id: hal-03940100

https://hal.science/hal-03940100

Submitted on 15 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
FOURIER REPRESENTATION OF THE DIFFUSION MRI SIGNAL
USING LAYER POTENTIALS

CHENGRAN FANG\*†, DEMIAN WASSERMANN\*, AND JING-REBECCA LI†

**Abstract.** The diffusion magnetic resonance imaging signal arising from biological tissues can be numerically simulated by solving the Bloch-Torrey partial differential equation. Numerical simulations can facilitate the investigation of the relationship between the diffusion MRI signals and cellular structures. With the rapid advance of available computing power, the diffusion MRI community has begun to employ numerical simulations for model formulation and validation, as well as for imaging sequence optimization. Existing simulation frameworks use the finite difference method, the finite element method, or the Matrix Formalism method to solve the Bloch-Torrey partial differential equation. We propose a new method based on the efficient evaluation of layer potentials. In this paper, the mathematical framework and the numerical implementation of the new method are described. We demonstrate the convergence of our method via numerical experiments and analyze the errors linked to various model and simulation parameters. Since our method provides a Fourier-type representation of the diffusion MRI signal, it can potentially facilitate new physical and biological signal interpretations in the future.

**Key words.** Bloch-Torrey equation, diffusion magnetic resonance imaging, simulation, layer potential, Fourier transform

**AMS subject classifications.** 35Q99, 65M32, 65Z05, 65E05, 31A10

**1. Introduction.** Diffusion magnetic resonance imaging (diffusion MRI) is a promising non-invasive imaging modality that can probe the tissue microstructure by encoding the motion of water molecules with magnetic gradient pulses [32, 13]. The goal of various imaging protocols is often the recovery of some biological parameters of interest, such as axon diameter and density [38, 30], dendrite structure [11, 29], effective diffusion coefficient [3]. However, in the past, researchers sometimes could not fully validate some of these proposed protocols due to the lack of ground truth. Therefore, recent works have started to include numerical simulations as a part of the validation process [28, 34]. Numerical simulations can facilitate the investigation of the effects of different pulse sequences and tissue features on the measured signal, and can be used for the development, testing, and optimization of novel MRI pulse sequences [20, 31]. For example, the simulation-based inference, which leverages the power of numerical simulations and statistical inference, is a novel direction for cytoarchitecture measurements [10].

Concerning the current numerical simulation methods for diffusion MRI applica-
tions, two popular groups of approaches are Monte-Carlo methods and Bloch-Torrey
PDE-based methods. Monte-Carlo methods use random walkers to mimic the diffu-
sion process in a geometrical configuration. Software packages include Camino [4],
DIFSIM [1], and a GPU-based Monte-Carlo simulator [22]. The Bloch Torrey PDE-
based methods solve the Bloch-Torrey partial differential equation, which describes
the evolution of the complex transverse water proton magnetization under the influ-
ence of diffusion-encoding magnetic field gradient pulses. The predominant numerical
methods to solve this PDE include the finite difference method [14], the finite element
method [17, 21, 2], and the Matrix Formalism method [6, 7, 18]. Although some recent
works in the diffusion MRI community [28, 34] still utilize Monte-Carlo simulations,

<footnote>\*INRIA Saclay, Equipe Parietal, 1 Rue Honoré d'Estienne d'Orves, 91120 Palaiseau, France</footnote>

<footnote>†INRIA Saclay, Equipe IDEFIX, UMA, ENSTA Paris, 828, Boulevard des Maréchaux, 91762 Palaiseau, France (jingrebecca.li@inria.fr).</footnote>
the Bloch-Torrey PDE-based methods have recently demonstrated their potential, including in high-performance computing settings [23, 24, 26, 17, 18] and in manifolds settings for thin-layer and thin-tube geometries [25]. In addition to numerical efficiency, some Bloch-Torrey PDE-based methods allow for a better understanding of the diffusion mechanism. Our previous works in Bloch-Torrey PDE-based neuron simulations demonstrate that diffusion MRI signals reflect the cellular organization of cortical gray matter, and these signals are sensitive to cell size and the presence of large neurons such as the spindle (von Economo) neurons [33, 36, 19, 5].

The Matrix Formalism method [6, 7], which decomposes the solution of the Bloch-Torrey PDE onto a Laplacian eigenbasis, provides another interesting perspective to the diffusion MRI signal. One can address many fundamental theoretical questions about the diffusion MRI signal thanks to the eigendecomposition. In some ways, the Matrix Formalism method inspired us to decompose the diffusion MRI signal into a Fourier type basis. Contrary to the Laplacian eigenbasis, the Fourier basis functions themselves do not depend on the geometrical confinement. This independence should allow for the comparison between various geometries and provides a new spectral perspective.

Based on potential theory from classical mathematics, we propose a new method that provides a Fourier type representation of the diffusion MRI signal. The main challenge to overcome involves the fundamental solution of the diffusion equation, also known as the heat kernel, which has a singularity in time. In theory, infinite Fourier modes are required to represent the heat kernel due to the singularity, while only finite Fourier modes are accessible for practical computation. This practical limitation may lead to the Gibbs phenomenon that could degrade the approximation accuracy [35]. In order to overcome this challenge, we follow the path of several previous works [8, 15, 16, 9] focusing on the evaluation of heat potentials. In particular, in [8], the authors proposed several fundamental ideas, such as 1) splitting the heat potential into a local in time part and a history part in order to overcome the singularity of the heat kernel; 2) approximating the local in time part by asymptotics; 3) leveraging the exponential decay of the history part to represent it using a few Fourier modes. These ideas are crucial to the Fourier type representation of the diffusion MRI signal that we derive in this paper.

Despite the intrinsic similarity between thermal conduction and diffusion process, in the literature, there have not been previous works about the representation of the diffusion MRI signal via potential theory, and certainly not by using a Fourier basis for layer potentials. As the first paper addressing this subject, we restrict ourselves to the 2D diffusion MRI setting with impermeable interfaces. We also restrict ourselves to simplified conditions on the diffusion-encoding gradient, specifically, we derive our method under the narrow pulse assumption, where the diffusion-encoding pulse duration is very short compared to the delay between the pulses. These two assumptions allow us to apply the theory developed for the diffusion kernel to the diffusion MRI application.

The main steps of our method are 1) transforming the Bloch-Torrey PDE to the diffusion equation using the narrow pulse assumption on the diffusion-encoding sequence; 2) formulating the solution of the diffusion equation using the single layer potential; 3) approximating the singular part of the single layer potential using an asymptotic expansion and solving the integral equation; 4) storing the non-singular part of the single layer potential using the Fourier coefficients, leveraging the fast decay in the spectrum; 5) computing the diffusion MRI signal using the above representation. We call our method the Fourier Potential Method (FPM).
The paper is organized as follows. Section 2 introduces the mathematical framework and the Bloch-Torrey PDE of diffusion MRI. Section 3 presents the Fourier Potential Method and error analysis. Section 4 contains numerical results, including convergence in the various simulation parameters. Section 5 contains conclusions and future work.

**2. Mathematical frame of diffusion MRI.** Suppose one would like to simulate the diffusion MRI signal due to spins inside a biological cell and assume that the spin exchange across the cell membrane is negligible under the simulation conditions. Let $\Omega$ be the domain that describes the geometry of the biological cell and let $\Gamma = \partial\Omega$ be the cell membrane.

**2.1. Bloch-Torrey PDE.** In diffusion MRI, a time-varying magnetic field gradient is applied to the tissue to encode water diffusion. Denoting the effective time profile of the diffusion-encoding magnetic field gradient by $f(t)$, and let the vector $g$ contain the amplitude and direction information of the magnetic field gradient, the complex transverse water proton magnetization in the rotating frame satisfies the Bloch-Torrey PDE:

$$ (2.1) \qquad \frac{\partial}{\partial t} M(\mathbf{x}, t) = -\gamma f(t) \mathbf{g} \cdot \mathbf{x} M(\mathbf{x}, t) + \nabla \cdot (\mathcal{D}_0 \nabla M(\mathbf{x}, t)), \quad \mathbf{x} \in \Omega, $$

where $\gamma = 267.513 \text{ rad} \, \mu s^{-1} T^{-1}$ is the gyromagnetic ratio of the water proton, $\gamma$ is the imaginary unit, $\mathcal{D}_0$ is the intrinsic diffusion coefficient in the neuron compartment $\Omega$. The magnetization is a function of position $\mathbf{x}$ and time $t$, and depends on the diffusion gradient vector $\mathbf{g}$ and the time profile $f(t)$.

A commonly used time profile (diffusion-encoding sequence) is the pulsed-gradient spin echo (PGSE) [32] sequence, with two rectangular pulses of duration $\delta$, separated by a time interval $\Delta - \delta$, for which the profile $f(t)$ is

$$ (2.2) \qquad f(t) = \begin{cases} 1, & 0 \le t \le \delta, \\ -1, & \Delta < t \le \Delta + \delta, \\ 0, & \text{otherwise.} \end{cases} $$

In the case that the rectangular pulses are narrow, i.e., $\delta \ll \Delta$, this allows the Bloch-Torrey PDE to be transformed to the diffusion equation. This assumption is called the narrow pulse approximation [32] and it is taken up in subsection 2.2.

The PDE in (2.1) needs boundary conditions. We assume negligible membrane permeability, meaning the zero Neumann boundary condition:

$$ \mathcal{D}_0 \nabla M(\mathbf{x}, t) \cdot \mathbf{n} = 0, \quad \mathbf{x} \in \partial\Omega, $$

where $\mathbf{n}$ is the unit outward pointing normal vector at $\mathbf{x}$. In addition, the PDE has the constant initial condition:

$$ M(\mathbf{x}, 0) = \rho, \quad \mathbf{x} \in \Omega, $$

where $\rho$ is the initial spin density.

The diffusion MRI signal is measured at echo time $t = T_E \ge \Delta + \delta$ for PGSE. This signal is the integral of $M(\mathbf{x}, T_E)$:

$$ s = \int_{\mathbf{x} \in \bigcup \{\Omega\}} M(\mathbf{x}, T_E) \, d\mathbf{x}, $$
where $\cup\{\Omega\}$ represents a set of biological compartments with impermeable membranes.

The signal $s$ is usually plotted against a quantity called the b-value [32, 27]. The b-value depends on $g$ and $f(t)$ and is defined as

$$b(g) = \gamma^2 \|g\|^2 \int_{0}^{T_E} du \left( \int_{0}^{u} f(s) ds \right)^2 .$$

For PGSE, the b-value is [32]:

$$(2.3) \qquad b(g, \delta, \Delta) = \gamma^2 \|g\|^2 \delta^2 (\Delta - \delta/3).$$

The reason for these definitions is that in a homogeneous medium, the signal attenuation is $e^{-D_0 b}$, where $D_0$ is the intrinsic diffusion coefficient.

**2.2. Narrow pulse approximation.** In this paper, we restrict ourselves to simplified conditions on the diffusion-encoding gradient, specifically, we derive our method under the narrow pulse assumption, where the pulse duration is very short compared to the delay between the pulses [32], i.e., $\delta \ll \Delta$. This will lead to the solution of a diffusion equation instead of the more complicated Bloch-Torrey PDE, as explained below.

Let us consider spins initially located at $x$. After the first pulse, the complex phase of these spins is $e^{-j\delta\gamma g \cdot x}$. This means the complex magnetization at $t = \delta$ due to a uniform distribution of initial spins with density $\rho$ can be written as:

$$M(x, \delta) \approx \rho e^{-j\delta\gamma g \cdot x}, \quad x \in \Omega.$$

Because the gradient magnetic field is turned off after the first pulse, the spins move but the phase of the spins does not change. Therefore, the magnetization between pulses satisfies the diffusion equation:

$$(2.4) \qquad \frac{\partial}{\partial t} M(x, t) = \nabla \cdot (D_0 \nabla M(x, t)), \quad x \in \Omega, t \in [\delta, \Delta],$$

subject to the zero Neumann boundary condition:

$$(2.5) \qquad D_0 \nabla M(x_0, t) \cdot n = 0, \quad x_0 \in \partial\Omega, t \in [\delta, \Delta],$$

where $n$ is the unit outward pointing normal vector at $x_0$, and initial condition at $t = \delta$ is:

$$(2.6) \qquad M(x, \delta) = \rho e^{-j\delta\gamma g \cdot x}, \quad x \in \Omega.$$

During the second pulse, at the point $x$, the additional accumulated complex phase is $e^{j\delta\gamma g \cdot x}$, so the magnetization at the position $x$ and time $TE$ is:

$$(2.7) \qquad M(x, TE) \approx M(x, \Delta) e^{j\delta\gamma g \cdot x}, \quad x \in \Omega.$$

We emphasize again that we used the assumption $\delta \ll \Delta$. The echo-time $TE$ is usually some time after the end of the second pulse (i.e. $TE \ge \Delta + \delta$).

The diffusion MRI signal $s$ is the total magnetization measured at the echo time:

$$(2.8) \qquad s = \int_{x \in \Omega} M(x, \Delta) e^{j\delta\gamma g \cdot x} dx.$$
**3. Method.** We derive our new method below.

**3.1. Solution of the diffusion equation and the diffusion MRI signal.**

Before we solve the diffusion equation using potential theory, we transform the initial and boundary conditions. We transform the diffusion equation in (2.4)–(2.6) such that it is subject to zero initial conditions and complex-valued non-zero Neumann boundary conditions. Define

$$ (3.1) \qquad \omega(\mathbf{x}, t) \equiv M(\mathbf{x}, t + \delta) - \rho e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 t} e^{-2\pi j \mathbf{q} \cdot \mathbf{x}}, \quad \mathbf{x} \in \Omega, t \in [0, \Delta - \delta], $$

where $\mathbf{q} = \delta\gamma\mathbf{g}/2\pi$. We will work on the quantity $\omega(\mathbf{x}, t)$ in (3.1), which satisfies the diffusion equation:

$$ (3.2) \qquad \frac{\partial}{\partial t} \omega(\mathbf{x}, t) = \nabla \cdot (\mathcal{D}_0 \nabla \omega(\mathbf{x}, t)), \quad \mathbf{x} \in \Omega, t \in [0, \Delta - \delta], $$

subject to non-homogeneous Neumann boundary conditions:

$$ (3.3) \qquad \mathcal{D}_0 \nabla \omega(\mathbf{x}_0, t) \cdot \mathbf{n} = \mathcal{D}_0 \mathcal{N}(\mathbf{x}_0, t, \mathbf{q}) \quad \mathbf{x}_0 \in \partial\Omega, t \in [0, \Delta - \delta], $$

and zero initial conditions:

$$ (3.4) \qquad \omega(\mathbf{x}, 0) = 0, \quad \mathbf{x} \in \Omega. $$

The Neumann forcing term is complex-valued, periodic in space in the direction $\mathbf{q}$, and decays exponentially in time:

$$ (3.5) \qquad \mathcal{N}(\mathbf{x}, t, \mathbf{q}) \equiv 2\pi\rho \mathbf{q} \cdot \mathbf{n} (je^{-2\pi j \mathbf{q} \cdot \mathbf{x}}) e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 t}. $$

The diffusion MRI signal $s$ can be reformulated in terms of $\omega$:

$$ (3.6) \qquad s = |\Omega| \rho e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 (\Delta - \delta)} + \int_{\mathbf{x} \in \Omega} \omega(\mathbf{x}, \Delta - \delta) e^{2\pi j \mathbf{q} \cdot \mathbf{x}} d\mathbf{x}. $$

In the above, the first term is explicit, the second term needs to be computed. We define a time dependent integral whose value at $t = \Delta - \delta$ gives second term:

$$ (3.7) \qquad \bar{\omega}(\mathbf{q}, t) \equiv \int_{\mathbf{x} \in \Omega} \omega(\mathbf{x}, t) e^{2\pi j \mathbf{q} \cdot \mathbf{x}} d\mathbf{x}, t \in [0, \Delta - \delta]. $$

The function $\bar{\omega}$ can be expanded by the Green's second identity:

$$ \begin{aligned} \bar{\omega}(\mathbf{q}, t) &= \frac{-1}{4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2} \left( \int_{\Omega} \nabla \cdot (\mathcal{D}_0 \nabla \omega(\mathbf{x}, t)) e^{2\pi j \mathbf{q} \cdot \mathbf{x}} d\mathbf{x} + B \right), \\ B &= \int_{\partial\Omega} 2\pi j \mathcal{D}_0 \mathbf{q} \cdot \mathbf{n} \omega(\mathbf{x}, t) e^{2\pi j \mathbf{q} \cdot \mathbf{x}} ds_{\mathbf{x}} - \int_{\partial\Omega} \mathcal{D}_0 \nabla \omega(\mathbf{x}, t) \cdot \mathbf{n} e^{2\pi j \mathbf{q} \cdot \mathbf{x}} ds_{\mathbf{x}}. \end{aligned} $$

Using the diffusion equation and the nonhomogeneous Neumann boundary conditions, we get an ordinary differential equation for $\bar{\omega}$:

$$ (3.8) \qquad \frac{d}{dt} \bar{\omega}(\mathbf{q}, t) = -4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 \bar{\omega}(\mathbf{q}, t) - 2\pi j \mathcal{D}_0 \int_{\partial\Omega} \mathbf{q} \cdot \mathbf{n} \omega(\mathbf{x}, t) e^{2\pi j \mathbf{q} \cdot \mathbf{x}} ds_{\mathbf{x}}, $$

which has an analytical solution:

$$ \begin{aligned} \bar{\omega}(\mathbf{q}, t) &= -\mathcal{D}_0 \int_{\partial\Omega} \int_0^t 2\pi j \mathbf{q} \cdot \mathbf{n} e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 (t-\tau)} \omega(\mathbf{x}, \tau) e^{2\pi j \mathbf{q} \cdot \mathbf{x}} d\tau ds_{\mathbf{x}}, \\ &= \mathcal{D}_0 \rho^{-1} \int_{\partial\Omega} \int_0^t \mathcal{N}^*(\mathbf{x}, t-\tau, \mathbf{q}) \omega(\mathbf{x}, \tau) d\tau ds_{\mathbf{x}}. \end{aligned} $$
The asterisk symbol * denotes the complex conjugation. It can be proved that (3.9) satisfies a recursive relationship in time:

$$ (3.10) \qquad \bar{\omega}(\mathbf{q}, t) = e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 \Delta t} \bar{\omega}(\mathbf{q}, t - \Delta t) + \mathcal{D}_0 \rho^{-1} \int_{\partial\Omega} \int_{t-\Delta t}^{t} \mathcal{N}^*(\mathbf{x}, t-\tau, \mathbf{q})\omega(\mathbf{x}, \tau) d\tau ds_{\mathbf{x}}. $$

Equation (3.7) and (3.9) are mathematically equivalent for evaluating the diffusion MRI signal (at $t = \Delta - \delta$). It can be seen that, while (3.7) requires the value of $\omega$ on the entire domain $\Omega$, (3.9) only needs the value of $\omega$ on the boundary, which is more computationally efficient. The recursion in time above also increases the computational efficiency. We will use the method of layer potentials to get the boundary values in the next section.

**3.2. The single layer potential representation.** The PDE in (3.2)–(3.4) has Neumann boundary conditions, zero initial conditions and zero forcing term, allowing us to represent the solution $\omega(\mathbf{x}, t)$ as a single layer potential, with a density function $\mu$ defined on $\partial\Omega$ [9]. In other words, $\omega(\mathbf{x}, t) = S[\mu](\mathbf{x}, t)$. The definition for the single layer potential is

$$ (3.11) \qquad \omega(\mathbf{x}, t) = S[\mu](\mathbf{x}, t) = \int_{0}^{t} \int_{\partial\Omega} \mathcal{D}_{0}G(\mathbf{x} - \mathbf{y}, t - \tau)\mu(\mathbf{y}, \tau)ds_{\mathbf{y}}d\tau, $$

where $G(\mathbf{x}, t)$ is the fundamental solution of the 2D diffusion equation in a box $[-L_1/2, L_1/2] \times [-L_2/2, L_2/2]$, with periodic boundary conditions. The fundamental solution $G(\mathbf{x}, t)$ has two equivalent representations [8]:

$$ (3.12) \qquad G_{\text{Gauss}}(\mathbf{x}, t) = (4\pi \mathcal{D}_0 t)^{-1} \sum_{\mathbf{z} \in \mathbb{Z}^2} e^{-\frac{\|\mathbf{x} - \mathbf{z} \odot \mathbf{L}\|^2}{4\mathcal{D}_0 t}}, $$

$$ (3.13) \qquad G_{\text{Fourier}}(\mathbf{x}, t) = \frac{1}{L_1 L_2} \sum_{\substack{\nu = \mathbf{z} \odot \mathbf{L} \\ \mathbf{z} \in \mathbb{Z}^2}} e^{-4\pi^2 \mathcal{D}_0 \|\nu\|^2 t} e^{2\pi j \nu \cdot \mathbf{x}}, $$

where $\odot$ and $\ominus$ are hadamard product and hadamard division, respectively, and $\mathbf{L} = [L_1, L_2]^T$. For the convenience of notation, in the following, we set $L_1 = L_2 = L$ and note by $\Delta\nu = \frac{1}{L}$. In this way, we rewrite (3.13) as

$$ (3.14) \qquad G_{\text{Fourier}}(\mathbf{x}, t) = \sum_{\substack{\nu = \mathbf{z} \odot \mathbf{L} \\ \mathbf{z} \in \mathbb{Z}^2}} e^{-4\pi^2 \mathcal{D}_0 \|\nu\|^2 t} e^{2\pi j \nu \cdot \mathbf{x}} \Delta \nu^2 $$

in order to recall its relationship with the Fourier transform. The imposition of periodic boundary conditions on the faces of the box allows us to use the discrete Fourier series.

The density function $\mu$ is chosen to be a causal function and is determined by imposing the Neumann boundary conditions on the geometry boundary $\partial\Omega$ [16]:

$$ \lim_{\mathbf{x} \to \mathbf{x}_0 \in \partial\Omega} \nabla S[\mu](\mathbf{x}, t) \cdot \mathbf{n} = \mathcal{N}(\mathbf{x}_0, t, \mathbf{q}), \quad \mathbf{x}_0 \in \partial\Omega, t \in [0, \Delta - \delta]. $$
Using the jump property of the trace of the double layer potential, the integral equation to be solved for $\mu$ is then the following:

$$ (3.15) \qquad \frac{1}{2}\mu(\mathbf{x}_0, t) + K[\mu](\mathbf{x}_0, t) = \mathcal{N}(\mathbf{x}_0, t, \mathbf{q}), \quad \mathbf{x}_0 \in \partial\Omega, t \in [0, \Delta - \delta], $$

with

$$ (3.16) \qquad K[\mu](\mathbf{x}_0, t) \equiv \int_0^t \int_{\partial\Omega} \mathcal{D}_0 \frac{\partial G}{\partial \mathbf{n}_{\mathbf{x}_0}}(\mathbf{x}_0 - \mathbf{y}, t - \tau) \mu(\mathbf{y}, \tau) ds_{\mathbf{y}} d\tau $$

being the principal value integral on the boundary. Solving the integral equation (3.15) for $\mu$ plays the pivotal role in our method. We present the detailed steps in the next sections.

### 3.3. Splitting the single layer potential into local and history parts.

The single layer potential $S[\mu]$ is split into a history part, $S_{\text{long}}[\mu]$, and a local in time part, $S_{\text{short}}[\mu]$. Since the local in time part $S_{\text{short}}[\mu]$ contains the singularity of the fundamental solution $G$, we approximate it by asymptotic formulas. The asymptotic trace formulas are only accurate in an interval near the singularity, so we limit their use to the interval $[t - \eta, t]$, with $\eta$ being a small quantity to be determined later. In other words,

$$ (3.17) \qquad S[\mu](\mathbf{x}, t) = S_{\text{short}}[\mu](\mathbf{x}, t) + S_{\text{long}}[\mu](\mathbf{x}, t), $$

with

$$ (3.18) \qquad S_{\text{short}}[\mu](\mathbf{x}, t) := \int_{t-\eta}^{t} \int_{\partial\Omega} \mathcal{D}_0 G_{\text{Gauss}}(\mathbf{x} - \mathbf{y}, t-\tau) \mu(\mathbf{y}, \tau) ds_{\mathbf{y}} d\tau, $$

$$ (3.19) \qquad S_{\text{long}}[\mu](\mathbf{x}, t) := \int_{0}^{t-\eta} \int_{\partial\Omega} \mathcal{D}_0 G_{\text{Fourier}}(\mathbf{x} - \mathbf{y}, t-\tau) \mu(\mathbf{y}, \tau) ds_{\mathbf{y}} d\tau. $$

Similarly, we decompose $K[\mu]$ into 2 parts:

$$ (3.19) \qquad K[\mu](\mathbf{x}_0, t) = K_{\text{short}}[\mu](\mathbf{x}_0, t) + K_{\text{long}}[\mu](\mathbf{x}_0, t), $$

with

$$ (3.20) \qquad K_{\text{short}}[\mu](\mathbf{x}_0, t) := \int_{t-\eta}^{t} \int_{\partial\Omega} \mathcal{D}_0 \frac{\partial G_{\text{Gauss}}}{\partial \mathbf{n}_{\mathbf{x}_0}}(\mathbf{x}_0 - \mathbf{y}, t-\tau) \mu(\mathbf{y}, \tau) ds_{\mathbf{y}} d\tau, $$

$$ (3.20) \qquad K_{\text{long}}[\mu](\mathbf{x}_0, t) := \int_{0}^{t-\eta} \int_{\partial\Omega} \mathcal{D}_0 \frac{\partial G_{\text{Fourier}}}{\partial \mathbf{n}_{\mathbf{x}_0}}(\mathbf{x}_0 - \mathbf{y}, t-\tau) \mu(\mathbf{y}, \tau) ds_{\mathbf{y}} d\tau. $$

Next, we compute or approximate the above history and local parts.

#### 3.3.1. Asymptotic trace formulas for the local part.
Based on the expressions derived in [8], the asymptotic trace formulas in two dimensions for the local parts, when $t > \eta$, are:

$$ (3.21) \qquad S_{\text{short}}[\mu](\mathbf{x}_0, t) = \sqrt{\frac{D_0 \eta}{\pi}} \mu(\mathbf{x}_0, t) + O(\eta^{3/2}), t > \eta $$

and

$$ (3.22) \qquad K_{\text{short}}[\mu](\mathbf{x}_0, t) = -\frac{\sqrt{D_0 \eta}}{2\sqrt{\pi}} \kappa(\mathbf{x}_0) \mu(\mathbf{x}_0, t) + O(\eta^{3/2}), t > \eta, $$
where $\kappa(\mathbf{x}_0)$ is the curvature at the point $\mathbf{x}_0 \in \partial\Omega$. The boundary $\partial\Omega$, which models the cell membrane, is a closed 2D plane curve. We assume it is twice differentiable. Let $\psi(\alpha) = (x(\alpha), y(\alpha))$ be a parametric representation of $\partial\Omega$. We choose a general parameter $\alpha$ such that $\psi(\alpha)$ is oriented counterclockwise. The curvature at the point $\mathbf{x}_0 = \psi(\alpha_0)$ is defined as

$$ (3.23) \qquad \kappa(\mathbf{x}_0) = \frac{x'y'' - y'x''}{(x'^2 + y'^2)^{3/2}} \Big|_{\alpha=\alpha_0}, $$

where primes refer to derivatives with respect to $\alpha$.

We also need to initialize values for $t \le \eta$. It has been derived in [9] that the expressions are:

$$ (3.24) \qquad S_{\text{short}}[\mu](\mathbf{x}_0, t) = \sqrt{\frac{D_0 t}{\pi}} \mu(\mathbf{x}_0, t) + O(t^{3/2}), t \le \eta, $$

and

$$ (3.25) \qquad K_{\text{short}}[\mu](\mathbf{x}_0, t) = - \frac{\sqrt{D_0 t}}{2\sqrt{\pi}} \kappa(\mathbf{x}_0) \mu(\mathbf{x}_0, t) + O(t^{3/2}), t \le \eta. $$

### 3.3.2. Fourier representation of history part.
For the smooth part of the single layer potential, a Fourier representation for the Dirichlet trace is proposed in [8]:

$$ (3.26) \qquad S_{\text{long}}[\mu](\mathbf{x}_0, t) = D_0 \sum_{\substack{\nu=\mathbf{z} \otimes \mathbf{L} \\ \mathbf{z} \in \mathbb{Z}^2}} \hat{f}(\nu, t) e^{2\pi j \nu \cdot \mathbf{x}_0} \Delta \nu^2, $$

and the Neumann trace is

$$ (3.27) \qquad K_{\text{long}}[\mu](\mathbf{x}_0, t) = D_0 \sum_{\substack{\nu=\mathbf{z} \otimes \mathbf{L} \\ \mathbf{z} \in \mathbb{Z}^2}} 2\pi j \nu \cdot \mathbf{n} \hat{f}(\nu, t) e^{2\pi j \nu \cdot \mathbf{x}_0} \Delta \nu^2, $$

where the Fourier coefficients are

$$ (3.28) \qquad \hat{f}(\nu, t) = \int_{0}^{t-\eta} \int_{\partial\Omega} e^{-4\pi^2 D_0 \|\nu\|^2 (t-\tau)} \mu(\mathbf{y}, \tau) e^{-2\pi j \nu \cdot \mathbf{y}} ds_{\mathbf{y}} d\tau. $$

To avoid history dependent time integration, we use the following recurrence formula for the Fourier coefficients

$$ (3.29) \qquad \hat{f}(\nu, t) = e^{-4\pi^2 D_0 \|\nu\|^2 \Delta t} \hat{f}(\nu, t - \Delta t) + \int_{t-\eta-\Delta t}^{t-\eta} \int_{\partial\Omega} e^{-4\pi^2 D_0 \|\nu\|^2 (t-\tau)} \mu(\mathbf{y}, \tau) e^{-2\pi j \nu \cdot \mathbf{y}} ds_{\mathbf{y}} d\tau, $$

so only local-in-time integrals are needed at each time step.

The above formulas hold when $t > \eta$. For $t \le \eta$, we initialize $S_{\text{long}}[\mu]$, $K_{\text{long}}[\mu]$, and $\hat{f}$ to be 0.
**3.4. Computation of the single layer density.** Based on the decomposition of the single layer potential and the approximation of the history and the local parts detailed previously, we can compute the density function $\mu$.

For $t \le \eta$, substituting (3.25) into (3.15) and solving the integral equation, we can get the approximation to the density

$$ (3.30) \qquad \mu(\mathbf{x}_0, t) = \frac{2\mathcal{N}(\mathbf{x}_0, t, \mathbf{q})}{1 - \sqrt{\frac{D_0 t}{\pi}} \kappa(\mathbf{x}_0)} + O(t^{3/2}), \quad \mathbf{x}_0 \in \partial\Omega, t \le \eta. $$

For $t \in (\eta, \Delta - \delta]$, the integral equation (3.15) can be rewritten as

$$ \frac{1}{2}\mu(\mathbf{x}_0, t) + K_{\text{short}}[\mu](\mathbf{x}_0, t) = \beta(\mathbf{x}_0, t), $$

where the right hand side is

$$ (3.31) \qquad \beta(\mathbf{x}_0, t) \equiv -K_{\text{long}}[\mu](\mathbf{x}_0, t) + \mathcal{N}(\mathbf{x}_0, t, \mathbf{q}). $$

We write the solution of the above integral equation as

$$ (3.32) \qquad \mu(\mathbf{x}_0, t) = 2(I + 2K_{\text{short}})^{-1} [\beta](\mathbf{x}_0, t), \quad \mathbf{x}_0 \in \partial\Omega, t \in (\eta, \Delta - \delta], $$

and expand the operator $(I + 2K_{\text{short}})^{-1}$ (corresponding to $K_{\text{short}}$ being a contraction) as

$$ (3.33) \qquad \mu(\mathbf{x}_0, t) = 2(I - 2K_{\text{short}} + 4K_{\text{short}}^2 + \dots + (-2)^n K_{\text{short}}^n + \dots) [\beta](\mathbf{x}_0, t). $$

We approximate $K_{\text{short}}^n[\beta]$ using (3.22) and we get

$$ (3.34) \qquad K_{\text{short}}^n[\beta](\mathbf{x}_0, t) = \frac{1}{(-2)^n} \left(\frac{D_0 \eta}{\pi}\right)^{n/2} \kappa^n(\mathbf{x}_0) \beta(\mathbf{x}_0, t) + O(\eta^{3/2}). $$

Then, we keep all terms of the operator expansion to obtain

$$ \begin{align} \mu(\mathbf{x}_0, t) &= 2\left(\beta(\mathbf{x}_0, t) - 2K_{\text{short}}[\beta](\mathbf{x}_0, t) + 4K_{\text{short}}^2[\beta](\mathbf{x}_0, t) + \dots\right) \nonumber \\ (3.35) \qquad \qquad \qquad \qquad &= 2\beta(\mathbf{x}_0, t)\left(1 + \left(\frac{D_0 \eta}{\pi}\right)^{\frac{1}{2}}\kappa(\mathbf{x}_0) + \frac{D_0 \eta}{\pi}\kappa^2(\mathbf{x}_0) + \dots\right) + O(\eta^{3/2}) \nonumber \\ \qquad \qquad \qquad \qquad &= 2\beta(\mathbf{x}_0, t)/\left(1 - \sqrt{\frac{D_0 \eta}{\pi}}\kappa(\mathbf{x}_0)\right) + O(\eta^{3/2}). \end{align} $$

**3.5. Computation of the single layer potential.** Once the density $\mu$ is obtained, we compute the single layer potential $S[\mu]$ in the following way.

When $t \le \eta$, the expression for the single layer potential is

$$ \begin{align*} S[\mu](\mathbf{x}_0, t) &= S_{\text{short}}[\mu](\mathbf{x}_0, t) \\ &= \sqrt{\frac{D_0 t}{\pi}} \frac{4\pi j \rho q \cdot n e^{-4\pi^2 D_0 \|q\|^2 t} e^{-2\pi j q \cdot x_0}}{1 - \sqrt{\frac{D_0 t}{\pi}} \kappa(\mathbf{x}_0)} + O(t^{3/2}), \quad \mathbf{x}_0 \in \partial\Omega, t \in [0, \eta]. \end{align*} $$

For $t \in (\eta, \Delta - \delta]$, the single layer potential has both a local part and a history part. The local part is

$$ S_{\text{short}}[\mu](\mathbf{x}_0, t) = \sqrt{\frac{D_0 \eta}{\pi}} \mu(\mathbf{x}_0, t) + O(\eta^{3/2}), \quad \mathbf{x}_0 \in \partial\Omega, t \in (\eta, \Delta - \delta]. $$
As for the history part $S_{long}[\mu]$, it can be approximated by the truncated Fourier series:
(3.36)

$$S_{long}[\mu](\mathbf{x}_0, t) = \mathcal{D}_0 \sum_{\mathbf{\nu}=-\mathbf{\nu}_{max}}^{\mathbf{\nu}_{max}} \hat{f}(\mathbf{\nu}, t) e^{2\pi j \mathbf{\nu} \cdot \mathbf{x}_0} \Delta \mathbf{\nu}^2 + E(\nu_{max}), \quad \mathbf{x}_0 \in \partial\Omega, t \in (\eta, \Delta - \delta].$$

In the above, we denote the error term due to truncating the infinite Fourier series up to $\nu_{max}$ by $E(\nu_{max})$. We do not have an analytical expression for $E(\nu_{max})$, but we will show later in the numerical results that it decays exponentially in $\nu_{max}$.

The addition of $S_{short}[\mu]$ and $S_{long}[\mu]$ gives the single layer potential $S[\mu]$ which is the solution of (3.2) on the boundary:

$$S[\mu](\mathbf{x}_0, t) = S_{short}[\mu](\mathbf{x}_0, t) + S_{long}[\mu](\mathbf{x}_0, t), \quad \mathbf{x}_0 \in \partial\Omega, t \in [0, \Delta - \delta].$$

At the current iteration step, the Fourier coefficients $\hat{f}$ that are still unknown, will be computed using the density function $\mu$ from the previous iterations, as explained in the following.

**3.5.1. Computation of the Fourier coefficients of the history part.** For $t \le \eta$, $\hat{f}$ is set to zero, as well as $K_{long}[\mu]$. For $t \in (\eta, 2\eta]$, $\hat{f}$ are computed using the density $\mu$ from the previous iterations:

$$(3.37) \qquad \hat{f}(\mathbf{\nu}, t) = e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{\nu}\|^2 \Delta t} \hat{f}(\mathbf{\nu}, t - \Delta t) + \underbrace{\int_{\partial\Omega} \int_{t-\eta-\Delta t}^{t-\eta} e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{\nu}\|^2 (t-\tau)} e^{-2\pi j \mathbf{\nu} \cdot \mathbf{y}} \mu(\mathbf{y}, \tau) d\tau ds_{\mathbf{y}}}_{\hat{f}_{\text{temp1}}(\mathbf{\nu}, t)}, \quad \mathbf{\nu} \in [-\nu_{max}, \nu_{max}]^2,$$

with

$$\hat{f}_{\text{temp1}}(\mathbf{\nu}, t) = \int_{\partial\Omega} 4\pi j \mathbf{q} \cdot \mathbf{n} e^{-2\pi j (\mathbf{q}+\mathbf{\nu}) \cdot \mathbf{y}} \underbrace{\int_{t-\eta-\Delta t}^{t-\eta} \frac{e^{-4\pi^2 \mathcal{D}_0 [\|\mathbf{\nu}\|^2(t-\tau) + \|\mathbf{q}\|^2\tau]}}{1 - \sqrt{\frac{\mathcal{D}_0 \tau}{\pi}} \kappa(\mathbf{y})}}_{p} d\tau \, ds_{\mathbf{y}}.$$

We apply the trapezoidal rule to the time integration $p$ to obtain

$$p = \begin{cases} -\frac{2\pi e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 t}}{\mathcal{D}_0 \kappa^2(\mathbf{y})} \left[ \kappa(\mathbf{y}) \sqrt{\frac{\mathcal{D}_0}{\pi}} (\sqrt{t-\eta} - \sqrt{t-\eta-\Delta t}) + \ln\left(\frac{1-\kappa(\mathbf{y})\sqrt{\frac{\mathcal{D}_0}{\pi}(t-\eta)}}{1-\kappa(\mathbf{y})\sqrt{\frac{\mathcal{D}_0}{\pi}(t-\eta-\Delta t)}}\right) \right], & \|\mathbf{\nu}\| = \|\mathbf{q}\|; \\ \\ e^{-4\pi^2 \mathcal{D}_0 [\|\mathbf{q}\|^2(t-\eta) + \|\mathbf{\nu}\|^2\eta]} \left[ \frac{1+e^{4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 - \|\mathbf{\nu}\|^2)\Delta t}}{\Delta t(4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 - \|\mathbf{\nu}\|^2)\Delta t - 1)} + \frac{e^{4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 - \|\mathbf{\nu}\|^2)\Delta t - 4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 - \|\mathbf{\nu}\|^2)\Delta t - 1}}{\Delta t(4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 - \|\mathbf{\nu}\|^2))^2 (1-\kappa(\mathbf{y})\sqrt{\frac{\mathcal{D}_0}{\pi}(t-\eta-\Delta t)})}, & \|\mathbf{\nu}\| \neq \|\mathbf{q}\|. \end{cases}$$

Once we compute the time integration $p$, the integration over the boundary $\partial\Omega$ can be approximated by uniform discretization in arc length.
Remaining on $t \in (\eta, 2\eta]$, next we compute the long time part $K_{long}[\mu]$ via the Fourier series

$$K_{long}[\mu](\mathbf{x}_0, t) = \mathcal{D}_0 \sum_{\substack{\nu=-\nu_{max} \\ \nu_{max}}}^{\nu_{max}} 2\pi j \mathbf{\nu} \cdot \mathbf{n} \hat{f}(\mathbf{\nu}, t) e^{2\pi j \mathbf{\nu} \cdot \mathbf{x}_0} \Delta \nu^2 + E(\nu_{max}).$$

With a slight abuse of notation, we use the same notation $E(\nu_{max})$ as in (3.36) for the error due to truncating the Fourier series at $\nu_{max}$.

Finally, the density function $\mu$ for (the current time) $t \in (\eta, 2\eta]$ is computed as:

$$\mu(\mathbf{x}_0, t) = \frac{2[\mathcal{N}(\mathbf{x}_0, t) - K_{long}[\mu](\mathbf{x}_0, t)]}{1 - \sqrt{\frac{\mathcal{D}_0 \eta}{\pi}} \kappa(\mathbf{x}_0)} + O(\eta^{3/2}).$$

On the rest of the time interval, $t \in (2\eta, \Delta - \delta]$, $\hat{f}$ still uses the density $\mu$ from previous iterations, but the formulas are different:

$$(3.38) \qquad \hat{f}(\mathbf{\nu}, t) = e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{\nu}\|^2 \Delta t} \hat{f}(\mathbf{\nu}, t - \Delta t) + \underbrace{\int_{\partial \Omega} \int_{t-\eta-\Delta t}^{t-\eta} e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{\nu}\|^2 (t-\tau)} e^{-2\pi j \mathbf{\nu} \cdot \mathbf{y}} \mu(\mathbf{y}, \tau) d\tau ds_{\mathbf{y}}}_{\hat{f}_{temp2}(\mathbf{\nu}, t)}, \quad \mathbf{\nu} \in [-\nu_{max}, \nu_{max}]^2.$$

In the above, the Fourier coefficients $\hat{f}(\mathbf{\nu}, t - \Delta t)$ at the previous time step are known, and the expression of $\mu(\mathbf{x}_0, \tau)$ for $\tau \in (\eta, \Delta - \delta - \eta]$ is

$$\mu(\mathbf{x}_0, \tau) = 2 \left( 1 - \sqrt{\frac{\mathcal{D}_0 \eta}{\pi}} \kappa(\mathbf{x}_0) \right)^{-1} [\mathcal{N}(\mathbf{x}_0, \tau) - K_{long}[\mu](\mathbf{x}_0, \tau)], \quad \mathbf{x}_0 \in \partial\Omega.$$

The integration on the right hand side of (3.38) is noted as $\hat{f}_{temp2}(\mathbf{\nu}, t)$ in which we substitute the expression of $\mu$ above. We split $\hat{f}_{temp2}$ into two parts and gather the terms that are independent of time

$$\hat{f}_{temp2}(\mathbf{\nu}, t) = \int_{\partial\Omega} 2 \left(1 - \sqrt{\frac{\mathcal{D}_0\eta}{\pi}} \kappa(\mathbf{y})\right)^{-1} e^{-2\pi j \mathbf{\nu} \cdot \mathbf{y}} \times \\ \underbrace{(2\pi j \mathbf{q} \cdot \mathbf{n} e^{-2\pi j \mathbf{q} \cdot \mathbf{y}} \int_{t-\eta-\Delta t}^{t-\eta} e^{-4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 \tau + \|\mathbf{\nu}\|^2 (t-\tau))} d\tau -}_{h_1} \\ \underbrace{\int_{t-\eta-\Delta t}^{t-\eta} K_{long}[\mu](\mathbf{y}, \tau) e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{\nu}\|^2 (t-\tau)} d\tau}_{h_2} ds_{\mathbf{y}}.$$

The time integration $h_1$ in the first part has an analytical expression

$$(3.39) \qquad h_1 = \begin{cases} \Delta t \cdot e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{\nu}\|^2 t} & \|\mathbf{q}\| = \|\mathbf{\nu}\| \\ e^{-4\pi^2 \mathcal{D}_0 [\|\mathbf{q}\|^2(t-\eta) + \|\mathbf{\nu}\|^2\eta]} \frac{e^{4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 - \|\mathbf{\nu}\|^2)\Delta t} - 1}{4\pi^2 \mathcal{D}_0 (\|\mathbf{q}\|^2 - \|\mathbf{\nu}\|^2)} & \|\mathbf{q}\| \neq \|\mathbf{\nu}\| \end{cases}$$
The time integration $h_2$ in the second part has to be calculated numerically. We apply the trapezoidal rule to $K_{long}[\mu](\mathbf{y}, \tau)$ and we get
(3.40)

$$h_2 = \begin{cases} \frac{\Delta t}{2} [K_{long}[\mu](\mathbf{y}, t - \eta) + K_{long}[\mu](\mathbf{y}, t - \eta - \Delta t)] & \|\nu\| = 0 \\ \left[ \frac{1 - e^{-4\pi^2 \mathcal{D}_0 \|\nu\|^2 \Delta t} (4\pi^2 \mathcal{D}_0 \|\nu\|^2 \Delta t + 1)}{(4\pi^2 \mathcal{D}_0 \|\nu\|^2)^2 \Delta t} K_{long}[\mu](\mathbf{y}, t - \eta - \Delta t) + \frac{e^{-4\pi^2 \mathcal{D}_0 \|\nu\|^2 \Delta t} + 4\pi^2 \mathcal{D}_0 \|\nu\|^2 \Delta t - 1}{(4\pi^2 \mathcal{D}_0 \|\nu\|^2)^2 \Delta t} K_{long}[\mu](\mathbf{y}, t - \eta) \right] e^{-4\pi^2 \mathcal{D}_0 \|\nu\|^2 \eta} & \|\nu\| \neq 0 \end{cases}$$

The values of $K_{long}[\mu]$ at time $t-\eta-\Delta t$ and $t-\eta$ have been computed in previous steps, thus the expressions for $h_1$ and $h_2$ can be computed in the current time step. Then we discretize uniformly in the arc length over the boundary to obtain $\hat{f}_{temp2}$ as well as $\hat{f}$.

Staying on $t \in (2\eta, \Delta - \delta]$, it is straightforward to recover the long time part $K_{long}[\mu]$ at time $t$ by applying the inverse discrete Fourier transform

$$K_{long}[\mu](\mathbf{x}_0, t) = \mathcal{D}_0 \sum_{\nu = -\nu_{max}}^{\nu_{max}} 2\pi j \nu \cdot \mathbf{n} \hat{f}(\nu, t) e^{2\pi j \nu \cdot \mathbf{x}_0} \Delta \nu^2 + E(\nu_{max}).$$

Again, with a slight abuse of notation, we use the same notation $E(\nu_{max})$ as in (3.36) for the error due to truncating the Fourier series at $\nu_{max}$. Finally, the density function $\mu$ at the current time $t$ is

$$\mu(\mathbf{x}_0, t) = \frac{2[\mathcal{N}(\mathbf{x}_0, t) - K_{long}[\mu](\mathbf{x}_0, t)]}{1 - \sqrt{\frac{\mathcal{D}_0 \eta}{\pi}} \kappa(\mathbf{x}_0)} + O(\eta^{3/2}),$$

which will be used for future iterations.

**3.6. Computation of the diffusion MRI signal.** After obtaining the single layer potential, the following procedure produces the diffusion MRI signal.

The diffusion MRI signal $s$ has the representation

$$s = |\Omega| \rho e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 (\Delta - \delta)} + \bar{\omega}(\mathbf{q}, \Delta - \delta).$$

The quantity $\bar{\omega}$ will be computed using the recursive relationship below (rewritten from (3.10)):

$$\bar{\omega}(\mathbf{q}, t) = e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 \Delta t} \bar{\omega}(\mathbf{q}, t - \Delta t) - \mathcal{D}_0 \int_{\partial\Omega} 2\pi j \mathbf{q} \cdot \mathbf{n} e^{2\pi j \mathbf{q} \cdot \mathbf{y}} \underbrace{\int_{t-\Delta t}^{t} e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 (t-\tau)} \omega(\mathbf{y}, \tau) d\tau}_{u} ds_y.$$

By applying the trapezoidal rule to the time integration $u$, we then get the expression

$$(3.41) \qquad u = \begin{cases} \frac{\Delta t}{2} [\omega(\mathbf{y}, t - \Delta t) + \omega(\mathbf{y}, t)] & \|\mathbf{q}\| = 0 \\ \frac{1 - e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 \Delta t} (4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 \Delta t + 1)}{(4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2)^2 \Delta t} \omega(\mathbf{y}, t - \Delta t) + \frac{e^{-4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 \Delta t} + 4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2 \Delta t - 1}{(4\pi^2 \mathcal{D}_0 \|\mathbf{q}\|^2)^2 \Delta t} \omega(\mathbf{y}, t) & \|\mathbf{q}\| \neq 0 \end{cases}$$
The variable $\omega$ is the single layer potential $S[\mu]$

$$ \omega(\mathbf{x}_0, t) = S[\mu](\mathbf{x}_0, t) = S_{short}[\mu](\mathbf{x}_0, t) + S_{long}[\mu](\mathbf{x}_0, t). $$

The short time part has an asymptotic expression

$$ S_{short}[\mu](\mathbf{x}_0, t) = 2\sqrt{\frac{\mathcal{D}_{0\eta}}{\pi}} \frac{\mathcal{N}(\mathbf{x}_0, t) - K_{long}[\mu](\mathbf{x}_0, t)}{1 - \sqrt{\frac{\mathcal{D}_{0\eta}}{\pi}} \kappa(\mathbf{x}_0)} + O(\eta^{3/2}) $$

with

$$ \mathcal{N}(\mathbf{x}_0, t) = 2\pi j \rho q \cdot n e^{-4\pi^2 \mathcal{D}_0 \|q\|^2 t} e^{-2\pi j q \cdot x_0}, $$

and

$$ K_{long}[\mu](\mathbf{x}_0, t) = \mathcal{D}_0 \sum_{\nu=-\nu_{max}}^{\nu_{max}} 2\pi j \nu \cdot n \hat{f}(\nu, t) e^{2\pi j \nu \cdot \mathbf{x}_0} \Delta \nu^2 + E(\nu_{max}). $$

The long time part is approximated by a Fourier series

$$ S_{long}[\mu](\mathbf{x}_0, t) = \mathcal{D}_0 \sum_{\nu=-\nu_{max}}^{\nu_{max}} \hat{f}(\nu, t) e^{2\pi j \nu \cdot \mathbf{x}_0} \Delta \nu^2 + E(\nu_{max}). $$

Finally, a uniform arc length discretization of the boundary allows the numerical computation of $\bar{\omega}$.

**4. Numerical results.** In this section, we study the convergence of the Fourier Potential method. The Matrix Formalism method [7] is capable of computing analytical signals for simple geometries, such as circles and spheres, using analytical expressions for the Laplace eigen-decomposition, so we use the Matrix Formalism signals as the reference signals. We note the diffusion MRI signal simulated by our method as $s$ and the analytical signal given by the Matrix Formalism method as $s_{ref}$.

The geometry on which we will conduct the convergence study is a circle of radius $r$ ($\phi = 2r, \kappa = 1/r$), where $\phi$ is the size of the geometry and $\kappa$ is the curvature. The default values for the physical parameters are below:

* $r = \{1, 2, 4\} \mu m, \kappa = \{1, 0.5, 0.25\} \mu m^{-1}$

* $\mathcal{D}_0 = 2 \times 10^{-3} \mu m^2 / \mu s$

* $\delta = 10^{-3} \mu s, \Delta = 5, 000 \mu s$

* $u_g = [1, 0]^T$

* $b = \{1, 4\} ms/\mu m^2$

We will study the dependence of the relative error, defined by

$$ (4.1) \qquad \varepsilon = \frac{|s - s_{ref}|}{s_{ref}}, $$

on the discretization parameters: spatial step $\Delta x$, the time step $\Delta t$, the maximum frequency $\nu_{max}$, the spectral step $\Delta \nu$, and the single layer local time interval $\eta$. As the convergence studies for various algorithm parameters are conducted, the default values for the fixed parameters are listed below:

* $\eta = 1 \mu s$

* $\nu_{max} = 10 \mu m^{-1}, \Delta \nu = 0.05 \mu m^{-1}$

* $\Delta x = 0.005 \mu m, \Delta t = 0.5 \mu s$
**4.1. The narrow pulse assumption error.** One important point to discuss here, before showing the convergence studies, is the choice of the duration of the diffusion-encoding gradient pulse, $\delta$. To satisfy the narrow pulse assumption, we need that $\delta \ll \Delta$. In Figure 1, we show the error due to the narrow pulse assumption, for a range of $\delta$ values. At $b = 1\,\text{ms}/\mu\text{m}^2$, the narrow pulse approximation error is around $10^{-2}$ at $\delta = 10^2\mu\text{s}$ for all three circle radii.

The validity of the narrow pulse assumption also depends on the separation between the two pulses $\Delta$. The narrow pulse assumption requires a small ratio of $\delta$ and $\Delta$. In Figure 2, we show the influence of this ratio on the relative error. For large $\Delta$ such as 20ms, the relative error is less than 5% with $\delta$ being 2ms.

Despite the fact that a relative error of a few percent is perfectly acceptable for diffusion MRI applications, for the sake of the numerical convergence study that follows, we have chosen much lower thresholds for the narrow pulse approximation error and picked an exceedingly small value of $\delta = 10^{-3}\mu\text{s}$, which is not achievable with current MRI scanners. The reason for this choice is we wanted the error from the narrow pulse assumption to be significantly smaller than the discretization errors of the numerical method as we refined the method parameters. In this way, the plateauing of the errors towards the narrow pulse approximation error occurs later in the refinement process. We note that at our choice of $\delta = 10^{-3}\mu\text{s}$, the narrow pulse approximation errors shown in Figure 1 range from $10^{-6}$ ($b = 1\,\text{ms}/\mu\text{m}^2$, $r = 4\mu\text{m}$) to $10^{-4}$ (higher b-values). These values will form the “floor” values for our convergence curves, to be shown next.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Influence of $\delta$ on the relative error. All discretization parameters are set to be the default. The sampled $\delta$'s are {0.0001, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 100} $\mu$s (from left to right).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Influence of $\delta/\Delta$ on the relative error. All discretization parameters are set to be the default. The geometry is a circle of radius 1µm and the b-value is 4ms/µm$^2$. The sampled ratios $\delta/\Delta$ are {0.1%, 0.5%, 1%, 2.5%, 5%, 7.5%, 10%} (from left to right).</figcaption></figure>

## 4.2. Duration of the local in time part of the single layer potential, $\eta$.
First, we study the duration $\eta$ of the local in time part of the single layer potential. The error term $O(\eta^{3/2})$ originates from the asymptotic trace formulas (3.21) and (3.22). In Figure 3 the curves show a clear convergence order of $\frac{3}{2}$ in $\eta$. In addition, the circle radius (curvature) and the b-value affect the errors: the errors are bigger for larger b-value and higher curvature.

We observe that minimum errors occur at $\eta = 1\mu s$. The values of the minimum errors coincide with the size of the narrow pulse approximation errors shown in Figure 1. At the smaller value of $\eta = 0.5\mu s$, the errors increased. The reason is that there is a tradeoff between two sources of error, one linked to $O(\eta^{3/2})$ and one to $E(\nu_{max})$. With a smaller $\eta$, the long time part $S_{long}[\mu]$ suffers more from the singularity of the heat kernel, thereby increasing the error $E(\nu_{max})$. After we decrease $\eta$ beyond a certain point, $E(\nu_{max})$ becomes the bottleneck for the accuracy, which will be studied next.

## 4.3. Maximum frequency.
A main feature of our method is that the history part of the single layer potential $S_{long}[\mu]$ has a spectral representation. The spectrum of the fundamental solution $G$ decreases exponentially with respect to the frequency $\nu$. As a result, the Fourier coefficients $\hat{f}$ are also subject to the exponential decay:

$$ (4.2) \qquad \hat{f}(\nu, t) = O(e^{-4\pi^2 D_0 \eta \|\nu\|^2}). $$

In order to numerically compute the spectrum of $S_{long}[\mu]$, we truncated it at $\nu_{max}$ and omitted all higher frequency components. The truncation gives rise to an error caused by the omitted Fourier modes, which we have denoted as $E(\nu_{max})$. Even though we
<figure><img src="image_4.png" /><figcaption>FIG. 3. Convergence curves regarding $\eta$. All discretization parameters except for $\eta$ are set to be the default. The sampled $\eta$'s are {0.5, 1, 2.5, 5, 12.5, 25, 50, 75, 100, 125, 250, 375, 500, 625, 750, 875, 1000} $\mu$s (from left to right). The slopes of the curves are around 3/2.</figcaption></figure>

do not have an analytical expression for $E(\nu_{max})$, considering the exponential decay of the Fourier coefficients, we could expect a rapid decrease of the truncation error.

We present the convergence curves in Figure 4. We note that the x-axis is linear, and the y-axis is logarithmic. Empirically, we observe that the error can be fitted by $c_1e^{-c_2\nu_{max}}$, where $c_1$ is a constant, and $c_2 = 1.17$. As expected, for the largest $\nu_{max}$, ($\nu_{max} > 9\mu m^{-1}$), the curves approach the errors due to the narrow pulse approximation.

**4.4. Spatial discretization.** Our method contains several boundary integrations. The geometries we used are circles, and we chose to have a piecewise linear approximation of $\partial\Omega$. This means the discretized geometries are regular polygons. Let us call the discretized segment length of the boundary $\Delta x$. On the other hand, the reference solution, the Matrix Formalism method, computes the Laplace eigenfunctions of exact circles.

Figure 5 illustrates the convergence curves in $\Delta x$. At the larger range of $\Delta x$, we observe exponential convergence in $\frac{1}{\Delta x}$, due to the exponential convergence of the trapezoidal rule for periodic functions (the integrand over a closed boundary being a periodic function). At the smaller range of $\Delta x$, we observe the plateauing towards the narrow pulse approximation errors.

In the middle range of $\Delta x$, we observe a convergence of $O(\Delta x^2)$, due to the approximation of the exact circle geometry by regular polygons. For better visualization of the convergence pattern, we plot the approximation error for the area of an exact circle by regular polygons. The area error $e$ is defined as the normalized difference
<figure><img src="image_5.png" /><figcaption>FIG. 4. Convergence curves regarding $\nu_{max}$. All discretization parameters except for $\nu_{max}$ are set to be the default. The sampled $\nu_{max}$'s are {0.1, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12} $\mu m^{-1}$ (from left to right).</figcaption></figure>

between the circle area and the area of a regular n-sided inscribed polygon $A_n$

$$ (4.3) \qquad e = \frac{\pi r^2 - A_n}{\pi r^2} \sim \frac{\Delta x^2}{6r^2} = \frac{(\kappa \Delta x)^2}{6}. $$

This explains the convergence order of $\Delta x^2$. Moreover, (4.3) also indicates the influence of curvature. High curvature geometries endure greater area errors and thus larger simulation errors.

**4.5. Temporal discretization.** Let the time step be $\Delta t$. We apply the trapezoidal rule to every time integration in our implementation, for instance, (3.37), (3.40), and (3.41). Theoretically, the trapezoidal integration error is $O(\Delta t^2)$ [12]. However, the local-in-time region size, $\eta$, which is an integer multiple of $\Delta t$, contributes an error from the asymptotic formula, as shown in Figure 3. This asymptotic formula error numerically dominates the $O(\Delta t^2)$ error from the trapezoidal integration. Thus, we do not show a plot of the trapezoidal rule convergence.

**4.6. Spectral discretization.** The spectral resolution $\Delta \nu$ is closely related to the size of the periodic box enclosing the geometry. Let the side length of the box be $L$. According to the Nyquist-Shannon sampling theorem, we should have

$$ \frac{1}{\Delta \nu} = L. $$

The inverse real relationship manifests itself in (3.12) and (3.13) as well. A necessary restriction on the box is that it must contain the entire domain $\Omega$, in our case, the
<figure><img src="image_6.png" /><figcaption>FIG. 5. Convergence curves regarding $\Delta x$. All discretization parameters except for $\Delta x$ are set to be the default. The sampled $\Delta x$'s are {0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5} µm (from left to right).</figcaption></figure>

domain being a circle, we get

$$ (4.4) \qquad \frac{1}{\Delta \nu} = L \geq \phi, $$

where we defined $\phi$ as twice the radius. In Figure 6, it is shown that the relative errors are greater than 100% when the box is smaller than the domain ($\frac{1}{\Delta \nu} < \phi$). As soon as the box contains the geometry, the errors reduce to the plateau values of the narrow approximation errors. Clearly, all simulations must satisfy the spectral discretization condition (4.4).

**4.7. Influence of q-vector.** Now we study the influence of the b-value/q-vector on the relative errors. We fix the diffusion time $\delta$ and $\Delta$, so the b-value is equivalent to the square of the magnitude of the q-vector. We chose to plot the relative error versus the magnitude of the q-vector because we explicitly formulate our method using q-vectors rather than b-values. The results are given in Figure 7. We note that the x-axis is logarithmic, and the y-axis is linear. The experiment results show that once the magnitude of the q-vector is large enough, the error increases logarithmically with the norm of q-vectors. For small q-vectors, the errors are within the range of the error floor ($10^{-4} - 10^{-6}$) imposed by the narrow pulse approximation. For larger $q$, the logarithmic dependence the error on $\|q\|$ requires further study to explain. We do not at this time have an explication for it.

**4.8. Extension to complex geometries.** In the previous sections, we used circles to study the convergence of our method. The Matrix Formalism method can compute the analytical solution on circles, which allowed us to show the convergence
<figure><img src="image_7.png" /><figcaption>FIG. 6. Convergence curves regarding $\Delta\nu$. Relative errors which are greater than 100 (2dB) are omitted for a better visualization. All discretization parameters except for $\Delta\nu$ are set to be the default. The sampled $\Delta\nu$'s are {0.025, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.75, 1, 2, 3, 4, 5} $\mu m^{-1}$ (from left to right).</figcaption></figure>

behavior of the Fourier Potential Method.

In fact, our method can simulate diffusion MRI signals on more complex 2D geometries. Here we present FPM simulation results on two realistic axons. The microscopy image (Figure 8) and the axon sections are obtained using the AxonDeepSeg segmentation framework [37]. With these irregular shapes, analytical solutions are not accessible, so we computed the reference signals by finite element simulations using the SpinDoctor toolbox [17]. We show, in Figure 9, the dMRI signals in 40 directions as well as the relative errors. Our method agrees with the finite element reference signals. For the middle b-value ($4ms/\mu m^2$), the relative error is less than 5%. One should note that the magnetization of the two adjacent axons is computed simultaneously by sharing the same Fourier basis. This feature is different from the Matrix Formalism method, which requires geometry-dependent bases.

The method that we derived in this paper is applicable to 2D geometries. To extend to 3 dimensions, the main changes to be made are the asymptotic trace formulas for the local part, i.e., (3.21) - (3.25), in particular, the curvature for 1D curves will need to be generalized to analogous quantities on 2D surfaces. As a consequence, the solution of the integral equation (3.35) will have a new formulation in 3D. Another change involves spatial integration on 2D surfaces instead of on 1D curves, the former being more numerically complicated than the latter.

**5. Conclusions.** In this paper, we derived a new representation of the diffusion MRI signal by solving the Bloch-Torrey PDE using potential theory. The decomposition of the single layer potential into singular and smooth parts allows the numerically
<figure><img src="image_8.png" /><figcaption>FIG. 7. Influence of $q$ on the relative error. All discretization parameters are set to be the default. The sampled $q$'s are {0.0225, 0.0318, 0.0450, 0.0551, 0.0637, 0.0712, 0.1007, 0.1424, 0.1743, 0.2013, 0.2251} $\mu\text{m}^{-1}$ and the corresponding $b$-values are {0.1, 0.2, 0.4, 0.6, 0.8, 1, 2, 4, 6, 8, 10} $\text{ms/}\mu\text{m}^2$ (from left to right).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. The microscopy image of axons from AxonDeepSeg. Two adjacent axons are selected.</figcaption></figure>

efficient storage of the smooth part in the Fourier basis. Time integrals that are in the form of certain exponentials allow us to use time recursion to avoid history-dependence. We numerically validated the convergence of our method and showed the error behavior in several simulation parameters.

One of the main features of our method is the availability of the spectrum of the smooth part of the magnetization field. The projection to the Fourier basis functions provides a unified spectrum space for different geometries. Since our method provides a Fourier like representation of the diffusion MRI signal, this can potentially facilitate new physical and biological signal interpretation in the future.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Comparison of FPM with FEM. (a) the normalized signals simulated by FPM and FEM. The gray areas illustrate the shapes of the two adjacent axons. The physical parameters are: $D_0 = 2 \times 10^{-3} \mu m^2/\mu s$, $\delta = 2 ms$, $\Delta = 100 ms$. The discretization parameters of the FPM are: $\eta = 50 \mu s$, $\nu_{max} = 2 \mu m^{-1}$, $\Delta\nu = 0.05 \mu m^{-1}$, $\Delta x = 0.01 \mu m$, and $\Delta t = 50 \mu s$. The signals are simulated in 40 directions evenly distributed on a unit circle. (b) the relative errors in percent.</figcaption></figure>

As the first paper addressing this subject, we restricted ourselves to the 2D diffusion MRI setting with impermeable interfaces. We also restrict ourselves to simplified conditions on the diffusion-encoding gradient, specifically, we derive our method under the narrow pulse assumption. Future work is planned to generalize beyond these two assumptions.

**Acknowledgments.** We would like 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. We also thank Leslie Greengard for introducing Jing-Rebecca Li to the methodology of the fast evaluation of heat layer potentials during a post-doc many years ago.

REFERENCES

[1] G. T. BALLS AND L. R. FRANK, *A simulation environment for diffusion weighted mr experiments in complex media*, Magn. Reson. Med., 62 (2009), pp. 771-778, http://dx.doi.org/10.1002/mrm.22033.

[2] L. BELTRACHINI, Z. A. TAYLOR, AND A. F. FRANGI, *A parametric finite element solution of the generalised bloch-torrey equation for arbitrary domains*, Journal of Magnetic Resonance, 259 (2015), pp. 126-134, https://doi.org/10.1016/j.jmr.2015.08.008, http://www.sciencedirect.com/science/article/pii/S1090780715001743.

[3] L. M. BURCAW, E. FIEREMANS, AND D. S. NOVIKOV, *Mesoscopic structure of neuronal tracts from time-dependent diffusion*, NeuroImage, 114 (2015), pp. 18-37, https://doi.org/10.1016/j.neuroimage.2015.03.061.

[4] P. COOK, Y. BAI, M. HALL, S. NEDJATI-GILANI, K. SEUNARINE, AND D. ALEXANDER, *Camino: Diffusion mri reconstruction and processing*, (2005).

[5] C. FANG, V.-D. NGUYEN, D. WASSERMANN, AND J.-R. LI, *Diffusion mri simulation of realistic neurons with spinductor and the neuron module*, NeuroImage, 222 (2020), p. 117198.

[6] D. GREBENKOV, *Laplacian eigenfunctions in nmr. i. a numerical tool*, Concepts in Magnetic Resonance Part A, 32A (2008), pp. 277-301, https://doi.org/10.1002/cmr.a.20117.

[7] D. S. GREBENKOV, *Laplacian eigenfunctions in nmr. ii. theoretical advances*, Concepts Magn. Reson., 34A (2009), pp. 264-296, https://doi.org/10.1002/cmr.a.20145.
[8] L. GREENGARD AND J. STRAIN, A fast algorithm for the evaluation of heat potentials, Comm. Pure Appl. Math., 43 (1990), pp. 949-963.

[9] H. HADDAR, J.-R. LI, AND S. SCHIAVI, Understanding the time-dependent effective diffusion coefficient measured by diffusion mri: the intracellular case, SIAM Journal on Applied Mathematics, 78 (2018), pp. 774-800.

[10] M. JALLAIS, P. L. RODRIGUES, A. GRAMFORT, AND D. WASSERMANN, Cytoarchitecture measurements in brain gray matter using likelihood-free inference, in International Conference on Information Processing in Medical Imaging, Springer, 2021, pp. 191-202.

[11] S. N. JESPERSEN, C. D. KROENKE, L. ASTERGAARD, J. J. ACKERMAN, AND D. A. YABLON-SKIY, Modeling dendrite density from magnetic resonance diffusion measurements, NeuroImage, 34 (2007), pp. 1473-1486, http://www.sciencedirect.com/science/article/pii/S1053811906010950.

[12] S. G. JOHNSON, Notes on the convergence of trapezoidal-rule quadrature, 2010.

[13] D. LE BIHAN, E. BRETON, D. LALLEMAND, P. GRENIER, E. CABANIS, AND M. LAVAL-JEANTET, MR imaging of intravoxel incoherent motions: application to diffusion and perfusion in neurologic disorders., Radiology, 161 (1986), pp. 401-407, http://radiology.rsna.org/content/161/2/401 abstract.

[14] J.-R. LI, D. CALHOUN, C. POUON, AND D. L. BIHAN, Numerical simulation of diffusion mri signals using an adaptive time-stepping method, Physics in Medicine and Biology, 59 (2014), p. 441, http://stacks.iop.org/0031-9155/59/i=2/a=441.

[15] J.-R. LI AND L. GREENGARD, On the numerical solution of the heat equation i: Fast solvers in free space, Journal of Computational Physics, 226 (2007), pp. 1891-1901, http://www.sciencedirect.com/science/article/B6WHY-4P4NPNB-1/2/38cdafb0fce301e05c730ed8fc503473.

[16] J.-R. LI AND L. GREENGARD, High order accurate methods for the evaluation of layer heat potentials, SIAM J. Sci. Comput., 31 (2009), pp. 3847-3860, https://doi.org/10.1137/080732389.

[17] J.-R. LI, V.-D. NGUYEN, T. N. TRAN, J. VALDMAN, C.-B. TRANG, K. V. NGUYEN, D. T. S. VU, H. A. TRAN, H. T. A. TRAN, AND T. M. P. NGUYEN, Spindoctor: A matlab toolbox for diffusion mri simulation, NeuroImage, 202 (2019), p. 116120, https://doi.org/10.1016/j.neuroimage.2019.116120, http://www.sciencedirect.com/science/article/pii/S1053811919307116.

[18] J.-R. LI, T. N. TRAN, AND V.-D. NGUYEN, Practical computation of the diffusion mri signal of realistic neurons based on laplace eigenfunctions, NMR in Biomedicine, 33 (2020), p. e4353.

[19] V. MENON, G. GUILLERMO, M. A. PINSK, V.-D. NGUYEN, J.-R. LI, W. CAI, AND D. WASSERMANN, Quantitative modeling links in vivo microstructural and macrofunctional organization of human and macaque insular cortex, and predicts cognitive control abilities, bioRxiv, (2019), https://doi.org/10.1101/662601, https://www.biorxiv.org/content/early/2019/06/06/662601, https://arxiv.org/abs/https://www.biorxiv.org/content/early/2019/06/06/662601 full.pdf.

[20] M. MERCREDI AND M. MARTIN, Toward faster inference of micron-scale axon diameters using monte carlo simulations, Magnetic Resonance Materials in Physics, Biology and Medicine, 31 (2018), pp. 511-530, https://doi.org/10.1007/s10334-018-0680-1, https://doi.org/10.1007/s10334-018-0680-1.

[21] D. V. NGUYEN, J.-R. LI, D. GREBENKOV, AND D. LE BIHAN, A finite elements method to solve the Bloch-Torrey equation applied to diffusion magnetic resonance imaging, Journal of Computational Physics, 263 (2014), pp. 283-302, http://www.sciencedirect.com/science/article/pii/S0021999114000308.

[22] K. V. NGUYEN, E. H. GARZON, AND J. VALETTE, Efficient gpu-based monte-carlo simulation of diffusion in real astrocytes reconstructed from confocal microscopy, Journal of Magnetic Resonance, (2018), https://doi.org/10.1016/j.jmr.2018.09.013, http://www.sciencedirect.com/science/article/pii/S1090780718302386.

[23] V. D. NGUYEN, A fenics-hpc framework for multi-compartment bloch-torrey models, in EC-COMAS Congress 2016, vol. 1, 2016, pp. 105-119.

[24] V.-D. NGUYEN, J. JANSSON, J. HOFFMAN, AND J.-R. LI, A partition of unity finite element method for computational diffusion mri, Journal of Computational Physics, (2018), https://doi.org/10.1016/j.jcp.2018.08.039, http://www.sciencedirect.com/science/article/pii/S0021999118305709.

[25] V.-D. NGUYEN, J. JANSSON, H. T. A. TRAN, J. HOFFMAN, AND J.-R. LI, Diffusion mri simulation in thin-layer and thin-tube media using a discretization on manifolds, Journal of Magnetic Resonance, 299 (2019), pp. 176-187, https://doi.org/https://doi.org/10.1016/j.jmr.2019.01.002, http://www.sciencedirect.com/science/article/pii/S1090780719300023.
[26] V.-D. NGUYEN, M. LEONI, T. DANCHEVA, J. JANSSON, J. HOFFMAN, D. WASSERMANN, AND J.-R. LI, Portable simulation framework for diffusion mri, Journal of Magnetic Resonance, 309 (2019), p. 106611, https://doi.org/https://doi.org/10.1016/j.jmr.2019.106611, http://www.sciencedirect.com/science/article/pii/S1090780719302502.

[27] D. S. NOVIKOV, E. FIEREMANS, S. N. JESPERSEN, AND V. G. KISELEV, Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation, NMR in Biomedicine, 32 (2019), p. e3998, https://doi.org/10.1002/nbm.3998, https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.3998.

[28] M. PALOMBO, A. IANUS, M. GUERRERI, D. NUNES, D. C. ALEXANDER, N. SHEMESH, AND H. ZHANG, Sandi: a compartment-based model for non-invasive apparent soma and neurite imaging by diffusion mri, NeuroImage, (2020), p. 116835.

[29] M. PALOMBO, C. LIGNEUL, AND J. VALETTE, Modeling diffusion of intracellular metabolites in the mouse brain up to very high diffusion-weighted Magnetic Resonance in Medicine, 77 (2017), pp. 343–350, https://doi.org/10.1002/mrm.26548.

[30] E. PANAGIOTAKI, T. SCHNEIDER, B. SIOW, M. G. HALL, M. F. LYTHGOE, AND D. C. ALEXANDER, Compartment models of the diffusion MR signal in brain white matter: A taxonomy and comparison, NeuroImage, 59 (2012), pp. 2241–2254, http://www.sciencedirect.com/science/article/pii/S1053811911011566.

[31] G. RENSONNET, B. SCHERRER, S. K. WARFIELD, B. MACQ, AND M. TAQUET, Assessing the validity of the approximation of diffusion-weighted-mri signals from crossing fascicles by sums of signals from single fascicles, Magnetic Resonance in Medicine, 79 (2018), pp. 2332–2345, https://doi.org/10.1002/mrm.26832, https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.26832, https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/mrm.26832.

[32] E. O. STEJSKAL AND J. E. TANNER, Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient, The Journal of Chemical Physics, 42 (1965), pp. 288–292, https://doi.org/10.1063/1.1695690.

[33] D. VAN NGUYEN, D. GREBENKOV, D. LE BIHAN, AND J.-R. LI, Numerical study of a cylinder model of the diffusion mri signal for neuronal dendrite trees, Journal of Magnetic Resonance, 252 (2015), pp. 103–113.

[34] J. VERAART, D. NUNES, U. RUDRAPATNA, E. FIEREMANS, D. K. JONES, D. S. NOVIKOV, AND N. SHEMESH, Noninvasive quantification of axon radii using diffusion mri, Elife, 9 (2020), p. e49855.

[35] A. VRETBLAD, Fourier analysis and its applications, vol. 223, Springer Science & Business Media, 2003.

[36] D. WASSERMANN, D. V. NGUYEN, G. GALLARDO, J.-R. LI, W. CAI, AND V. MENON, Sensing Von Economo Neurons in the Insula with Multi-shell Diffusion MRI. International Society for Magnetic Resonance in Medicine, 2018, https://hal.inria.fr/hal-01807704. Poster.

[37] A. ZAIMI, M. WABARTHA, V. HERMAN, P.-L. ANTONSANTI, C. S. PERONE, AND J. COHEN-ADAD, AxonDeepSeg: automatic axon and myelin segmentation from microscopy data using convolutional neural networks, Scientific Reports, 8 (2018), p. 3816, https://doi.org/10.1038/s41598-018-22181-4, https://www.nature.com/articles/s41598-018-22181-4 (accessed 2022-04-22). Number: 1 Publisher: Nature Publishing Group.

[38] H. ZHANG, P. L. HUBBARD, G. J. PARKER, AND D. C. ALEXANDER, Axon diameter mapping in the presence of orientation dispersion with diffusion mri, NeuroImage, 56 (2011), pp. 1301–1315, http://www.sciencedirect.com/science/article/pii/S1053811911001376.