Convergence of methods for coupling of microscopic and mesoscopic reaction-diffusion simulations
aa r X i v : . [ q - b i o . S C ] A p r Convergence of methods for coupling of microscopic and mesoscopicreaction-diffusion simulations
Mark B. Flegg a , Stefan Hellander b , Radek Erban a a Mathematical Institute, University of Oxford, 24-29 St Giles’ Oxford OX1 3LB, United Kingdom;e-mails: fl[email protected], [email protected] b Department of Information Technology, Uppsala Universitet, Box 480, 751 06 Uppsala, Sweden;e-mail: [email protected]
Abstract
In this paper, three multiscale methods for coupling of mesoscopic (compartment-based)and microscopic (molecular-based) stochastic reaction-diffusion simulations are investigated.Two of the three methods that will be discussed in detail have been previously reportedin the literature; the two-regime method (TRM) and the compartment-placement method(CPM). The third method that is introduced and analysed in this paper is the ghost cellmethod (GCM). Presented is a comparison of sources of error. The convergent propertiesof this error are studied as the time step ∆ t (for updating the molecular-based part of themodel) approaches zero. It is found that the error behaviour depends on another fundamentalcomputational parameter h , the compartment size in the mesoscopic part of the model. Twoimportant limiting cases, which appear in applications, are considered:(i) ∆ t → h is fixed;(ii) ∆ t → h → √ ∆ t/h is fixed.The error for previously developed approaches (the TRM and CPM) converges to zero onlyin the limiting case (ii), but not in case (i). It is shown that the error of the GCM convergesin the limiting case (i). Thus the GCM is superior to previous coupling techniques if themesoscopic description is much coarser than the microscopic part of the model. Keywords:
Multiscale simulation, reaction-diffusion, particle-based model.
1. Introduction
Multiscale stochastic reaction-diffusion methods which use models with different levels ofdetail in different parts of the computational domain are applicable to a number of biologicalsystems, including modelling of intracellular calcium dynamics [12], MAPK pathway [20]and actin dynamics [9]. In these applications, a detailed modelling approach (which requiressimulation of trajectories and reactive collisions of individual biomolecules) is only neededin a small part of the computational domain. The main idea of multiscale methods is thensimple to formulate [10]: we use a detailed modelling approach in the small subdomainof interest and a coarser model in the rest of the computational domain. In this paper,
Preprint submitted to Journal of Computational Physics October 15, 2018 etailed molecular-based (microscopic) models will be given in terms of Brownian dynamics[3, 25]. The remainder of the computational domain will be divided into compartments and amesoscopic (compartment-based) model will be used, i.e. we will simulate the time evolutionof the numbers of molecules in the corresponding compartments [5, 19].There have been a number of approaches developed for coupling different reaction-diffusion models. They include coupling of mesoscopic (compartment-based) models withcoarser (mean-field) PDE-based descriptions [13, 2, 26, 23], coupling of microcopic (molecular-based) models with mean-field PDEs [15, 18, 14], and coupling of microscopic and mesoscopicmodels [10, 11, 20, 21] A successful multiscale algorithm requires an accurate implementationof inter-regime transfer of molecules. In this paper, we will study convergence properties oftwo algorithms for coupling microscopic and mesoscopic descriptions which were previouslypublished in the literature: the two-regime method (TRM) [10, 11] and the compartment-placement method (CPM) [20]. One of the conclusions of our analysis is that these algorithmsdo not converge in the limit of small time steps and a fixed compartment size. Thus, we alsopropose another approach, the ghost cell method (GCM) which is suitable for this parameterregime.We will consider a reaction-diffusion model in the computational domain Ω ⊂ R N forboth N = 1 and N = 3. We will divide Ω into two parts, open sets Ω M and Ω C , whichsatisfy Ω M ∪ Ω C = Ω and Ω M ∩ Ω C = ∅ , (1)where an overbar denotes the closure of the corresponding set. The microscopic simulationtechnique is used in Ω M . Each molecule, j , in Ω M is considered to be a point particle atsome location in space, X j ( t ), at time t , which is updated according to discretized Brownianmotion, i.e. X j ( t + ∆ t ) = X j ( t ) + p D j ∆ t ζ , (2)where D j is the diffusion constant of the j -th molecule, ∆ t is a small prescribed time stepand ζ is a vector containing zero mean, unit variance normally distributed random numbers.In this paper, we will study the convergence of multiscale methods in the limit ∆ t → M , we have to specify what willbe done in the remainder of the domain, Ω C , where the mesoscopic model is used. In thispaper, we distinguish the following two cases:(i) the mesoscopic model is kept fixed in the limit ∆ t → t approaches zero.The resolution of the mesoscopic model (compartment size) will be denoted by h . Of par-ticular interest is the error that is caused as a direct result of the coupling and thus we willuse the parameter h as a measure of the compartment size at/on the interface between thetwo modelling subdomains. In the case of regular cubic compartments of volume h , theparameter h is simply the length of an edge of each cube. We will also consider unstructuredmeshes where the compartment size h will be suitably generalized. Using h , the cases (i)–(ii)can be formulated as follows:(i) ∆ t → h is fixed; 2ii) ∆ t → h → √ ∆ t/h is fixed.Both limits (i) and (ii) are important in applications. We will see that the error at theinterface ∂ Ω M ∩ ∂ Ω C of previously developed methods [10, 11, 20] converges to zero in thelimit (ii). This limit requires the refinement of the mesoscopic model. However, the standardmesoscopic model converges in the limit h →
2. The two-regime method (TRM)
The two-regime method (TRM) [10, 11] couples microscopic and mesoscopic subdomains bycareful selection of the jump rate over the interface from the mesoscopic compartments andcareful placement of these molecules into the microscopic domain. To date, the TRM hasbeen used with mesoscopic subdomains with regular meshes [12, 9]. The advantage of usingthis technique is that accuracy can be gained in ‘regions of interest’ Ω M without the needto run computationally expensive microscopic simulations over the entire domain Ω. In thissection we will briefly cover the two different simulation paradigms and then discuss howthese paradigms are combined using the TRM. The defining characteristic of ‘microscopic’ simulation techniques for diffusion is that eachmolecule in the system is simulated individually on a continuous domain. In particular, thesetechniques follow the trajectory of each Brownian molecule to a resolution dependent on thetime steps that are used. For illustrative purposes we consider here a time-driven microscopicalgorithm. That is, an algorithm with a defined constant time step. Furthermore, we willnot be considering volume exclusion effects in this manuscript. Each molecule, j , is thereforeconsidered to be a point particle at some location in space, X j ( t ), at time t . The Browniandiffusion of these molecules is modelled by (2). Reactions may take place between thesediffusing molecules at a particular time step if the reactants are within a given reactionradius of each other [24, 22].Molecule interactions with boundaries depend on the type of boundary: boundaries canbe reflective, adsorbing or reactive (partially adsorbing) [6]. Considering that √ D ∆ t is muchsmaller than the local radius of curvature of the boundary, then the boundary is locally flaton the scale of relative motion of the molecules in one time step. In the case of absorbing3oundaries, molecules are removed from the system when they are updated to a positionoutside of the boundary. Since we simulate Brownian motion using a finite time step, wehave to take into account that a molecule can interact with the boundary during the timestep [ t, t + ∆ t ] even if its computed position at time t + ∆ t is inside the simulation domain.The probability, P m that this molecule-boundary interaction occured within the time interval( t, t + ∆ t ] is dependent on the diffusion constant and the initial and final normal distancesfrom the boundary the molecule is found (∆ x i and ∆ x f respectively) P m (∆ x i , ∆ x f , D, ∆ t ) = exp (cid:18) − ∆ x i ∆ x f D ∆ t (cid:19) . (3)This probability will also be important when it comes to coupling of microscopic simulationswith mesoscopic simulations via an interface in the two-regime method [10, 11]. Mesoscopic approaches to reaction-diffusion processes are simulated on a lattice. For thepurposes of the TRM we will describe how this is done for a regular cubic lattice. Thedistance between each node is h . In a mesoscopic model, molecules can be thought to existonly at lattice nodes rather than existing in continuous space. The state of the simulationat any moment of time is defined by a set of numbers describing the copy numbers N i,j ofmolecules of the i -th type at the j -th lattice point. Considering the diffusion of (non-reacting)molecules, the expected state of the system E( N i,j ) is described by the equation: d E( N i,j ) dt = D i X k ( q k,j E( N i,k ) − q j,k E( N i,j )) , (4)where q k,j is the propensity per molecule to go from the k -th compartment to the j − thcompartment. It is possible to show that for a regular lattice with spacing h , q k,j = ( D i /h , if k and j are adjacent lattice points,0 , if k and j are not adjacent lattice points, (5)results in the recovery of the discretized form of the diffusion partial differential equation andcan therefore be used to approximate a diffusion process on the lattice correct to order h .The simulation of a mesoscopic reaction-diffusion process usually makes use of event-drivenalgorithms, such as the Gillespie algorithm [17] or the Gibson-Bruck algorithm [16]. Weshall conceptualize the mesoscopic simulation by considering that when a molecule is at aparticular lattice point, rather than existing at the node, it is somewhere at random insidethe compartment belonging to the node defined by the lattice dual mesh [5]. That is, for aregular cubic lattice with node spacing h , each molecule which is at a particular lattice pointis thought to exist inside the cubic compartment of side length h for which the lattice pointis at the center. It is important to note that the state of the molecule has no specific locationbut rather is thought to exist in a probabilistic sense uniformly over its compartment.4 .3. Interfacing microscopic and mesoscopic simulations Interfacing microscopic and mesoscopic simulations of reaction-diffusion processes using theTRM has previously been derived for mesoscopic regimes that use regular cubic lattices[10, 11]. The TRM is proposed by partitioning the domain Ω into subdomains (1) separted bythe interface I = ∂ Ω M ∩ ∂ Ω C . In both subdomains, molecules behave as they would normallyaccording to the rules of that particular regime. We describe the TRM with an event-drivenmesoscopic simulation in Ω C and a time-driven microscopic simulation with constant timestep ∆ t in Ω M . Reactions do not cause any issue within the domain because they occurlocally. We focus, therefore, on the correct manner in which molecules may migrate over theinterface I . It is assumed that the TRM is simulated such that √ D ∆ t ∼ h ≪
1. A diagramof the numerical TRM scheme using a regular cubic lattice can be seen in two dimensionsin Figure 1. A detailed TRM algorithm may be found in the reference [10]. In order that amolecular migration over the interface is smooth with optimally small error, the propensityΓ per molecule to cross the interface I from each adjacent compartment is dependent on theparameters h and ∆ t . For a regular cubic mesoscopic lattice,Γ( h, ∆ t ) = 2 r Dπ ∆ th , (6)where D is the diffusion constant of the migrating molecule. The TRM considers thatmicroscopic molecules in Ω M cease to be microscopic molecules, in principle, when theymigrate over the interface. Molecules are therefore absorbed by the interface I from Ω M and placed in the closest compartment in Ω C . Equation (3) is used to absorb all moleculeswhich interacted with the interface. If this is not used then molecules effectively migrate intoΩ C and back out again without changing from a microscopic molecule to a mesoscopic one.This is crucial for coupling of the two regimes as outlined in the derivation in the reference[10]. Furthermore, molecules must be precisely placed in Ω M when migrating from Ω C . Inparticular, the perpendicular distance x the molecule is placed from the interface into Ω M isfound by sampling from the distribution f ( x ) f ( x ) = r π D ∆ t erfc (cid:18) x √ D ∆ t (cid:19) , (7)where erfc ( x ) = 2 /π R ∞ x exp( − t ) dt is the complementary error function. In higher dimen-sions, the initial position of molecules migrating into Ω M can be chosen to be uniformlydistributed tangentially to the interface in the region of the originating compartment [11].Then the error associated with the TRM is O ( h ). We shall investigate the error associatedwith the TRM in 1D in a later section of this manuscript and compare it with the GCMmethod introduced in Section 4.
3. Compartment-placement method (CPM)
In this section, we will discuss how mesoscopic simulation is implemented on an irregularlattice [5]. We will then present a brief description of the CPM [20].5 C Ω M I h Mesoscopic Microscopic
D/h Jump propensity per particleCopy numbers
Γ(h,Δt)
Jump propensity per particleAll molecules absorbed by compartments x Figure 1:
Graphical representation of the TRM on a regular square lattice.
Mesoscopic simulations on Cartesian meshes are convenient in the sense that they are mem-ory lenient. However, complex geometries and surfaces with high curvature, are easier toresolve accurately with an unstructured mesh. Living cells can have different shapes andeukaryotes have a complicated internal structure with two-dimensional membranes and a one-dimensional cytoskeleton [1]. The geometrical flexibility of unstructured meshes is thereforean advantage when considering simulations of realistic biological problems.Consider a domain Ω. The domain is covered by a primal mesh, such that the bound-ary ∂ Ω is covered with non-overlapping triangles and the domain Ω is covered with non-overlapping tetrahedra (resp. triangles in 2D). A dual mesh is constructed from the primalmesh, see Figure 2, from the bisectors of the tetrahedra (resp. triangles) that use the nodesas vertices. The diffusion of molecules is now modelled as discrete jumps between the nodesof the dual mesh. The rate q i,j at which a molecules jump from voxel V i to V j is given bythe diffusion constant of the molecule and the finite element discretization of the Laplacianon the primal mesh. For details on how the unstructured meshes and the diffusion matrixare generated the reader is referred to [5]. The algorithm for the CPM is presented in a similar way to the TRM. The algorithm pro-gresses asynchronously by updates in the mesoscopic simulation and microscopic simulationseparately [20]. The jump rates from compartments that are on the interface I betweenregimes are calculated from the underlying mesh over the entire domain. That is, the jumprates are calculated by computing the mesoscopic jump rates between interfacial compart-6 ompartments from dual meshNodes (a) (b) Figure 2:
Schematic of CPM computational domain. (a)
The primal mesh is indicated with red dashed lines.The nodes are connected to form triangles. The bisectors are then drawn in to create the dual mesh (bluedotted lines). Compartments are drawn from the dual mesh with one node at the center of each compartment.One example compartment is shown in blue. (b)
The domain is split into mesoscopic Ω C and microscopic Ω M domains. Jumps between compartmentsand from the compartments into Ω M are calculated using a finite element discretization of the Laplacian.The copy numbers of molecules in each compartment in Ω C are stored whilst in Ω M each molecule has itsown position in continuous space. ments and “compartments” that are adjacent to the interface in the microscopic domain Ω M (see Figure 2).Molecules that start in a compartment in Ω C and, at the end of the time step, haveended up in Ω M are initialized uniformly inside the “compartment” which they jump into,and is the process from which the CPM has been named. Molecules in Ω M migrate backto Ω C via microscopic domain diffusion (2). When a molecule appears inside one of themesoscopic compartments from Brownian motion, it is encorporated into that compartmentby increasing the copy number inside this compartment.The CPM has been determined using heuristics. Molecules that are in compartmentsobey mesoscopic rules for diffusive migration. This includes molecules that are on interfacialcompartments. They jump to compartments in Ω M as though they were still in Ω C . Whenthis occurs, initialization of the molecules must take place. The molecules are initiateduniformly over the compartment in which they are placed. Molecules are not placed atthe node at the center of this compartment because this would unphysically concentratemolecules at this point and reactions would occur between possible reactants apon migrationover the interface. Conversely, molecules that are in Ω M obey microscopic rules for diffusivemigration (Brownian trajectory). When this Brownian trajectory leads to a compartment, itcan no longer be described using the microscopic description and is added to the compartment7G.1] Initialize lattice over whole domain Ω and construct dual mesh (compartments).Generate interface I on the edges of compartments to separate Ω M from Ω C . Choose∆ t and set time t = 0. Determine q k,j using finite element method between allcompartments [5].[G.2] Initialize the initial state of the system by placing molecules in compartments inΩ C and placing molecules in continuous space in Ω M . Count and store numbersof molecules in ghost cells, those compartments in Ω M which are adjacent to theinterface I .[G.3] Determine the time τ for the next event (reaction or diffusive) in Ω C or diffusivejumps to and from ghost cells and Ω C .[G.4] If t + τ < ∆ t + ∆ t ⌊ t/ ∆ t ⌋ then change the state of the system to reflect the nextevent corresponding to τ and update time t := t + τ . If this event is a diffusive jumpfrom ghost cell to Ω C choose a molecule at random within the relavant ghost cellto migrate. If this event is a diffusive jump from Ω C to a ghost cell then initializethis molecule with uniform probability over the ghost cell.[G.5] If t + τ ≥ ∆ t + ∆ t ⌊ t/ ∆ t ⌋ then update the positions of all molecules in Ω M using(2). Check for reactions in Ω M [11]. All molecules incident on the interface I arereflected. Update time t := ∆ t + ∆ t ⌊ t/ ∆ t ⌋ [G.6] Repeat steps [G.3]–[G.5] until the desired end of the simulation. Table 1: The ghost cell method algorithm. in which it lands. As we shall see, this heuristic approach can lead to inaccuracies. Theinaccuracies can be minimized if h ∼ D ∆ t (that is, if the size of the compartment isapproximately the size of a microscopic molecular jump).
4. The ghost cell method (GCM)
Here we will consider a new method for interfacing mesoscopic and microscopic simulations.This method uses different assumptions to the TRM and CPM and is therefore implementeddifferently. We call this method the ghost cell method (GCM) since microscopic moleculesin Ω M feel the presence of a pseudo-compartment allowing for instantaneous jumping fromΩ M to Ω C in the same way that molecules within compartments jump instantaneously. Thesteps of the GCM are given in Table 1.The key assumption that is used in the TRM and CPM is that molecules in Ω M mi-grate via diffusion (2) over the interface I whereby they become parts of the correspondingcompartment. In the GCM, this assumption is relaxed. Instead, molecules migrate over theinterface using the compartment-based approach in both directions. Microscopic moleculesin Ω M near the interface feel the presence of a layer of “ghost” cells (compartments). In thestep [G.2] in Table 1, we calculate the numbers of molecules in these “ghost” cells. They are8 RM/CPM paradigm GCM paradigm
Absorbedat interfaceTransfered at set propensity Reflectedat interfaceAbsorbed by ghost cell
Transfered at set propensity (a) (b)
Figure 3:
A diagram illustrating the fundamental differences between (a) the TRM/CPM paradigm and (b) the GCM paradigm. used in the step [G.4] to create a fully compartment-based simulation of transition acrossthe interface I .To justify the GCM, let us consider a simulation of diffusion in a domain Ω for which amesoscopic method was implemented. Then consider the same domain where a microscopicsimulation is implemented. Let the molecules of the microscopic simulation be “binned”according to compartments of the mesoscopic simulation. The expected number of moleculesbinned into each compartment should match that of the mesoscopic simulation to withinthe precision of the mesoscopic method. This is because both simulations are accuraterepresentations of the same phenomena, diffusion. This is the philosophy behind the GCM.Molecules which are binned into ghost compartments near the interface may jump intocompartments in Ω C via the rates prescribed by the mesoscopic algorithm. If both regimesare correct individually then the flux over the interface I is the same as though a mesoscopicalgorithm was used over the whole domain. To ensure that microscopic molecules do notmigrate to Ω C via diffusion (2), they are reflected at the interface I in the step [G.5]. Figure 3demonstrates the principle differences between a TRM/CPM and a GCM description of theinterface. In Appendix A we provide a mathematical analysis of the GCM in one dimensionto demonstrate that the expected concentration and flux of molecules over the interface arematched. The theoretical error associated with the GCM scales as √ ∆ t which is on thesame order as that of the TRM. Unlike the TRM, this error, as we will see in the later partof this manuscript, is reduced to zero by reducing √ ∆ t/h .The ghost cell method is implemented using the algorithm in Table 1. This algorithm isgiven for an event-driven mesoscopic simulation and a time-driven microscopic simulation,however it can also be extended to event-driven microscopic simulations.9 . Numerical results and discussion In this section, we present numerical examples comparing the TRM, CPM and GCM. First,we demonstrate how the error associated with the interface I is dependent on choices ofmesh spacing h in the mesoscopic subdomain at the interface and the time step chosenfor the microscopic subdomain ∆ t for both the TRM and the GCM using one dimensionalsimulations. We use a simple diffusion test problem to compare the diffusive flow over the interface withan exact solution which can be analytically obtained. We use the domain Ω = (0 ,
1) andsubdomains Ω C = (0 , . , Ω M = (0 . , I = { . } . Weinitially position N = 5 × molecules according to the distribution g ( x ) = 2 x , x ∈ Ω. Weconstruct regular spaced compartments of width h = 0 . C and “bin” the moleculesgenerated in Ω C into these compartments. We allow these molecules to diffuse throughoutthe domain Ω with a diffusion constant D = 1 using the TRM or GCM until time t = 1. Atthe boundary x = 0 molecules are absorbed and placed at x = 1. At the boundary x = 1molecules are reflected. In this way, N g ( x ) is the steady state distribution of this systemand 0 . N is the steady state number of molecules in Ω C . We define a measure of the error E to this test problem for each simulation scheme E = P j N j (1) − . N N , (8)where N j is the copy number of molecules in the j -th compartment evaluated at t = 1 andthe sum is taken over all compartments in Ω C .In order to see the effect of the compartment spacing near the interface h on the error E for both the TRM and GCM we start with the set of regular compartments(0 , h ) , ( h , h ) , . . . , (0 . − h , . , which have nodes (centres of compartments) at h /
2, 3 h /
2, . . . , 0 . − h /
2. Then we usethe following lattice refinement technique designed specifically so that the position of theinterface does not change (see Figure 4):[R.1] Delete the two nodes closest to the interface.[R.2] Introduce into the space between the new node closest to the interface and the interface(a distance of ∆ x ) three nodes placed consecutively a distance of 2∆ x/ C Ω M I Before refinementAfter refinement Δx ghost cell size matches interface compartment Figure 4:
Diagram of one iteration of the lattice refinement [R.1]–[R.3]. to the ghost cell in the GCM. The refinement algorithm [R.1]–[R.3] is repeated m times suchthat the size of the final compartment in Ω C (and ghost cell), h m , is h m = h (cid:18) (cid:19) m . (9)A diagram representing one iteration of the refinement technique [R.1]–[R.3] is shown inFigure 4. The error is computed for various final compartment sizes h m ( m = 0 , , . . . , t k ( k = 0 , , . . . ,
10) where∆ t k = 2 k ∆ t , (10)and ∆ t = 5 × − .Figures 5 and 6 show how the absolute error k E k given by (8) depends on both parameters h m (compartment size on the interface) and ∆ t for the TRM and GCM algorithms respec-tively. The error due to the interface in the TRM includes a shift of h m / C [10]. This is because of the “initialization”of molecules from Ω M into Ω C . Unlike the initialization of molecules from Ω C into Ω M ,molecules that are transported in the reverse direction cannot be placed carefully accordingto a continuous distribution but must necessarily be placed in the nearest compartment.This initialization has an expected position of h m / h m / C with a continuous distribution, for symmetry reasons one would expect this to be done witha distribution of f ( x ) given by (7). The average distance, therefore, that a molecule wouldideally be placed into Ω C is R ∞ xf ( x ) dx = √ πD ∆ t/
2. Therefore, the error that is due to11 igure 5:
Surface plot of the absolute error k E k defined by (8) as a function of compartment size at theinterface h m and time step ∆ t for the one dimensional test problem using the TRM. unphysical shifting of molecules is proportional to the expected shift of molecules as theyare transferred from Ω M to Ω C . That is E ∝ h m − √ πD ∆ t .In Figure 5 a dotted red line showing h m = √ πD ∆ t approximately follows the path ofthe minimum absolute error. The discrepancy between the actual minimum absolute errorand the dotted red line in Figure 5 can be attributed to higher order error that is inherent inthe mesoscopic approximation to the diffusion equation. To show that E ∝ h m − √ πD ∆ t ,Figure 7 is a plot of error E versus h m − √ πD ∆ t . The plot is generated by using variousvalues of h m (see legend) and then plotting a number of points while changing ∆ t . Whilst itis clear that the graph is approximately linear, the higher order mesoscopic error is clearlyseen in the form of a vertical displacement of this curve about the origin. The effect thatthe higher order mesoscopic error has on the interface is difficult to quantify because it willdepend on the particular molecular system. Therefore, the best choice of parameters thatcan be chosen for the TRM is h m ∼ √ πD ∆ t .In Figure 6, we see that the error of the GCM depends on ∆ t and specifically on itsrelative size compared to h (the analysis of the GCM is provided in Appendix A). Rapidlyincreasing error (quickly saturating the color bar in Figure 6) is observed when h m ∼ √ πD ∆ t .The higher order mesoscopic error artefact can also be seen in Figure 6 since this artefactis independent of the coupling mechanism (see the larger absolute error for large values of12 igure 6: Surface plot of the absolute error k E k defined by (8) as a function of compartment size at theinterface h m and time step ∆ t for the one dimensional test problem using the GCM. h ). The GCM is therefore most accurate for very small values of ∆ t . Whilst in practicemaking ∆ t small may significantly increase the computing time, small ∆ t is often requiredfor accurate microscopic simulation (for example, capturing reactions with high resolution)and in such cases the GCM is more appropriate than the TRM. In this section we will demonstrate how, when using an unstructured mesh, the error asso-ciated with the GCM coupling converges as ∆ t → h ∼ √ D ∆ t where h is the average size of boundary compartments. Boththe error associated with the CPM and GCM are due to imbalances in the flux of moleculesover the interface. We implement the CPM and GCM in three spatial dimensions using atetrahedral primal mesh as described in Section 3. The implementation builds on the freelyavailable software URDME [4].We consider a cube with side length L = 1. The cube is first discretized with an unstruc-tured mesh and then divided into a mesoscopic region Ω C , and a microscopic region Ω M ,where Ω M is the set of all voxels with a vertex ( x, y, z ) such that x < . C = Ω \ Ω M .Here Ω is the set of all voxels. The partitioning is illustrated in Figure 8 for two differentmesh sizes.We start each simulation with N = 2 · molecules whose initial positions aresampled from a uniform distribution. The diffusion constant of the molecules is D = 1, and13 h − √ π D ∆ t E rr o r , E h = 0.1h = 0.051h = 0.026h = 0.013h = 0.0068h = 0.0035 Figure 7:
Scatter plot of the error E versus h m − √ πD ∆ t for the TRM. The different color points representdifferent values of the compartment size at the interface h (see legend) and in each instance ∆ t is variedfrom × − to × − . we simulate the system for time t=0 .
1. Since we start with a uniform distribution and themolecules only diffuse and do not react, we expect the distribution to be uniform at thefinal time. As the interface is parallel with the y − z -plane, we expect that the distributionsof molecules in the y - and z -directions are uniform, but that we get a small error in thedistribution of molecules in the x -direction. We now divide the x -axis into 10 bins of equallength, and then count the number of molecules in each bin at the final time. Mesoscopicmolecules are binned by first sampling a continuous position from a uniform distribution onthe voxel. We expect N /
10 molecules in each bin, and can therefore estimate the error E by E = P i =1 | N i − N / | N . (11)In Figure 9 we have computed k E k for different mesh sizes and time steps. As expected, theerror decreases as we refine the mesh and decrease the time step.In the CPM method, mesoscopic (resp. microscopic) molecules stay mesoscopic (resp.microscopic) during a time step. This implies that the time step should be chosen sufficientlysmall such that a molecule does not diffuse across several voxels. On the other hand, if the14a) (b) Figure 8:
The cube [0 , is partitioned into a mesoscopic region (grey) and a microscopic region (white). (a) a coarse mesh ; (b) a fine mesh . −5 −4 −1 ∆ t E Figure 9:
The error k E k of the GCM method for different mesh sizes and time steps. The error decreaseswith decreasing time step, as expected. −5 −4 −3 −2 −2 −1 ∆ t E CPM methodGCM methodh /(2D) Figure 10:
Comparison of the error k E k produced by the GCM and CPM for interfacing microscopic andmesoscopic simulations as a function of the time step in the microscopic simulation domain ∆ t . The lengthscale h is defined to be the cubic root of the average volume of a voxel. time step is too small the distribution of molecules in space will be biased towards themicroscopic region. This can be seen by considering a microscopic molecule diffusing intothe mesoscopic regime. If the time step is small, it is likely that it will be close to themicroscopic regime at the end of the time step, but if it ends up on the mesoscopic side itwill nevertheless be considered uniformly distributed in the voxel at the end of the time step.Thus, the time step should not be chosen too small relative to the size of the voxels, or theerror due to the spatial splitting will become large.Since the GCM converges with decreasing time step, but performs worse for larger timesteps, one could suspect that there is a regime where the CPM in [20] performs better thanthe GCM. At some point, however, the error of the GCM will become small and outperformthe CPM in [20]. The errors of the different methods are compared in Figure 10 for a meshwith 49101 voxels. Indeed, we see that the CPM method performs better for time stepsdown to almost ∆ t = 10 − , at which point the error of the GCM method becomes smaller.
6. Summary
In this paper we have compared two existing mesoscopic-microscopic coupling techniques forstochastic simulations of reaction-diffusion processes with a new convergent method calledthe ghost cell method (GCM). Here we will summarize the specific sources of error of theTRM, CPM and GCM, when they converge, how they may be optimized for accuracy andnotes on their computational efficiency. 16 .1. Summary of the two-regime method
The TRM couples molecules by considering that mesoscopic compartments contain moleculesthat are evenly distributed in a probabilistic sense. As these molecules diffuse over the in-terface they are placed according to the distribution f ( x ) given by (7). Molecules migratingin reverse from the microscopic regime to the mesoscopic regime must be absorbed by theinterfacial compartments and be indistinguishable from other molecules in these compart-ments by virtue of this paradigm. As such, instead of migrating an average distance over theinterface proportional to √ ∆ t it becomes evenly distributed over the compartment with anexpected location of h/ h − √ πD ∆ t ) / h − √ πD ∆ t . The nature of this error isthat, if the expected net flux of molecules over the interface is 0, then no error due to thepresence of the interface will be experienced. This is important to note, since, this is not thecase for both the CPM and GCM methods. Furthermore, since the error is proportional to h − √ πD ∆ t it clearly converges in the limiting case (ii) described in the introduction butnot in the limiting case (i).Whilst the TRM can give controllably accurate results, it can be computationally morecostly to implement. This is because perfect absorption of molecules is required on theboundary and this means that each molecule in the molecular-based domain needs to bechecked for interaction with the boundary in a given time step using (3). The CPM is a coupling mechanism that, whilst heuristically derived, can produce accurateresults under some circumstances and do so with minimum computational cost. Moleculesare placed within a pseudo-compartment in the molecular-based domain via diffusion fromthe compartment-based domain. In reverse molecules are placed in compartments fromthe molecular-based domain via diffusion of these molecules in the continuous domain overthe interface. The antisymmetry that is seen in the methods of migration, mesoscopicto microscopic and microscopic to mesoscopic, result in a boundary layer in the expecteddistribution of molecules at the interface. This boundary layer is caused by the fact thatmolecules diffusing from the molecular region to the compartmental region are considereduniformly distributed on the compartment at the end of the time step. If ∆ t is smallcompared to h , this will be a poor approximation. It should thus be noted that the error ofthe CPM does not converge in the limiting case (i) described in the introduction, howeverthe error appears to converge according to limiting case (ii).The CPM is computationally efficient. Its only inefficiency is that, unlike the TRM, itrequires the knowledge of a pseudo-compartment in the microscopic domain. The TRM istherefore more appropriate than the CPM for coupling completely independent simulationalgorithms, since the CPM requires its own custom algorithm to be implemented fully.17 .3. Summary of the ghost cell method The GCM couples molecules according to a discrete domain on each side of the interface.Molecules that are in the microscopic domain are binned according to a ghost compart-ment/cell and jump into the mesoscopic domain using jump rates derived using the meso-scopic approach. In such a way, symmetry is conserved in the method of migration frommesoscopic to microscopic and from microscopic to mesoscopic, unlike the CPM. It is im-portant that the molecules are binned correctly for this coupling to work accurately. Theerror, therefore, can be attributed to molecules that are in the ghost cell when they shouldnot be, or not in the ghost cell when they should be. Therefore, if the compartment size h at the interface (and of the ghost cell) is much larger than the resolution of the particletracking in the microscopic domain, the correct number of molecules will be in the ghostcell. The error therefore converges in the limit of small ∆ t so long as h is sufficiently coarse.Furthermore, it is possible to show that, unlike the TRM, this source of error is not due to adisplacement of molecules but an unballanced flux of molecules (see Appendix A) and willtherefore appear even if the expected net flux over the interface is 0. The GCM, however,is convergent in the limiting case (i) but not (ii) from the introduction giving the GCM aunique advantage over both the CPM and TRM.The GCM is computationally efficient for small ∆ t since the jump rates from the micro-scopic domain to the mesoscopic domain are determined by the ghost cell size and not thetime step (like the TRM for example). Acknowledgements:
The research leading to these results has received funding from theEuropean Research Council under the
European Community ’s Seventh Framework Pro-gramme (
FP7/2007-2013 )/ ERC grant agreement
No. 239870. This publication was basedon work supported in part by Award No KUK-C1-013-04, made by King Abdullah Uni-versity of Science and Technology (KAUST). Stefan Hellander has been supported by theNational Institute of Health under Award no. 1R01EB014877-01 and the Swedish ResearchCouncil. Radek Erban would also like to thank Brasenose College, University of Oxford, fora Nicholas Kurti Junior Fellowship; the Royal Society for a University Research Fellowship;and the Leverhulme Trust for a Philip Leverhulme Prize.
Appendix A. Mathematical justification for the ghost cell method
Here we present an analysis of the GCM in one-dimension. We show that error of theGCM that is produced on the interface between mesoscopic and microscopic subdomainsconverges in the case (i). Specifically, we see convergence of the interface-derived error asΛ = √ D ∆ t/h →
0. This property of convergence is unique to the GCM when compared withother reported coupling mechanisms. In showing that the interface-derived error vanishes inthe small time step limit, we will show that rapid variation within the boundary layer of theinterface vanishes as Λ →
0, leaving behind a linear approximation of the true distributionof molecules. Since the error at the interface will be of order h it is accurate to the sameorder as the mesoscopic algorithm itself. 18ithout loss of generality, consider an interface at x = 0 on an infinite one-dimensionaldomain. To the left of this interface ( x <
0) a compartment-based model is used with fixedcompartment size h . To the right of the interface ( x >
0) a molecular-based algorithm isused and is updated at fixed time increments of ∆ t , i.e. Ω C = ( −∞ ,
0) and Ω M = (0 , ∞ ).We denote the compartments in Ω C by C k = ( − kh, − kh + h ), where k = 1 , , . . . . Then theinterface compartment is C = ( − h, C M = (0 , h ).Molecules in Ω C are described only by their compartment. Their compartment changeswith an exponentially distributed random time with a rate that is given by D/h . This rateis conditional on initial and final states being compartments. The rate given by D/h ischosen in such a way that the expected number of molecules in each compartment matchesthat of a discretized diffusion equation (see (4) and (5)). These rates, however, breakdownin the case of the TRM because the initial and final states of jump across the interfaceare not compartments but rather the final state is a molecule in Ω M . The jump across theboundary for the TRM is given by (6). Molecules in Ω M have one thing in common withthose in compartments. A domain that is modeled microscopically and then binned intocompartments shows the same expected behaviour as a domain modeled with compartmentsto leading and first order accuracy in h in the limit as ∆ t →
0. Therefore, in an attempt tointerface the two regimes together it may be appropriate to bin molecules in Ω M into a ghostcell/compartment C M near the interface. The molecules that are in C M will then have thesame properties as a compartment from the perspective of the interface compartment C . Tothis end, any molecule in C M may spontaneously change state from the molecular domainto C . We expect that since the interface compartment-bound molecules see a compartmentstate for molecules in C M , the change of state from C to a random position within C M willoccur with a normal inter-compartmental rate. In the following analysis we show that thisis the case.We shall test the hypothesis by matching the master equations for C and the probabilitydistribution in Ω M in such a way that no rapid variation in probability to find molecules,¯ p ( x, t ), is apparent at the interface. We shall assume that the rate for molecules to jumpinto Ω M from C is Γ + and are placed in an initial position from the interface given bythe probability distribution f ( x ). Molecules in Ω M spontaneously jump into C with a rateΓ − g ( x ). Functions f ( x ) and g ( x ) are normalized such that they have a unit integral overΩ M . We shall show thatΓ + = Γ − = Dh and g ( x ) = f ( x ) = h , for x ∈ C M ;0 , otherwise. (A.1)We will find it convenient for the sake of notation to introduce the parameters α ± = h Γ ± D , and Λ = √ D ∆ th . To show (A.1) we focus on the purely diffusive problem, since bulk reactions have no effect onboundary conditions. In order to limit the flux of molecules jumping into C , all moleculesin Ω M that hit the interface by Brownian motion are reflected back to Ω M .19e denote the probability of finding a molecule in compartment C k , k = 1 , , . . . , by p k ( t ) h (so that p k ( t ) approximates the probability density function at the node within thiscompartment). If we denote by p ( x, t ) the probability density function of the discrete-time molecular-based algorithm, then the transmission/reflection rules give us the followingmaster equation p ( t + ∆ t ) = (cid:0) − (1 + α + )Λ (cid:1) p ( t ) + Λ p ( t ) + α − Λ Z ∞ g ( x ) p ( x, t ) d x, (A.2) p ( x, t + ∆ t ) = Z ∞ p ( y, t ) √ πD ∆ t (cid:20) exp (cid:18) − ( x − y ) D ∆ t (cid:19) + exp (cid:18) − ( x + y ) D ∆ t (cid:19)(cid:21) d y + α + D ∆ tp ( t ) f ( x ) h − α − D ∆ tp ( x, t ) g ( x ) h . (A.3)In the vicinity of x = 0 there is a boundary layer of width O ( h ) so long as f ( x ) and g ( x ) vanish for x ≫ h [6]. We rescale (A.2) and (A.3) using the (dimensionless) boundarylayer coordinate ξ = x/h . We also denote p inner ( ξ, t ) = p ( hξ, t ), f inner ( ξ ) = h f ( hξ ) and g inner ( ξ ) = h g ( hξ ). The rescalings of f and g by h are done to keep the integrals of thesefunctions equal to 1. Thus, in the boundary layer coordinates, (A.2) and (A.3) become p ( t + ∆ t ) = (cid:0) − (1 + α + )Λ (cid:1) p ( t ) + Λ p ( t )+ α − Λ Z ∞ g inner ( ξ ) p inner ( ξ, t ) d ξ, (A.4) p inner ( ξ, t + ∆ t ) = Λ − Z ∞ p inner ( η, t ) (cid:2) K (cid:0) Λ − ( η − ξ ) (cid:1) + K (cid:0) Λ − ( η + ξ ) (cid:1)(cid:3) d η + α + Λ p ( t ) f inner ( ξ ) − α − Λ p inner ( ξ, t ) g inner ( ξ ) , (A.5)where K ( x ) = (4 π ) − / exp( − x / D ∆ t to h and we wish to show that as Λ → h , which is the same size of the errorassociated with the mesoscopic discretization in Ω C . In order to join these models smoothlywe require in Ω C that p ( t ) = p ( − h/ , t ) = p (0 , t ) − h p x (0 , t ) + O ( h ) + O (Λ) , (A.6) p ( t ) = p ( − h/ , t ) = p (0 , t ) − h p x (0 , t ) + O ( h ) + O (Λ) , (A.7) p ( t + ∆ t ) = p ( t ) + O (∆ t ) , (A.8)while, for the molecular-based side, we want variation from the linear approximation in theboundary layer to be limited to O (Λ) up to order h accuracy, so that p inner ( ξ, t ) = p (0 , t ) + h ξ p x (0 , t ) + O ( h ) + O (Λ) , (A.9) p inner ( ξ, t + ∆ t ) = p inner ( ξ, t ) + O (∆ t ) . (A.10)20he prescription of a consistent probability density p (0 , t ) and derivative p x (0 , t ) for bothsides of the interface, along with linear approximations sufficiently close to the interfaceequates to continuity and differentiability over the interface which is the matching conditionthat we are attempting to achieve.Substituting (A.6)–(A.10) into (A.4) and (A.5) and equating terms of the same order in h and leading order in Λ gives the following conditions that must be placed on g ( x ), f ( x ), α + and α − :(i) O ( h Λ ) terms from equation (A.4) give condition: α + = α − Z ∞ g inner ( ξ ) , d ξ. (A.11)Condition (A.11) states how the relative rates for molecules to transition to and fromΩ M and Ω C must be dependent on the relative sizes of C and C M .(ii) O ( h Λ ) terms from equation (A.4) give condition:2 = α + + α − Z ∞ ξg inner ( ξ ) , d ξ. (A.12)Condition (A.12) states how the rates for molecules to transition to and from Ω M andΩ C depend on the average distance molecules are placed from the interface when placedwithin C M . This is the same condition given in Ω C for jumps between the compartments.(iii) O ( h Λ ) terms from equation (A.5) give condition: α + f inner ( ξ ) = α − g inner ( ξ ) . (A.13)Condition (A.13) states that molecules must be placed into Ω M with the same probabilityweighting that they are taken out and placed back into Ω C .(iv) O ( h Λ ) terms from equation (A.5) are automatically satisfied.The GCM that is presented in this manuscript uses parameters (A.1) which satisfy the threeconditions (A.11)–(A.13) listed above. Such a scheme, therefore, has an error that is nogreater than the error of the mesoscopic scheme in the limit Λ →
0, in other words, in thelimit ∆ t → h remains constant. References [1] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter.
Molecular Biologyof the Cell . Garland Science, New York, 2007.[2] F. Alexander, A. Garcia, and D. Tartakovsky. Algorithm refinement for stochasticpartial differential equations: I. linear diffusion.
Journal of Computational Physics ,182(1):47–66, 2002. 213] S. Andrews and D. Bray. Stochastic simulation of chemical reactions with spatial reso-lution and single molecule detail.
Physical Biology , 1:137–151, 2004.[4] B. Drawert, S. Engblom, and A. Hellander. URDME: A modular framework for stochas-tic simulation of reaction-transport processes in complex geometries.
BMC Syst. Biol. ,6:76, 2012.[5] S. Engblom, L. Ferm, A. Hellander, and P. L¨otstedt. Simulation of stochastic reaction-diffusion processes on unstructured meshes.
SIAM Journal on Scientific Computing ,31:1774–1797, 2009.[6] R. Erban and S. J. Chapman. Reactive boundary conditions for stochastic simulationsof reaction-diffusion processes.
Physical Biology , 4(1):16–28, 2007.[7] R. Erban and S. J. Chapman. Stochastic modelling of reaction-diffusion processes:algorithms for bimolecular reactions.
Physical Biology , 6(4):046001, 2009.[8] R. Erban, S. J. Chapman, and P. Maini. A practical guide to stochastic simulationsof reaction-diffusion processes. 35 pages, available as http://arxiv.org/abs/0704.1908,2007.[9] R. Erban, M. Flegg, and G. Papoian. Multiscale stochastic reaction-diffusion modelling:application to actin dynamics in filopodia. to appear in the
Bulletin of MathematicalBiology , 2013.[10] M. Flegg, J. Chapman, and R. Erban. The two-regime method for optimizing stochasticreaction-diffusion simulations.
Journal of the Royal Society Interface , 9(70):859–868,2012.[11] M. Flegg, J. Chapman, L. Zheng, and R. Erban. Analysis of the two-regime methodon square meshes. submitted to
SIAM Journal on Scientific Computing , 23 pages,available as http://arxiv.org/abs/1304.5487 2013.[12] M. Flegg, S. R¨udiger, and R. Erban. Diffusive spatio-temporal noise in a first-passagetime model for intracellular calcium release. to appear in the
Journal of ChemicalPhysics , 2013.[13] E. Flekkøy, J. Feder, and G. Wagner. Coupling particles and fields in a diffusive hybridmodel.
Physical Review E , 64:066302, 2001.[14] B. Franz, M. Flegg, J. Chapman, and R. Erban. Multiscale reaction-diffusion algo-rithms: PDE-assisted Brownian dynamics. available as http://arxiv.org/abs/1206.5860,to appear in SIAM Journal on Applied Mathematics, 2013.[15] T. Geyer, C. Gorba, and V. Helms. Interfacing Brownian dynamics simulations.
Journalof Chemical Physics , 120(10):4573–4580, 2004.2216] M. Gibson and J. Bruck. Efficient exact stochastic simulation of chemical systems withmany species and many channels.
Journal of Physical Chemistry A , 104:1876–1889,2000.[17] D. Gillespie. Exact stochastic simulation of coupled chemical reactions.
Journal ofPhysical Chemistry , 81(25):2340–2361, 1977.[18] C. Gorba, T. Geyer, and V. Helms. Brownian dynamics simulations of simplified cy-tochrome c molecules in the presence of a charged surface.
Journal of Chemical Physics ,121(1):457–464, 2004.[19] J. Hattne, D. Fange, and J. Elf. Stochastic reaction-diffusion simulation with MesoRD.
Bioinformatics , 21(12):2923–2924, 2005.[20] A. Hellander, S. Hellander, and P. L¨otstedt. Coupled mesoscopic and microscopic sim-ulation of stochastic reaction-diffusion processes in mixed dimensions.
Multiscale Mod-eling and Simulation , 10(2):585–611, 2012.[21] M. Klann, A. Ganguly, and H. Koeppl. Hybrid spatial Gillespie and particle trackingsimulation.
Bioinformatics , 28(18):i549–i555, 2012.[22] J. Lipkova, K. Zygalakis, J. Chapman, and R. Erban. Analysis of Brownian dynamicssimulations of reversible bimolecular reactions.
SIAM Journal on Applied Mathematics ,71(3):714–730, 2011.[23] E. Moro. Hybrid method for simulating front propagation in reaction-diffusion systems.
Physical Review E , 69:060101, 2004.[24] S.A. Rice.
Diffusion Limited Reactions . Amsterdam: Elsevier, 1 edition, 1985.[25] J. van Zon and P. ten Wolde. Green’s-function reaction dynamics: a particle-basedapproach for simulating biochemical networks in time and space.
Journal of ChemicalPhysics , 123:234910, 2005.[26] G. Wagner and E. Flekkøy. Hybrid computations with flux exchange.