Lattice Boltzmann method for simulation of diffusion magnetic resonance imaging physics in multiphase tissue models

We report an implementation of the lattice Boltzmann method (LBM) to integrate the Bloch-Torrey equation, which describes the evolution of the transverse magnetization vector and the fate of the signal of diffusion magnetic resonance imaging (dMRI). Motivated by the need to interpret dMRI experiments in biological tissues, and to offset the small time-step limitation of classical LBM, a hybrid LBM scheme is introduced and implemented to solve the Bloch-Torrey equation. A membrane boundary condition is presented which is able to accurately represent the effects of thin curvilinear membranes typically found in biological tissues. As implemented, the hybrid LBM scheme accommodates piece-wise uniform transport, dMRI parameters, periodic and mirroring outer boundary conditions, and finite membrane permeabilities on non-boundary-conforming inner boundaries. By comparing with analytical solutions of limiting cases, we demonstrate that the hybrid LBM scheme is more accurate than the classical LBM scheme. The proposed explicit LBM scheme maintains second-order spatial accuracy, stability, and first-order temporal accuracy for a wide range of parameters. The parallel implementation of the hybrid LBM code in a multi-CPU computer system, as well as on GPUs, is straightforward and efficient. Along with offering certain advantages over finite element or Monte Carlo schemes, the proposed hybrid LBM constitutes a flexible scheme that can by easily adapted to model more complex interfacial conditions and physics in heterogeneous multiphase tissue models and to accommodate sophisticated dMRI sequences.

Figure 1

Two-pulse Stejskal-Tanner PGSE sequence with rectangular bipolar gradient waveforms for diffusion MRI.

Figure 2

(a) Periodic parallel fiber model of heterogeneous tissue, with intracellular (in) and extracellular (ex) compartments and membrane of infinitesimal thickness. (b) Two-dimensional periodic computational domain and representative elementary volume (REV) of size.

Figure 3

Stencils for lattice Boltzmann scheme. (a) Two-dimensional, five-speed (D2Q5) and (b) three-dimensional, seven-speed (D3Q7).

Figure 4

The particle probability distribution functions near the membrane interface and relationship to the local lattice. Adapted from Li et al. [70].

Figure 5

Grid convergence rate and spatial truncation error of (a) classical and (b) hybrid LBM schemes, for the impulse response test case, for a 2D square domain ($L$ = 200 $\mu \text{m}$) with $\tau$ = 0.625, $b=1000{\mathrm{s}/\mathrm{mm}}^2, \delta$ = 3 ms, $\mathrm{\Delta }$ = 6 ms, and TE = 12 ms. The trend line marked with $n=2$ corresponds to quadratic convergence.

Figure 6

Comparison of spatial (a) and temporal (b) convergence rates for hybrid and classical LBM schemes for the uniform gradient test case at $t=20$ ms, with $G=23.16$ mT/m and $\tau$ = 0.60.

Figure 7

(a) Convergence of packed disk geometry and angled membrane geometry for different angles of $\phi$ and (b) schematic of angled membrane geometry.

Figure 8

Field map of 3D magnetization on a cross section at each grid point in the $z$ direction (exploded view). Due to the cylinder alignment and periodic boundary conditions there is no $z$ dependence in the signal, even though the gradient was applied obliquely.

Figure 9

(a) Comparison of LBM scheme with analytical solutions (dotted lines) for impermeable and permeable slabs and an impermeable disk for increasing $q$ values. (b) Comparison of the LBM scheme with the analytical solution of a periodic membrane from Ref. [21] for increasing gradient strength as the signal enters the localization regime wherein the signal transitions to a $\sim exp(-q^{1/3})$ dependence.

Figure 10

(a) Demonstration of the LBM scheme's ability to match the analytical short time solution (solid and dashed lines) for permeable and impermeable slabs as well as an impermeable disk when the signal exhibits a dependence on the surface-to-volume ratio. $D=2.3 \mu {\text{m}}^2/\text{ms}, q=50 {\text{mm}}^{-1},\delta$ = 50 $\mu \text{s}$, TE = $\mathrm{\Delta } +$ 1.0 ms, $\delta x$ = 0.5 $\mu \text{m}$, and $\delta t$ = 10 $\mu \text{s}$. (b) LBM scheme in the long time limit for a periodic geometry and a geometry with short-range disorder. In the long time limit, the LBM scheme matches the analytical solution presented in Ref. [25] for the short-range disorder as well as the periodic membrane solution given by Ref. [76], which is valid for all times. $D=2.3 \mu {\text{m}}^2/\text{ms},\kappa =50 \mu \text{m}/\text{s}, b=100 \text{s}/{\text{mm}}^2,\delta$ = 50 $\mu \text{s}$, TE = $\mathrm{\Delta }+\delta ,\delta x$ = 0.5 $\mu \text{m}$, and $\delta t$ = 10 $\mu \text{s}$.

Figure 11

Performance of hybrid LBM code after domain decomposition and parallelization with one-to-one mapping between domains, MPI processes and computer cores. (a) Speedup for $N=400$ and $N=1000$ with two to 48 cores. (b) Speedup (open symbols) and parallel efficiency (solid symbols) for $N=1000$ with two to 168 cores. Domain and timing parameters were $D=2.0 \mu {\text{m}}^2$/ms, $T_2$ = 100 ms, $q=40$ ${\text{mm}}^{-1}$, TE = 24 ms, $\mathrm{\Delta }$ = 20 ms, and $\delta$ = 4 ms, $\tau$ = 0.58, $\delta t$ = 0.010 ms, and $\delta x$ = 1.0 $\mu \text{m}$.

Figure 12

(a) Comparison of the convergence behavior for four versions of the membrane boundary condition. The full boundary condition, the boundary condition with ${\mathrm{\Delta }}_m=0.5$, with $\theta =0$ and with both ${\mathrm{\Delta }}_m=0.5$ and $\theta =0$. (b) L2 error between the membrane boundary conditions for $\theta =0$ (dashed lines) and both ${\mathrm{\Delta }}_m=0.5$ and $\theta =0$ (solid lines) for different angles and increasing $b$ values.

Figure 13

(a) Biphasic geometry domain of skeletal muscle cross section. (b) dMRI signal map (arbitrary units) of biphasic geometry at $t=\text{TE}$ with $\mathrm{\Delta }=100$ ms with the gradient applied in the horizontal direction. (c) Skeletonized geometry domain of skeletal muscle cross section. (d) dMRI signal map (arbitrary units) of skeletonized geometry at $t=\text{TE}$ with $\mathrm{\Delta }$=100 ms with the gradient applied in the horizontal direction. For both geometries, $\delta x$ = 0.333 $\mu \text{m}$ and $\delta t$ = 8.33 $\mu \text{s}$.

Figure 14

Time-dependent radial diffusivity of biphasic and skeletonized geometries with associated power-law fits demonstrating agreement in the fitted exponents to experimentally observed measurements [25]. Fitting of the power law was performed using all data points except the two shortest diffusion times.