Lattice Boltzmann method for simulation of diffusion magnetic resonance imaging physics in multiphase tissue models
Noel M. Naughton Caroline G. Tennyson John G. Georgiadis
 
	



Fast uncertainty quantification of reservoir simulation with variational U-Net
    Larry Jin,   Hannah Lu,  Gege Wen

Energy Resources Engineering

Stanford University

{zjin, hannahlu, gegewen}@stanford.edu
    December 30, 2023
================================================================================================================================


empty

	

	We report the first 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. 








§ INTRODUCTION


Random molecular motion in the presence of tailored magnetic field gradients imparts a phase dispersion in the nuclear spin transverse magnetization. The resulting signal loss has been employed to quantify the statistics of that motion and probe microscopic diffusion barriers in heterogeneous media <cit.>. Both diffusion-weighted nuclear magnetic resonance (NMR) and diffusion-weighted magnetic resonance imaging (dMRI) exploit this phenomenon to non-invasively probe the microscopic structure of porous media such as sedimentary rocks <cit.> and biological tissues <cit.>. In biological tissues, dMRI has successfully been used to sensitize the measured NMR signal to microstructural restrictions to the free diffusion of water within the tissue. This measurement has successfully been made in a variety of different tissues such as the brain <cit.>, skeletal muscle <cit.>, cardiac muscle <cit.>, breast tissue <cit.>, liver <cit.>, and cancerous tumors <cit.>. The present work is motivated by the need to interpret the signal measured during water diffusion through heterogeneous biological tissues in terms of the underlying microstructure.

The spin ensemble physics of dMRI are mathematically described by the Bloch-Torrey partial differential equation, which is a semi-classical model describing the evolution of the bulk magnetization of the spin ensemble in space and time <cit.>. The Bloch-Torrey equation is a linear diffusion-reaction equation, with a reaction term that is a function of space and time. This equation can fully accommodate dMRI physics by modeling the effect of externally applied magnetic field gradients (whose timing defines the dMRI sequence) and the diffusion and bulk flow of spins. Spin ensemble physics undergoing diffusion can also be described though the use of a diffusion propagator subjected to magnetic gradients <cit.>. Many reduced dMRI models have been employed based on limiting cases <cit.>, a-priori Brownian motion statistics <cit.>, or effective medium models <cit.>. Reviews of the different approaches to microstructure modeling in dMRI may be found in <cit.> with an emphasis on neural microstructure. 

Notable progress has been made in developing analytical models that describe the evolution of the signal <cit.>; however, our focus here is on numerical models of dMRI based on realistic representations of individual cell geometry and tissue microstructure, which are contained in tissue-based representative elementary volumes (REV). Such numerical models are useful in at least two ways: (i) simulating the signal of a dMRI experiment within a given microscopic reconstruction of the tissue microstructure (a tissue model), or (ii) developing and validating more accurate reduced-order dMRI models. In particular, numerical models can be used to simulate the evolution of the dMRI signal in complex tissue geometries, such as domains based on histological images of tissue <cit.>, for which no such analytical models exist. 

Previous numerical schemes for the solution of the Bloch-Torrey equation include Monte Carlo <cit.>, finite difference <cit.>, and finite element <cit.> methods, with Monte Carlo methods being the most widely employed. The majority of schemes employ a forward-Euler temporal discretization <cit.> of the Bloch-Torrey reaction term, which is first order in time. Higher order temporal discretization schemes using finite elements have been introduced recently, such as an explicit Runge-Kutta <cit.> and a second-order implicit scheme based on Crank-Nicolson <cit.>. Based on boundary-conforming finite elements, these schemes have definite advantages over finite differences in terms of describing complex boundaries. This advantage is shared with Monte Carlo methods, which are very simple to code but require careful optimization in order to run efficiently <cit.>.

The objective of the present work is to revisit and develop a novel simulation method for the numerical integration of the Bloch-Torrey equation in a specified tissue-based continuum REV with an applied linear gradient based on the lattice Boltzmann method (LBM). As a mesoscopic method based on the discrete Boltzmann equation, LBM is particularly efficient for simulating transport processes in complex heterogeneous biological tissue, whereby each lattice node can be assigned unique physics or transport properties. LBM is competitive relative to other computational methods because it involves uncomplicated algorithms, handles complex boundary conditions efficiently and accurately <cit.>, and is naturally amenable to parallelization <cit.>. Some of the challenges of LBM are the constraint it imposes on the time step (typical of explicit schemes) and the requirement to derive special boundary conditions for the probability distribution functions in order to preserve the consistency and accuracy of the numerical scheme <cit.>.

This work introduces a hybrid implementation of the LBM to integrate the Bloch-Torrey equation in heterogeneous tissue models, which carries the following advantages:


  * Obviates the problem of the classical LBM implementations, which require small temporal steps when applied to reaction-diffusion problems with dominant reaction terms.


  * Retains the second-order spatial accuracy in multi-compartmental domains containing complex permeable interfaces.
 

  * The numerical algorithm can be easily parallelized and executed efficiently with high parallel efficiency in multi-core computer systems. Based on spatial domain decomposition, the expectation is that the kinetic nature of the LBM and the locality of the operations involved result in execution times that scale linearly with the number of cores.

The overarching aim of this work is to support the claim that the proposed scheme is accurate, fast, and can accommodate complex geometries of relevance to more general tissue models.




§ METHODS



 §.§ Diffusion-weighted imaging


The governing equation describing hydrogen proton (1H) spin dynamics in the presence of diffusion during an MRI experiment is the Bloch-Torrey equation <cit.>.  Neglecting coherent (advective) fluid transport, the Bloch-Torrey differential equation can be formulated in a coordinate frame rotating at a fixed Larmor frequency (determined by the MR scanner permanent magnetic field) as follows: 

    ∂M/∂ t = 
    -i γ [ x·G(t)] M - 
    M/T_2 ( x ) + 
    ∇· (D ( x ) ∇M ) ;   M ( x ,t)= 
    (  M ( x ,t) ) + 
    i  (  M ( x ,t) ) ,
			 
where M(x,t) is a complex variable representing the bulk (transverse) magnetization of the spins, i is the imaginary unit, γ is the gyromagnetic ratio for 1H, x is the spin position vector, G(t) is the time-varying magnetic field gradient vector used to encode diffusion, T_2 is the spin-spin relaxation time, and D is the diffusion coefficient. Except for γ, all variables listed above are local, in that they represent the ensemble average of spin behavior at a given spatial location. The problem described by Eq. (<ref>) is supplemented with an initial condition M(x,0) and appropriate boundary conditions, which will be discussed in section <ref>. Unless explicitly stated, the initial condition throughout is M(x,0) = 1+i0. Decomposing the transverse magnetization M(x,t), as shown in Eq. (<ref>), the Bloch-Torrey equation yields two coupled reaction-diffusion equations for (M(x,t)) and (M(x,t)), respectively. The coupling occurs through the first term of the right-hand side of Eq. (<ref>), which depends on the specific dMRI sequence. 




A typical dMRI sequence is the Stejskal-Tanner Pulse-Gradient-Spin-Echo (PGSE) sequence <cit.>. Although more sophisticated sequences are in use, we will employ PGSE here since it is adequate to represent MR physics and study local diffusion without image formation. This sequence involves a bipolar magnetic gradient pulse (diffusion gradients), with the gradient vector G(t) controlled by the operator, cf. Figure <ref>. The resulting evolution of M(x,t) generates a time-varying magnetic flux which constitutes the dMRI signal and is acquired upon the appearance of a “spin echo” at time t = TE. Assuming spatially uniform spin density, the dMRI signal S is obtained by integrating |M| over the domain, here defined as a representative elementary volume (REV). The gradient magnitude is typically constant in space but varies in time, so it is convenient to express it in separable form, G(t)=G_0  f(t). By judiciously choosing a set of vectors G_0, each oriented along a specific, non-collinear direction, the signal can be sensitized to probe the dynamics of diffusion along these directions. In addition to the gradient orientation, three parameters describe the PGSE sequence: the diffusion gradient amplitude |G_0|, the gradient pulse duration δ, and the delay Δ between the gradient pulses, cf. Figure <ref>. These parameters can be grouped to define a diffusion decay factor, b=γ^2 |G_0|^2δ^2(Δ - δ/3). Additionally, in dMRI it can be useful to define a parameter q as q=(γG_0 δ /2π)^2. By taking measurements with multiple q-vectors, the average diffusion propagator can be reconstructed <cit.>.



Here heterogeneous tissues are considered to be fibrous inclusions encased in an extracellular matrix. The fibers are surrounded by thin permeable membranes cf. Figure <ref>. Such domains are commonly found in biological tissues and much effort has been devoted to understanding their influence on the dMRI signal in order to use dMRI as a probe of the tissue structure. While such work is important to the larger goal of relating the dMRI signal to the underlying tissue microstructure, the focus of this paper is firmly on the forward problem of solving the Bloch-Torrey equation in tissues with thin permeable membranes. Here the Bloch-Torrey equation is solved in a representative elementary volume (REV) containing intracellular (in) and extracellular (ex) subdomains, cf. Figure <ref>. We will consider problems with isotropic diffusion and piece-wise uniform T_2 and D. Referring to the two subdomains in Figure <ref>, for example, there are two diffusion coefficients, D_in and D_ex, for intracellular and extracellular compartments, respectively. This property notation will be suppressed in the following, until it is explicitly reinstated. 



 §.§ The lattice Boltzmann method


Historically, numerical methods of analyzing dMRI physics have been split between particle-tracking based Monte Carlo methods and continuum-based finite difference and finite element methods. Lattice Boltzmann methods (LBM) are mesoscale methods that exist between microscopic Monte Carlo methods, which consider the dynamics of individual particles and macroscopic methods like finite elements, which directly discretize the continuum-based Bloch-Torrey equation. In contrast, the LBM is based on simplified kinetic models that incorporate the necessary microscopic physics to allow the averaged properties to obey desired macroscopic equations. Setting aside the Bloch-Torrey equation for a moment, LBM solves a discretized version of the Boltzmann distribution through consideration of a discrete-velocity distribution function (g_i) that describes the distribution of particle velocities in the different lattice directions <cit.>. This discretization of the Boltzmann equation leads to the lattice Boltzmann equation, 

    g_i(x + e_i·δ t, t+δ t) = g_i(x,t) + Ω_i(x,t) ,

which describes how the particles g_i(x,t) move to the neighboring node g_i(x + e_i·δ t, t+δ t) with velocity e_i after being influenced by the collision operator Ω_i(x,t), which models the collision and subsequent redistribution of fictitious particles at each node. Through proper selection of Ω_i(x,t), often by adopting the Bhatnagar-Gross-Krook (BGK) form of the collision operator <cit.>, it is possible to recover the Navier-Stokes equations of fluid mechanics <cit.>, leading to the LBM's success in modeling a variety of different fluid mechanical domains <cit.>. In the case of the Bloch-Torrey equation, there is no bulk fluid velocity, and the Bloch-Torrey equation can be viewed as a reaction-diffusion equation, for which LBM schemes have previously been presented <cit.> (see Appendix A). Reviews of the LBM for fluid flow may be found in <cit.> while reviews focused on LBM solutions of heat and mass transfer problems, which have a similar formulation to the Bloch-Torrey equation, are available in <cit.>.

We note that that although traditional methods of modeling dMRI often consider the evolution of a diffusion propagator, the connection between such a propagator and the discrete-velocity distribution function discussed here is a superficial one. LBM considers a fictitious particle distribution, which is a numerical technique that enables solutions to both reaction-diffusion equations as well as problems that include advection, in which case the lack of a diffusion propagator analogue is clear. There is no direct connection to the actual diffusing spin packets, though such a connection may be an interesting future avenue of inquiry. 





 §.§ Hybrid lattice Boltzmann method


An order of magnitude analysis of the terms of the right hand side of Eq. (<ref>) is presented in Appendix A. The phase of the magnetization vector M(x,t) exhibits fast oscillations when the term x·G(t) becomes large. As the truncation error analysis in Appendix A shows, the implementation of the classical LBM method to solve the reaction-diffusion equation (<ref>) introduces a truncation error term that grows with the square of the REV length size. To remove this dependence on domain size, a hybrid lattice Boltzmann method involving the factorization of the operator in terms of a reaction (slow) and a diffusion (fast) operator is introduced here where, for each timestep,

    M(x,t) = exp( -i γ [ x·G(t)] δ t^' - 1/T_2δ t^') M^' (x,t) ,

with M^'(x,t) as an intermediate function. The exponent in Eq. (<ref>), with δ t^' denoting the reaction time step (employed in the discretization of the reaction term), has been reported first in <cit.> and has since been used in most schemes to integrate the Bloch-Torrey equation <cit.>. As shown in Appendix A, this functional form of the exponent is appropriate for G(t) piece-wise constant in time, like in the case of the PGSE gradient pulse sequence, cf. Figure <ref>. For sequences with gradient pulses of different time-dependence, treatment can be generalized. Eq. (<ref>) is recovered from Eq. (<ref>), accurate to first order in δ t^' (Appendix A), if M^' obeys the following diffusion equation

    ∂M^'/∂ t= ∇· (D ∇M^') .




Eq. (<ref>) is then integrated with the classical LBM algorithm over a diffusion time step δ t, as shown below. The proposed hybrid LBM scheme for the integration of Eq. (<ref>) is essentially a time splitting scheme: 

    . [       Diffusion split during [t,t+δ t]:  M(x,t)→M(x,t+δ t)  from Eq. <ref>,;                               then, initialize M^' (x,t) = M(x,t+δ t) , and; Reaction split during [t,t+δ t^']:  M^' (x,t) →M(x,t+δ t^')  from Eq. <ref> ]}

As written, the practical implementation of the time splitting scheme (<ref>) requires that the diffusion and reaction time steps are identical (δ t=δ t^'), but this is not necessary. In fact, one of the advantages of separating diffusion and reaction steps is that different time steps or time-splitting schemes can be used depending on the stiffness of the Eqs. (<ref>) - (<ref>) <cit.>. For example, using δ t=k δ t^', with an integer k > 1, solve Eq. (<ref>) for 1 step, and then solve Eq. (<ref>) for k steps so that the timing is consistent. As the analysis in Appendix A indicates, δ t=δ t^' is a choice that is consistent with the range of physical and numerical parameters pertinent to the present work. In the following, when we refer to the time step, δ t=δ t^' is assumed. 

The integration of Eq. (<ref>) during the diffusion split is performed with the classical LBM algorithm presented in Eq. <ref>. The most common version of this algorithm is based on a single relaxation parameter (Bhatnagar–Gross–Krook model) and can be expressed as: 

    g_i(x + e_i·δ t, t+δ t) - g_i(x,t) = -1/τ[g_i(x,t)-g_i^eq(x,t)]

where g_i is the particle probability distribution function defined on a discrete lattice, i denotes the lattice direction, x is the space coordinate on the lattice, δ x is the lattice spacing (grid size), e_i is the lattice (speed) vector, δ t is the diffusion time step, g_i^eq is the discretization of the fully-mixed equilibrium state for g_i, and τ is the dimensionless relaxation time. In this case, like the magnetization M, the function g_i is a complex variable (g_i = (  g_i  + i  (  g_i )). Since the zero-th moment of g_i is equal to M, the magnetization vector components are recovered by taking the sum of these functions over the lattice directions. Following the Chapman-Enskog analysis of Eq. (<ref>), the diffusion equation (<ref>) can be recovered, accurate to O(δ t, δ x^2) if the relaxation time parameter τ is defined as 

    τ = 1/2 + δ t/ε_D(δ x)^2D ,
 
where ε_D is a positive constant related to the weighting factors ω_i. Because advection is neglected, the equilibrium distribution function for the LBM scheme (irrespective of whether the reaction term is included) is given by

    g_i^eq(x,t)=ω_i M(x,t) .

Here, ω_i are weighting factors chosen to allow Eq. (<ref>) to be recovered <cit.>. To be consistent with the LBM lattice topology, ω_i must satisfy the isotropy and symmetry conditions

    ∑_i ω_i e_i = 0;   ∑_i ω_i e_iae_ib= ε_D δ_ab;   and  ∑_i ω_i = 1 ,

where e_ia and e_ib are the spatial components of e_i <cit.>.

For clarity, the theoretical development reported in this section is confined to 2-D isotropic diffusion, thus a 2-D square lattice, 5-speed model (D2Q5), as shown in Figure <ref>, has sufficient symmetries for a consistent spatial discretization of Eq. (<ref>). In fact, an analysis of the 2D advection-diffusion equation <cit.> indicates that the D2Q5 stencil produces more accurate and robust results than D2Q9, which is the 9-speed stencil. The extension of this scheme to 3D is straight-forward, as is demonstrated in Section <ref>, so is not considered in detail here.  
 
The D2Q5 lattice speed vectors in Eq. (<ref>) are given by

    e_i =
    {[           (0,0)           (i=0); (± 1,0),(0,± 1)   (i = 1,2,3,4) ]. .

For a D2Q5 lattice, ε_D=1/3, and the weighting factors for the equilibrium distribution are

    ω_i =
    {[           ε_D         (i=0);         ε_D/2 (i = 1,2,3,4) ]. .

Under the classical LBM scheme, the evolution equation, Eq. (<ref>), is integrated in two steps over δ t, a collision step (g_i →ĝ_i) followed by a streaming step (ĝ_i →g̅_i). The collision step is 

    ĝ_i(x,t) = g_i(x,t)-1/τ[g_i(x,t)-g_i^eq(x,t)] ,

where g_i is the initial particle distribution at the beginning of the time step, g_i^eq is the equilibrium particle distribution Eq. (<ref>), and ĝ_i is the particle distribution function following the collision step, which is the input to the streaming step. The streaming step is  

    g̅_i (x+δ x e_i ,t+δ t) = ĝ_i(x,t) .

From here the reaction step is initialized as g_i^'(x_n,t^k) = g̅_i (x_n,t^k+1), where t^k=k δ t^'. The magnetization vector M^'(x,t) is computed by the zero-th moment of the particle probability distribution function

    M^'(x,t)= ∑_ig_i^'(x,t).
 
As proven in <cit.>, the above LBM scheme is unconditionally stable for τ>1/2, which is always satisfied given that the diffusion coefficient D in Eq. (<ref>) is positive.	

Turning to the reaction step, the integration of Eq. (<ref>) during the reaction split accounts for the effects of the diffusion gradient pulse and T_2 relaxation on the magnetization. We start with a discretized version of Eq. (<ref>)

    M(x_n ,t^k+δ t^')=exp(-i γ [x_n ·G(t^k)]δ t^')exp(-δ t^'/T_2)M^'(x_n ,t^k)

where x_n denotes the coordinate location, and  G(t^k)=G_0 f(t^k) is the gradient vector at time t^k. By substituting Eq. (<ref>) into Eq. (<ref>), the distribution function after the completion of the reaction step at t^k+1=t^k+δ t^' becomes	

    g_i(x_n ,t^k+1)=exp(-i γ[x_n·G_0f(t^k)]δ t^')exp(-δ t^'/T_2)g_i^'(x_n,t^k) .

The hybrid LBM scheme is summarized below in terms of the sequence of the particle distribution functions (g_i→ĝ_i→g̅_i→g_i^') computed at each step

    . [                                   Collision at t^k:                          g_i(x_n,t^k) →ĝ_i(x_n,t^k)                                          Eq. <ref>;;                                   Streaming at t^k:              ĝ_i(x_n,t^k)→g̅_i (x_n+δ x e_i ,t^k+1)                                          Eq. <ref>;; 2cthen, initialize g_i^'(x_n,t^k)=g̅_i (x_n,t^k+1).;                                    Reaction at t^k:                    g_i^'(x_n,t^k) →g_i(x_n  ,t^k+1)                                           Eq. <ref> ]}.


The phase in the first term on the right hand side of Eq. (<ref>), which is a function of space and time, couples the computation of the distribution functions. Assume, for example, that a diffusion gradient G_0 is applied in the (x,y) plane, and there are [N+1]×[N+1] lattice nodes in the 2D REV shown in Figure <ref>. The phase in Eq. (<ref>) at time t^k and location x_n becomes

    -Δφ^k_n = γ [ x_n ·G_0]f(t^k)δ t^',    with  n∈ [0,1,2,… N] × [0,1,2,… N] .

In this implementation of the LBM code, the distribution function g_i is then updated in terms of its real (g_i) and imaginary (g_i) components, as follows

    ( g_i(x_n,t^k+1)) = [ ( g_i(x_n,t^k)) cos(Δφ^k_n) - ( g_i(x_n,t^k)) sin(Δφ^k_n)] exp(-δ t^'/T_2) 
    
    ( g_i(x_n,t^k+1)) = [ ( g_i(x_n,t^k)) cos(Δφ^k_n) + ( g_i(x_n,t^k)) sin(Δφ^k_n)] exp(-δ t^'/T_2)

For simplicity, the bold font notation for g_i and the other sequence members is henceforth suppressed. 



 §.§ Boundary Conditions


In the following, we show how the boundary conditions for the presented LBM scheme are expressed in terms of the distribution functions for external boundaries, which are periodic, and internal boundaries consisting of the mathematical membranes separating the “in” and “ex” subdomains.



  §.§.§ Membrane boundary condition




Biological cells are delineated by thin semi-permeable membranes that are weakly diamagnetic, which means that they do not disturb the magnetization significantly. The boundary condition at a thin permeable membrane involves the conservation of mass flux of water spins carrying the scalar M, without any loss in the membrane. Assuming equal spin density on both sides of the membrane implies conservation of magnetization flux across the membrane. Letting n denote the unit vector normal to the membrane and pointing towards the extra-cellular space, and introducing the membrane permeability κ, this conservation principle imposes the following boundary condition  

    D_exn·∇M_ex=D_inn·∇M_in= κ(M_ex-M_in) ,

where M_in and M_ex denote the values on the intra- and extra-cellular sides of the membrane, respectively.  This interfacial condition, which is a mixed boundary condition, needs to be reformulated in terms of the particle distribution functions in order to be integrated in the LBM scheme. 

The second-order Dirichlet and Neumann boundary conditions presented by Li et al. <cit.> are used in this work. These boundary conditions are based on the idea of `bounce back' from the membrane. Additionally, their use of spatial interpolation allows preservation of the membrane geometry and their application to curved geometries. These boundary conditions are presented in detail in <cit.>. Li et al. <cit.> derived second-order accurate Dirichlet and Neumann boundary conditions for general curved boundaries, and later Li et al. <cit.> extended these boundary conditions to develop an interfacial treatment for conjugate heat and mass transfer. Guo et al. <cit.> further developed boundary conditions for jumps in concentration or flux at the interface. The extension presented here is the case of the membrane boundary condition, Eq. (<ref>), used in place of the continuity equation (M_in=M_ex). 

The variable distance from the lattice point in the intracellular region to the point at which the membrane cuts the lattice link is denoted by Δ_m δ x. So, Δ_m expresses the dimensionless distance between the internal node closest to the membrane and the membrane (0<Δ_m<1). Here we keep the LBM collision → streaming nomenclature with the distributions denoted by g_i→ĝ_i→ g_i^'. Each membrane boundary condition is enforced at the end of the LBM collision step, so extra subscripts are necessary to distinguish particle distributions based on the direction the particles move. Figure <ref> illustrates the nomenclature used in this section. Four lattice velocities are defined (e_α, e_α̅, e_β, and e_β̅), with α and β denoting lattice velocities moving in the direction of the intracellular to extracellular domain while α̅ and β̅ denote lattice velocities in the opposite direction.

A detailed derivation of this reformulation is presented in Appendix B. The final equations of the membrane boundary condition are

    g^'_α̅(x_i,t)=     
    A_1^iĝ_α(x_i,t)+
    A_2^iĝ_α(x_ii,t)+
    A_3^iĝ_α̅(x_i,t)+
       
    B_1^iĝ_α̅(x_e,t)+ 
    B_2^iĝ_α̅(x_ee,t)+ 
    B_3^iĝ_α(x_e,t)+ 
       
    C_1^iĝ_β(x^'_i,t)+ 
    C_2^iĝ_β(x^'_ii,t)+ 
    C_3^iĝ_β̅(x^'_i,t)+
       
    D_1^iĝ_β̅(x^'_e,t)+ 
    D_2^iĝ_β̅(x^'_ee,t)+ 
    D_3^iĝ_β(x^'_e,t)

and

    g^'_α(x_e,t)=    
    A_1^eĝ_α̅(x_e,t)+A_2^eĝ_α̅(x_ee,t)+A_3^eĝ_α(x_e,t)+
       
    B_1^eĝ_α(x_i,t)+ 
    B_2^eĝ_α(x_ii,t)+B_3^eĝ_α̅(x_i,t)+
       
    C_1^eĝ_β̅(x^'_e,t)+ 
    C_2^eĝ_β̅(x^'_ee,t)+ 
    C_3^eĝ_β(x^'_e,t)+
       
    D_1^eĝ_β(x^'_i,t)+ 
    D_2^eĝ_β(x^'_ii,t)+ 
    D_3^eĝ_β̅(x^'_i,t)  .

where the coefficients are given in Eqs. (<ref>) and (<ref>). It is noted that in the case of infinite permeability, these boundary conditions match those presented in Li et al. <cit.> for conjugate heat and mass transfer. If the membrane is impermeable, the boundary conditions reduce to those originally presented by Li et al. <cit.>.




  §.§.§ Half-lattice link membrane boundary condition

While the above boundary condition is valid for general membrane geometries, in the case of a straight membrane that is perpendicular to the lattice direction (θ = 0), then Φ_n,α̅ = Φ_n,in and Φ_n,α = Φ_n,ex, yielding a simplified version of the coefficients 

    A_j^i= (c_d4c_d4^∗ c_njϵ_Dδ x+(c_d4^∗ c_n4 c_dj+c_n4^∗ c_d4 c_nj)κ δ t)/F 
       B_j^i= (c_d4c_n4 (c_nj^∗-c_dj^∗) κ δ t )/F 
       C_j^i =  D_j^i =  0

and

    A_j^e= (c_d4c_d4^∗ c_nj^∗ϵ_Dδ x+(c_d4c_n4^∗ c_dj^∗+c_n4 c_d4^∗ c_nj^∗ )κ δ t)/F 
       B_j^e= (c_d4^∗ c_n4^∗(c_nj -c_dj) κ δ t)/F 
       C_j^e =  D_j^e =  0

with 

    F = c_d4 c_d4^∗ϵ_Dδ x+(c_d4 c_n4^∗+c_n4c_d4^∗)κ δ t .

In such a simplified case, it is worthwhile to place the membrane at the half-way point between nodes, setting Δ=0.5. In this case, and for the coefficients chosen in Eq. (<ref>) and Eq. (<ref>), Eqs. (<ref>) - (<ref>) further reduce to  

    A_1^i = A_1^e = P/1+P  and   
    B_1^i = B_1^e = 1/1+P

with 

    P=ϵ_D/2κδ x/δ t

and all other coefficients equal to zero, allowing the boundary condition to be expressed as

    g_α̅^'(x_i,t)= 1/1+P ĝ_α̅(x_e,t)+P/1+P ĝ_α(x_i,t)
    
       g_α^'(x_e,t)= P/1+P ĝ_α(x_e,t)+1/1+P ĝ_α̅(x_i,t)
    .

The boundary condition is applied after the collision step of the LBM algorithm and in place of the streaming step. It is expressed only in terms of distribution functions at nodes x_i and x_e, reducing the complexity and computational cost of the boundary condition as only the nearest neighboring nodes are necessary and the coefficients (P/1+P and 1/1+P) are constants that can be precomputed to increase efficiency. We conclude this section with a physical interpretation of the factors in Eq. (<ref>) by considering the particle distribution functions involved. The presence of the membrane splits the population of the particles moving towards the membrane (from either side) into a portion 1/1+P that cross and a portion P/1+P that is reflected back. In the limit of infinite permeability (κ→∞  and P→0), Eq. (<ref>) reduces to Eq. (<ref>). Conversely, for the limit of impermeability (κ→0  and P→∞), Eq. (<ref>) reduces to the standard bounce-back condition for a homogeneous Neumann boundary <cit.>. 
	


  §.§.§ Modified periodic boundary condition

	
	The typical method to terminate the prescribed external boundary conditions for the REV in 2D is to consider an infinite periodic solution domain exhibiting a spatial translation symmetry along x and y, as shown in Figure <ref>. Given the the signal phase's linear spatial dependence, owing to the dMRI gradient term [x·G_0 f(t)] in Eq. (<ref>), conventional periodic conditions do not apply. Let us consider a spatial period L (length of the REV), and the two (vertical) boundaries marked “Left” and “Right” in Figure <ref>. As demonstrated in <cit.>, the magnetization on these boundaries obeys the following constraint
	
    M(x_Left,t)=exp[i φ(L,t)]M(x_Right,t);  
       φ(L,t)=γ [(x_Right-x_Left)·G_0]∫_0^tf(t^')dt^'.

	This means that the periodic boundary condition results in a phase difference that is proportional to the component of G_0 that is perpendicular to the boundaries, the spatial period L, and an integral factor that varies with time. If G_0 is applied only along x, then the periodic boundary condition along the two (horizontal) boundaries marked “Up” and “Down” in Figure <ref> reduces to the conventional form 
	
    M(x_Down,t)=M(x_Up,t) .

	To further simplify the presentation, we set t^k=k δ t, and keep only the x-coordinate dependence below, with x_n=n δ x, where n∈[0,1,2, … N]. Since the spatial period is L=N δ x, the phase in Eq. (<ref>) can be discretized by approximating the integral by a sum (low order approximation is consistent with Eq. (<ref>))
	
    φ(L,t^k)=γ L |G_0|∑_m=0^k-1f(t^m) δ t.

	
	As with the membrane boundary condition, all external boundary conditions are applied between the collision and streaming step. The modified periodic boundary condition has been adapted to the LBM scheme as follows. Two external “buffer” lattice rows are introduced at n∈[-1,N+1] in order to complete the streaming step at n=0 and n=N. The following assignments are applied to these rows after the collision step, 
	
    ĝ_i(x_-1,t^k)=exp[i φ(L,t^k)] ĝ_i(x_N-1,t^k) 
    	
       ĝ_i(x_N+1,t^k)=exp[-i φ(L,t^k)] ĝ_i(x_1,t^k) .

	After the streaming, Eq. (<ref>), and reaction initialization step, the correct phase difference, Eq. (<ref>), is maintained for the distributions at the nodes on the “Left” and “Right” boundaries at n=0 and n=N, respectively,
	
    g_i^'(x_0,t^k)=exp[i φ(L,t^k)] g_i^'(x_N,t^k) .

	During the reaction step, each distribution in Eq. (<ref>) gains phase according to Eqs. (<ref>) and (<ref>), 
	
    g_i(x_0,t^k+1)=exp(i Δφ^k_0)exp(-δ t^' / T_2)  g_i^'(x_N,t^k) 
    	
       g_i(x_N,t^k+1)=exp(i Δφ^k_N)exp(-δ t^' / T_2)  g_i^'(x_0,t^k) .

	Considering Eqs. (<ref>) and (<ref>), it is straightforward to show that
	
    exp(-i Δφ^k_0)exp[i φ(L,t^k)]exp(i Δφ^k_N) = 
    	exp[i φ(L,t^k+1)].

	By incorporating Eq. (<ref>), Eqs. (<ref>) and (<ref>) yield 
	
    g_i(x_0,t^k+1)=exp[i φ(L,t^k+1)] g_i(x_N,t^k+1) ,

	which is consistent with the modified periodic boundary condition Eq. (<ref>). This implies that, by making the assignments from Eq. (<ref>) to the nodes on buffer rows after the collision step, the correct phase difference, Eq. (<ref>), is maintained for the distributions at the appropriate boundary nodes and at the completion of each time step. For the “Up” and “Down” boundaries, the conventional periodic condition given by Eq. (<ref>) is satisfied if we make the following assignment after the collision step: 
	
    ĝ_i(x_Down,t^k)= ĝ_i(x_Up,t^k) .

	
	

  §.§.§ Mirroring boundary condition

	The modified periodic boundary condition allows implementation of the periodic boundary condition, however, it can also be adapted to implement a mirroring boundary condition <cit.>. A mirroring boundary condition reflects the domain across a boundary, effectively doubling the analyzed REV. Such a boundary condition is particularly useful when considering non-uniform geometries that have cells cross the REV's boundaries, such as when considering domains derived from tissue micrographs. The mirroring boundary condition avoids the possibility of geometrical discontinuities without having to manually edit the image to make both sides of the domain agree, as is necessary if a periodic boundary condition is imposed <cit.>. In the mirroring boundary condition, the buffer node (x_m) geometry is equal to the geometry of the node on the boundary of the domain (x_n). The mirroring boundary condition is similar to the modified boundary condition described in Eq. (<ref>), however, there are two notable differences. The first relates to the exchange of the lattice velocities. Because the node is mirrored instead of translated, the lattice directions are also mirrored. Recalling the notation used for the membrane boundary condition, this means that the lattice velocities at the buffer node in direction e_i are computed using the lattice velocities in direction e_ from the source nodes. The second notable aspect of the mirroring boundary condition is related to this swapping of the lattice velocities. The lattice velocity distribution at each node is influenced by the gradient direction. Under the mirroring boundary condition, the mirrored node is effectively subject to a gradient that is also mirrored, and thus in a different direction than the gradient direction at all other nodes in the domain. To account for this, Eq. (<ref>) is modified to `unwind' the magnetization at the boundary back to the origin (x_o), defined as the location where G_0 ·x=0,
	
    g_i(x_o,t)=
    	g_(x_n,t)
    	exp(i γ (x_n-x_o)·G_0 ∫_0^tf(t^')dt^').

	Here, the 0^th order moment of the lattice velocity distribution M corresponds to the expected value if no diffusion-sensitizing gradient were gradient. However, the lattice velocity distribution will still exhibit the effect of the gradient direction. If the gradient direction is perpendicular to the boundary's edge, the direction of the gradient for the mirrored boundary node is a reflection of the original gradient direction. In this case, the complex conjugate (g_^∗) describes the lattice velocity distribution in the gradient direction necessary to describe the mirrored node. Thus, the magnetization can be `rewound' to the buffer node, 
	
    g_i(x_m,t)=
    	g^∗_(x_o,t)
    	exp(i γ (x_o-x_m)·G_0 ∫_0^tf(t^')dt^').

	Eqs. (<ref>) and (<ref>) can be combined to describe the mirroring boundary condition when applied after the collision step as
	
    (ĝ_i(x_m,t)) =    [ ĝ_(x_n,t) exp( i γ( x_n+x_m-2x_o) ·G_0 ∫_0^t) ] 
    (ĝ_i(x_m,t)) =    [ ĝ^∗_(x_n,t) exp( i γ( x_n+x_m-2x_o) ·G_0 ∫_0^t) ] .

	In the case of the gradient applied parallel to the boundary edge, no reflection of the gradient direction is necessary, leading to the boundary condition being described by
	
    ĝ_i(x_m,t) =  ĝ_(x_n,t) exp( i γ( x_n-x_m ) ·G_0 ∫_0^t) .

	For the more general case of the gradient not aligning with one of the boundary edges, a more complicated rotation of g_i(x_o,t) is necessary, which is not considered here. 
	


 §.§ Parallelization of LBM hybrid scheme

Considering the length scales probed in a typical dMRI experiment (∼1-100  μ m) relative to the typical size of a dMRI voxel  (∼1  mm^3), the simulation of the entire magnetization field in a voxel would result in a very large computational problem. Although the quasi-periodic structure of many tissues can be exploited to reduce the size of this problem somewhat, it is still necessary to develop an efficient code adapted to LBM to solve this problem directly. To address this, a parallel implementation of the hybrid LBM scheme based on domain decomposition <cit.> and a Fortran code with message passage interface (MPI) is presented. We employ a 3D version of the hybrid LBM code (D3Q7 stencil) on an REV that is a thin rectangular prism with grid size N × N × 1. The REV is decomposed into multiple non-overlapping domains that are assigned to separate MPI processes. In this case, the REV is partitioned into strips, and each strip is handled by one MPI process, which in turn is assigned to one core. For every time step, each core executes the hybrid LBM algorithm over the assigned domain and exchanges boundary information with adjacent strips while the two strips at the edges of the REV exchange boundary information with each other subject to the applied external boundary conditions. Two quantities must be minimized to maximize parallelization efficiency: the number of neighboring domains each individual MPI process communicates with and the amount of information passed between each domain. Here we focused on minimizing the number of neighboring domains, while holding the length of the boundary between strips fixed. The performance of the parallel implementation is quantified in terms of the following ratios (p denotes the number of cores)

    Speedup = Exceution time with 2 cores/Exceution time with p cores; 
        Parallel efficiency = Speedup/p.

Due to the particulars of the MPI implementation, it is not possible to run the code with only one core, hence the definition of speedup ratios in Eq. (<ref>).





§ RESULTS

In this section, we first analyze the accuracy and convergence of the proposed hybrid LBM scheme. This is done by first comparing the hybrid LBM scheme, summarized by Eq. (<ref>), to that of the classical LBM scheme, which is described by replacing Eq. (<ref>) with Eq. (<ref>). All computations involve the numerical integration of the Bloch-Torrey equation (Eq. <ref>) to simulate the evolution of dMRI signal under the PGSE sequence without imaging gradients (Figure <ref>). These results are then compared with analytical solutions of the Bloch-Torrey equation. We then assess the accuracy of the hybrid LBM scheme applied to solve the Bloch-Torrey equation as a function of spatial resolution and conclude by examining the accuracy of the proposed membrane boundary conditions. 

Following this analysis, we demonstrate the ability of the hybrid LBM scheme to match solutions of the Bloch-Torrey equation in a number of limiting cases for which analytical solutions exist. We also demonstrate the ability of the hybrid LBM scheme to efficiently scale in a parallel implementation of the computer code as well as the straightforward extension of the scheme to 3D. We conclude this section with an analysis of the error introduced by making various assumption about the orientation and location of the membrane and provide a demonstration of the utility of the hybrid LBM scheme in analyzing complex multiphase domains, such as those typical in biological domains as represented by micrographs. Based on the discussion following Eq. (<ref>), we set the time step for the reaction operator identical to that used for the diffusion operator, δ t^'= δ t. 





 §.§ Comparison of truncation error for classical and hybrid LBM schemes


Here, we study the truncation error relative to analytical solutions of Eq. (<ref>) obtained first with an impulse initial condition and second with a uniform initial condition. The predictions of the classical and hybrid LBM schemes were compared to the exact solution of the Bloch-Torrey equation in a periodic domain with uniform diffusion coefficient and subject to an initial Dirac delta distribution M(x,t=0)=δ(x) <cit.>. A square REV with L=200 μ m was used, and the remaining physical parameters were b = 1000 s/mm2, δ = 3 ms, Δ = 6 ms and TE = 12 ms. The LBM simulations were performed on a N× N grid for N = 50, 100, 200, and 400 lattice points in each direction, and τ = 0.625 was kept fixed. The normalized L2 norms of the errors are plotted in Figure <ref> for four time points: t_1 = δ/2 = 1.5 ms, t_2 = (Δ + δ)/2 = 4.5 ms, t_3 = Δ + δ/2 = 7.5 ms, and t_4 = TE = 12 ms. A trend line corresponding to L2 error ∼δ x^n for n = 2 is also included for reference. The results indicate that the spatial convergence of both schemes is second order in space, with the exception of the classical scheme at t_4. In that case, the scheme reaches the asymptotic convergence regime only for grid sizes smaller than δ x = 2 μ m.
 
In order to highlight the difference between the classical and hybrid LBM scheme as the size L of the domain increases, the numerical solution of the above homogeneous problem was repeated with a uniform initial condition M(x,0)=1+i0. For this case, the solution is trivial: M(x,t)=exp(-bD)  exp(-t/T_2). The L2 norms of the errors are plotted in Figure <ref> as a function of lattice spacing and time step. For the classical method, the error is given for two domain sizes (L=20 μ m, and L=100 μ m) to illustrate how the size of the domain affects accuracy. This effect is not manifested for the hybrid scheme. Although both schemes are second order in space, the accuracy of the hybrid scheme is independent of L and higher than that of the classical scheme. Figure <ref> indicates that both schemes are first-order accurate in time, as per the discussion in Appendix A. 

To analyze the accuracy of the membrane boundary condition, a periodic domain with a permeable membrane at an angle φ was constructed as illustrated in Figure <ref>. The domain was defined with a distance between the membranes being approximately 40 μ m. In order to allow a direct analysis of the convergence of the scheme, this distance was slightly adjusted such that the number of nodes in each direction of the REV was an integer, thus allowing consistent refinement of the domain with no other changes to the REV geometry. The hybrid LBM scheme results were compared with the analytical short pulse approximation solutions for periodic, permeable membranes <cit.>. To allow comparison with the short pulse approximation, a gradient duration of 1 μ s was used, or, if such a gradient duration was less than the prescribed time step, δ = δ t was used instead. Along with the straight, angled membrane, a packed disk domain similar to Figure <ref> with a diameter of 40 μ m and intracellular volume fraction of 0.50 was also simulated to examine the convergence of the scheme for curved boundaries. There is currently no analytical solution for such a domain so the scheme was compared with a highly refined solution. Results for both the angled membrane and the packed disked are shown in Figure <ref>. For both domains, simulation parameters were D = 2.3 μ m^2/ms, κ = 50 μ m/s, Δ = 20 ms, TE = 25 ms and b-value = 1000 s/mm^2. The results for the angled membrane domain are second order for all angles while the packed disk is slightly less than second order.




Given that the first term in the factor  exp(-i γ[x·G(t)] δ t^'-(T_2)^-1δ t^') in Eq. (<ref>) changes only the phase of M, while the second term decays in time, the numerical stability of the hybrid LBM scheme is controlled by the stability of the classical LBM scheme for Eq. (<ref>). LBM is unconditionally stable for the diffusion equation <cit.>, however, this conclusion does not cover the effect of the membrane boundary condition on the stability of the numerical scheme. In their stability analysis of the Dirichlet and Neumann boundary conditions used to develop the membrane boundary condition, Li et al. <cit.> show that if -1<c_d1<1, the boundary conditions are stable. Further, they show that the lower bound of stability can be extended below -1 depending on the chosen relaxation coefficient τ and Δ_m. In the implementation presented here c_d1 = 2(Δ_m-1), so the lowest possible value of c_d1 occurs for Δ_m=0 when c_d1=-2. Due to the coupling between domains, this value of c_d1 must be in the stability regimes for both Δ_m=0 and Δ^∗_m=1. In this case, the boundary conditions will be stable if τ≳  0.6. In the opposite direction, the largest possible value of c_d1 is c_d1=1 at Δ_m=1, which is within the original stability regime and so always stable. For the simplified case of Δ_m fixed at Δ_m=0.5 then c_d1=2(Δ-1)=-1, so the LBM scheme with membrane boundary conditions located halfway between boundaries is always stable, however, for membrane boundaries not located halfway between lattice nodes, there is a restriction on τ in order to maintain numerical stability.



 §.§ Three-dimensional version of hybrid LBM scheme


Previous results thus far have been confined to 2D using a D2Q5 stencil, however, the extension of the hybrid LBM scheme to 3D is straightforward. To demonstrate, a D3Q7 stencil (Figure <ref>) was used to simulate the dMRI signal by solving a 3D periodic array of permeable cylinders with circular cross section whose axis is aligned with the z-coordinate (Figure <ref>a). The cylinder diameter is 55 μ m and the packing fraction is 0.65, while the physical parameters are D_in = 1.5 μ m^2/ms, D_ex = 2 μ m^2/ms, T_2,in = 30 ms, T_2,ex = 10 ms and κ = 10 μ m/s. The REV is a rectangular prism with grid size 55 × 55 × 5. The temporal and spatial steps are δ t = 0.025 ms, and δ x = 1.0 μ m. Modified periodic boundary conditions were implemented on all external boundaries, and a PGSE sequence with b = 1000 s/mm2, TE = 24 ms, Δ = 20 ms, and δ = 4 ms was used. The gradient G_0 was applied along an oblique direction with directional cosines (1/√(3),1/√(3),1/√(3)). The field map of the magnetization at t = TE is shown in Figure <ref>. Because the cylinder axis is aligned with the z-direction and the boundary conditions are periodic, there is no z-dependence in the signal, even though the gradient G_0 has a nonzero z-component. This spatial symmetry of the result is a consequence of the symmetry of the inner boundaries and outer boundary conditions and does not constitute a limitation of the general 3D implementation of the LBM scheme.






 §.§ Comparison with analytical solutions of the Bloch-Torrey equation





There are a number of analytical solutions to the Bloch-Torrey equation that involve permeable membranes which the hybrid LBM scheme can be compared with. Here we present four benchmarks where the hybrid LBM scheme is compared with known analytical solutions, demonstrating the LBM scheme's ability to accurately match such analytical solutions in a variety of different cases. The hybrid scheme is first compared with the analytical solution for a periodic slab geometry with both permeable and impermeable membranes <cit.> as well as for an impermeable cylinder <cit.>. All three cases consider the effect of increasing gradient strength and assume the short gradient pulse approximation (δ→ 0). Simulations were performed for diameters of 5.0 μ m with D = 2.3 μ m^2/ms, Δ = 100 ms, and TE = 110 ms. For the permeable slab case, κ = 50 μ m/s. For the disk, δ x = 0.1 μ m, δ t = 0.333 μ s, and δ = 1.0 μ s. For the slabs, δ x = 0.1 μ m, δ t = 0.714 μ s, and δ = 5.0 μ s. Figure <ref> shows that the hybrid scheme is able to successfully match the phase cancellations of the signal, leading to the observed diffraction patterns for all three cases <cit.>. 

For strong gradients, the signal enters a so-called localization region where the signal departs from the Gaussian behavior of the signal and instead illustrates a ∼exp(-q^1/3) dependence on the gradient strength <cit.>. To verify that the LBM scheme is able to recreate this behavior, simulations for an impermeable slab with a diameter of 160 μ m were performed with D = 2.3 μ m^2/ms, Δ = 60 ms, δ = 60 ms, TE = 120 ms, δ x = 0.5 μ m, and δ t = 12.2 μ s. Figure <ref> shows the hybrid LBM scheme correctly replicates the transition to the localization regime as the gradient increases. 

Solutions of Bloch-Torrey equation demonstrates time-dependent behavior in both the long- and short-time limits. In the short-time limit (Δ << D/2a), the signal demonstrates a dependence on the surface-to-volume ratio that the hybrid LBM scheme is able to match for both a impermeable and permeable (κ=50 μ m/s) slab with a diameter of 10 μ m <cit.> as well as an impermeable 10 μ m diameter disk <cit.>. The signal behavior was examined for Δ between 0.1 and 5.2 ms. The effective diffusion coefficient was computed using the low b-value representation of the diffusion coefficient: D_eff = -ln(S)/b <cit.>.

In the long time limit, dMRI signal exhibits a power law dependence t^-γ, with the exponent related to the organization of the membranes <cit.>. For periodic membranes the effective diffusion coefficient exhibits a t^-1 dependence while for short-range disorder the signal exhibits a t^-1/2 dependence. Slabs with both periodic membranes and short-range disorder were simulated with an average diameter of 10 μ m allowing comparison of the LBM scheme with the 1D results from <cit.>. Figure <ref> shows that the hybrid LBM scheme is able to accurately exhibit the expected long-time behavior for both domains.



 §.§ Parallelization of hybrid LBM scheme




To investigate speedup due to parallelization, a parallel code was implemented using Fortran 90 with Intel’s IFORT v14.0.2 compiler and MVAPICH2 v2.1. Simulations were run on SDSC's Comet cluster <cit.>, which consists of 1944 nodes with 2 x 12 core CPU processors (Intel Xeon E5-2680 v3 2.5 Ghz), and 128 GB DDR4 DRAM running CentOS 6.7. Two homogeneous REVs with square cross-sections and side lengths of 0.4 mm and 1.0 mm were simulated. These domains correspond to N × N × 1 grids with sizes N = 400 and N = 1000, respectively. LBM algorithms require high memory throughput and, as such, are often limited by memory bandwidth <cit.>. To examine the scaling of the algorithm beyond this known limitation, simulations were preformed using up to 48 cores and adjusting the number of cores per node so the entire problem could be held in cache when possible. All simulations were repeated 5 times and the average execution time was used to determine scaling performance. The results are shown in Figure <ref> and indicate that for N=1000 the speedup is linear over the considered range. For N=400, there is linear speedup up to 24 cores before the performance begins to degrade. For the current partitioning of the REV, the cost of message passing between cores scales as ∼N while the number of operations per core (p) scales as ∼N^2/p. Thus, the cost of message passing relative to operation count per core scales as ∼p/N, implying that performance should degrade at a lower number of cores for coarser grids, explaining the performance degradation for N=400. Simulations were also performed for N=1000 using up to 168 cores. The results, plotted in Figure <ref>, indicate that the code exhibits ideal (linear) speedup, with a minimum parallel efficiency of 77% as defined in Eq. (<ref>).






 §.§ Effect of simplified membrane boundary condition assumptions


While Figure <ref> demonstrates the second order accuracy of the proposed membrane boundary condition, it is worthwhile to consider the accuracy of simplified versions of the membrane boundary condition, such as when Δ_m=0.5 or θ =0. These boundary conditions are both easier to implement as well as less computationally expensive. If these simplified versions provide sufficient accuracy, it may be preferable to use them, particularly for complex domains when computing Δ_m and θ is not straightforward. To begin, the convergence of a 40 μ m diameter packed disk with an intracellular volume fraction of 0.50 is examined for four different versions of the membrane boundary condition. They are 1) the full membrane boundary condition, 2) the boundary condition with Δ_m=0.5, 3) with θ=0, and 4) with both Δ_m = 0.5 and θ=0. A convergence study was performed with D = 2.3 μ m^2/ms, κ = 50 μ m/s, Δ = 20 ms, δ = 5 ms, TE = 30 ms and b-value = 1000 s/mm^2. Results are shown in Figure <ref>. The membrane boundary condition for Δ_m=0.5 converges with the full boundary condition as the grid is refined while the two membrane boundary conditions that assume θ=0 demonstrate zeroth-order accuracy as they converge to a different value than the full boundary condition, though the result is within 2% of the result from the full membrane boundary condition. These results indicate that the assumption of θ = 0 is the more limiting of the two assumptions. 

To better understand the error introduced by the assumption of θ=0, the angled domain of Figure <ref> was reexamined using simplified membrane boundary conditions that assume either θ=0 or both θ=0 and Δ_m=0.5. Multiple angles φ were examined for increasing b-values. The angle φ was varied between 0 and π/2 for b-values between 100 and 2000 s/mm^2. Simulation parameters were the same as for the packed disk as well as δ x = 0.5 μ m and δ t = 12.25 μ s. These results were  compared with the full membrane boundary condition to quantify the L2 error introduced by these assumptions. Figure <ref> shows that the error increases with gradient strength (b-value). It should be noted that the magnitude of the dMRI signal decreases with b-value, and the error for all b-values was always less than 1% of the original signal value at t=0. The error is greatest for angles of φ=π/4, which, due to symmetries of the geometry, is the greatest deviation possible from θ=0. 

The combined results of Figure <ref> suggest that the assumption of θ=0 should be avoided when possible. However, avoiding such an assumption is not always possible. Particularly when dealing with complex domains, such as those patterned off of realistic biological tissues, it may not be straightforward to compute the local angle θ. While Figure <ref> demonstrates that for high b-values this assumption will lead to errors in the simulated dMRI signal, for b-values <1000 s/mm^2 the error introduced by this assumption is on the order of <4% and <1% for b-values <500 s/mm^2. Considering typical SNR values of dMRI measurements are often in the range of 20-50, this suggests that, for b-values <1000 s/mm^2, the error introduced by the assumption of θ=0 will be less than the noise in the dMRI signal. Further, any comparisons of simulations with realistic tissue structures will require some approximations of the tissue shape that will also introduce errors. Here, the assumption of θ=0 greatly simplifies the analysis as computing θ based on images in non-trivial. Further work is necessary to better understand how the error introduced from the assumption of θ=0 is influenced by changes in the b-value or small variations in the tissue structure, however, for low to moderate b-values, these errors do not necessarily preclude its use in analyzing complex tissue geometries.



 §.§ Simulations of histology-informed domains







One area where dMRI has found success in measuring tissue microstructure is in skeletal muscle <cit.>. Additionally, these measurements are often made with relatively low b-values, identifying it as an area where the simplified half-link membrane boundary condition can be applied to analyze complex, curved domains derived from tissue micrographs of skeletal muscle cross-sections. Here a micrograph was obtained from digital images available in the open literature <cit.>. To derive realistic REVs, morphological image processing of the micrograph was performed. This involved thresholding with ImageJ <cit.> and segmenting manually to produce a biphasic geometry consisting of extra- and intra-cellular domains (Figure <ref>). Skeletal muscle is tightly packed and often assumed to fill the entire domain. To analyze this case, another domain was created. Using a custom Matlab script, a watershed transform was performed on this image to dilate each cell so that the extracellular space was reduced to a skeletonized outline with the interface between the cells defined by a single permeable membrane (Figure <ref>). 

LBM simulations of the dMRI signal in the REVs of Figure <ref> and Figure <ref> were performed using the hybrid LBM scheme over a uniform, non-boundary conforming grid for diffusion times in the range 10-1000 ms. Membrane boundary conditions using the half-link membrane boundary condition were applied at the interfaces between cells while mirror boundary conditions were applied at the edges of the domains. Other dMRI sequence parameters were b = 400 s/mm^2, δ = 10 ms, and TE = Δ+δ. For long diffusion times, a stimulated echo (STEAM) sequence is often employed. To approximate this sequence, a generalized diffusion-weighted sequence, which emulates the STEAM sequence by ignoring the effects of T_1 and T_2 relaxation <cit.>, was employed. For the biphasic domain, the intracellular and extracellular diffusion coefficients were equal (D_in=D_ex = 1.5 μ m^2/ms) while the skeletonized domain only had one diffusion coefficient (D = 1.5 μ m^2/ms). The membrane permeability was set to κ = 50 μ m/s. We note that if one desired to directly compare these two domains, one would need to adjust the membrane permeability in the skeletonized domain to account for two membranes sandwiched together.  Mirroring boundary condition's were applied on all sides. Figures <ref> and <ref> give the corresponding field maps for Δ = 100 ms.

Novikov et al. <cit.> showed that the radial diffusivity of skeletal muscle, that is, the average diffusion coefficient in the cross-sectional plane, demonstrates a diffusion time dependence that is ∼t^-1/2. Fitting a power law to the LBM scheme's results shows that both the biphasic and skeletonized domain exhibit time dependence that is consistent with this experimentally observed behavior (Figure <ref>). This match with experimentally observed results suggests that the hybrid LBM scheme may be a useful tool to examine how the dMRI signal evolves in biologically realistic domains and under different assumptions, For example, it can facilitate the quantification of the effect of the extracellular space on the signal.  




§ DISCUSSION


	Reporting on the first implementation of the lattice Boltzmann method to solve the Bloch-Torrey equation (<ref>), we proposed the hybrid LBM scheme summarized by Eq. (<ref>). The fundamental difference with the classical version of LBM, which is given by replacing Eq. (<ref>) by Eq. (<ref>), is the factorization of the operator, Eq. (<ref>), which splits the reaction and diffusion temporal discretization. This splitting addresses the stiffness of the problem, which is characterized by the disparity between reaction and diffusion time scales, as discussed in the first paragraph of Appendix A. Based on comparisons of two cases with analytical solutions, we demonstrated that the hybrid scheme is more accurate than the classical LBM scheme. Both schemes are second-order accurate in space and first-order in time. When  τ and D are kept constant, Eq. (<ref>) implies that δ t ∼δ x^2. In other words, we need to decrease δ t while maintaining δ t/δ x^2 constant in order to increase the approximation accuracy. As Eq. (<ref>) indicates, the time step limitation for the classical LBM is more restrictive than that for the hybrid LBM, which is given by Eq. (<ref>). This implies that the accurate integration of Eq. (<ref>) with classical LBM requires an order of magnitude smaller time step than with the hybrid scheme for L ∼ 100 μ m, with the concomitant requirement that the grid size has to decrease by two orders of magnitude (δ x^2 ∼δ t). This result explains why the truncation error of the classical scheme is higher for the larger domain size, cf. Figure <ref>. 

	The hybrid LBM method shares the same clear advantage as LBM methods in terms of using uniform grids to discretize curved boundaries between various compartments in the REV while retaining second-order spatial accuracy and stability. Consistent with the kinetic nature of the LBM scheme, the membrane boundary conditions (Eqs. (<ref>) and (<ref>)) connect the particle distribution functions on either side of the interface directly with the membrane permeability κ and avoid the need to approximate transmembrane derivatives, as is the case with finite difference or finite element schemes. The full membrane boundary condition is capable of accurately maintaining the second-order spatial accuracy of the scheme when considering the effects of curvilinear boundaries that intersect the lattice at arbitrary angles. Additionally, the presented membrane boundary condition is valid in other heat and mass transfer conditions as well as when advection is considered. Simplifying the membrane boundary condition to assume that the membrane is perpendicular to the lattice (θ=0) introduces error that increases for increasing gradient strength. While this assumption introduces error into the simulation, its simplicity of implementation as well as not requiring the angle of the membrane-lattice intersection make it appealing for use with more complicated tissue geometries such as those derived from histological images. 

	Like the classical LBM scheme, the hybrid version based on the time splitting scheme described by Eqs. (<ref>) and (<ref>) is unconditionally stable for τ>0.6 (because both time splits are stable), so the step sizes are determined by approximation accuracy considerations. Concerning the choice of δ t and δ t^', and in view of Eq. (<ref>), let us consider the limitations placed on the diffusion time step δ t for realistic values of the diffusion coefficients in tissue. For typical values τ∼0.6 and D∼2 μ m^2/ms, Eq. (<ref>) requires δ t/δ x^2∼0.02 ms/μ m^2. For a spatial resolution of δ x∼1 μ m in an REV with L∼100 μ m, this requires δ t∼0.02 ms. Given that this diffusion step also satisfies the requirement for the reaction step in Eq. (<ref>), we have employed δ t=δ t^' in the present study. Formal optimization of the hybrid LBM would involve a study of the variation of the truncation error as a function of τ, like in <cit.>, and an investigation of varying the sequence or the step size of the diffusion and reaction splits (δ t≠δ t^') <cit.>, but both are outside the scope of the present study. 

	While improved explicit <cit.> and implicit <cit.> temporal integration schemes have been proposed for the solution of the Bloch-Torrey equation with the finite element method, the disadvantage of LBM schemes relative to higher-order temporal schemes is offset by the amenability of the former to parallelization. Parallel computing is where LBM schemes have a performance advantage over finite difference or finite element schemes. Our LBM code for solving the Bloch-Torrey equation on a 1000 × 1000 × 1 grid exhibits a parallel efficiency of 77% at 168 cores, vs. an efficiency of ∼55% for the finite element scheme with ∼3.5 million degrees of freedom at 256 cores <cit.>. Our speedup is optimal (linear) for the whole range, cf. Figure <ref>. Our parallelization algorithm relies on an uncomplicated domain decomposition scheme and one-to-one mapping of the MPI processes to CPU cores. Further gains in performance are anticipated by accounting for special computer architecture or by employing GPU cores <cit.>.

	In terms of memory allocation, the LBM scheme requires a total allocation of  44× N^3 for the numerical integration of the Bloch-Torrey equation for N× N × N lattices (in 3D). Based on these estimates we can describe the complexity of the hybrid LBM scheme as 44× N_t× N^3, which is ∼10^10, given discretization parameters discussed previously. Analysis of the computing performance of Monte Carlo methods to integrate the Bloch-Torrey equation in two-compartment tissue models morphologically similar to the ones used in this study indicates that a minimum complexity of ∼10^9 is needed in order to avoid statistical error <cit.>. Moreover, Yeh et al. <cit.> repeat each simulation 10 times, thus raising this minimum to ∼10^10. Consideration of the effect of thin cell membranes on the computing performance increases the complexity significantly with lattice-free Monte Carlo methods being required to adequately model the effect of curvilinear membranes. It is here that the advantages of a LBM scheme become evident in comparison to more widely adopted Monte Carlo methods. By being able to accurately resolve the effects of curvilinear membranes, LBM allows retention of the computational efficiently gained from using a structured grid. Further, the imposition of the external periodic conditions is straightforward for LBM schemes, in contrast with the finite element scheme <cit.> where they had to be approximated by introducing an artificial permeability to mimic diffusion at the external boundaries. A comprehensive comparison of the relative advantages of each of these three approaches will be the subject of future work.

	The hybrid LBM scheme is able to accurately match analytical solutions of the Bloch-Torrey equation in both the short and long-time limit as well as for increasing gradient strength. It is also able to match experimentally observed diffusion-time dependence when applied to a domain derived from histological images of skeletal muscle. The extension from 2D (D2Q5 stencil) to 3D (D3Q7 stencil) is straightforward owing to the simplicity of the spatial discretization and implementation of the boundary conditions. These results demonstrate the ability of the LBM scheme to be applied in a variety of cases. In particular, its ability to straightforwardly integrate a histologically-derived domain illustrates the flexibility of the LBM scheme to incorporate complex tissue domains. Further, the hybrid LBM scheme is not limited to PGSE sequences or linear gradients. It is capable of simulating arbitrary pulse sequences as well as non-linear magnetic gradients by modifying the term in brackets in Eq. <ref>, though the inclusion of non-linear gradients would introduce a domain size dependence in the error term of the hybrid LBM scheme (see Eq. (<ref>)). Future applications of the hybrid-LBM scheme in the field of dMRI include investigation and development of reduced analytical models <cit.>, analysis of increasing complex tissue structures, particularly those derived from histological images, and even analysis of the inverse problem of estimating microstructure properties from dMRI measurements <cit.>. 

	We conclude with several comments regarding possible extensions of the hybrid LBM scheme developed here. The presented hybrid LBM scheme is very general and can accommodate other transport phenomena, multiple tissue compartments, and other MRI sequences. The scheme can be readily extended to incorporate more complex physics and be applied in the study of a variety of biological tissues such as neural tissue, cardiac and skeletal muscle, liver, and cancer tumors. It can also accommodate more complex dMRI gradient waveforms and sequences involving other MRI contrast mechanisms (e.g. perfusion, magnetic susceptibility, elastography, etc.), or imaging gradients (slice selection, phase encoding, or readout). Starting with its first application for modeling unrestricted diffusion <cit.>, LBM has since accommodated anisotropic diffusion and advection <cit.>, coupled diffusion <cit.>, coupled reaction-diffusion between multiple species <cit.>, finite cell membrane permeability <cit.>, phase field models <cit.>, and interstitial flow <cit.>. Such processes are pertinent to biophysics problems involving transport and evolution of large biomolecules in blood-perfused cellular systems as well as physics other than diffusion. For example, LBM was applied to model protein diffusion inside mammalian cells <cit.> and to model cardiac electrophysiology <cit.>. As an additional indication of its versatility, LBM has been recently employed to integrate fractional order diffusion-advection-reaction equations <cit.> and has also shown promise in incorporating fluid-structure interactions <cit.>. This is not to say that other numerical methods could not have been employed for the phenomena mentioned above, but rather, owing to the local nature of the operations, the same LBM formulation can be easily adapted to accurately and efficiently simulate a vast range of physical phenomena. The combination of the numerical accuracy, efficiency, and ability to incorporate additional physical phenomena make lattice Boltzmann schemes an exciting alternative to currently used Monte Carlo and finite element based schemes in the simulation of diffusion-weighted MRI.




§ CONCLUSIONS


	Motivated by the need to interpret the dMRI signal from biological tissue, we have laid the foundation and performed the analysis of a hybrid implementation of the LBM to integrate the Bloch-Torrey equation in heterogeneous tissue models. In its current implementation, the hybrid LBM scheme accommodates finite membrane permeability, piece-wise uniform diffusion coefficients, a wide range of dMRI parameters, periodic and mirroring boundary conditions, and interphase conditions accounting for flux continuity. By splitting the reaction and diffusion time steps, the algorithm maintains the explicit nature of the (classical) LBM implementation. We have shown via truncation error analysis and numerical tests that this splitting obviates the requirement of small temporal steps introduced by the strong reaction term in the Bloch-Torrey equation. 

	Another attractive feature of the classical LBM scheme is also maintained here: the phase boundaries are discretized on uniform 2D and 3D lattices while still maintaining the ability to accurately solve for the effects of curvilinear permeable membranes located at arbitrary angles to the lattice. We have shown that the hybrid scheme retains second-order spatial accuracy and stability for a wide range of membrane orientations and typical dMRI parameter values. Further, the presented hybrid LBM scheme accurately replicates the behavior of analytical solutions in a variety of limiting cases, illustrating the robustness of the presented scheme. Our results indicate that the associated LBM code is very compact and can be easily parallelized and executed efficiently on a general multi-core computer with a excellent scaling for up to 168 CPU cores and a parallel efficiency above 77%. Additionally, the LBM scheme is based on a uniform grid mesh, which, when combined with the efficient parallelization of the scheme, allows for straightforward application of the scheme to large, realistic tissue structures such as those derived from histological images. 

	Overall, the proposed lattice Boltzmann method present an exciting development in the numerical simulation of diffusion-weighted nuclear magnetic resonance physics. LBM schemes allow accurate treatments of general curvilinear membranes while retaining the computational efficiency and advantages of a structured grid-based scheme. Additionally, the ability of LBM schemes to incorporate a wide variety of additional physical phenomena such as advection, susceptibility, fluid-structure interaction illustrate the flexible and extensible nature of these schemes. 




§ ACKNOWLEDGMENTS

	We are grateful to Mr A. Z. Wang for histology image processing. We acknowledge the financial support by NSF (grants CBET-1236451 & CMMI-1437113, and a Graduate Research Fellowship to NMN) and NIH (grants HL090455 and EB018107).  Partial support by the R.A. Pritzker chair fund is also acknowledged. Parallel computations were performed on the Extreme Science and Engineering Discovery Environment’s (XSEDE) Comet cluster located at the San Diego Supercomputer Center, which is supported by National Science Foundation grant number ACI-1548562. The contributions of two anonymous reviewers are also gratefully acknowledged.
	
APPENDICES 




§ APPENDIX A


The first part of this appendix presents a truncation error analysis of the time-splitting method implemented in the hybrid LBM scheme. There are two methods to derive the macroscopic equation from the evolution of the particle probability distribution function: multiple time scales (Chapman-Enskog expansion), and asymptotic analysis <cit.>; here we employ a combination of them. 

First, we present a scaling analysis of the Bloch-Torrey differential equation (<ref>), rewritten for a general gradient pulse G(t)=G_0 f(t) and for piece-wise uniform diffusion coefficients as follows

    ∂M/∂ t =  -i γ[x·G_0 f (t)]M  -  M/T_2  +   D∇·∇M        γδ x |G_0|∼ 10^-2/ms        1/T_2∼ 10^-2/ms     D/δ x^2∼1/ms.

The above equation is a homogeneous reaction-diffusion differential equation so the relative order of magnitude of the various terms does not depend on the magnitude of M. Eq. (<ref>) is defined in t∈[0,TE] and x∈ REV. Using a δ x of ∼1 μ m, typical whole-body MRI scanner parameters, and typical properties of biological tissue (T_2 = 100 ms and D = 1.0 μ m^2/ms), an order of magnitude analysis of the terms in the right hand side of Eq. (<ref>) reveals the disparity between the reaction (first and second term) and diffusion (third term) time scales. The reaction rates are more than two orders of magnitude slower than the diffusion rate. This is the motivation for the splitting scheme associated with the hybrid LBM, which is a concept that has been explored in prior studies of such models <cit.>. We also note that as δ x decreases, the disparity between the diffusion and reaction rates increases, allowing the difference in scales to be maintained even for large gradient strengths provided a sufficiently small δ x.

Second, we can rename the linear reaction operator in Eq. (<ref>) as follows

    ∂M/∂ t=R(x,t) M + D∇·∇M ,  with  R(x,t)=-i γ[x·G_0 f(t)]-1/T_2) ,

and generalize Eq. (<ref>) to express the evolution of magnetization starting past a reference time instant t^k ∈[0,t_E]

    M(x,t)=exp[∫_t^k^tR(x,t^'') dt^'']M^'(x,t)   for   t≥ t^k .

Note that M(x,t^k)=M^'(x,t^k). Differentiate Eq. (<ref>) with respect to time to obtain

    ∂/∂ tM(x,t)=R(x,t)M(x,t)+ E(x,t^k;δ t^')∂/∂ tM^'(x,t)

where

    E(x,t^k;δ t^')=exp[∫_t^k^t^k+δ t^'R(x,t^'') dt^''],   and  δ t^'≥ 0 .

Let us suppress all dependent variables everywhere except in the expression E(x,t^k;δ t^'), and require that M^' obeys the diffusion equation, within a certain truncation error TE_D

    ∂/∂ tM^'=D∇·∇M^'+TE_D

so Eq. (<ref>) becomes

    ∂/∂ tM=R  M+   E(x,t^k;δ t^')  [ D∇·∇M^'+TE_D] .

Using Eq. (<ref>), we can show that 

    E(x,t^k;δ t^') ∇·∇M^'= ∇·∇M+2iγ F(t^k;δ t^') G_0·∇M - γ^2F^2(t^k;δ t^') |G_0|^2 M

where 

    F(t^k;δ t^')=∫_t^k^t^k+δ t^'f(t")dt" .

This allows casting Eq. (<ref>) in the form of Eq. (<ref>) 

    ∂/∂ tM=  R M+ D∇·∇M+[2iγ D F(t^k;δ t^') G_0·∇M- γ^2D F^2(t^k;δ t^') |G_0|^2M]   
    
    + [E(x,t^k;δ t^') TE_D] .

The terms contained in the two square brackets constitute the truncation error of the hybrid LBM scheme proposed here. The terms in the first bracket containing Eq. (<ref>) correspond to the error introduced in the treatment of the reaction part of Eq. (<ref>) according to Eq. (<ref>), and their magnitude can be estimated by assessing the magnitude of the integral in Eq. (<ref>). By recognizing that f(t)∼ O(1), we can easily see that F(t^k;δ t^')∼ O(δ t^'),  so that the formal order of magnitude of the two reaction truncation error terms is γ D|G_0| δ t^' and  γ^2D|G_0|^2 δ t^'^2. In order for Eq. (<ref>) to be consistent with Eq. (<ref>), both these terms have to be much smaller than the smallest term in Eq. (<ref>), which is the reaction term according to the order of magnitude analysis of Eq. (<ref>) (for clarity, we will only consider the part of the reaction term related to the gradient, for ultra-short T_2, its effect can be straightforwardly incorporated). This requirement implies D  δ t^' / δ x^2≪ 1 and γ|G_0| D  δt^'^2 /δ x ≪ 1, leading to constraints on the time step of 

    δ t^'≪δ x^2/D ∼ 1  ms
       δ t^'≪√(δ x/γ|G_0| D) ∼ 10  ms .


Since the diffusion problem Eq. (<ref>) is integrated with the classical LBM scheme, we can estimate TE_D from a truncation error analysis of that scheme. This error can be obtained by modifying, according to our Eqs. (<ref>) and (<ref>), the expression (A23) obtained by the Chapman-Enskog expansion in the Appendix of <cit.>, and by separating the effect of the diffusion time step δ t from the lattice grid size δ x

    TE_D=3 δ t  τ^2-τ+1/6/τ-1/2 ∂^2M^'/∂ t^2 + h.o.t. = 3ϵ_D δ x^2/D(τ^2-τ+1/6) ∂^2M^'/∂ t^2 + h.o.t.

The expressions in Eq. (<ref>), which are equivalent via the use of Eq. (<ref>), recover the known fact that the truncation error of the classical LBM is first-order in time and second-order in space. Returning to the last term in Eq. (<ref>), we can see that this is also the contribution of TE_D to the hybrid LBM error if we require that E(x,t^k;δ t^')∼ O(1). Referring to Eqs. (<ref>) and (<ref>), we note that the diffusion gradient term in E(x,t^k;δ t^') is periodic and so always O(1), so the leading contribution to the overall truncation error of the diffusion term is O(δ t,δ x^2), and of the reaction term is O(δ t^'). 

 
The second part of this Appendix addresses the truncation error of the classical LBM scheme applied in the solution of the Bloch-Torrey (<ref>), i.e. without the time splitting scheme Eq. (<ref>)). This involves a modification of the collision step of the LBM scheme as described in section <ref>, which unlike Eq. (<ref>), now reads

    ĝ_i(x,t)=g_i(x,t)-1/τ[g_i(x,t)-g_i^eq(x,t)]+δ t ω_i R(x,t) M(x,t)

where M(x,t) is computed by summing over g_i, as shown in Eq. (<ref>). Including the reaction term in Eq. (<ref>), results in a different version of Eq. (<ref>):

    ∂/∂ tM=R  M+ D∇·∇M+TE_D+TE_R .

The reaction truncation term, TE_R, can be evaluated by starting from the relevant truncation error expression (A23) in the Appendix of <cit.> (after correcting an error; the reaction term is only O(δ t) and not O(δ t^2)).

    TE_R=τ δ t∂/∂ t[R  M]= τ δ t [∂R/∂ tM+R∂M/∂ t]

and using Eq. (<ref>) 

    TE_R=τ δ t [∂R/∂ tM+R(R M + D∇·∇M )] .

Referring to the order of magnitude analysis performed for Eq. (<ref>), for the classical LBM scheme, the length scale of the domain (L) should be used in the reaction term instead of δ x because the strength of the reaction term R at each node in Eq. (<ref>) is determined from the nodes location (Eq. (<ref>)), rather than using a periodic function as in the hybrid splitting scheme. In this case, the leading order term in TE_R is the second term in Eq. (<ref>), 

    τ δ t R R ∼ O{τ δ t (γ L|G_0|)^2}

Again, for consistency, this truncation error term has to be much smaller than the smallest term in Eq. (<ref>), which now becomes the diffusion term. Since τ∼ O(1), this requirement implies that  

    δ t (γ L|G_0|)^2≪ D/δ x^2 →δ t≪D/(γ L|G_0| δ x)^2  ∼ 10^-1 ms

for a typical domain length L=100 μ m. This reaction error term scales with domain size as L^2, explaining the domain dependent results for the classical LBM scheme observed in Figure <ref>. 






§ APPENDIX B

In this appendix, the derivation of the membrane boundary condition is presented. For a general curved interface, Li et al. <cit.> proposed Dirichlet and Neumann boundary conditions describing the treatment of the nodes closest to the interface. The Dirchlet boundary conditions can be written as

    g^'_α̅(x_i,t) = c_d1ĝ_α(x_i,t) + c_d2ĝ_α(x_ii,t) + c_d3ĝ_α̅(x_i,t) + c_d4ϵ_D Φ_d,in


    g^'_α(x_e,t) = c_d1^∗ĝ_α(x_e,t) + c_d2^∗ĝ_α(x_ee,t) + c_d3^∗ĝ_α̅(x_e,t) + c_d4^∗ϵ_D Φ_d,ex

while the Neumann boundary conditions are

    g^'_α̅(x_i,t) = c_n1ĝ_α(x_i,t) + c_n2ĝ_α(x_ii,t) + c_n3ĝ_α̅(x_i,t) + c_n4δ t/δ xΦ_nα̅


    g^'_α(x_e,t) = c_n1^∗ĝ_α(x_e,t) + c_n2^∗ĝ_α(x_ee,t) + c_n3^∗ĝ_α̅(x_e,t) + c_n4^∗δ t/δ xΦ_nα.

Here the coefficients c_d1-c_d4 and c_n1-c_n4 are coefficients related to the membrane lattice link distance Δ_m while c_d1^∗-c_d4^∗ and c_n1^∗-c_n4^∗ relate to Δ_m^∗ = 1-Δ_m. Φ_d,in and Φ_d,ex are the concentrations at the membrane on the intracellular and extracellular sides respectively while Φ_nα and Φ_nα are the fluxes in the lattice direction. To maintain second-order accuracy, the coefficients for the Neumann boundary condition must be 

    c_n1 = 1,   
    c_n2 = -2Δ_m - 1/2Δ_m+1,   
    c_n3 = 2Δ_m - 1/2Δ_m+1,   and  
    c_n4 = 2/2Δ_m+1.

For the Dirichlet case, the second-order accuracy is preserved under certain relationships between the coefficients. For definiteness, here we use 

    c_d1 = 2(Δ_m-1),   
    c_d2 = -(2Δ_m - 1)^2/2Δ_m+1,   
    c_d3 = 2(2Δ_m - 1)/2Δ_m+1,   and  
    c_d4 = 3-2Δ_m/2Δ_m+1.

In both cases, c_d1^∗-c_d4^∗ and c_n1^∗-c_n4^∗ are the same coefficients as those in Eq. (<ref>) and Eq. (<ref>) but with Δ_m^∗ substituted for Δ_m. 

Eq. (<ref>) and Eq. (<ref>) describe the flux along the lattice direction, however, to implement the boundary condition, this flux must be related to the flux normal to the membrane. In the two-dimensional case, the following relationships exists

    Φ_nα̅ = {   1/c^'_d4[
    	(c^'_n1-c^'_d1) ĝ_β(x_i^',t) + (c^'_n2-c^'_d2) ĝ_β(x_ii^',t) +   (c^'_n3-c^'_d3) ĝ_β̅(x_i^',t)  ]sinθ
     -    1/c_d4[
    	(c_n1-c_d1) ĝ_α(x_i,t) + 
    	(c_n2-c_d2) ĝ_α(x_ii,t) + 
    	(c_n3-c_d3) ĝ_α̅(x_i,t) ]sinθ
     +    c^'_n4/c^'_d4δ t/δ xΦ_n,in}/ 
    	[
    	c_n4/c_d4δ t/δ xsinθ + c^'_n4/c^'_d4δ t/δ xcosθ]

and

    Φ_nα = {   1/c^'_d4[
    (c^'_n1-c^'_d1) ĝ_β̅(x_e^',t) + (c^'_n2-c^'_d2) ĝ_β̅(x_ee^',t) +   (c^'_n3-c^'_d3) ĝ_β(x_e^',t)  ]sinθ
     -    1/c^∗_d4[
    (c^∗_n1-c^∗_d1) ĝ_α̅(x_e,t)  + 
    (c^∗_n2-c^∗_d2) ĝ_α̅(x_ee,t) + 
    (c^∗_n3-c^∗_d3) ĝ_α(x_e,t) ]sinθ
     +    c^'_n4/c^'_d4δ t/δ xΦ_n,ex}/ 
    [
    c^∗_n4/c^∗_d4δ t/δ xsinθ + c^'_n4/c^'_d4δ t/δ xcosθ] ,

where c_d1^'-c_d4^' and c_n1^'-c_n4^' are the coefficients from Eq. (<ref>) and Eq. (<ref>) for Δ_m=0, and Φ_n,in and Φ_n,ex are the fluxes normal to the membrane <cit.>. The subscripts `e' and `ee' in x denote lattice nodes immediately adjacent to the membrane in the extracellular domain, while the subscripts `i' and `ii' denote corresponding adjacent nodes in the intracellular domain. The superscripted x^' refers to the extrapolated values from nodes within a respective domain for 
x^'_e= x^'_i= x_m, 
x^'_ii= x^'_i+e_β̅ δ t, and 
x^'_ee= x^'_e+e_βδ t,
where e_β and e_β̅ are in directions orthogonal to the lattice direction (cf. Figure <ref>).

To define the membrane boundary condition, we begin by considering the interfacial conditions, Eq. (<ref>), which result in two relations that can be expressed in terms of Dirichlet (Φ_d =M) and Neumann (Φ_n =D n·∇M) boundary conditions at either side of the interface,

    Φ_n,in = κ(Φ_d,ex -Φ_d,in)

and

    Φ_n,ex = D_ex∂M_ex/∂ n = 
    D_in∂M_in/∂ n = -Φ_n,in.


We consider the distribution functions representing particles towards the membrane in extra-and intra-cellular domains, denoted by g_α̅^'(x_i, t) and g_α^'(x_e, t), respectively (the reaction initialization step of g_i^'(x_n,t)=g̅_i (x_n,t+δ t) is implied). Substituting Eq. (<ref>) and Eq. (<ref>) into Eqs. (<ref> – <ref>), combining with Eq. (<ref>) and Eq. (<ref>), and rearranging yields

    g^'_α̅(x_i,t)=     
    A_1^iĝ_α(x_i,t)+
    A_2^iĝ_α(x_ii,t)+
    A_3^iĝ_α̅(x_i,t)+
       
    B_1^iĝ_α̅(x_e,t)+ 
    B_2^iĝ_α̅(x_ee,t)+ 
    B_3^iĝ_α(x_e,t)+ 
       
    C_1^iĝ_β(x^'_i,t)+ 
    C_2^iĝ_β(x^'_ii,t)+ 
    C_3^iĝ_β̅(x^'_i,t)+
       
    D_1^iĝ_β̅(x^'_e,t)+ 
    D_2^iĝ_β̅(x^'_ee,t)+ 
    D_3^iĝ_β(x^'_e,t)

and

    g^'_α(x_e,t)=    
    A_1^eĝ_α̅(x_e,t)+A_2^eĝ_α̅(x_ee,t)+A_3^eĝ_α(x_e,t)+
       
    B_1^eĝ_α(x_i,t)+ 
    B_2^eĝ_α(x_ii,t)+B_3^eĝ_α̅(x_i,t)+
       
    C_1^eĝ_β̅(x^'_e,t)+ 
    C_2^eĝ_β̅(x^'_ee,t)+ 
    C_3^eĝ_β(x^'_e,t)+
       
    D_1^eĝ_β(x^'_i,t)+ 
    D_2^eĝ_β(x^'_ii,t)+ 
    D_3^eĝ_β̅(x^'_i,t)  .

The coefficients in Eq. (<ref>) and Eq. (<ref>) (which are the same as Eq. (<ref>) and Eq. (<ref>) in the main text) are  

    A_j^i = (c_d4c_n4^'cosθ(c_n4^' c_d4^∗ c_njϵ_Dδ xcosθ+c_n4^' c_n4^∗ c_njκ δ t+c_d4^' c_n4^∗ c_nj ϵ_Dδ xsinθ)   
    
    + c_n4(cosθc_n4^' c_d4^∗ c_dj (c_n4^' κ δ t
    + c_d4^'ϵ_Dδ xsinθ)   
    
    + c_d4^' c_n4^∗ c_djsinθ(2c_n4^' κ δ t+ c_d4^'ϵ_D δ xsinθ)   ))/F 
       B_j^i =  (c_d4c_n4c^'_n4^2 κ δ tcosθ (c_nj^∗-c_dj^∗))/F 
       C_j^i =  (c_d4c_n4sinθ(c_nj^'-c_dj^')
    (c_n4^' c_d4^∗ ϵ_D δ xcosθ+ c_n4^' c_n4^∗ κ δ t+ c_d4^' c_n4^∗ ϵ_D δ xsinθ))/F 
       D_j^i =  (c_d4c_n4c_n4^∗ c_n4^' κ δ tsinθ (c_nj^'-c_dj^'))/F

and

    A_j^e= (c_d4c_n4^'cosθ(c_n4^' c_d4^∗ c_nj^∗ϵ_Dδ xcosθ+c_n4^' c_n4^∗ c_dj^∗κ δ t+c_d4^' c_n4^∗ c_dj^∗ ϵ_Dδ xsinθ)   
    
    + c_n4(cosθc_n4^' c_d4^∗ c_nj^∗ (c_n4^' κ δ t+ c_d4^'ϵ_Dδ xsinθ)   
    
    +c_d4^' c_n4^∗ c_nj^∗sinθ(2c_n4^' κ δ t+ c_d4^'ϵ_D δ xsinθ)   ))/F 
    
       B_j^e=(c_d4^∗ c_n4^∗c^'_n4^2 κ δ tcosθ (c_nj -c_dj))/F 
       C_j^e= (c_d4^∗ c_n4^∗sinθ(c_nj^'-c_dj^')(c_n4^' c_d4 ϵ_D δ xcosθ+c_n4^' c_n4 κ δ t+ c_d4^' c_n4 ϵ_D δ xsinθ))/F 
       D_j^e= (c_d4^∗ c_n4c_n4^∗ c_n4^' κ δ tsinθ (c_nj^'-c_dj^'))/F

with 

    F = (c_d4 c_n4^'cosθ(cosθc_n4^' c_d4^∗ ϵ_D dx+c_n4^' c_n4^∗ κ δ t+c_d4^' c_n4^∗ ϵ_Dδ xsinθ)    
      + c_n4(cosθc_n4^' c_d4^∗ (c_n4^' κ δ t+  c_d4^'ϵ_Dδ xsinθ)    
     + c_d4^' c_n4^∗sinθ(2c_n4^' κ δ t+ c_d4^'ϵ_D δ xsinθ)   ))  .

Here θ is the angle between the normal to the membrane and lattice direction e_α̅ while the coefficients  c_dj  and c_nj  are the coefficients of the Dirichlet and Neumann boundary conditions for a lattice membrane distance of Δ_m, c_dj^∗ and c_nj^∗ are the same coefficients for the  extracellular fraction Δ_m^∗=(1-Δ_m), and c_dj^' and c_nj^' are the coefficients for Δ_m=0. 



REFERENCES 




