The blending region hybrid framework for the simulation of stochastic reaction-diffusion processes
Christian A. Yates, Adam George, Armand Jordana, Cameron A. Smith, Andrew B. Duncan, Konstantinos C. Zygalakis
aa r X i v : . [ q - b i o . Q M ] S e p The blending region hybrid framework for the simulation ofstochastic reaction-diffusion processes
Christian A. Yates , ∗ , Adam George , Armand Jordana , Cameron A. Smith , Andrew B.Duncan , Konstantinos C. Zygalakis ∗ ∗ E-mail: [email protected]; [email protected]
Abstract
The simulation of stochastic reaction-diffusion systems using fine-grained representations canbecome computationally prohibitive when particle numbers become large. If particle numbersare sufficiently high then it may be possible to ignore stochastic fluctuations and use a moreefficient coarse-grained simulation approach. Nevertheless, for multiscale systems which exhibitsignificant spatial variation in concentration, a coarse-grained approach may not be appropri-ate throughout the simulation domain. Such scenarios suggest a hybrid paradigm in which acomputationally cheap, coarse-grained model is coupled to a more expensive, but more detailedfine-grained model enabling the accurate simulation of the fine-scale dynamics at a reasonablecomputational cost.In this paper, in order to couple two representations of reaction-diffusion at distinct spatialscales, we allow them to overlap in a “blending region”. Both modelling paradigms provide avalid representation of the particle density in this region. From one end of the blending regionto the other, control of the implementation of diffusion is passed from one modelling paradigmto another through the use of complementary “blending functions” which scale up or down thecontribution of each model to the overall diffusion. We establish the reliability of our novelhybrid paradigm by demonstrating its simulation on four exemplar reaction-diffusion scenarios.Key index words: hybrid modelling, stochastic reaction-diffusion, multiscale modelling, par-tial differential equation, hybrid modelling framework
Many biological and physical systems are inherently multiscale in nature [59, 72, 46, 38, 9, 22].The modelling of such systems therefore requires multiscale representations which, by theirnature, are not well captured using a single modelling paradigm. There is a trade-off between,on the one hand, ensuring that models are sufficiently detailed that they accurately captureknown biological and physical phenomena of interest and, on the other, achieving model outputsin a timely manner.The appropriate representation of travelling waves of cells in developmental or maintenancecontexts is a classic example of a multiscale phenomenon for which the trade-off between cheap-but-coarse and expensive-but-accurate modelling paradigms is evident. For a pulled wave-frontthe wave speed is determined by the low-density dynamics at the front of the wave [39]. It istherefore important to represent cell movement and proliferation dynamics at the front usingan appropriately detailed model. A model that is too coarse may neglect important features ofthe real process. Behind the wave, cell density is higher making a fine-grained representationmore computationally expensive. Since the fine details are less important in this region we cansubstitute the more detailed model for a cheaper, coarser representation. Coupling modellingregimes at different scales is an open question to which a variety of solutions have previouslybeen proposed [76, 45, 65, 58, 29, 23, 57, 44, 8, 21, 19, 56, 20, 9, 30, 40, 61, 24, 25, 28, 1, 2, 52,43, 18, 66, 73]. For more details on the different types of hybrid methods available we direct theinterested reader to [60].In this paper we focus on the three main modelling paradigms used for representing reaction-diffusion systems. At the coarsest scale (which we refer to as the macroscopic scale) we representthe concentration of reactant species by partial differential equations (PDEs) [35, 37, 36, 68,51, 50, 32]. For validity, these models typically require high concentrations since assumptionsunderlying the use of PDEs break down for low copy numbers. Continuum models such as thesecan usually be simulated extremely efficiently using a wide variety of well-established numericalmethods, however, they lack the realism of finer-scale models.At the next level down, the mesoscopic scale, reactant species are represented as individualparticles and are compartmentalised into contiguous, non-overlapping subdivisions of space [12,11, 34, 16, 4, 74, 75, 46]. Particles are assumed to be well-mixed within a compartment andcan interact with others in their compartment. These models can capture stochasticity in thebehaviour of the particles and can be simulated efficiently when copy numbers are low. However,when particle numbers become large, simulations can become prohibitively slow in comparisonto macroscale representations. They also lack the accuracy of more fine-grained models sincethe individual particle identities and positions are not retained.The finest representation we consider is Brownian-dynamics models at the microscopic scale[3, 42, 14, 13]. In these models the trajectories of all particles are simulated (typically usinga discrete fixed time-step paradigm) in continuous space [63, 3, 71, 64]. For a system of N particles, an appropriate simulation algorithm must generate ZN Gaussian random variables(where Z is the dimension of the system) in order to update the particle positions. For simu-lations incorporating pairwise interactions, N pairwise distances must also be updated at eachtime-step . Consequently, these methods can be extremely computationally intensive. They do,however, provide a comprehensive and accurate individual representation capable of incorpo-rating stochasticity into particle positions and interaction times. More details on the specificimplementation of each of these three modelling paradigms will be given in the next section.In general, the aim of a hybrid method is to exploit the complementary advantages and negate Note that by a careful partitioning of space the number of comparisons can be reduced dramatically to almost O ( N ) when particles are only compared with others in their local neighbourhood [55]. the complementary weaknesses of models at different scales. Using a coarse, cheap representationin a region of space in which particle density is high allows for significant computational savings incomparison to the purely fine-scale simulation. Conversely, implementing a fine-scale individual-based representation in regions in which low-copy number effects are of paramount importancecan give significant improvements in accuracy in comparison to coarser models. Consequently,one way to achieve accurate simulations that are also computationally tractable is to combinethe models’ strengths in a hybrid representation.In this paper we propose a novel hybrid method for coupling PDEs at the macroscale tocompartment-based models at the mesoscale and a related novel hybrid method for couplingcompartment-based models at the mesoscale to Brownian-dynamics models at the microscale.In each case, the coarser regime is coupled to the finer regime through an overlap region. In thisoverlap region, which from now on we will refer to as the blending region , both representationsof the reaction-diffusion dynamics are valid. In the blending region the strength of diffusionfor each model is determined by a spatially-varying blending function which is prescribed to beunity on one end of the overlap region and zero on the other. The blending functions for thetwo models are complementary so that the sum of the two blending functions at any point inthe domain is equal to unity. These functions control the relative contribution of each modelto the diffusion dynamics. This approach is reminiscent of that taken by [10] in a non-spatialcontext. In [10] two different non-spatial models for stochastic chemical kinetics were coupledin copy-number space through a blending region in which both models co-existed.The remainder of the paper is organised as follows. In Section 2 we describe the individualreaction-diffusion models that we couple together and provide a brief justification for why themodels can be considered “equivalent” and hence are suitable candidates for coupling. In Section3 we present the mechanics of the two hybrid blending methods and prove their effectiveness,in Section 4, by simulating a number of test scenarios and determining whether any bias isintroduced by the blending methods. We conclude in Section 5 with a short summary of ourfindings and suggestions for extensions to this work. Within this section, we describe the three different modelling scales that we will couple in orderto create two distinct spatially-coupled hybrid methods. In Section 2.1 we describe a generalmacroscale PDE for reaction-diffusion systems with a single species, as well as different numericalapproaches for its solution. Section 2.2 contains a discussion of mesoscale compartment-basedmodels and their simulation, while in Section 2.3 we introduce the microscale individual-baseddynamics. In Section 2.4 we briefly discuss how each of these representations of reaction-diffusionprocesses at different scales might be considered to be equivalent in an appropriate limit.
Partial differential equations, the macroscale models we employ in this paper, can be consideredto be appropriate representations of the mean behaviour of particles at high concentrations.The primary advantage of the PDE representation is that there exists a wide range of well-established and well-understood tools for their numerical simulation. In rare, simple cases,PDEs are amenable to mathematical analysis. However, they typically fail to model low copynumber behaviour.A generic PDE which describes the spatio-temporal evolution of the concentration of a singlespecies, c ( x , t ), at position x and time t takes the form: ∂c∂t ( x , t ) = ∇ · ( D ( x ) ∇ c ( x , t )) + R ( c ( x , t ) , x , t ) , x ∈ R Z , t ∈ [0 , T ] , (1)where consistent initial and boundary conditions need also to be specified. Here reactionsare represented by the function R , Z is the dimension of space and T is the final time towhich we wish to evolve the solution. Note that the spatially varying diffusion coefficient,represented by D ( x ), sits inside the first derivative, but not the second. As noted by [69],there is no canonical choice of operator describing spatially dependent diffusion. In physicalapplications the form of the macroscopic diffusion equation should be dictated by the underlyingmicroscopic or mesoscopic process. Since the spatial dependence of the diffusion coefficient inour hybrid methods is introduced purely as a modelling convenience and does not correspondto any microscopic or mesoscopic ground truth, we are effectively free to choose the form ofthe diffusion operator. We adopt the form considered by [5] (see equation (1)). We choose thetransition rates in the corresponding compartment-based representation (see Section 2.2) andthe drift and diffusion coefficients of the corresponding microscopic position evolution equation(see Section 2.3) so that diffusion in the overlap regions of the hybrid methods satisfies the sameform of non-constant coefficient diffusion equation.For the majority of this paper we focus on the following one-dimensional PDE in the regionΩ = [ a, b ]: ∂c∂t = ∂∂x (cid:18) D ( x ) ∂c∂x (cid:19) + R ( c ( x, t )) , (2)with constant flux boundary conditions D ( a ) ∂c∂x (cid:12)(cid:12)(cid:12) x = a = J a , D ( b ) ∂c∂x (cid:12)(cid:12)(cid:12) x = b = J b . (3)For a discussion of the implementation of the numerical solution of the PDEs employed in thispaper please refer to Appendix A. Note that there is no explicit spatial dependence in thereaction term in equation (2). Compartment-based methods are coarse-grained stochastic representations. The spatial domainis typically divided into compartments, each of size h , in which particles are assumed to be well-mixed. The reaction-diffusion dynamics are characterised by a set of possible events. Eventsare either reactions, in which particles can interact with others within their own compartmentaccording to some prespecified reaction rates, or jumps to adjacent compartments with rateswhich depend on the macroscopic diffusion coefficient, D ( x ), and the compartment size, h .Specifically, in order to capture diffusion which corresponds to the macroscopic equation (1) wemust choose the rates of jumping to be different depending on the direction of the jump (seeequations (39) and (40) for more detail).Throughout this paper we refer to models at this scale as mesoscopic or compartment-based .For a discussion of the implementation of the numerical simulation of the compartment-basedmodels employed in this paper please refer to Appendix B. Individual-based methods require the recording and updating of large numbers of particles’positions. Relative positions for each pair of particles must also be maintained at every stepif higher-order reactions (higher than first-order) or volume-exclusion are to be modelled. Forlarge particle numbers, N , the O ( N ) computational complexity means that individual-basedsimulation algorithms can become extremely expensive .In what follows we employ a fixed-time-step algorithm, although we note that continuous-time algorithms for Brownian reaction-diffusion dynamics are also available [71]. The evolutionof particle i ’s position, y i ( t ), between times t and t + ∆ t b in the case of space-dependent diffusion(corresponding to PDE (2) and compartment-based jump-rates given by equations (39) and (40))can be simulated according to the following discrete-time update equation y i ( t + ∆ t b ) = y i ( t ) + ∆ t b d D ( x )d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = y i ( t ) + q D ( y i ( t ))∆ t b ξ i , (4)where ξ i ∼ N (0 ,
1) is a Gaussian random variable with mean 0 and variance 1. If required,reactions can be implemented according to a variety of different algorithms [71, 63]. In thispaper, we employ the λ - ρ method [14]. If two eligible particles come within a reaction radius, ρ ,of each other they interact with a given rate, λ , according to the appropriate reaction pathway.We refer to these models at this scale as off-lattice , microscopic or individual-based modelsin what follows. In attempting to couple together representations of the same phenomenon at different scales weneed to ensure that, under certain assumptions, they are representations of the same process.Pioneering work in establishing the connection between stochastic and deterministic models wasundertaken by [27], [70] and [41]. In this section we concisely summarise the ways in which themodels outlined above can be considered to be equivalent and direct the interested reader toresources which contain more detailed arguments.In order to transition from the mesoscale to the macroscale, we can first use the reaction-diffusion master equation to derive the deterministic mean-field representation of the compartment-based particle numbers [14, 4, 48]. It should be noted that for second- and higher-order reactions,the mean-field equations are only approximations of the true mean behaviour of the stochas-tic system [15]. Taking the diffusive limit of the mean-field equations gives a correspondingreaction-diffusion PDE.The Fokker-Planck equation can be used to connect a microscale stochastic differential equa-tion (SDE) model of diffusion to a macroscale model describing the evolution of the probability As previously noted some of this complexity can be offset by a careful partitioning of space allowing particlesto be compared only with others in their local neighbourhood [55]. density of a particle’s position [15, 54]. For example, the canonical diffusion equation is themacroscopic Fokker-Planck equation corresponding to non-interacting particles undergoing sim-ple Brownian motion.Although we do not use this macroscopic-microscopic coupling directly in this work, weemploy it indirectly in order to link the microscopic and mesoscopic descriptions together throughtheir connection to the same PDE. Alternatively, first-passage time theory can be applied to aparticle which moves subject to a given SDE in order to derive jump rates between neighbouringcompartments in a compartment-based representation [53, 74]. Connections between the modelsat microscale and mesoscale are stated more rigorously by [33].
In this section we discuss the two main algorithms of this paper. In particular, in Section 3.1we present the central unifying idea behind both of our hybrid methods. The methods can bothbe understood as operator-splitting algorithms in which, in a central overlap region between thetwo regimes, diffusion is dealt with by both regimes using spatially varying diffusion coefficients.We discuss how to couple the methods discussed in Section 2 in order to accommodate thissplit-diffusion paradigm. In Section 3.2 we give the specific details of how to convert mass fromone modelling regime to another to ensure both models are synchronised and valid representa-tions of the particle density in the blending region. We then present, in Section 3.3, a genericalgorithm for coupling the PDE with the compartment-based approach, as well as a similarlygeneral algorithm for coupling the compartment-based approach with Brownian dynamics. Weemphasise that the generic methods we present for coupling two regimes are independent of thenumerical implementations chosen to simulate each regime. However, for ease of use and repro-ducibility we have provided details of the numerical implementations we chose in AppendicesA-C.
In order to illustrate the conceptual framework behind our algorithms we consider the followingconstant coefficient diffusion PDE in Ω = [ a, b ]: ∂c∂t = ∂∂x (cid:18) D ∂c∂x (cid:19) , (5)with the following zero-flux boundary conditions: D ∂c∂x (cid:12)(cid:12)(cid:12) x = a = D ∂c∂x (cid:12)(cid:12)(cid:12) x = b = 0 . (6)Divide the domain, Ω, into three subdomains Ω = [ a, I ] , Ω = [ I , I ] , Ω = [ I , b ] and writethe constant diffusion coefficient D = D ( x ) + D ( x ) where D ( x ) = D, a ≤ x < I ,f ( x ) , I ≤ x < I , , I ≤ x ≤ b, (7)and D ( x ) = , a ≤ x < I ,f ( x ) , I ≤ x < I ,D, I ≤ x ≤ b, (8)where f and f are monotonically decreasing/increasing functions, respectively, with f ( x ) = D − f ( x ) and f ( I ) = f ( I ) = D and f ( I ) = f ( I ) = 0 in order to ensure continuity of D and D across Ω.Equation (5) can now be written as ∂c∂t = ∂∂x (cid:18) D ( x ) ∂c∂x (cid:19)| {z } + ∂∂x (cid:18) D ( x ) ∂c∂x (cid:19)| {z } , (9)with corresponding boundary conditions( D ( x )+ D ( x )) ∂c∂x (cid:12)(cid:12)(cid:12) x = a = D ∂c∂x (cid:12)(cid:12)(cid:12) x = a = 0 and ( D ( x )+ D ( x )) ∂c∂x (cid:12)(cid:12)(cid:12) x = b = D ∂c∂x (cid:12)(cid:12)(cid:12) x = b = 0 . (10)In addition we specify the initial condition c ( x,
0) = c ( x ). It is straightforward to show that,because D ( x ) = 0 in [ I , b ], the operator indicated by 1 in equation (9) does not influence theconcentration of c in that region. In a similar way, because D ( x ) = 0 in [ a, I ], the operatorindicated by 2 in equation (9) does not influence the concentration of c in that region. Now let φ τ , φ τ be the flow maps associated with the propagation of the operators 1 and 2 in equation(9) until time τ . Specifically this means that the solution of the following equations ∂c (1) ∂t = ∂∂x D ( x ) ∂c (1) ∂x ! , D ( a ) ∂c (1) ∂x (cid:12)(cid:12)(cid:12) x = a = D ( I ) ∂c (1) ∂x (cid:12)(cid:12)(cid:12) x = I = 0 , (11a) ∂c (2) ∂t = ∂∂x D ( x ) ∂c (2) ∂x ! , D ( I ) ∂c (2) ∂x (cid:12)(cid:12)(cid:12) x = I = D ( b ) ∂c (2) ∂x (cid:12)(cid:12)(cid:12) x = b = 0 , (11b)subject to initial conditions c ( i ) ( x,
0) = c ( i )0 ( x ) can be written as c ( i ) ( x, τ ) = φ ( i ) τ ( c ( i )0 )( x ) , for i =1 ,
2, respectively .The idea behind splitting methods is that one can now obtain an approximation for thesolution of equation (9) at time τ by using an appropriate composition of the flow maps φ (1) τ and φ (2) τ . In particular, the simplest splitting method is given by c ( x, τ ) ≈ ( φ (1) τ ◦ φ (2) τ )( c )( x ) , (12)where we note that the ordering of the composition is unimportant.At a first glance this seems like an unnecessarily complicated approach for obtaining anapproximation for the solution of equation (5). However, choosing the flow maps φ τ and φ τ to represent propagation operators for two different model types allows us to seamlessly blend Note that due to the choice of blending functions the boundary condition at I in (11a) and at I in (11b)are automatically satisfied. the distinct numerical update rules of the different modelling regimes described in Section 2.For example, when coupling the PDE to the compartment based model, φ τ might represent anupdate operator for the numerical solution of the PDE up to time τ , whilst φ τ might representsteps of the position-jump Markov processes described in Section 2.2 up until time τ .Due to the properties of the diffusion functions D i ( x ), the two models only co-exist in theblending region [ I , I ]. Therefore, in applying the operator splitting update illustrated in equa-tion (12), we only need to worry about how the concentration of the numerical solution of thePDE in the blending region translates to particle numbers for the compartment-based approachand vice versa. We must ensure that any PDE solution update in the blending region imple-mented by operator φ τ is also reflected in the compartment-based solution. Equivalently, anyupdate to the compartment-based solution in the blending region implemented via φ τ must bereflected in the PDE solution. In a similar way, when coupling the compartment-based model toBrownian dynamics, one need only worry about how the particle numbers for the compartment-based approach in the blending region impact on the particle positions of the off-lattice Browniandynamics and vice versa. Outside the two blending regimes the two representations are effec-tively decoupled in terms up their update operators. In this section we illustrate how to couple two distinct representations of reaction-diffusionprocesses in the blending region. First we tackle a PDE-compartment-based hybrid pairing,followed by a coupling between compartment-based and Brownian-based particle dynamics.
Conversion between PDE and compartment-based model:
We assume that the nu-merical solution of the PDE is calculated on the discrete mesh (see figure 7 in Appendix A foran illustration) of size ∆ x in [ a, I ] and that compartment-based dynamics are simulated withcompartments of size h in [ I , b ]. It is natural to assume that h ≥ ∆ x , as a fine discretisation ofthe PDE mesh is required in order to minimise the error between the numerical solution and theexact solution it approximates. Note, however, that this is not a limitation of our algorithm andthat h ≤ ∆ x would also be possible. There are n = I − I ∆ x PDE solution voxels in the overlapregion [ I , I ] and n = I − I h compartments in the same region, where n , n ∈ N . For ease ofcomputation we assume that n = γn , with γ ∈ N so that there are an integer number of PDEsolution voxels per compartment. There are also n p = ( I − a ) / ∆ x PDE solution voxels in thepurely PDE region, [ a, I ], and n c = ( b − I ) /h compartments in the purely compartment-basedregime [ I , b ]. The numerical solution of the PDE in voxel i is labelled q i for i = 1 , . . . , n p + n and the number of particles in compartment i is labelled C i for i = 1 , . . . , n + n c .In each time interval of length ∆ t p we assume, without loss of generality, that the PDEsolution is updated first and the compartment-based solution second. After the propagation ofthe discrete PDE solution operator in the time interval [ t, t +∆ t p ], assume that the concentrationsin PDE voxels of the blending region have changed. Consequently it is necessary to modify the Note that we describe the coupling between the two regimes in the blending region using the terminologyof the finite volume PDE discretisation that we employ in our numerical examples (see Section 4). However,we also note that finite volume voxels can be substituted for finite difference or finite element mesh points in astraightforward manner. corresponding compartment-based description in the blending region [ I , I ] before propagatingthe compartment-based model in the region [ I , b ]. More precisely, for compartment i in theblending region, set C i = γ X j =1 q n p + γ ( i − j ∆ x. (13)Because we are required to synchronise the representations of the solutions in the two regimesaccording to equation (13), the number of particles contained in the i -th compartment in theblending region is no longer an integer. Nevertheless, when it comes to performing the stochasticsimulation algorithm we work with these non-integer values to calculate the time until the nextevent. This could potentially be an issue when the copy numbers in a compartment are low, butarguably this would imply that we were using the PDE description to represent concentrationsin a region of the domain for which this is not appropriate. A similar synchronisation is imple-mented once the compartment-based model has been propagated and the number of particlesin the blending region has changed. In particular, if δC i corresponds to the integer change inparticle numbers in the compartment i in the blending region, then one adds uniformly δC i /γ ∆ x to the PDE solution in each of the PDE voxels, i.e q n p + γi + j = q n p + γi + j + δC i γ ∆ x , j = 1 , · · · , γ. (14)Reactions in the blending region are always implemented according to the compartment-based paradigm. If reactions occur then particle numbers in compartments are updated and thecorresponding change is also implemented in the appropriate PDE voxels, as in equation (14). Conversion between compartment based and individual particle models
Withoutloss of generality assume that the compartment-based model is employed in [ a, I ] and theBrownian-based model is employed in [ I , b ] with the two models being simultaneously employedin the blending region [ I , I ]. Compartment-based dynamics are simulated with compartmentsof size h in [ a, I ]. There are n c = ( I − a ) /h compartments in the purely compartment-basedregion, [ a, I ], and n = ( I − I ) /h compartments in the overlap region [ I , I ]. The number ofparticles in compartment i is, as before, labelled C i for i = 1 , . . . , n c + n . Brownian particlesare simulated off-lattice with positions updated according to the dicretised SDE (4) in [ I , b ].In each time interval of length ∆ t b assume, without loss of generality, that the compartment-based solution is updated first, followed by the Brownian-based dynamics. During the prop-agation of the compartment-based solution it is likely that the numbers of particles in thecompartments of the blending region have changed. Consequently we need to alter the po-sitions of Brownian particles in the blending region. If a particle jumps from compartment i to a neighbouring compartment j in the hybrid region, then we select a Brownian particleuniformly at random from amongst the particles which currently reside in compartment i andmove it a distance ± h with the sign of the displacement corresponding to the direction of thecompartment-based particle’s jump i.e. y k = y k ± h, (15)where k indexes the randomly selected Brownian particle from compartment i .0 (a) (b) Figure 1.
Schematic representations of (a) the PDE-compartment hybrid and (b) thecompartment-Brownian hybrid. In panel (a) the green curve in the green region [ a, I ]represents the PDE solution in the purely PDE region of the domain. The red curve and thered boxes represent equivalent PDE- and compartment-based representations of the mass inthe red blending region. The blue boxes in the blue region of the domain represent the numberof particles in each compartment in the purely compartment region of the domain. In panel(b) the blue boxes in the blue region of the domain represent the number of particles in eachcompartment in the purely compartment region of the domain. The red boxes and the redcircles represent equivalent compartment- and Brownian-based representations of the mass inthe red blending region. The yellow circles in the yellow region of the domain representindividual particles in the purely Brownian region of the domain. Note that we have giveneach Brownian particle a different height to aid clarity of visualisation, but in reality allparticles lie on the x -axis in these one-dimensional simulations.If a particle in compartment n c + 1 (the first compartment in the blending region) jumpsleftwards out of the blending region (according to the compartment-based jump rates) and intothe purely compartment-based region then a Brownian particle in the compartment n c + 1 isselected uniformly at random and removed from the simulation (as well as particle numbers inthe affected compartments being updated). Conversely, if a compartment-based particle jumpsto the right, out of the last compartment in the purely compartment-based regime into the firstcompartment in the blending region, then a Brownian particle is added with its position chosenuniformly at random in this compartment, [ I , I + h ] (as well as particle numbers in the affectedcompartments being updated). Note that the jump rates in the compartment-based model,which implement diffusion corresponding to equation (2), are such that, with our chosen blendingdiffusion coefficients, the rate of jumping to the right out of the final compartment is zero, sothat no compartment-based particles can erroneously jump into the purely-Brownian regime.Similarly, the diffusion coefficient of the Brownian particles at the pure-compartment/blendingregion interface is zero. Technically, with our finite time-step implementation of diffusion itmight be possible for Brownian particles to erroneously jump over the interface into the purely1compartment-based regime . On the rare occasions that a Brownian particle is chosen to jumpover the interface (as an artefact of the numerical implementation) we simply reflect it back intothe blending region. Since the diffusion coefficient is low in the boxes close to the interface thisvery rarely happens, and when it does the error caused by reflecting the particle is minimal.Once the particle-based method has been propagated, it is usually necessary to update thenumber of particles in the compartments of the blending region, C i for i = n c + 1 , . . . , n c + n .Rather than tracking every Brownian-particle movement to see whether it has crossed overa compartment boundary, instead we simply sum the number of Brownian-based particles ineach compartment at the end of the Brownian update to find the numbers of particles in eachcompartment of the blending region: C n c + i = N X k =1 I y k ∈ [ I +( i − h,I + ih ] , for i = 1 , . . . , n , (16)where I y ∈ [ I +( i − h,I + ih ] is the indicator function which takes the value 1 if the Brownian particlelies in the ( n c + i )th compartment and 0 otherwise.Reactions in the blending region (similarly to the PDE-compartment hybrid method) arealways implemented using the compartment-based paradigm. If a reaction occurs in the hybridregion then the appropriate Brownian particles are added (with positions chosen uniformlyat random across the corresponding compartment) or removed (with the particle(s) selecteduniformly at random from amongst those in the compartment). Having established the conversion rules in the previous section we are now in the position topresent two hybrid algorithms. In particular, Algorithm 1 is the algorithm that couples diffusionin the PDE and compartment-based models, while Algorithm 2 is the algorithm that couplesdiffusion in the compartment-based models with Brownian-based dynamics. We have presentedboth of these algorithms with maximum generality in order to emphasise that the specific simu-lation methodologies are not important. In the next section we implement these algorithms witha finite volume PDE solver, the spatial Gillespie algorithm for compartment-based dynamics andthe λ − ρ Brownian reaction-diffusion paradigm for the Brownian dynamics, in order to provideconcrete examples of their implementation. Algorithms for the implementation of these threemethods are given in Appendices A, B and C respectively.
In this section we demonstrate that our proposed algorithms correctly simulate four test prob-lems of increasing complexity. The first two problem are simulations of pure diffusion withdifferent initial conditions, demonstrating that the fluxes over the interface of the hybrid modelare consistent with the expected behaviour of the finer-scale representation in each hybrid model. Whilst there do exist integrators for diffusion processes which can guarantee that this situation does nothappen [6], implementing such an approach is beyond the scope of the article. Algorithm 1:
Coupling a PDE solution with a compartment-based approach
Input:
PDE mesh size – ∆ x ; compartment size – h ; time-step for the solution of the PDE– ∆ t p ; left and right ends of the blending region – I , I ; initial concentration forthe PDE – c init ; initial particle numbers – C init ; final time – T . Set t = 0. while t < T do Simulate diffusion due to the PDE in [ a, I ] (and reactions due to the PDE in [ a, I ])between t and t + ∆ t p for diffusion coefficient given by D ( x ) using an appropriatenumerical solver. Update the compartment-based particle numbers in [ I , I ] according to equation (13). Simulate diffusion and reactions due to the compartment-based approach in [ I , b ]between t and t + ∆ t p for diffusion coefficient given by D ( x ) using an appropriatestochastic simulation algorithm, taking as an initial condition the updated particlenumbers from line 4. Update the PDE solution in [ I , I ] according to equation (14). Set t = t + ∆ t p . end The third problem, one of morphogen gradient formation, evidences the successful implementa-tion of reactions in our hybrid algorithms. Finally, in the fourth test problem we implement asecond-order reaction system in three dimensions, demonstrating the applicability of the methodto more complicated scenarios.For each of the first three test problems, the one-dimensional domain we employ is Ω =[ a, b ] = [0 , I = 1 / I = 2 /
3. The remainder of the parameter values for examples1 and 2 are specified in table 1, for example 3 in table 2 and for example 4 in table 3. Theblending functions for these three problems (and by simple extension for the fourth problem)are defined as the simple linear functions f ( x ) = 2 − x, (17) f ( x ) = 3 x − , (18)which scale the contribution of each method to the diffusion coefficient linearly between 0 and1 across the blending region. These, in conjunction with equations (7) and (8), define thediffusion coefficients for both regimes across the whole domain. For each of the first threeexample and for both couplings we will quantify the qualitative comparisons (provided by densitycomparison snapshots) with error plots displaying the evolution of the difference between theaveraged profiles of our hybrid methods the mean-field PDE (see equations (19)-(22)). In thefourth example (for which the PDE is not an exact description of the mean behaviour of theindividual-based methods) we will compare the averaged profiles of our hybrid methods withthe averaged profiles of the finer scale ‘ground truth’ (e.g. mesoscale or miscroscale) simulations(see equations (31)-(36)).3 Algorithm 2:
Coupling a compartment-based approach with Brownian dynamics.
Input:
Compartment size – h ; time-step for the solution update of the Browniandynamics – ∆ t b ; left and right ends of the blending region – I , I ; initial particlenumbers – C init ; initial Brownian particle positions – y ; final time – T . Set t = 0. while t < T do Simulate diffusion and reactions due to the compartment-based approach in [ a, I ]between t and t + ∆ t b for diffusion coefficient given by D ( x ) using an appropriatestochastic simulation algorithm. Update the Brownian particle positions in [ I , I ] according to equation (15) (ifappropriate). Simulate diffusion due to the Brownian particle dynamics in [ I , b ] (and reactions dueto the Brownian particle dynamics in [ I , b ]) between t and t + ∆ t b for diffusioncoefficient given by D ( x ) using an appropriate numerical solver, taking as an initialcondition the updated concentration from line 4. Update the compartment-based particle numbers in [ I , I ] according to equation (16). Set t = t + ∆ t b . end Parameter Value Description N ,
1] Spatial domain D K
20 Number of compartments h /
30 Compartment width∆ x /
300 PDE voxel width∆ t p − PDE time-step∆ t b − Brownian time-stepM 500 Number of repeats
Table 1.
Table of Parameter values used for the pure diffusion simulation of test problems 1and 2.
The first test of our hybrid algorithms is to determine whether, when simulating diffusion, theyare capable of maintaining the uniform steady state distribution across the domain withoutintroducing any bias. We initialise particles uniformly across the domain and implement zero-flux boundary conditions.In figure 2 (as well as for figures 3-4) the top three figures are for the PDE-compartmentcoupling and the bottom three figures for the compartment-Brownian coupling. The left-mostpanels display the density profile of the hybrid methods at time t = 0 . x P a r t i c l e d e n s i t y (a) x P a r t i c l e d e n s i t y (b) Time −0.0050.0000.005 R e l a t i v e e rr o r PDE subdomainBlended subdomainCompartment subdomain (c) x P a r t i c l e d e n s i t y (d) x P a r t i c l e d e n s i t y (e) Time −0.0050.0000.005 R e l a t i v e e rr o r Compartment subdomainBlended subdomainParticle subdomain (f)
Figure 2.
Density and error plots for test problem 1 - pure diffusion with a uniform initialcondition. Panels (a)-(c) are for the PDE-compartment hybrid method and (d)-(f) for thecompartment-Brownian hybrid method. Panels (a) and (d) are snapshots at time 0.1, and (b)and (e) at time 1. In panels (a) and (b) the green line is the PDE part of the hybrid method,the red bars represent the number of particles in each compartment in the blending region andthe blue bars represent the number of particles in each compartment in the purelycompartment-based region. In panels (d) and (e) the blue bars represent the number ofparticles in each compartment in the purely compartment-based region, the red bars representthe number of particles in each compartment in the blending region and the yellow bars thenumber of particles (appropriately binned for visualisation purposes) in the purely Brownianregion. In all four density comparison panels the black dashed line represents the analyticalsolution of the diffusion equation with given initial condition. Vertical red lines mark theposition of the interfaces. Panel (c) shows the relative error (described in the main text)between the density given by PDE-compartment hybrid method and the density given by theanalytical solution of the diffusion equation with the same initial condition. Similarly panel (f)shows the relative error (described in the main text) between the density given bycompartment-Brownian hybrid method and the density given by the analytical solution of thediffusion equation with the same initial condition. Results shown are for N = 1000 particlesand are averaged over 500 repeats. All other parameters are given within table 1.the density profile at t = 1. In both left and middle panels the mean-behaviour of the stochasticmodel simulated across the whole of the domain is displayed as a black, dashed line for compari-son. The right-most panels display the evolution through time of the relative mass error of eachregion of the domain: [ a, I ], [ I , I ] and [ I , b ]. For the PDE-compartment coupling the relativemass error (RME) is the difference between the average (over 500 repeats - unless otherwisestated) number of particles in the given region in the hybrid method and the correspondingnumber in the same region in the analytical solution of the PDE, u ( x, t ), divided by the number5of particles in the relevant region of the analytical solution of the PDE (to normalise): RM E P ( t ) = R Ω P ¯ c ( x, t ) dx − R Ω P u ( x, t ) dx R Ω P u ( x, t ) dx , (19) RM E H ( t ) = P i ¯ C i ( t ) I c i ∈ Ω H − R Ω H u ( x, t ) dx R Ω H u ( x, t ) dx , (20) RM E C ( t ) = P i ¯ C i ( t ) I c i ∈ Ω C − R Ω C u ( x, t ) dx R Ω C u ( x, t ) dx , (21)where Ω P = [ a, I ] is the purely PDE region of the domain, Ω H = [ I , I ] is the blending regionand Ω C = [ I , b ] is the purely compartment region of the domain. The averaged solution of thePDE component of the hybrid method at position x at time t is denoted ¯ c ( x, t ) and the averagedcompartment particle numbers in voxel i of the hybrid method are denoted ¯ C i . The positions c i are the centres of the compartments.For the compartment-Brownian coupling the relative mass error is the difference between theaverage (over 500 repeats - unless otherwise stated) number of particles in each region given bythe hybrid method and the number of particles in the analytical solution of the mean-field PDEmodel in the corresponding region, divided by the number of particles in the relevant regionof the analytical solution of the PDE (to normalise). In the pure compartment and blendingregions these are given by equations (21) and (20) respectively, with the altered definition ofΩ C = [ a, I ] for equation (21). For the purely Brownian region the RME is given by RM E B ( t ) = ¯ B − R Ω B u ( x, t ) dx R Ω B u ( x, t ) dx , (22)where Ω B = [ I , b ] and ¯ B represents the average number of Brownian particles in the purelyBrownian regime.Figure 2 demonstrates that both of our hybrid blending methods pass this most-basic testof maintaining a uniform distribution across the domain. The interfaces between the differentmodelling regimes are effectively undetectable. Qualitatively, the density plots all show goodagreement between the hybrid methods and the analytical solution to the mean-field diffusionequation. This is confirmed by the relative error plots (panels 2(c) and 2(f)) which demonstratelow errors which fluctuate around zero with no discernible long-term bias. The second test problem is designed to determine whether the hybrid methods can cope withhigh levels of flux across their interfaces. As with the previous example, we model pure diffusionwith no reactions, but this time with a different initial condition. All the particles are distributeduniformly within [ a, I ] and the system is allowed to equilibrate. The results of these simulationsfor both the PDE-compartment hybrid and the compartment-Brownian hybrid are given in figure3. In figure 3 we have initialised the particles uniformly in the left-hand-most third of thedomain, corresponding to the purely PDE region in the PDE-compartment hybrid and the6 x P a r t i c l e d e n s i t y (a) x P a r t i c l e d e n s i t y (b) Time −0.010.000.01 R e l a t i v e e rr o r PDE subdomainBlended subdomainCompartment subdomain (c) x P a r t i c l e d e n s i t y (d) x P a r t i c l e d e n s i t y (e) Time −0.010.000.01 R e l a t i v e e rr o r Compartment subdomainBlended subdomainParticle subdomain (f)
Figure 3.
Density and error plots for test problem 2 - pure diffusion with a step functioninitial condition in [ a, I ]. Descriptions, including definitions of relative errors are as in figure 2.purely compartment-based region in the compartment-Brownian hybrid . As in test problem1, both of our hybrid methods correctly match the evolution of the density of the mean-fielddiffusion equation, as evidenced quantitatively by the relative error plots 3(c) and 3(f). The formation of a morphogen gradient from a uniform initial condition constitutes the thirdtest of our hybrid simulation algorithms. Particles are allowed to diffuse freely throughout thedomain and degrade at a rate µ . To counteract the degradation and ensure a non-trivial steadystate, particles are introduced at the left-hand boundary, x = a = 0, with flux DJ , and a zero-flux boundary condition is implemented at x = b = 1. Since the reactions we have introduced arefirst order, the continuum mean-field model corresponding to the described set up is governedby the following PDE: ∂c∂t = D ∂ c∂x − µc, for x ∈ (0 ,
1) and t ∈ (0 , T ) , (23)with boundary conditions ∂c∂x (0 , t ) = − J, ∂c∂x (1 , t ) = 0 , t ∈ (0 , T ) , (24)and initial condition c ( x,
0) = c , for x ∈ [0 , , where c = DJµ . (25) We see similarly agreeable results when particles are initialised in the third of the domain [ I , b ] correspondingto the purely compartment or purely Brownian regions respectively. N (0) 1000 Initial number of particlesΩ [0 ,
1] Spatial domain D J ,
000 Rate of influx at the left boundary µ
10 Rate of particle decay K
20 Number of compartments h /
30 Compartment width∆ x /
300 PDE voxel width∆ t p − PDE time-step∆ t b − Brownian time-stepM 1000 Number of repeats
Table 2.
Table of parameters for the morphogen gradient simulation (Test problem 3).The initial condition is chosen so that we begin with the same number of particles as there willbe at steady state, but distributed uniformly across the domain . The parameters we employfor the simulations shown in figure 4 are given in table 2. Specifically, influx parameter, J ,and degradation parameter, µ , are chosen to ensure an average of 1000 particles populating thedomain throughout the simulation.Figure 4 illustrates the solutions of our two hybrid methods and those of the correspondingmean-field model (given by equation (23)). As with the previous two test problems, qualita-tive density profiles are in close agreement and quantitative error plots show low error and nosustained bias about zero. The final scenario we will use to demonstrate the accuracy of our hybrid methods is a systemof diffusing particles interacting through the following pair of chemical reactions:2 A κ −→ ∅ , ∅ κ −→ A, (26)which occur within the cuboidal domain Ω ⊆ R of volume V , where Ω = [0 , × [0 , × [0 , I at x = 10 / I at x = 20 /
3. The Note that we have chosen this initial condition to ensure the PDE-compartment algorithm functions appro-priately. Whilst the compartment-to-Brownian algorithm can deal naturally with low particle numbers, as notedearlier, there is the potential for low particle numbers to break the PDE-compartment algorithm. Potentially,when particle numbers are low in the blending region, fractional particle numbers in a compartment could causea particle to be chosen to jump out of one compartment even though there is not sufficient mass for this to occur.The solution to this problem, as will be proposed in the discussion, is to introduce adaptive blending regimes,which ensure the PDE representation is only employed in regions of the domain where particle concentrations aresufficiently high to justify its use. x P a r t i c l e d e n s i t y (a) x P a r t i c l e d e n s i t y (b) Time −0.010.000.01 R e l a t i v e e rr o r PDE subdomainBlended subdomainCompartment subdomain (c) x P a r t i c l e d e n s i t y (d) x P a r t i c l e d e n s i t y (e) Time −0.010.000.01 R e l a t i v e e rr o r Compartment subdomainBlended subdomainParticle subdomain (f)
Figure 4.
Density and error plots for test problem 3 - morphogen gradient formation with auniform initial condition. Descriptions, including definitions of relative errors are as in figure 2except that panels (a) and (d) are density profiles evaluated at t = 0 .
01 rather than at t = 0 . h x × h y × h z . The blending region is itself a cuboidal region in which both thecoarse and fine models co-exist as equivalent representations of the mass in that region. Forthis translationally invariant example the blending functions are simply a function of x . Thismeans that only diffusion parallel to the x -axis is impacted in the blending region. Of course,for differently shaped domains and interfaces, the blending functions may be functions of allthree coordinates chosen to scale-diffusion as required providing f ( x, y, z ) + f ( x, y, z ) = D forall ( x, y, z ) ∈ B , the blending region, and both f ( I ) = f ( I ) = D and f ( I ) = f ( I ) = 0are satisfied, where I ∈ R and I ∈ R are surfaces specifying the interfaces which form theboundaries of the blending region.The mean-field PDE that corresponds to the reaction system (26) under the Poisson momentclosure assumption is given by ∂c∂t = D ∇ c − κ c + κ , (27)with corresponding zero-flux boundary conditions on each of the domain’s boundaries: ∂c∂n (cid:12)(cid:12)(cid:12)(cid:12) ∂ Ω = 0 . (28)For the simulations whose results are displayed in figure 6, we initialise the particles accordingto a linear gradient so that the initial density decreases in the positive x -direction. Explicitlyparticle density profiles are initialised according to the following density profile: c ( x, y, z ) = 183 − x , for ( x, y, z ) ∈ [0 , × [0 , × [0 , , (29)9giving N = 465 particles initially. The PDE part of the hybrid simulation can be initialisedexactly according to this profile. For the regions of the domain modelled by stochastic com-ponents of the hybrid method (e.g. in compartment-based regions or Brownian-based regions)the density profile is normalised and used as a probability density function (pdf) to assign po-sitions to the appropriate number of particles corresponding to that region of the domain. Inthe blending regions, particles are initialised according to the finer-scale simulation method andthe coarse scale density is matched appropriately. For example in the compartment-Brownianhybrid method we initialise, on average, one third of the particles with y and z coordinateschosen uniformly at random in [0 ,
1] and x -coordinates chosen from the pdf P ( x ) = ≤ x < / , − x for 10 / ≤ x < / , / ≤ x < . (30)Once the positions of the Brownian particles have been specified, the particles can then bebinned into compartments to determine the compartment-based initial condition in that region. (a) (b) Figure 5.
Schematic representations of the two-dimensional (a) PDE-compartment hybridand (b) compartment-Brownian hybrid. In panel (a) the green surface in the green regionrepresents the PDE solution in the purely PDE region of the domain. The red surface and thered columns represent equivalent PDE- and compartment-based representations of the mass inthe red blending region. The blue columns in the blue region of the domain represent thenumber of particles in each compartment in the purely compartment region of the domain. Inpanel (b) the blue boxes in the blue region of the domain represent the number of particles ineach compartment in the purely compartment region of the domain. The red boxes and thered circles represent equivalent compartment- and Brownian-based representations of the massin the red blending region. The yellow circles in the yellow region of the domain representindividual particles in the purely Brownian region of the domain.The hybrid method in three dimensions proceeds in an entirely analogous way to the one-0dimensional algorithms described above with full three-dimensional simulation of the PDE,compartment-based method and Brownian-based method. As before, in the blending region thetwo different descriptions are kept in sync with each other at every time step. Figure 5 providesschematic representations of the two coupling methods in two dimensions (in order to illustratehow the method generalises from one dimension).Parameter Value Description N (0) 465 Initial number of particlesΩ [0 , × [0 , × [0 ,
1] Spatial domain D κ . κ . ρ .
06 Particle interaction radius∆ t b − Brownian time-step P λ . × − Probability of reaction when inside the interaction radius K
20 Number of compartments h x / h y h z x /
300 PDE voxel width∆ t p − PDE time-stepM 500 Number of repeats
Table 3.
Table of parameters for the bimolecular reaction simulation (26) (test problem 4).The layout for figure 6 is the same as for figures 2-4. The only difference is the calculationof the relative mass error. For this final example, which includes second-order reactions, thesolution of the mean-field PDE model will no longer match the mean behaviour of either thecompartment-based model or the Brownian-based model. Consequently, in order to calculatethe relative mass error, we use the average behaviour of the finest-scale model in each hybridrepresentation (e.g. the compartment-based representation in the PDE-compartment modeland Brownian-based representation in the compartment-Brownian model) simulated across thewhole domain as the ground truth. For the PDE-compartment coupling the relative mass erroris, then, the difference between the average (over 500 repeats) number of particles in the givenregion in the hybrid method and the corresponding average (over 500 repeats) number in thesame region in the purely compartment-based simulation, divided by the number of particles in1the relevant region of the purely compartment-based simulation (to normalise):
RM E P ( t ) = R Ω P ¯ c ( x, y, z, t ) dx − P i,j,k ¯ F i,j,k ( t ) I c i,j,k ∈ Ω P P i,j,k ¯ F i,j,k ( t ) I c i,j,k ∈ Ω P , (31) RM E H ( t ) = P i,j,k ¯ C i,j,k ( t ) I c i,j,k ∈ Ω H − P i,j,k ¯ F i,j,k ( t ) I c i,j,k ∈ Ω H P i,j,k ¯ F i,j,k ( t ) I c i,j,k ∈ Ω H , (32) RM E C ( t ) = P i,j,k ¯ C i,j,k ( t ) I c i,j,k ∈ Ω C − P i,j,k ¯ F i,j,k ( t ) I c i,j,k ∈ Ω C P i,j,k ¯ F i,j,k ( t ) I c i,j,k ∈ Ω C , (33)where, as before, Ω P is the purely PDE region of the domain, Ω H is the blending region and Ω C is the purely compartment region of the domain. The averaged solution of the PDE componentof the hybrid method at position ( x, y, z ) at time t is denoted ¯ c ( x, y, z, t ), the averaged compart-ment particle numbers in compartment ( i, j, k ) of the hybrid method are denoted ¯ C i , j, k andthe averaged compartment particle numbers in compartment ( i, j, k ) of the fully compartment-based ‘ground truth’ simulation are denoted ¯ F i,j,k . The positions c i,j,k are the centres of thecompartments indexed ( i, j, k ).For the compartment-Brownian coupling the relative mass error is the difference between theaverage (over 500 repeats) number of particles in each region given by the hybrid method and theaverage number of particles in the same region in the purely Brownian-based simulation, dividedby the number of particles in the relevant region of the purely Brownian-based simulation (tonormalise): RM E C ( t ) = P i,j,k ¯ C i,j,k ( t ) I c i,j,k ∈ Ω H − ¯ E C ( t )¯ E C ( t ) , (34) RM E H ( t ) = P i,j,k ¯ C i,j,k ( t ) I c i,j,k ∈ Ω H − ¯ E H ( t )¯ E H ( t ) , (35) RM E B ( t ) = ¯ B ( t ) − ¯ E B ( t )¯ E B ( t ) , (36)where, as before, Ω B is the purely Brownian region of the domain and ¯ B ( t ) represents themean number of Brownian particles in the purely Brownian region of the hybrid method and¯ E C ( t ), ¯ E H ( t ) and ¯ E B ( t ) represent the mean number of Brownian particles in Ω C , Ω H and Ω B respectively at time t in the fully Brownian ‘ground truth’ simulations.There are some special points to note about the models which incorporate second-orderreactions. Firstly, as noted above, the solution of mean-field PDE, which we will employ in thePDE region of the PDE-compartment hybrid method, will not correspond to the mean behaviourof the compartment-based method. This is a result of the moment-closure approximation whichmust be used in order to derive a closed PDE for the mean behaviour. As a consequence,we might expect some disparity between the solution of the hybrid method and the solutionof the fully-compartment-based simulation that we take to be the ground truth in the PDE-compartment-based hybrid. Fortunately, for our compartment-Brownian hybrid method, [14]provide a method for matching reaction rates in compartment-based simulations to those inBrownian-based simulations, which we make use of.2We must also be careful to choose our parameters carefully in the compartment-Brownianmethod. If compartment-sizes are too small in the compartment-based method then particles canbecome too sparsely distributed and second-order reactions lost. [14] provide a way to alter thereaction rate (depending on the compartment size) to maintain the same overall reaction rate asa well mixed system. This correction, however, only holds down to a certain compartment size,beyond which second-order reactions are irrevocably lost. It is worth noting that [34] postulatedthe convergent reaction-diffusion master equation representation (in which particles can interactwith others in neighbouring boxes), which is consistent with the spatially continuous Doi modelfor reaction-diffusion even as box sizes become small. [31] numerically approximate mesoscopicreaction rates that are consistent with the popular Smoluchowski Brownian dynamics model upto a given lower limit on mesh size.In the Brownian-based method we need to ensure that the time-step is chosen to be suffi-ciently small that particles do not jump ‘too far’ between position updates. If particles jumplarge distances in each time-step then it is possible that particles which should have been giventhe opportunity to react with each other may not come into close enough proximity and somesecond-order reactions may be lost. Choosing the reaction radius of particles to be large mayhelp to mitigate this somewhat, but brings its own problems. The size of the interaction radiusis calculated by considering particles in free-space [14]. In reality, in our simulations, particlesare often close to boundaries. The proportion of the particle’s interaction radius that overlapsthe exterior of the domain is not able to interact with particles inside the domain, so the rate ofsecond-order reactions is again effectively reduced. Since, for a given reaction rate, the size ofthe interaction radius increases with the time-step, reducing the time-step is often sufficient tosolve both of these problems. We note that, with the exception of the PDE not matching themean behaviour of the compartment-based method, these issues are all inherent to the individ-ual modelling paradigms we have chosen to couple, and are not specific to the hybrid methodswe have developed. With sensible simulation parameter choices these issues can be overcome.The results of our simulations are plotted in figure 6. In Figures 6(d) and 6(e), whichcompare densities for the the compartment-Brownian hybrid paradigm we have good qualitativeagreement with the ground truth (the ubiquitously Brownian-based model). These qualitativeresults are further corroborated in figure 6(f) in which the low and unbiased relative mass errorover time are demonstrated.The density plots in figures 6(a) and 6(b) for the PDE-compartment hybrid coupling alsoappear to demonstrate good qualitative agreement. However, when considering the relativemass error in the different regions, in figure 6(c), we observe that, although low, the relativemass errors appear to be biased. This, as discussed above, should not be a surprise since themean-field PDE does not capture the mean behaviour of the compartment-based model, whichwe assume to be the ground truth for the relative mass error calculations. The overall massexpected in the fully compartment-based model at equilibrium would exceed that predicted bythe mean-field PDE. In agreement with this expectation we find that the total mass in all threeregions of the domain is less than it would be in the fully compartment-based simulations withthe problem being particularly acute in the PDE region. A simple comparison of the expecteddensities at time t = 1 shows that the maximum magnitude of the PDE relative error withrespect to the compartment based model is roughly 3 × − , demonstrating that the size of therelative error we find between our hybrid method and the solution of the fully compartment-3 x P a r t i c l e d e n s i t y (a) x P a r t i c l e d e n s i t y (b) Time −0.020.000.02 R e l a t i v e e rr o r PDE subdomainBlended subdomainCompartment subdomain (c) x P a r t i c l e d e n s i t y (d) x P a r t i c l e d e n s i t y (e) Time −0.020.000.02 R e l a t i v e e rr o r Compartment subdomainBlended subdomainParticle subdomain (f)
Figure 6.
Density and error plots for test problem 4 with an initial condition which exhibits aconstant gradient. Descriptions, excluding definitions of relative errors are as in figure 2. Allparameters are given in table 3.based simulations is of an appropriate order or magnitude, as it is similar to the difference in theconcentration when comparing the equilibrium profile of the full PDE to the fully compartmentbased method, adjusted for the specific voxel size.
When modelling multiscale phenomena it is often the case that concentrations vary spatiallyto such a degree that in one region of the domain a coarse, computationally inexpensive modelcan be tolerated, whereas, in another region of the domain a more accurate, but more expensiverepresentation is required.In this paper we have proposed a general hybrid blending mechanism which facilitates thespatial coupling of two reaction-diffusion modelling paradigms at different levels of detail inorder to accommodate the modelling of such multiscale phenomena. Our method employsa blending region and a corresponding blending function. The blending function scales upor down (respectively) the relative contribution to diffusion of a coarse or fine (respectively)representation of the reaction-diffusion process across the blending region such that diffusion ishandled to a different degree by each modelling representation.Specifically, we have developed an algorithm which couples a PDE representation of areaction-diffusion process to a compartment-based representation and, separately, an algorithmwhich couples a compartment-based representation at the coarse scale to a Brownian-based rep-resentation at the fine scale. Other algorithms exist to achieve such couplings [76, 19, 21, 24].Some of these algorithms are conceptually complex - relying variously on artificially introduced‘psuedo-compartments’, ‘ghost cells’ and ‘overlap regions’ - technically challenging to implement,4and strongly parameter dependent - working only in specific parameter regimes. We believe ourblending method provides a conceptually simple and easily implementable coupling methodol-ogy - requiring only an intuitively defined blending function to couple the two regimes together.This methodology might be readily employed to couple other modelling regimes (for examplePDE and Brownian modelling regimes) to form novel hybrid methods under a unified frameworkor implemented simply by non-experts for physical and biological applications.We have demonstrated, through four representative examples, that both of our coupling al-gorithms are able to handle a wide range of reaction-diffusion processes from simple diffusionthrough to reaction-diffusion processes incorporating first- and second-order reactions. The hy-brid methods are capable of representing these processes accurately (low error) and without bias(in the situation for which there is no discrepancy between mean-field behaviour of the coupledmodels) or with the expected bias (when such a discrepancy exists). Due to the computationalsavings afforded by coupling a cheap coarse model with an expensive fine-scale model, we canscale up particle numbers in our simulations in order to demonstrate that the hybrid algorithmsperform arbitrarily well in comparison to the full finest-scale model. For this reason we do notprovide explicit time comparisons of our methods, but rather focus on their accuracy.There are several directions in which we intend to extend this work, but which are notappropriate for inclusion in this initial proof-of-principle paper. Firstly, and perhaps moststraight-forwardly we would like to extend these hybrid methods to deal with more complexdomain geometries. Although we have demonstrated that our blending hybrid methods cancope with three dimensional reaction-diffusion processes, in real biological scenarios boundariesare likely to be curved and there is the potential for the requirement that interfaces betweencoarse and fine regimes are non-planar.Secondly, the dynamic nature of many biological processes mean that concentrations changesignificantly over time. If we are to ensure that the coarse modelling regime represents regions ofhigh concentration and the fine modelling regime regions of low concentration, then it is necessaryfor interfaces that border the blending region and the blending region itself to be dynamic. Themain challenge associated with dynamic interfaces is the conversion of one particle type intoanother. Fortunately, this challenge has been overcome previously by a number of differenthybrid methods, whose dynamic interface methodologies we might readily adapt to our hybridparadigm in a follow up work [56, 29, 65]. Related concerns are the need for the creation orremoval of multiple interfaces in scenarios in which particle concentrations oscillate in spaceand time. Similarly, reaction-diffusion simulations in which more than just a single speciesare interacting may require different interfaces for each of the different species. This raisespotentially difficult questions about how to carry-out reactions between species represented bydistinct modelling paradigms in the same region of space.A final direction in which we would like to extend this work is by considering entirely newhybridisation methods. For example, rather than having the two distinct modelling paradigmsrepresenting the same particles (as we have in the blending region) requiring both regimes tobe updated when one changes, it might be practicable to have the two modelling paradigms co-existing across the whole of the domain, but representing different proportions of the particlesdepending on the concentration. Such a method would remove the requirement for interfacesbetween the regions of the domain, effectively doing away with many of the concerns related todynamically and spatially changing concentrations raised earlier in this section.5Since biological and physical experiments can be carried out at increasingly high levels ofdetail, we are gaining more intricate and specific information about a wide variety of multiscaleprocesses. In order to test experimentally generated hypotheses about such processes we needto have modelling frameworks which are capable of replicating experimental behaviour to a highdegree of accuracy. The blending hybrid methods presented in this paper provide a straight-forward way to couple modelling paradigms with different levels of detail, which will facilitatemore accurate and more efficient multiscale modelling. Consequently, we expect that both ourown future work and the work of others, building on just such hybrid paradigms, will enablebiochemical simulations which go beyond what is tractable with current approaches.
AppendicesA Numerical simulation of the PDE
We now provide more details on the specifics of the macroscopic model that we employ through-out the paper, including an algorithm for its implementation. There exist a number of well-developed, efficient numerical methods for the solution of such reaction-diffusion PDEs [62,47, 17, 7]. Typically to implement these algorithms we discretise the PDE on a spatial mesh.This results in a system of ordinary differential equations (ODEs). These ODEs can then beintegrated forwards in time using standard numerical techniques.For the PDE (2) we start by dividing [ a, b ] into M voxels each of size ∆ x = ( b − a ) /M andwe define x j = ∆ x ( j − / x j is the centre of the voxel j (see figure 7). Typically thegrid spacing of the PDE solution method is very fine (much finer than the discretisation of spacein the compartment-based method) in order to mimic the true continuous-space PDE solutionas closely as possible. We discretise the PDE using the finite volume method over the grid infigure 7. For time integration we use the simple θ -method [47]. Figure 7.
Schematic illustrating a spatial discretisation of the one-dimensional domain [ a, b ],which is used to simulate equation (2) numerically. ∆ x represents the size of the voxels and x i for i = 0 , . . . , M − q j ( t ) = 1∆ x Z x j +1 / x j − / c ( t, x ) d x, (37)which corresponds to the average concentration per voxel, and D j = D ( x j ), where we note that j is not necessarily integer valued. By integration of PDE (2) over the finite volume voxels we6then obtain the semi-discrete approximation, d q dt = A q + b + R ( q ) , (38)where A = 1∆ x − D / D / D / − ( D / + D / ) D / . . . . . . . . . D M − / − ( D M − / + D M − / ) D M − / D M − / − D M − / , b = ∆ x − ( − J a , , · · · , , J b ) and q ( t ) = ( q ( t ) , q ( t ) , · · · , q TM − ( t )). We now solve the semi-discrete approximation using the θ -method . The complete method is described in Algorithm3. Algorithm 3:
An algorithm for the numerical solution of equation (2) using a first-orderfinite volume method
Input:
PDE mesh size – ∆ x ; time-step for the solution of the PDE – ∆ t p ; left and rightends of the domain – a, b ; initial concentration for the PDE – c init ; finalintegration time – T ; value of θ . Set t = 0, calculate ˆ c such that ˆ c j = 1∆ x Z x j +1 / x j − / c init ( x ) d x, and set q = ˆ c , n = 0. while t < T do q n +1 − q n ∆ t p = (1 − θ ) A q n + θA q n +1 + b + (1 − θ ) R ( q n ) + θR ( q n +1 ) Set t = t + ∆ t p , n = n + 1. end B Simulation of the compartment-based method
We now provide more details on the specifics of the mesoscopic model that we employ throughoutthe paper, including an algorithm for its implementation. More precisely, we have used an event-driven approach in order to simulate our compartment-based dynamics. The most commonly We employed a fully implicit method (i.e. θ = 1) for the test problems with linear reaction terms. For thetest problem with non-linear reaction terms we integrated the dynamics explicitly (i.e. θ = 0). a, b ] into K compartments, each of size h = ( b − a ) /K . Inorder to replicate the density dependent diffusion specified in the macroscopic model describedby equation (2) we require that the rates at which particles jump to the left and the right arenot equal in regions in which the diffusion coefficient is non-constant. Specifically, we mustevaluate the jump rates based on the diffusion coefficient at compartment boundaries [49]. Thisis visualised in figure 8 where we denote the jumping rate of a particle in compartment i intocompartment i + 1 with d + i , while we denote the jumping rate of a compartment i particle intocompartment i − d − i . Without loss of generality, assuming that the left-hand boundaryof the compartment-based regime is at a , the left jump rate for compartment i is given by d − i = D ( a + ( i − h ) h , for i = 2 , . . . , K. (39)and the right jump rate is given by d + i = D ( a + ih ) h , for i = 1 , . . . , K − . (40)Jump rates d − and d + K at the boundaries can be chosen in order to replicate the chosen boundaryconditions [67]. In the case of zero-flux boundary conditions (equivalent to setting J a = J b = 0in equations (3)), these jump rates are simply chosen to be d − = d + K = 0. Reaction propensityfunctions are specified to bring about the desired reaction rate . Once all the event rates havebeen specified then one can simulate the system using Gillespie’s direct method [26].We next provide a detailed implementation for the spatial Gillespie algorithm over a timeinterval of size ∆ t . This algorithm is designed to replace line 5 of Algorithm 1 and line 3 ofalgorithm 2. Without loss of generality assume the compartment-based region occupies [ a, I ], as It should be noted that the rate of second- and higher-order reactions depends, non-trivially, on the com-partment size, h , and that the desired rate of such higher-order reactions may not be implementable for someparticularly small compartment sizes [14]. Figure 8.
Schematic illustrating the jump rates, d ± i between compartments when simulatingthe mesoscopic reaction-diffusion paradigm on the domain [ a, b ] with compartment-size h .in the compartment-Brownian hybrid method. However, we note the caveat that for the PDE-compartment hybrid method the compartment-based region would occupy [ I , b ]. As alreadynoted, left and right jumping rates from compartment i are different and given in equations (39)9and (40), respectively. Algorithm 4:
Simulating the spatial Gillespie algorithm for a time-interval ∆ t . Input:
Compartment size – h ; time-step for the solution of the PDE or Brownian updatetime-step – ∆ t ; left and right ends of the blending region – I , I ; the totalnumber of compartments – K ; number of reactions – M ; initial particle numbers– C init ; compartment-based diffusion coefficient – D ( x ). Set t = 0. while t < ∆ t do Using equations (39) and (40), respectively, calculate the propensity functionscorresponding to the left jumps, α i = d − i , and right jump, α i + K = d + i , respectively,for i = 1 , . . . , K . Calculate also the propensity functions for each reaction, m , ineach compartment, α ( m +1) K + i for m = 1 , . . . , M and i = 1 , . . . , K . Calculate the sum of the propensity functions α = K X i =1 α i + α i + K + M X m =1 α ( m +1) K + i ! . (41) Draw the time, τ , until the next reaction from an exponential distribution withparameter α : τ = − α ln( u ) , (42)where u is drawn from a uniform distribution with support (0 , Update the time: t = t + τ . if t > ∆ t then Break end Choose the j th reaction to fire. Each reaction, j , (for j = 1 , . . . , (2 + M ) K ) is chosenwith probability α j /α (proportional to its propensity function). Implement the particle movement or reaction specified by reaction j by updating thecorresponding particle numbers. end Note that the Gillespie algorithm steps forwards in discrete time-steps. However, the time-steps themselves are drawn form a continuous distribution so that the solution time-points of theGillespie algorithm do not match up with those of the fixed time-step algorithms for PDE andBrownian-based simulation. Consequently, our technique to couple the two simulation method-ologies is to simulate the compartment-based dynamics until such a time as ∆ t is exceeded forthe first time. Since a PDE or Brownian update step is due at time ∆ t we do not implement thefinal Gillespie reaction whose time-step took us over the ∆ t time limit. Instead we implement aPDE or Brownian update step accordingly and correspondingly update the propensity functionsready to begin Algorithm 4 again.0 C Simulation of Brownian dynamics
In this appendix we provide a provide a detailed implementation algorithm for Brownian-baseddynamics, which is designed to replace line 5 in Algorithm 2.
Algorithm 5:
Simulating a Brownian update step of length ∆ t b . Input:
Brownian update time-step – ∆ t b ; left and right ends of the blending region – I , I ; number of reactions – M ; number of zeroth- first- and second-orderreactions – Z , F and S (respectively); probability of zeroth-, first- andsecond-order reactions in time-step ∆ t – P z , P f and P s (respectively); initialposition of particle i for i = 1 , . . . , N – y i ; Brownian-based diffusion coefficient – D ( x ). for i = 1 : N do Use equation (4) to update particle positions, implementing reflective boundaryconditions for any particles that cross the interface at I or the boundary at b . end if there are zeroth-order reactions then for z=1,. . . , Z do if u z < P z then Position a new particle of the appropriate type uniformly across the purelyBrownian domain. end end end if there are first-order reactions then for each appropriate particle do for f=1,. . . , F do if u f < P f then implement reaction f by removing reacting particle and/or placing productparticle(s) as appropriate. end end end end if there are second-order reactions then Calculate the distances between each pair of particles capable of reacting with eachother. for each appropriate reaction pair do for s = 1, . . . , S do if particle pair are within ρ s of each other then if u s < P s then implement reaction s by removing reacting particles and/or placingproduct particles as appropriate. In the above u z , u f and u s are drawn from a uniform distribution with support (0 , P z and P f , respec-tively, are calculated simply by multiplying the rate of reaction with the time-step, ∆ t b . Thecalculation of P s for second-order reaction, s , is somewhat more complicated and depends onthe choice of reaction radius, ρ s . For more details on this and the placement of new particlesafter reaction see [14]. Note that Brownian reactions are only implemented in the region [ I , b ],since outside this region reactions are implemented using the compartment-based regime.In theory, the fact that the diffusion coefficient of the Brownian-based particles falls to zeroat I should mean that particles cannot cross the interface there, rendering the implementationof reflecting boundary conditions at I in step 2 redundant. However, in practice, the finitetime-step we use to update the Brownian particles means that, with low probability, particlescan jump across the interface and must consequently be reflected back. Acknowledgements
Part of this work was conceived during the authors’ stay at the Newton Institute for the programStochastic Dynamical Systems in Biology: Numerical Methods and Applications. This work wassupported by EPSRC grant no EP/K032208/1. This work was partially supported by a grantfrom the Simons Foundation. CAS is supported by a scholarship from the EPSRC Centre forDoctoral Training in Statistical Applied Mathematics at Bath (SAMBa), under the projectEP/L015684/1. AG was partially supported by a summer research placement from the Institutefor Mathematical Innovation at the University of Bath.
References
1. F. Alexander, A. Garcia, and D. Tartakovsky. Algorithm refinement for stochastic partialdifferential equations: I. Linear diffusion.
J. Comput. Phys. , 182(1):47–66, 2002.2. F. Alexander, A. Garcia, and D. Tartakovsky. Algorithm refinement for stochastic partialdifferential equations: II. Correlated systems.
J. Comput. Phys. , 207(2):769–787, 2005.23. S. Andrews and D. Bray. Stochastic simulation of chemical reactions with spatial resolu-tion and single molecule detail.
Phys. Biol. , 1(3-4):137–151, 2004.4. R. Baker, C. Yates, and R. Erban. From microscopic to macroscopic descriptions of cellmigration on growing domains.
Bull. Math. Biol. , 72(3):719–762, 2010.5. D. Benson, P. Maini, and J. Sherratt. Analysis of pattern formation in reaction diffu-sion models with spatially inhomogenous diffusion coefficients.
Math. Comput. Model. ,17(12):29–34, 1993.6. A. Beskos and G. Roberts. Exact simulation of diffusions.
The Annals of Applied Proba-bility , 15(4):2422–2444, 2005.7. S. Brenner and C. Carstensen.
Finite element methods , chapter 1. Encyclopedia ofComputational Mechanics. John Wiley & Sons, Ltd, 2004.8. K.-H. Chiam, C. Tan, V. Bhargava, and G. Rajagopal. Hybrid simulations of stochasticreaction-diffusion processes for modeling intracellular signaling pathways.
Phys. Rev. E ,74(5):051910, 2006.9. U. Dobramysl, S. Rüdiger, and R. Erban. Particle-based multiscale modeling of intracel-lular calcium dynamics.
Multiscale. Model. Sim. , 14(3):997–1016, 2015.10. A. Duncan, R. Erban, and K. Zygalakis. Hybrid framework for the simulation of stochasticchemical kinetics.
J. Comput. Phys. , 326:398–419, 2016.11. J. Elf and M. Ehrenberg. Spontaneous separation of bi-stable biochemical systems intospatial domains of opposite phases.
Syst. Biol. , 1(2):230–236, 2004.12. S. Engblom, L. Ferm, A. Hellander, and P. Lötstedt. Simulation of stochastic reaction-diffusion processes on unstructured meshes.
SIAM J. Sci. Comput. , 31(3):1774–1797,2009.13. R. Erban and S. Chapman. Reactive boundary conditions for stochastic simulations ofreaction–diffusion processes.
Phys. Biol. , 4(1):16–28, 2007.14. R. Erban and S. Chapman. Stochastic modelling of reaction–diffusion processes: algo-rithms for bimolecular reactions.
Phys. Biol. , 6(4):1–18, 2009.15. R. Erban, S. Chapman, and P. Maini. A practical guide to stochastic simulations ofreaction-diffusion processes. arXiv preprint arXiv:0704.1908 , 2007.16. R. Erban, M. Flegg, and G. Papoian. Multiscale stochastic reaction–diffusion modeling:application to actin dynamics in filopodia.
Bull. Math. Biol. , 76(4):799–818, 2014.17. R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods.
Handbook of numericalanalysis , 7:713–1018, 2000.18. L. Ferm, A. Hellander, and P. Lötstedt. An adaptive algorithm for simulation of stochasticreaction-diffusion processes.
J. Comput. Phys. , 229(2):343–360, 2010.319. M. Flegg, S. Chapman, and R. Erban. The two-regime method for optimizing stochasticreaction–diffusion simulations.
J. Roy. Soc. Interface , 9(70):859–868, 2012.20. M. Flegg, S. Chapman, L. Zheng, and R. Erban. Analysis of the two-regime method onsquare meshes. (SIAM) J. Sci. Comput. , 36(3):B561–B588, 2014.21. M. Flegg, S. Hellander, and R. Erban. Convergence of methods for coupling of microscopicand mesoscopic reaction-diffusion simulations.
J. Comput. Phys. , 289(C):1–17, 2015.22. M. Flegg, S. Rüdiger, and R. Erban. Diffusive spatio-temporal noise in a first-passagetime model for intracellular calcium release.
J. Chem. Phys. , 138(15):154103, 2013.23. E. Flekkøy, J. Feder, and G. Wagner. Coupling particles and fields in a diffusive hybridmodel.
Phys. Rev. E , 64(6):066302, 2001.24. B. Franz, M. Flegg, S. Chapman, and R. Erban. Multiscale reaction-diffusion algorithms:PDE-assisted Brownian dynamics.
SIAM J. Appl. Math. , 73(3):1224–1247, 2013.25. T. Geyer, C. Gorba, and V. Helms. Interfacing brownian dynamics simulations.
J. Chem.Phys. , 120(10):4573–4580, 2004.26. D. Gillespie. Exact stochastic simulation of coupled chemical reactions.
J. Phys. Chem. ,81(25):2340–2361, 1977.27. D. Gillespie. Deterministic limit of stochastic chemical kinetics.
J. Phys. Chem. B ,113(6):1640–1644, 2009.28. C. Gorba, T. Geyer, and V. Helms. Brownian dynamics simulations of simplified cy-tochrome c molecules in the presence of a charged surface.
J. Chem. Phys. , 121(1):457–464, 2004.29. J. Harrison and C. Yates. A hybrid algorithm for coupling PDE and compartment-baseddynamics.
J. Roy. Soc. Interface , 13(122):20160335, 2016.30. A. Hellander, S. Hellander, and P. Lotstedt. Coupled mesoscopic and microscopic simu-lation of stochastic reaction-diffusion processes in mixed dimensions.
Multiscale. Model.Sim. , 10(2):585–611, 2012.31. S. Hellander, A. Hellander, and L. Petzold. Reaction rates for mesoscopic reaction-diffusion kinetics.
Phys. Rev. E , 91(2):023312, 2015.32. T. Hillen and K. Painter. A user’s guide to PDE models for chemotaxis.
J. Math. Biol. ,58(1):183–217, 2009.33. S. Isaacson. Relationship between the reaction–diffusion master equation and particletracking models.
J. Phys. A.-Math. Theor. , 41(6):065003, 2008.34. S. Isaacson. A convergent reaction-diffusion master equation.
J. Chem. Phys. ,139(5):054101, 2013.435. E. Keller and L. Segel. Initiation of slime mold aggregation viewed as an instability.
J.Theor. Biol. , 26(3):399–415, 1970.36. E. Keller and L. Segel. Model for chemotaxis.
J. Theor. Biol. , 30(2):225–234, 1971.37. E. Keller and L. Segel. Traveling bands of chemotactic bacteria: a theoretical analysis.
J. Theor. Biol. , 30(2):235–248, 1971.38. S. Khan, Y. Zou, A. Amjad, A. Gardezi, C. Smith, C. Winters, and T. Reese. Seques-tration of CaMKII in dendritic spines in silico . J. Comput. Neurosci. , 31(3):581–594,2011.39. J. King and R. O’Dea. Pushed and pulled fronts in a discrete reaction–diffusion equation.
J. Eng. Math. , 102(1):89–116, 2017.40. M. Klann, A. Ganguly, and H. Koeppl. Hybrid spatial gillespie and particle trackingsimulation.
Bioinformatics , 28(18):i549–i555, 2012.41. T. Kurtz. The relationship between stochastic and deterministic models for chemicalreactions.
J. Chem. Phys. , 57(7):2976–2978, 1972.42. J. Lipková, K. Zygalakis, S. Chapman, and R. Erban. Analysis of brownian dynamicssimulations of reversible bimolecular reactions.
SIAM J. Appl. Math. , 71(3):714–730, 2011.43. W.-C. Lo and S. Mao. A hybrid stochastic method with adaptive time step control forreaction–diffusion systems.
J. Comput. Phys. , 379:392–402, 2019.44. W.-C. Lo, L. Zheng, and Q. Nie. A hybrid continuous-discrete method for stochasticreaction–diffusion processes.
R. Soc. Open Sci. , 3(9):160485, 2016.45. E. Moro. Hybrid method for simulating front propagation in reaction-diffusion systems.
Phys. Rev. E , 69(6):060101, 2004.46. R. Mort, R. Ross, K. Hainey, O. Harrison, M. Keighren, G. Landini, R. Baker, K. Painter,I. Jackson, and C. Yates. Reconciling diverse mammalian pigmentation patterns with afundamental mathematical model.
Nat. Commun. , 7:10288, 2016.47. K. Morton and D. Mayers.
Numerical Solution of Partial Differential Equations . Cam-bridge University Press, 2005.48. H. Othmer, S. Dunbar, and W. Alt. Models of dispersal in biological systems.
J. Math.Biol. , 26(3):263–298, 1988.49. H. Othmer and A. Stevens. Aggregation, blowup, and collapse: the ABC’s of taxis inreinforced random walks.
SIAM J. Appl. Math. , 57(4):1044–1081, 1997.50. K. Painter and T. Hillen. Volume-filling and quorum-sensing in models for chemosensitivemovement.
Can. Appl. Math. Quart. , 10(4):501–543, 2002.551. K. Painter, P. Maini, and H. Othmer. Stripe formation in juvenile
Pomacanthus ex-plained by a generalized Turing mechanism with chemotaxis.
Proc. Natl. Acad. Sci. USA ,96(10):5549–5554, 1999.52. M. Plapp and A. Karma. Multiscale random-walk algorithm for simulating interfacialpattern formation.
Phys. Rev. Lett. , 84(8):1740, 2000.53. S. Redner.
A Guide to First-Passage Processes . Cambridge University Press, 2001.54. H. Risken.
The Fokker-Planck Equation. Methods of Solution and Applications . 1989.55. M. Robinson and M. Bruna. Particle-based and meshless methods with aboria.
SoftwareX ,6:172–178, 2017.56. M. Robinson, M. Flegg, and R. Erban. Adaptive two-regime method: application to frontpropagation.
J. Chem. Phys. , 140(12):124109, 2014.57. D. Rossinelli, B. Bayati, and P. Koumoutsakos. Accelerated stochastic and hybrid methodsfor spatial simulations of reaction–diffusion systems.
Chem. Phys. Lett. , 451(1):136–140,2008.58. T. Schulze, P. Smereka, and W. E. Coupling kinetic Monte-Carlo and continuum modelswith application to epitaxial growth.
J. Comput. Phys. , 189(1):197–211, 2003.59. J. Sherratt. An analysis of vegetation stripe formation in semi-arid landscapes.
J. Math.Biol. , 51(2):183–197, 2005.60. C. Smith and C. Yates. Spatially extended hybrid methods: a review.
J. Roy. Soc.Interface , 15(139), 2018.61. C. Smith and C. Yates. The auxiliary region method: A hybrid method for couplinga PDE to Brownian-based dynamics for reaction-diffusion systems.
R. Soc. Open Sci. ,5(8):180920, 2018.62. G. Smith.
Numerical solution of partial differential equations: finite difference methods .Oxford University Press, 1985.63. M. Smoluchowski. Versuch einer mathematischen Theorie der Koagulationskinetik kolloi-der Lösungen.
Z. Phys. Chem. , 92(129-168):9, 1917.64. T. Sokolowski, J. Paijmans, L. Bossen, T. Miedema, M. Wehrens, N. Becker, K. Kaizu,K. Takahashi, M. Dogterom, and P. Ten Wolde. eGFRD in all dimensions.
J. Chem.Phys. , 150(5):054108, 2019.65. F. Spill, P. Guerrero, T. Alarcon, P. Maini, and H. Byrne. Hybrid approaches for multiple-species stochastic reaction–diffusion models.
J. Comput. Phys. , 299:429–445, 2015.66. R. Strehl and S. Ilie. Hybrid stochastic simulation of reaction-diffusion systems with slowand fast dynamics.
J. Chem. Phys. , 143(23):234108, 2015.667. P. Taylor, R. Baker, and C. Yates. Deriving appropriate boundary conditions, and ac-celerating position-jump simulations, of diffusion using non-local jumping.
Phys. Biol. ,12(1):016006, 2015.68. A. Turing. The chemical basis of morphogenesis.
Phil. Trans. R. Soc. B. , 237(641):37–72,1952.69. N. van Kampen. Diffusion in inhomogeneous media.
J. Phys. Chem. Solids. , 49(6):673–677, 1988.70. N. van Kampen.
Stochastic processes in physics and chemistry . North-Holland, 2007.71. J. van Zon and P. ten Wolde. Green’s-function reaction dynamics: a particle-basedapproach for simulating biochemical networks in time and space.
J. Chem. Phys. ,123(23):234910, 2005.72. V. Volpert and S. Petrovskii. Reaction–diffusion waves in biology.
Phys. Life Rev. ,6(4):267–310, 2009.73. S. Winkelmann and C. Schütte. Hybrid models for chemical reaction networks: Multiscaletheory and application to gene regulatory systems.
J. Chem. Phys. , 147(11):114115, 2017.74. C. Yates and R. Baker. Importance of the Voronoi domain partition for position-jump reaction-diffusion processes on non-uniform rectilinear lattices.
Phys. Rev. E ,88(5):054701, 2013.75. C. Yates, R. Baker, R. Erban, and P. Maini. Going from microscopic to macroscopic onnon-uniform growing domains.
Phys. Rev. E , 86(2):021921, 2012.76. C. Yates and M. Flegg. The pseudo-compartment method for coupling partial differ-ential equation and compartment-based models of diffusion.