Introduction of diffusion simulations

Numerical simulations provide the flexibility to validate biophysical models and optimize NMR sequences under a wide range of acquisition parameters and tissue properties, thanks to the total control over phan- tom properties and the possibility of creating substrates mimicking specific biological tissue microstructure. Numerical simulations could be performed in several ways, based on (1) matrix formalisms, (2) finite- difference and finite-element methods, and (3) Monte Carlo (MC) simulations.

The matrix formalism relies on finding solutions for the eigenmodes of the diffusion propagator, and is computationally fast due to its implementation using matrix multiplications. Practically, these eigenmodes were only found analytically for diffusion within simple shapes. To extend the applicability of the matrix formalism for a complicated geometry, the eigenmodes of corresponding Laplace equations can potentially be estimated numerically.

To numerically solve the mesoscopic Bloch-Torrey equations for complex micro-geometries with permeable membranes, finite-difference methods were proposed to simulate the diffusion. Also, Bloch-Torrey equations can be solved numerically by finite-element methods for various geometries with impermeable or permeable membranes.

MC simulations are computationally expensive since they collect the ensemble average of numerous particles' random walks, resulting in statistically meaningful parameter estimates. Nonetheless, MC simulations are by far the most popular since they offer great flexibility to simulate diffusion along with other contrast mechanisms in all kinds of micro-geometries.

How to set up Monte Carlo simulations?

$\textbf{Microstructure geometry}$
To mimic a realistic microstructure geometry, we can arrange, pack and combine multiple simple shapes, in either 2d or 3d, with a length scale similar to the biological tissue of interest. For simulations including a non-porous medium, simulation results depend dramatically on the packing geometry. The size of a numerical phantom geometry should be much larger than the diffusion length (depending on the total simulated diffusion time $t$ and intrinsic diffusivity $D_0$).

$\textbf{Particle number }N_p$
Based on the central limit theorem, MC simulations necessitate including as many particles as possible to achieve an accurate parameter estimation of the ensemble average, i.e. error $\propto1/\sqrt{N_p}$. A clear cut-off for a minimally required number of particles has not been proposed so far, but a rule of thumb often used is that $\approx10^5-10^6$ particles diffusing from random initial positions generally provide reliable simulation results.

$\textbf{Step size}$
The step size is limited by the smallest microstructural length scales in a numerical phantom. The time for each step $\delta t$, as well as the step size
$$\delta x=\sqrt{2dD_0\delta t}\,,\quad\quad (1)$$
can be limited by the diffusion gradient waveform spectrum. For MC simulations of pulsed gradient, $\delta t$ is limited by the gradient pulse width $\delta$ via $\delta t\leq\delta$, whereas for oscillating gradient, $\delta t$ is limited by the gradient frequency $f$ via $\delta t\ll f$. For the simulation of water exchange, MT and surface relaxation, extra conditions related to the permeation probability need to be satisfied.

$\textbf{Step number }N_\text{step}$
The step number is proportional to the diffusion time. Empirically, the first few thousand steps provide unreliable results, caused by the discretization of a realistic continuous diffusion process into separate steps. This can be evaluated by simulating the kurtosis, which due to discretization, results in $K(t)\sim-1/N_\text{step}$, suggesting a minimal $N_\text{step}>1000$ for an error $<0.1$%.

$\textbf{Boundary condition (substrate edge)}$
A numerical phantom's finite boundary box, if treated inappropriately, may influence the simulation result. To further expand the finite numerical phantom, periodic or mirroring (hard wall) boundary conditions are often used.

$\textbf{Impermeable and permeable membranes}$
If a particle encounters an impermeable membrane within a step, the particle is reflected back by the membrane and experiences a perfectly elastic collision. The displacements before and after the collision are then summed up to the step size given by Eq. (1).

While encountering a permeable membrane, the particle has a probability $1-P_\text{EX}$ to be reflected elastically by the membrane, and a probability $P_\text{EX}$ to penetrate through, where
$$P_\text{EX}\simeq\frac{\kappa_0}{D_0}\delta x\cdot C_d\,,\quad\quad(2)$$
with permeability $\kappa_0$ and $C_d=1,\pi/4$, and $2/3$ for dimensionality $d=1,2,3$.

$\textbf{Particle-membrane interaction}$
The most commonly implemented interaction is elastic collision (specular reflection), which works well for impermeable membranes and for permeable and/or absorbing membranes. Other kinds of interaction include diffuse reflection and equal-step-length random leap for non-permeable membranes. For diffuse reflection, the particle reflected by a membrane consumes the rest of its step length in a random direction over the same side of the membrane. For equal-step-length random leap, the particle encountering a membrane cancels the original step and randomly chooses another direction to leap until the chosen step does not cross any membranes.

$\textbf{Relieve the computational bottleneck}$
For simulations within a medium composed of multiple compartments, the most computationally expensive calculation is to identify which compartment a particle resides in at a given time. A way to solve this problem is to build a lookup table. Such algorithm reduces the search to a local sub-volume scale and thereby significantly speeds up simulation run times.

$\textbf{Other contrast}$
Diffusion simulations can be performed by incorporating other MR contrasts, such as T1 and T2 relaxation, water exchange, magnetization transfer/surface relaxation, susceptibility effect/BOLD effect, and IVIM.

$\textbf{Simulation results}$
The simulated diffusion displacements $\Delta x$ yields the diffusivity and kurtosis along $\hat{n}$ via
$$D(t,\hat{n})=\frac{\langle(\Delta x\cdot\hat{n})^2\rangle}{2t}\,,$$
and
$$K(t,\hat{n})=\frac{\langle(\Delta x\cdot\hat{n})^4\rangle}{\langle(\Delta x\cdot\hat{n})^2\rangle^2}-3\,,$$
and the simulated trajectories $x(t)$ yields the diffusional signals by $S=\langle e^{i\phi}\rangle$, where $\phi=-\gamma\int g(t)\cdot x(t)dt$ with the gyromagnetic ratio $\gamma$ and diffusional gradient $g(t)$.

Oh, No! Common mistakes

$\textbf{Diffusion propagator: Sphere point picking}$
The single-step diffusion propagator is modeled by a binomial distribution in 1d, a circle in 2d, and a spherical surface in 3d. A potential mistake in 3d simulations is to sample the step orientation over $\theta=\pi v$ and $\phi=2\pi u$of a spherical coordinate, rather than over $\cos\theta=1-2v$ and $\phi=2\pi u$ with two uniformly distributed random numbers $u,v$ between 0 and 1.

$\textbf{Particle-membrane interaction and exchange}$
To speed up MC simulations, diffuse reflection and equal-step-length random leap are often used to replace the elastic collision. When the step number is sufficiently large and the step size is reasonable small (compared with the length scale in the microstructure), the two non-elastic interactions results agree with results from elastic collision with impermeable membranes. However, for permeable and absorbing membranes, the two non-elastic interactions could lead to bias in simulated diffusion metrics/signals due to the particle density re-distribution around the membranes.

$\textbf{Membrane permeation: permeability based}$
The original form of the permeation probability $P$ through a membrane of permeability $\kappa$ depends on the distance $\delta s$ between the random walker and the encountered membrane when $\delta s$ is smaller than the step size $\delta x$:
$$\frac{P}{1-P}\propto\frac{\kappa\delta s}{D_0}\,.$$
To reduce the computational load, we would like to approximate $\delta s$ with $\delta x$, by averaging over all possible $\delta s$ up to $\delta x$. This approximation leads to the probability in Eq. (2), which is applicable when $P\ll 1$ is satisfied, i.e.
$$\kappa_0\ll\sqrt{\frac{D_0}{2d\delta t}}\cdot\frac{1}{C_d}\,,$$
indicating that, for a large $\kappa_0$, a sufficient small time step $\delta t$ should be applied. In this case, $\kappa\approx\kappa_0$.

Actually, the choice of the small time step is surprisingly strict; for example, when $P\sim5$%, Eq. (2) is biased by roughly 5%.

$\textbf{Membrane permeation: diffusivity based}$
Instead of assigning a nominal permeability $\kappa$ for a permeable membrane, the permeation probability could be defined based on the spin concentration ($C_1, C_2$) and intrinsic diffusivity ($D_1, D_2$) over the left and right side of the membrane:
$$P_{1\to2}=\frac{C_2\sqrt{D_2}}{C_1\sqrt{D_1}}\,, \quad \quad P_{2\to1}=1\,.$$
It seems that this method introduces neither adjustable parameters for membrane permeability nor particle density transition (offset) over the membrane; however, this is true only for infinitely small time step. For a finite time step $\delta t$, an extra $\delta t$-dependent permeability $\kappa_0'$ may pronounce in simulations:
$$\kappa_0'=\sqrt{\frac{D_2}{2d\delta t}}\cdot\frac{1}{C_d}\,,$$
complicating the choice of time step, especially in 2d and 3d geometries.

Outlook: To do list

So far, most of the diffusion simulation studies, including non-MC studies, are implemented in geometries containing artificial geometrical shapes. Benefitting from recent advances in histology analysis techniques, it is possible now to perform simulations in realistic microstructure reconstructed from microscopy images of biological tissues. Furthermore, advances in computational resources, including the use of clusters and GPUs, enable ones to perform diffusion simulations along with complicated MR contrast within a reasonable amount of time. Here, we provide a to do list for the near future of this field:

1. Where is my extra-cellular space? Realistic microstructure generation
2. Fast and Fabulous: Parallel computation
3. Diffusion simulation and Bloch simulator in MRI system

Acknowledgements

I would like to thank Els Fieremans, Dmitry S Novikov, and Marco Palombo for the discussion.