Coupling particle-based reaction-diffusion simulations with reservoirs mediated by reaction-diffusion PDEs
Margarita Kostré, Christof Schütte, Frank Noé, Mauricio J. del Razo
CCoupling particle-based reaction-diffusion simulations with reservoirs mediated byreaction-diffusion PDEs.
Margarita Kostr´e † , Christof Sch¨utte † , ‡ , Frank No´e ‡ ,a ) and Mauricio J. del Razo ‡ ,a ) Abstract:
Open biochemical systems of interacting molecules are ubiquitous in life-related pro-cesses. However, established computational methodologies, like molecular dynamics, are still mostlyconstrained to closed systems and timescales too small to be relevant for life processes. Alternatively,particle-based reaction-diffusion models are currently the most accurate and computationally feasi-ble approach at these scales. Their efficiency lies in modeling entire molecules as particles that candiffuse and interact with each other. In this work, we develop modeling and numerical schemes forparticle-based reaction-diffusion in an open setting, where the reservoirs are mediated by reaction-diffusion PDEs. We derive two important theoretical results. The first one is the mean-field for opensystems of diffusing particles; the second one is the mean-field for a particle-based reaction-diffusionsystem with second-order reactions. We employ these two results to develop a numerical schemethat consistently couples particle-based reaction-diffusion processes with reaction-diffusion PDEs.This allows modeling open biochemical systems in contact with reservoirs that are time-dependentand spatially inhomogeneous, as in many relevant real-world applications. † Zuse Institute Berlin, Germany. ‡ Freie Universit¨at Berlin, Department of Mathematicsand Computer Science, Germany a ) Corresponding authors. E-mails: [email protected] [email protected]
I. INTRODUCTION
Complex systems of interacting particles/agents thatexchange energy and matter with a large reservoir areextremely common, from a city exchanging infected cit-izens at airports during a pandemic to a living cell ex-changing chemicals with its environment. In the contextof molecular biology, most biochemical reaction systems,either inside living cells or composed by them, interactwith some form of reservoir [1]. They also consume chem-ical energy for their survival, produce waste and dissipateheat; they operate in an open non-equilibrium setting. Interms of physical chemistry, every living system must bean open system —a closed system has no life [1]. It isthus fundamental to develop models of biochemical reac-tion systems capable of exchanging materials and energywith their environment. Although this is the guiding mo-tivation for this work, the results here presented can alsobe applied in other areas such as agent-based modeling.Molecular dynamics (MD) is theoretically capable ofmodeling biochemical systems accurately at cellular andsub-cellular scales. However, in timescales relevant tolife processes, even a system with one or two molecules isalready large enough to render any MD simulation com-putationally unfeasible. Moreover, although there is on-going research on open MD systems [2–4], establishedcomputational protocols rely by design on the simulatedMD system being closed [3]. Consequently, the most ac-curate and computationally reliable simulations for opensystems at these scales are based on stochastic particle-based reaction-diffusion (PBRD) models. Novel methodsare emerging capable of coupling PBRD with MD, inte- grating the accuracy of MD into efficient PBRD simula-tions [5–7]. This, along with other recent research, pointsout the need for hybrid reaction-diffusion methods thatare accurate, open and consistent across multiple scales.In this paper, we concentrate on multiscale modelsfor PBRD simulations of open systems, where we focuson coupling particle-based simulations to macroscopicparticle/chemical reservoirs. These reservoirs are givenby a mean concentration of chemical species that canvary in time and space. We model these reservoirs asreaction-diffusion partial differential equations (PDEs).The goal of this work is to achieve a mathematically con-sistent coupling between the PBRD simulations and thereaction-diffusion PDEs. This entails two major chal-lenges: • Determine how to consistently couple a particle-based simulation to a constant concentration chem-ical reservoir. Note the particle-based simulationwill be in the grand canonical ensemble since thenumber of particles is not constant. This will besolved by calculating the mean-field for an opensystem of diffusing particles. • Calculate the relations between the reaction ratesand diffusion coefficients in the particle-based simu-lation and the macroscopic reaction-diffusion PDE.This can be especially cumbersome for second-orderreactions. We will solve this by calculating themean-field for reaction-diffusion systems with upto second-order reactions.Unlike the case of homogeneous reaction theory [8–10], the connection between microscopic, mesoscopic andmacroscopic scales for reaction-diffusion phenomena isstill a matter of recent research [11–15]. One of the maindifficulties to establish this connection is to relate themacro and microscopic rates for second-order reactions.We will present results on the relation between these ratesat first order. The results for the most general case arecurrently a work in progress [16]. a r X i v : . [ q - b i o . Q M ] M a y Other very relevant work in this area is [17], where theauthors introduce a method to couple Brownian dynam-ics simulations with mean-field PDEs. However, theydo not implement their methods for the nontrivial caseof second-order reactions, and they conceptually do notconsider the PDE domain as a reservoir. The work [18]builds and improves on these ideas. The authors includea second-order reaction example. However, the relationbetween the particle-based and PDE reaction rate doesnot seem correct. It is based on [19], where the diffusionis included in the rate relation, which is only correct ifdiffusion is averaged out. In this work, we show that thediffusion coefficient does not play a role in the relationbetween the microscopic and the macroscopic reaction-diffusion PDE reaction rates.The theory and numerical implementation details areexplained in this paper; the code is available in GitHubunder an MIT license [20].
II. REACTION-DIFFUSION MODELS
Reaction-diffusion processes take different forms de-pending on the number of particles involved and thescale of interest (Fig. 1). When the number of par-ticles is small, the stochastic fluctuations due to diffu-sion and chancy reactions need to be taken into accountwith a probabilistic particle-based approach. However, ifthe number of particles is large, small fluctuations in thenumber of particles become negligible, and deterministicconcentration dynamics in the form of reaction-diffusionPDEs become a more suitable alternative.
Particle-basedreaction di ff usionReaction Di ff usionPDEChemical MasterEquationReaction rate equations (law of mass action) spatial resolutionlow (well-mixed) high nu m b e r o f p a r ti c l e s Reaction di ff usionmaster equationCompartmentODEs l a r g e s m a ll Figure 1. Some models of reaction-diffusion processes orga-nized by their spatial scaling and the number of particles.Only the most relevant models for this work are shown here.
In the case of well-mixed systems, the spatial com-ponent of these models is averaged out yielding eitherthe well-known rate equations or the chemical masterequation, depending on the number of particles (Fig.1). In this work, we concentrate on spatially inhomoge-neous systems at different scales. We develop theoreticaland simulations techniques that bridge PBRD simula-tions with reaction-diffusion PDEs, in particular in the context of open systems. We give below a brief overviewof the two relevant reaction-diffusion models.
A. Particle-based reaction-diffusion
The particle-based approach to model reaction-diffusion follows the approach used in the ReaDDy2 soft-ware [21] (other well-known software packages are [22–24]). It consists on simulating each molecule as an spher-ical particle. The Brownian diffusion of each molecule ismodeled using overdamped Langevin dynamics, dx ( t ) = √ Ddw ( t ) , (1)where x ( t ) is the position of the molecule at time t and w ( t ) is a collection of independent Wiener processes (oneper coordinate). Reactions are modeled depending onthe type of reaction (Fig. 2): • Zeroth order reactions, ∅ κ − (cid:42) A : a new A moleculeis placed uniformly in the whole domain with rate κ . • First order reactions, A κ ab −− (cid:42) B : molecule A istransformed into B with rate κ ab . • Second order reactions, A + B κ − (cid:42) C : if the relativedistance between A and B is less than the reac-tion radius σ , the particles react with rate α . Thenew molecule C is placed in an averaged positionbetween A and B . Figure 2. Illustration of some of the possible reactions ina particle-based reaction-diffusion simulation. a. Diagramof the Doi model for bimolecular reactions. When the twoparticles are closer than a distance of σ , they react with rate α . Note it is not the same as the macroscopic rate κ . b. An example of a unimolecular reaction (first-order) reaction,where A simply transforms into B . c. The backward reactionof the binding given by the Doi model. The products shouldbe placed uniformly at a distance δr such that δr ≤ σ tosatisfy detailed balance [25]. The mathematical model for zeroth and first-order reac-tions is the same as in the well-mixed case, where thereis no spatial dependence. Higher-order reactions can bedecomposed into several second-order reactions, so thetheory for second-order reactions results the most rele-vant. The mathematical theory we use to model second-order reactions is based on diffusion-influenced reactions,specifically in the Doi volume reactivity model (Fig. 2a)[26, 27].The Doi model consists of an isolated pair of molecules A and B , where A is fixed at the origin, B is a distance r from A , and it undergoes Brownian diffusion (Eq. (1))with diffusion coefficient D . Further, B can only reactwith A with rate α if within the reaction radius r ≤ σ . The probability of finding B at distance r at time t ,providing it started at a distance r , is f ( r, t | r ), whichobeys the following Fokker-Planck equation ∂ t f ( r, t | r ) = D ∇ f ( r, t | r ) − χ r ≤ σ ( r ) f ( r, t | r ) , where χ r ≤ σ ( r ) is the indicator function (1 if r < σ and 0otherwise). Note as the steady state is reached ( t → ∞ ), f ( r, t | r ) goes to zero since the particle B will reactwith probability 1. Following generalizations of diffusion-influenced reactions to reversible reactions [25, 28–35],we can also model the reversible reaction C − (cid:42) A + B ,by placing the reactants uniformly within a distance σ of each other. This choice ensures detailed balance issatisfied. Alternative solutions exists when there is aninteraction potential involved [25] or when numerical ef-ficiency is a priority [36].Note the parameters required for a particle-basedreaction-diffusion simulation are the microscopic ratesand the reaction radius for second-order reactions. B. Reaction diffusion PDEs
When a systems has a very large number of particles,its chemical kinetics are better described by deterministicconcentration-based approach. Consider c ( t ) the vectorof concentrations of N chemical species, which are in-volved in M different reactions. The kinetics are givenin terms of the following PDE, ∂ t c = D ∇ c + R ( c ) , where here D is a diagonal matrix with the diffusion coef-ficient of each species and R ( c ) encapsulates all the M re-actions. This equation without the diffusion term wouldbe of the form of the well-known law of mass action [37].As an example, consider the predator-prey dynamics A κ − (cid:42) A, A + B κ − (cid:42) B, B κ − (cid:42) ∅ , (2)where A represent the preys and B the predators. Theconcentration dynamics are given by the Lotka-Volterraequations ∂ t c A = D A ∇ c A + κ c A − κ c A c B ,∂ t c B = D B ∇ c B + κ c A c B − κ c B , where c A and c B represent the concentration of predatorsand preys, respectively. Naturally, initial and boundaryconditions need to be provided to close the system.Reaction-diffusion PDEs, like the one just presented,can be solved using standard numerical methods, likefinite difference schemes, finite elements and spectralmethods. In this work, we use finite differences [38–40]since they are simple and fit the purpose of this work. Wemainly use the Crank-Nicolson method combined with anoperator splitting approach if necessary, see Appendix Afor brief implementation details. III. COUPLING PARTICLE-BASED MODELSTO CHEMICAL RESERVOIRS
In this section, we derive two results that are funda-mental to couple PBRD simulations with reservoirs me-diated by reaction-diffusion PDEs. These results addressthe two main challenges mentioned in the introduction,and they are summarized in Fig. 4.The first result determines how to consistently cou-ple a particle-based model to a constant concentrationreservoir. It does so by matching the mean-field limitof a particle-based diffusion process to its correspond-ing macroscopic PDE description. In this setting, theparticle-based model is in contact with a chemical reser-voir, so we call it a grand canonical diffusion process. Ina simulation context, this result can be easily extendedto reservoirs with spatially and time-dependent concen-trations given by a reaction-diffusion PDE.The second result shows how reaction-diffusion PDEscan be recovered as the mean-field of particle-basedreaction-diffusion processes. This is essential to developconsistent coupling numerical schemes since it determinesthe relation between the microscopic and macroscopic pa-rameters.
A. Mean-field of grand canonical diffusionprocesses
We want to consistently couple PBRD simulations toPDE mediated chemical reservoirs. To parametrize thecoupling, we need to match the mean-field dynamics ofparticle-based models with its corresponding macroscopicPDE behavior.We begin with a one-dimensional system with an arbi-trary number of noninteracting diffusing particles and acoupling to a constant concentration reservoir on one end;this is a grand canonical diffusion process. In the macro-scopic setting, this corresponds to the diffusion PDE ∂ t c ( x, t ) = D ∇ c ( x, t ) , (3) ∂ x c ( x, t ) | x =0 = 0 c ( R, t ) = c R , with c R a constant. The boundary condition at x = 0is not really relevant for our purpose, but we assumeNeumann for simplicity.In the particle-based setting, this corresponds to parti-cles diffusing independently following standard Brownianmotion. The reservoir on one end can absorb and intro-duce new particles into the system with a certain rate.To obtain the mean-field of the particle-based system, itwill be convenient to discretize the domain [0 , R ] into N cells of size δx . The number of particles in the i th cellis denoted by n i . The state of the systems is given by P ( n , . . . , n N , t ), which corresponds to the probability ofhaving { n , . . . , n N } particles in the cells { . . . , N } attime t . We further assume cell N is in contact with areservoir in cell N + 1 of volume V R with a constant con-centration of particles c R at all times (Fig. 3b.). As theparticles diffuse independently, the jump rates of eachparticle from cell i to neighboring cell j are simply thediffusion jump rates [41], q i,j = Dδx , i (cid:54) = j. (4)We can now write a master equation for the dynamics of P ( n , . . . , n N , t ) [41], dP ( n , . . . , n N , t ) dt = − P ( n , . . . , n N , t ) N (cid:88) i =1 [ q i,i +1 + q i,i − ] n i + N − (cid:88) i =1 (cid:20) P ( n i + 1 , n i +1 − q i,i +1 ( n i + 1)+ (5) P ( n i − , n i +1 + 1) q i +1 ,i ( n i +1 + 1) (cid:21) + P ( . . . , n N + 1) q N,N +1 ( n N + 1) + P ( . . . , n N − q N +1 ,N n R , where we used P ( n i + 1 , n i +1 −
1) = P ( . . . , n i + 1 , n i +1 − . . . ) to simplify notation. We refer to this equation asthe grand canonical master equation since it describes anopen system that can exchange particles with its environ-ment. The terms in the first sum of Eq. (5) correspond totransitions that leave the current state; the second sumcorresponds to transitions into the current state; and thelast two terms corresponds to the transitions into the cur-rent state due to interactions with the reservoir. In orderto recover the Neumann boundary condition at x = 0,we set q , = 0, see Fig. 3 and Eq. (5) for reference.Note the jump rate of particles from the reservoir intothe system γ = q N +1 ,N is not yet known. We will referto γ as the injection rate.We will next study the continuous limit of the mean-field of Eq. (5) and obtain the value of the injection rate γ . The mean-field of the Eq. (5) is given by (cid:88) { ¯ n } n i dP ( n , ..., n N , t ) dt := d (cid:104) n i (cid:105) dt , (6)where (cid:104) n i (cid:105) is the expected number of particles at cell i . After some algebra, we obtain the following equation[41, 42] d (cid:104) n i (cid:105) dt = (cid:104) n i +1 (cid:105) q i +1 ,i − (cid:104) n i (cid:105) [ q i,i +1 + q i,i − ] + (cid:104) n i − (cid:105) q i − ,i . (7) Figure 3. Diagram of the grand canonical master equation foran open system. It allows for an arbitrary number of diffusingparticles, and it is coupled to a material reservoir. This isthe result of discretizing the one-dimensional particle-baseddiffusion processes in contact with a constant concentrationreservoir.
Writing everything in terms of concentrations c i = (cid:104) n i (cid:105) /δx and substituting the rates, we can take the limit δx →
0. Note that the concentration remains boundedbecause (cid:104) n i (cid:105) goes to zero at the same rate as the volumeshrinks. This yields [41] ∂ t c ( x, t ) = D ∇ c ( x, t ) , (8)which is not surprisingly the diffusion equation for theconcentration c ( x, t ). It is also straightforward to check,we recover the Neumann boundary condition at x = 0by using a ghost cell approach [38]. However we needto be careful at the boundary in contact with the reser-voir. We denote F i := (cid:104) n i (cid:105) /δx the mean concentrationat cell i . We write Eq. 7 for i = N and substitute thecorresponding rates dF N dt = q N +1 ,N c R + Dδx ( − F N + F N − ) . We add a ghost cell N + 1 that represents the reservoir.As the concentration in the material bath is constant, weset c R = F N +1 and rewrite this equation as dF i dt = Dδx ( F N +1 − F N + F N − )+ q N +1 ,N c R − Dδx F N +1 . The first term corresponds to the discretized diffusionequation (Eq. 3) that we want to recover in the continu-ous limit. Therefore, the additional terms must be zero.As c R = F N +1 , this implies that the injection rate is γ = q N +1 ,N = Dδx . (9)This is the jump rate of particles from the reservoir intothe system, such that the macroscopic Eq. (3) withits boundary conditions are recovered in the continuousmean field limit. Not surprisingly, it matches the diffu-sion jump rate.This result establishes the connection between theparticle-based and the concentration-based approach foropen systems. It can be used to implement particle-based a.b. Figure 4. Diagram summarizing the two main results of Sec-tion III. a. Mean field limit of the particle-based diffusionopen system. If the injection rate of particles from the reser-voir is set to γ = D/δx in the master equation, the meanfield yields a constant concentration boundary condition. b. Mean field limit of PBRD. It relates the miscroscopic parame-ters of the Doi model, α and σ , with the macroscopic reactionrate κ . simulations in contact with material reservoirs with con-stant concentrations, see Section IV. It is also straight-forward to extend it to more complicated systems, suchas time and space-dependent reservoirs.Note we assumed that the concentration c R remainsconstant, even when extracting particles from the sys-tem. This is only possible if we make the number of par-ticles in the reservoir and its volume both infinite whilekeeping the concentration constant. We implicitly makethis assumption when taking the continuous limit.To implement a simulation of this process, we can dis-cretize time at first order, so the jump rates of each parti-cle (Eq. (4)) become jump probabilities p i,j = ∆ tD/δx , i (cid:54) = j [32]. We can then use these, along with the jumpprobabilities of reservoir particles, to implement the sim-ulation. However, as the probability of one particle jump-ing either left or right is at most one, then 2∆ tD/δx ≤ t to,∆ t ≤ δx / (2 D ) . (10)This will be important for the implemenation of the Hy-brid scheme in Section IV.The results of this section are based on [41], and theycan also be obtained using the law of large numbers. Weadvise the reader to consult [41] for additional details. B. Mean-field of particle-based reaction-diffusionprocesses
In this section, we obtain the mean-field dynamics ofparticle-based reaction-diffusion processes. This resultprovides a mathematical bridge between the microscopicparticle-based approach and the macroscopic determin-istic approach for reaction-diffusion processes. It estab-lishes the relation between the microscopic and macro-scopic parameters. The diffusion coefficients and the re-action rates for up to first order reactions remain the same for particle-based models and macroscopic models[41]. However, this is not the case for second-order reac-tions (bimolecular reactions). We thus focus on derivingthis result for bimolecular reactions, A + B → C. (11)We denote α as the microscopic reaction rate based onthe Doi model (Section II A). We want to determine therelation between the particle-based simulations of this re-action with the corresponding macroscopic deterministicreaction-diffusion PDE ∂ t c A = D A ∇ c A − κc A c B , (12) ∂ t c B = D B ∇ c B − κc A c B , where c A and c B denote the deterministic concentrationsof A and B , and κ is the macroscopic reaction rate. Thiswill naturally provide a connection between α and κ . Weobtain this result for two dimensions, but it extends nat-urally to higher dimensions.We begin with the particle-based description. Onceagain, it is convenient to discretize the two dimensionaldomain, Ω, in n square cells denoted by V , V , ...., V n of length and height h . Denote by X Ai ( t ) and X Bi ( t ) thenumber of A and B particles in cell i at time t . Followingthe Doi model for particle-based reactions [21, 26], weassume A, B react with rate α if they are are closer than adistance σ . Using the reaction diffusion master equationformalism [15], the change in number of A particles incell i can be written using the Kurtz representation [8] X Ai ( t ) = X Ai (0) − n (cid:88) j =1 U j (cid:18)(cid:90) t αφ ij X Ai ( s ) X Bj ( s ) ds (cid:19) + D ti ( X A ) , (13)where X A = { X A , . . . X AN } and the term D ti denotes thediscrete change of number of particles at cell i due todiffusion after a time t . This operator can be explicitlywritten in terms of the discrete Laplace operator andthe elapsed time t . Each U j denotes a unit rate Poissonprocess, where the rate function (cid:82) t λ ( s ) ds depends onthe propensity λ ( s ) = αφ ij X Ai ( s ) X Bj ( s ) of the reaction.We can analogously formulate an equation for X Bi anddo the same process as below.Note the reaction propensity not only depends on themicroscopic rate α but also on the quantity φ ij . Thiswill help us quantify how likely the reaction is to happen,given that the A particles are in cell V i and B particlesare in cell V j . This is defined as φ ij = |R ∩ V ij || V ij | , (14)where R is the reactive region in the 4-dimensional spacedefined by positions x and y such that | x − y | ≤ σ , V ij isthe hypercube formed by V i × V j with volume h and thebars denote we are taking the volume of these regions.This quantity is simply the ratio between the reactivevolume contained within V ij and the total volume of V ij .If A particles are well-mixed in V i and B particles arewell-mixed in V j , it gives the probability of A and B particles being close enough to react, see [15] for moredetails on this quantity.Rewriting Eq. (13) in terms of concentrations, we ob-tain C Ai ( t ) = C Ai (0) − | V i | n (cid:88) j =1 U j (cid:18)(cid:90) t αφ ij | V ij | C Ai ( s ) C Bj ( s ) ds (cid:19) + D ti ( C A ) , (15)where we used the linearity of the diffusion operator.Note the rate of the unit Poisson process depends on theconcentrations, which are random variables themselves.These are therefore doubly-stochastic processes, and theyare called mixed Poisson processes or more generally Coxprocesses [43, 44]. The expected value of a mixed Poissonprocess corresponds to the expected value of the randomrate function [43]. Therefore the expected value of eachof the Poisson processes U j from Eq. (15) is (cid:90) t αφ ij | V ij |(cid:104) C Ai ( s ) C Bj ( s ) (cid:105) ds. In general, (cid:104) C Ai C Bj (cid:105) = (cid:104) C Ai (cid:105)(cid:104) C Bj (cid:105) + Cov( C Ai , C Bj ) sincethey are not independent. Assuming the number of par-ticles is large enough, as the particles are only correlatedthrough reactions limited to a small reaction volume, thecovariance is negligible in comparison to the product ofthe means. We can use these results to calculate theexpected value of equation (15) (cid:104) C Ai ( t ) (cid:105) = (cid:104) C Ai (0) (cid:105) − n (cid:88) j =1 | V j | (cid:90) t αφ ij (cid:104) C Ai ( s ) (cid:105)(cid:104) C Bj ( s ) (cid:105) ds + D ti ( (cid:104) C A (cid:105) ) , (16)where we used the linearity of the discrete diffusion op-erator to pass the expectation into the argument. Wefurther derive with respect to t both sides of the equa-tion, yielding d (cid:104) C Ai (cid:105) dt = D i ( (cid:104) C A (cid:105) ) − n (cid:88) j =1 αφ ij | V j |(cid:104) C Ai (cid:105)(cid:104) C Bj (cid:105) , (17)where D i is now the time-continuous diffusion operatorin terms of rates instead of jump probabilities.Equation (17) is still discrete in space; we now takethe continuous limit h →
0. In this limit, the first termconverges to the continuous diffusion operator D , see [41,45] for details. The limiting behavior of the second termis not trivial, since we need to solvelim h → n (cid:88) j =1 αφ ij | V j |(cid:104) C Ai (cid:105)(cid:104) C Bj (cid:105) , (18) where φ ij , as defined in Eq. (14), can be expressed asthe following integral φ ij = 1 | V ij | (cid:90) V i (cid:90) V j χ | x − y | <σ dydx, (19)where χ is the indicator function. If our domain is atwo-dimensional cube, then x and y , each correspond to2-dimensional vectors, and we have an integral over ahypercube. Substituting Eq. (19) into Eq. (18), weobtainlim h → αh (cid:104) C Ai (cid:105) n (cid:88) j =1 (cid:104) C Bj (cid:105) (cid:90) V i (cid:90) V j χ | x − y | <σ dydx , where we used that | V i | = h =, | V ij | = h . Rearrangingthe terms, we obtain= lim h → αh (cid:104) C Ai (cid:105) (cid:90) V i n (cid:88) j =1 (cid:104) C Bj (cid:105) (cid:90) V j χ | x − y | <σ dydx , As σ is small, we can assume the value (cid:104) C Bj (cid:105) does notchange much in the domain of interest, | x − y | < σ ,around cell i . We thus approximate it by its centralvalue (cid:104) C Bi (cid:105) . Using this approximation, we can combinethe sum over j and the integral over V j into one integralover the whole domain Ω, yielding= lim h → (cid:18) αh (cid:104) C Ai (cid:105)(cid:104) C Bi (cid:105) (cid:90) V i (cid:90) Ω χ | x − y | <σ dydx (cid:19) . We further approximate the integral of the indicatorfunction over Ω by an integral on the whole 2-dimensionalspace. This approximation is exact everywhere except ina small region close to the boundaries of Ω. Thus, ityields the area of a circle of radius σ .= lim h → (cid:18) αh (cid:104) C Ai (cid:105)(cid:104) C Bi (cid:105) πσ (cid:90) V i dx (cid:19) , The remaining integral is simply h , so= lim h → (cid:0) απσ (cid:104) C Ai (cid:105)(cid:104) C Bi (cid:105) (cid:1) = απσ c A ( x i ) c B ( x i ) , where c A = (cid:104) C A (cid:105) and c B = (cid:104) C B (cid:105) are continuous func-tions in space and time. We need to be careful with thislimit since the indexing can change as h → . Here weassumed the mean concentrations (cid:104) C Ai (cid:105) and (cid:104) C Bi (cid:105) are al-ways centered in a fixed value x i as h →
0. Note thisresult is only an approximation, albeit a very accurateone for σ (cid:28) Ω. With this result, the continuous spacelimit of Eq. (17) is then the familiar reaction diffusionPDE ∂c A ∂t = D ( c A ) − απσ c A c B . Note the same equation can be obtained for the kinet-ics of c B . The operator D corresponds to the well-knowndiffusion operator and comparing to Eq. (12), the macro-scopic rate κ is simply κ = απσ . (20)This relation holds for two dimensions, and it is of theform κ = αV react , where V react is the volume of the re-active volume. In three dimensions, we can analogouslyshow this relation has the same form with the correspond-ing reactive volume κ = α πσ / . The result we just derived shows the reaction-diffusionPDE is the mean-field of the particle-based simulationbased on the Doi model, assuming the number of par-ticles is sufficiently large (otherwise covariances need tobe taken into account). It also provides a connection be-tween the microscopic bimolecular reaction rate α andthe reaction radius σ with the macroscopic reaction rate κ . As a side note, it is interesting that the same result isobtained when σ (cid:28) h . Also note the diffusion coefficientdoes not play a role in this relation.We completely neglected the covariances when takingthe expectation of the mixed Poisson process from Eq.(15). Taking into account the covariances would yield amodified mean-field behavior valid at mesoscopic scales.The authors are currently working on a formal and moregeneral version of this result that takes covariances intoaccount [16]. IV. HYBRID SCHEME
In this section, we use the results from Section IIIto derive a hybrid scheme to couple particle-based sim-ulations with reservoirs mediated by reaction-diffusionPDEs.The essence of the algorithm is illustrated in Fig. 5.The domain is split into the particle domain and the con-centration domain. The dynamics on the concentrationdomain, which functions as the reservoir, are modeledby a reaction-diffusion PDE, which is solved either ana-lytically or with standard finite difference methods. Asthis is well documented in the literature [38, 40] (see Ap-pendix A), we concentrate on the particle domain wherethe dynamics are governed by three processes: injection,reaction and diffusion. We describe next these processesand show how they are combined into one algorithm.
A. Injection
We use the result from Section III A, mainly Eq. (9),to obtain a consistent scheme to inject particles from thereservoir into the particle system. Along the edge be-tween the particle and the concentration domain, we cre-ate squared boundary cells of edge length δx (Fig. 5a).Every time iteration and for every boundary cell i in the concentration domain, we calculate the average concen-tration c i within the boundary cell and convert it intonumber of particles N i (multiplying by cell volume V i ).The resulting values will be in general non-integers. Theinteger part corresponds to the number of virtual par-ticles, and the fractional part to one fractional virtualparticle. We call them virtual particles because they donot belong to the particle domain. Given this informa-tion, the injection of particles for a time-step ∆ t followsthe following procedure. Injection procedure
When called for a time interval τ , this procedure follows these steps:1. Let each virtual particle jump into the correspond-ing neighboring boundary cell in the particle do-main with the injection rate γ (Eq. (9)). This cor-responds to jumping with probability 1 − exp( − γτ ).2. Let the fractional virtual particles jump in the sameway with a rate γ scaled by the corresponding frac-tion value. For example, if N i = 15 .
87, then therate for the fractional virtual particle is 0 . γ .3. For every successful jump event, place a new par-ticle uniformly in the corresponding boundary cellof the particle domain. B. Reaction
Most reactions are or can be decomposed into uni-molecular or bimolecular reactions, so we only focus onthese types of reactions. We employ the methodology in-troduced in Section II A to implement the particle-basedsimulation in the particle domain. We further use theresult from Section III B to establish a relation betweenthe parameters of the PBRD simulation and the param-eters of the reaction-diffusion PDE. The diffusion coeffi-cients and unimolecular reaction rates remain the sameon both models, while the bimolecular reaction rate fromthe PDE can be obtained from PBRD parameters follow-ing Eq. (20).
Zeroth and first-order reactions
These reactionsdepend uniquely on zero or one particle, e.g. creation ofparticles and conformational changes, respectively. Thediffusion coefficients of zeroth and first-order reactionrates are the same in the PBRD simulation as in thePDE. Given a reaction rate k , on the particle-based sim-ulation, the probability of a reaction to happen within atime interval ∆ t is given by1 − exp( − k ∆ t ) . (21)The boundary coupling remains the same as in the injec-tion process, as shown in Fig. 5a. Second order reactions
These reactions depend ontwo particles, so they are also called bimolecular reac-tions. Given a microscopic bimolecular reaction rate α ,we calculate analogously the probability of a reactionwithin a time interval ∆ t as 1 − exp( − α ∆ t ), but we only a. b. Figure 5. Illustrations of the boundary coupling in the hybrid scheme. The domain is divided by an interface into the particleand the concentration domain. a. Boundary coupling for one diffusing species. In the particle domain, particles diffuse freelyfollowing Brownian motion. If a particle diffuses into the concentration domain (the reservoir), it is eliminated. Along theboundary cells in the concentration domain, we convert the concentration into particle number, generally a non-integer value.The integer part is the number of virtual particles, each can jump with rate γ into the boundary cells in the particle domain.The fractional part corresponds to a fractional virtual particle, whose jump rate is scaled by this fraction. Note virtual particlesare only drawn for illustration purposes, and they do not have a specific position within the boundary cell. The same procedureapplies for multiple species with up to unimolecular reactions. b. Boundary coupling for a system with three species (
A, B, C )involving a bimolecular reaction A + B → C . The coupling is analogous to the one in Fig. 5a and follows the particle-baseddynamics described in Section II A. If a red particle (A) is close enough to a blue one (B), they can react, and the product(C) is placed in the average position between the two particles. Unlike the coupling from Fig. 5a, the positions of the virtualparticles are sampled uniformly within the boundary cell. This allows for bimolecular reactions to occur within all boundarycells across the coupling boundary, which makes the coupling accurate and robust. If the position of the reaction product iswithin the concentration domain, it is eliminated. This coupling can be applied to a general system with an arbitrary numberof species with up to second-order reactions. do so if the reactants are within a distance σ of eachother. The microscopic bimolecular reaction rate needsto be consistent with the PDE reaction rate κ . FollowingEq. (20), the relation should be κ = απσ .In order to incorporate bimolecular reactions, we needto modify the boundary coupling (Fig. 5b). Consideringa reactant in the particle domain can react with anotherreactant in the concentration domain, we risk losing ac-curacy in the boundary region if we do not allow for thesereactions to happen. We can solve this by applying thefollowing steps at a given time step (Fig. 5b):1. Uniformly sample the positions of the virtual parti-cles within the boundary cells in the concentrationdomain. The fractional virtual particle is placedwith a probability equal to the fraction value.2. Allow the particles in the particle domain, alongwith the particles and virtual particles in all theboundary cells to react.3. If a reaction happens, sample the location of theproduct, usually given by the average position ofthe reactants. If the product location is in theparticle domain, place the new particle, otherwiseeliminate it from the simulation.Note the number and position of virtual particles areresampled at the beginning of each time step. Order of reactions
There are several possible waysto deal with the order of reactions at a given time step.One possibility is to choose the reaction event by draw-ing the next reaction event uniformly from all possible events while avoiding conflicting events to happen simul-taneously. In this work, we apply this approach in con-junction with a Strang splitting [46], which leads to thereaction procedure. An alternative approach is to weightthe probability of possible reaction events with their re-spective reaction probability, as done in ReaDDy 2 [21].
Reaction procedure
When called for a time interval τ , this procedure follows these steps:1. Sample all possible zeroth and first-order reactionshappening within τ /
2. Select uniformly which re-action happens and apply them while avoiding con-flicting events.2. Sample all possible second-order reactions happen-ing within τ . Select uniformly which reaction hap-pens and apply them while avoiding conflictingevents.3. Sample, select and apply again zeroth and first-order reactions happening within τ / τ . We’ve found this approach the most stable atthe coupling boundary. C. Diffusion
All the particles in the particle domain, including theadded ones during the injection and reaction steps, dif-fuse following standard Brownian motion. This can besimulated by applying the Euler-Maruyama scheme [47]to the dynamics of each molecule (Eq. (1)), x n +1 = x n + √ D ∆ tξ n , (22)where x n is the position of one molecule at the n th timeiteration. The time is t = n ∆ t and ξ n is a vector witheach entry sampled at every time step from a normal dis-tribution with mean zero and variance one. The diffusioncoefficient D will be different for every chemical species.The dimensions of ξ corresponds to the dimensionality ofthe problem (one, two or three).If a particle diffuses from the particle domain into theconcentration domain, it must be removed from the sim-ulation. The desired boundary conditions along the otherboundaries of the particle domain should be set. D. Algorithm
The process we are modeling is composed of severalcoupled processes: injection, reaction and diffusion, eachof which can follow a different timescale. From a numer-ical standpoint, it is more stable, robust and accurateto integrate them using a Strang splitting [41, 46]. Inone time step τ , the Strang splitting integration couldbe as follows: integrate injection and reactions for τ / τ and integrate again injection andreactions for τ /
2. Unlike a straightforward integration,the Strang splitting allows for some newly created parti-cles, due to injection or reactions, to diffuse and maybeeven react again within the same time step, improvingthe modeling accuracy of the coupling between the pro-cesses within one time step.In this algorithm, we implement two Strang splittings.The main Strang splitting integrates the injection, reac-tion and diffusion processes. The secondary one is imple-mented within the reaction procedure to smoothly inte-grate zeroth and first-order reactions with second-orderreactions. The scheme is as follows:
Main input variables : time iterations N , time step size∆ t , boundary cell width δx , diffusion coefficients of allspecies involved and reaction parameters for all reactionsconsidered.For every time step t ∈ { , ∆ t, ..., N ∆ t } :1. For every species X and every boundary cell i in theconcentration domain, calculate its average concen-tration c Xi and convert it into number of particles N Xi = c Xi V i . Then, get the corresponding numberof virtual and fractional virtual particles at eachboundary cell.2. For the species involving bimolecular reactions,sample uniformly the locations of the virtual par-ticles within the boundary cells. 3. Inject particles from the concentration domain intothe particle domain for half a time step, ∆ t/
2, fol-lowing the injection procedure. Note that the in-jection rate depends on the diffusion coefficient, soit will be different for different chemical species.4. Determine and apply reactions in the particle do-main occurring within half a time step, ∆ t/
2, fol-lowing the reaction procedure. This consists of asecondary Strang splitting, where zeroth and first-order reactions are applied for ∆ t/
4, second-orderreactions for ∆ t/ t/
4. Note that for bimolecularreactions, the particles can react with the virtualparticles in the boundary cells. If there are only ze-roth and first-order reactions, the secondary Strangsplitting is not necessary.5. Diffuse all particles in the particle domain for a fulltime step ∆ t using the Euler-Maruyama scheme ofEq. (22). Note the diffusion coefficient is differentfor different species.6. Determine and apply reactions for another half atime step, ∆ t/
2, in the same way as in step 4.7. Inject particles for another half a time step, ∆ t/ • Particles that are injected in step 3 cannot be usedfor reactions in step 4. However, they can be usedfor reactions in step 6. • Particles that reacted in step 4 or were generatedduring a reaction in step 4, cannot react again inanother reaction in step 4. However, they can beused for reactions in step 6. • Particles that reacted in step 6 or were generatedduring a reaction in step 6, cannot react again inanother reaction in step 6. They can only be usedfor reactions in the next time step.Following Eq. (10), we further recommend to choosethe time step ∆ t and boundary cell width δx to satisfythe relation ∆ t = δx / (2 D ). This maximizes the discretejumping probabilities [32, 41], and it is the most accurateif new particles are placed uniformly on the boundarycells of the particle domain. In general, ∆ t ≤ δx / (2 D )must be satisfied. Note we can choose several values of δx , one per species with a different diffusion coefficient.0 V. NUMERICAL RESULTS
In this section, we show the numerical results for fourexamples. The first example is a test case for diffu-sion processes to verify the coupling with a time- andspatially-dependent reservoir. The second example ver-ifies the hybrid scheme by implementing a system withfirst-order reactions. The third example implements thehybrid scheme for the Lotka-Volterra system, which in-cludes second-order reactions. This verifies the couplingfor a more realistic and complex reaction system with bi-molecular reactions. The last example couples a particle-based Lotka-Volterra system, where the reservoir is notmodeled by a PDE but by a simple constant function.This illustrates that the coupling scheme works for casesbeyond reservoirs mediated by reaction-diffusion PDEs.
A. Diffusion of one species
Consider the diffusion PDE with open boundaries forthe concentration c of a chemical species, ∂ t c ( x, t ) = D ∇ c ( x, t ) with initial condition c ( x,
0) = c δ ( x − x )and x = 2 ,
0. The solution of this equation in d dimen-sions is a Gaussian c ( x, t ) = c e −| x − | / Dt (4 πDt ) d/ . (23)We apply the scheme from Section IV to this problemin one and two dimensions d = 1 ,
2, and show its solutionin Fig. 6. We define the particle domain by x ∈ ( ∞ , x ∈ (0 , ∞ ), and wechoose δx = 0 .
05. We use the average values of c ( x, t )in the boundary cell/cells delimited by x = (0 , δx ] asthe concentration values for the particle-based simulationreservoir at each time step. The results show a compari-son between the reference analytic result and the averageconcentration obtained from several particle-based simu-lations. We can observe an excellent match, illustratinga successful coupling of a particle-based simulation to amaterial bath mediated by a PDE.In the pure diffusion case, we obtain very good resultseven when averaging over a relatively small number ofsimulations. In more complex cases, we will average overa larger number of simulations, and we will verify thescheme using the Jensen-Shannon (JS) divergence. B. Proliferation of one species (first-order)
In this example, we model the diffusion of one species A , along with the first-order reaction A − (cid:42) A. The corresponding reaction diffusion PDE is ∂ t c = D ∇ c + κ c, (24) a. Particle domain Concentration domain b. Figure 6. Diffusion coupling results in one and two dimen-sions. a. Diffusion coupling results in one dimension at fourdifferent times. The analytic solution of the diffusion PDEwith D = 1 is shown in orange (Eq. (23) with n = 1). Theaverage concentration of 200 particle-based simulations in theparticle domain is shown as a blue histogram. The parametersused were δx = 0 .
05 and ∆ t = δx / (2 D ). b. Coupling resultsfor the two dimensional extension with the same parameters,and the colorbar values denoting concentration. The solutionis showed at two times. The top row shows the analytic solu-tion (Eq. (23) with n = 2). The bottom row shows the hybridsimulation with the interface at x = 0. In the particle domain(left half), it shows a histogram of the average concentrationof 200 particle-based simulations. In the concentration do-main (right half), it shows the reservoir, which correspondsto the reference value given by the PDE. with c ( x, t ) the concentration of A . We use a diffusioncoefficient of D = 0 . κ = 0 .
1. In Fig. 7a, we solve this equation using a finitedifference scheme on the domain x ∈ [0 , × [0 ,
12] withan indicator function as initial condition and Neumannboundary conditions. The indicator function consists ofa concentration of 50 in a 2 × . , Figure 7. Coupling results for diffusion with a proliferation first-order reaction at three different times t = 4 , ,
9. The colorbar indicates the value of the concentration. a. Reference solution using the finite difference scheme in a domain of 12 × . × .
12 and a time step of 0 . b. Solution of the hybrid system with the interface at x = 6. The lefthalf corresponds to the particle domain, and it shows the average over 3000 particle-based simulations, each using a time stepof ∆ t = 0 .
01. The boundary cell width δx is chosen to satisfy ∆ t = δx / (2 D ). The right half is the same as in the referencesolution, and it serves as the material reservoir for the particle-based simulations following the scheme from Section IV. c. JSdivergence between the reference concentration and the averaged concentration of the hybrid simulations at the same threetimes. The x-axis is the number of averaged hybrid simulations. Each point is calculated using 500 bootstrapped samples. and 0 elsewhere. This is be our reference solution. Figure7b shows the solution using the hybrid scheme from Sec-tion IV. In the particle-based simulation, particles thatproliferate are placed in the same location as the sourceparticle.In Fig. 7c., we verify the scheme using the JS di-vergence, which is a measure of the difference betweendistributions or histograms. We compare the referencesolution against averages of the hybrid simulation. Asthe number of hybrid simulations used to produce theaverage grows, the JS divergence becomes closer to zero,showing the expected convergence.
C. Lotka-Volterra dynamics
The Lotka-Volterra dynamics are a chemical kineticssystem, which can be understood in terms of preys A andpredators B . If a predator meets one prey, the prey canbe eaten and the predator multiplies. Moreover, the preycan multiply and predators can die, both independently.This is an important system since it captures the complexdynamics that could appear in most relevant reaction-diffusion applications. For instance, the models used inepidemiology for the spread of infectious diseases, suchas the SIR model, follow similar dynamics. The kinetics are condensed in the following reactions A − (cid:42) AA + B − (cid:42) B (25) B − (cid:42) ∅ If the system is not well-mixed and the numbers of preda-tors and preys is large, the macroscopic dynamics can bedescribed by the PDE ∂ t c A = D A ∇ c A + κ c A − κ c A c B ,∂ t c B = D B ∇ c B + κ c A c B − κ c B , (26)where the macroscopic rates, κ , κ and κ , correspond tothe proliferation rate of preys, the rate at which preys areeaten and the rate at which predators die, respectively.The domain is x ∈ [0 , × [0 , × ,
5) and 0 elsewhere. For the predators, itconsists of a concentration of 10 in a 1 × , κ = 0 . κ = 0 . σ = 0 .
01, they react with α = 0 . κ = απσ ,following Eq. (20). The preys diffuse with D A = 0 . D B = 0 .
1. We use Neumann(reflective) boundary conditions in all the boundaries.2
Figure 8. Solutions for the Lotka-Volterra dynamics with diffusion at three different times t = 4 , , ×
10 domain. Thecolor bar indicates the value of the concentration. a. Reference solution of the preys using a finite difference scheme with atime step of 0 .
002 and grid size of 0 . b. Solution of the preys in the hybrid simulation with the interface at x = 5. Theleft half of the domain consists of the average over 3000 particle-based simulation using a time step of ∆ t = 0 . δx A ) and one for the predators ( δx B ), each satisfies the relation ∆ t = δx k / (2 D k )with k = A or B . The right half is the same as the reference solution, and it is used as the reservoir for the particle-basedsimulation. The coupling used is the one described in Section IV. c. Reference solution of the predators ( B ) corresponding tothe same simulation as in Fig. 8a. d. Solution of the hybrid simulation for the predators corresponding to the same simulationas in Fig. 8b. e. & f. JS divergence, for preys and predators respectively, calculated between the reference concentration andthe averaged concentration of the hybrid simulations at the three times. The x-axis is the number of hybrid simulations usedto calculate the average, and each point is calculated using 500 bootstrapped samples.
In Fig. 8, we show the reference simulation andthe averaged hybrid simulation solution results for bothprey and predators. The coupling produces an excel-lent match. We further verify the results using the JSdivergence, which shows convergence as the number ofaveraged simulations is increased.
D. Reservoirs constant in time
We focus again on the Lotka-Volterra dynamics de-scribed by Eq. (25). However, instead of coupling theparticle-based simulation to a material reservoir medi-ated by a PDE, we couple it to a reservoir with a constantconcentration in time (not in space). This is helpful when3
Figure 9. Coupling results for a Lotka-Volterra system coupled to a constant in time prey reservoir at three different times t = 4 , ,
9. The color bar indicates the value of the concentration. a. Reference solution of preys using a finite difference schemewith 60 ×
30 grid cells and a time step of 0 . b. Solution of preys in the hybrid simulation, consisting of the average over3000 particle-based simulation. Each particle-based simulations used a time step of ∆ t = 0 .
01 and was coupled to the reservoirfollowing the scheme from Section IV. We use two boundary cell widths, one for the preys ( δx A ) and one for the predators( δx B ), each satisfies the relation ∆ t = δx k / (2 D k ) with k = A or B . c. Reference solution of the predators corresponding to thesame simulation as in Fig. 9a. d. Solution of the hybrid simulation for the predators corresponding to the same simulation asin Fig. 9b. e. and f. JS divergence, for preys and predators respectively, calculated between the reference concentration andthe averaged concentration of the hybrid simulations at the three plotted times. The x-axis is the number of averaged hybridsimulations, and each point is calculated using 500 bootstrapped samples. modeling a reservoir with a specific spatial distribution.The domain is x ∈ [0 , × [0 , σ = 0 . κ = 0 . κ = απσ and ∆ t = 0 .
01. The initial condi-tion is zero preys and an indicator function with a con- centration of predators of 30 in a rectangle of 2 × , πx/ VI. DISCUSSION
The main goal of this work was to improve currentmodeling techniques for biochemical open systems, asthey are extremely relevant to model life-related pro-cesses. In this paper, we contributed to this goal bydeveloping models and numerical schemes, which arecapable of consistently coupling particle-based reaction-diffusion processes with reservoirs mediated by reaction-diffusion PDEs, i.e. with time-dependent and spatiallynon-homogeneous reservoirs.The coupling was rendered possible by the two theo-retical results from Section III. The first result derivedthe mean-field of a particle-based diffusion model in con-tact with a constant concentration reservoir. This re-sulted in a diffusion PDE with a constant concentrationboundary condition and elucidated the relation betweenthe reservoir dynamics in the two models. The secondresult derived the mean-field limit of reaction-diffusionprocesses with a bimolecular reaction. We recovered thecorresponding reaction-diffusion PDE, and we obtained aprecise connection between the microscopic and macro-scopic parameters, specifically the bimolecular reaction rates and the reaction radius. Section IV further employsthese two theoretical results to build the coupling numer-ical scheme, which can be used in any reaction-diffusionsystem with up to second-order reactions.In Section V, we implement the coupling schemefor four representative examples: pure diffusion of onespecies, proliferation of one species, Lotka-Volterra dy-namics and Lotka-Volterra dynamics coupled to a reser-voir constant in time. In order to verify the scheme, wecompare the average of several particle-based simulationswith the theoretical mean-field given by the reaction-diffusion PDE. The difference between the averaged andthe reference solution is quantified with the JS diver-gence. We obtain excellent results for all the examples.We finally point out the result from Section III B ne-glected the covariances when taking the expectation ofthe mixed Poisson process. This is a valid approxima-tion if there is a very large number of particles. How-ever, for mesoscopic scales in the number of particles, thecovariances would yield a non-negligible modified mean-field behavior, which would possibly yield alternative al-gorithms. The authors are currently researching thesetopics [16].
VII. ACKNOWLEDGEMENTS
We gratefully acknowledge support by the DeutscheForschungsgemeinschaft (grants SFB1114, projects C03and A04), the Berlin Mathematics Research CenterMATH+, Project AA1-6, and the European researchcouncil (ERC starting grant307494 “pcCell”). [1] H. Qian, “Phosphorylation energy hypothesis: openchemical systems and their biological functions,”
Annu.Rev. Phys. Chem. , vol. 58, pp. 113–142, 2007.[2] A. Agarwal, J. Zhu, C. Hartmann, H. Wang, andL. Delle Site, “Molecular dynamics in a grand ensem-ble: Bergmann–lebowitz model and adaptive resolutionsimulation,”
New J. Phys. , vol. 17, no. 8, p. 083042, 2015.[3] L. Delle Site and R. Klein, “Liouville-type equations forthe n-particle distribution functions of an open system,”
Phys. Rev. Res. , 2019.[4] L. Delle Site, C. Krekeler, J. Whittaker, A. Agarwal,R. Klein, and F. H¨ofling, “Molecular dynamics of opensystems: Construction of a mean-field particle reservoir,”
Adv. Theory Simul. , vol. 2, no. 5, p. 1900014, 2019.[5] M. Dibak, M. J. del Razo, D. De Sancho, C. Sch¨utte, andF. No´e, “MSM/RD: Coupling markov state models ofmolecular kinetics with reaction-diffusion simulations,”
J. Chem. Phys. , vol. 148, no. 21, p. 214107, 2018.[6] L. Sbail`o and F. No´e, “An efficient multi-scale green’sfunction reaction dynamics scheme,”
The Journal ofchemical physics , vol. 147, no. 18, p. 184106, 2017.[7] A. Vijaykumar, P. G. Bolhuis, and P. R. ten Wolde,“Combining molecular dynamics with mesoscopic green’sfunction reaction dynamics simulations,”
The Journal ofchemical physics , vol. 143, no. 21, p. 214102, 2015. [8] D. F. Anderson and T. G. Kurtz,
Stochastic analysis ofbiochemical systems . Springer, 2015, vol. 1.[9] T. G. v. Kurtz, “The relationship between stochastic anddeterministic models for chemical reactions,”
J. Chem.Phys. , vol. 57, no. 7, pp. 2976–2978, 1972.[10] H. Qian, “Nonlinear stochastic dynamics of mesoscopichomogeneous biochemical reaction systems—an analyti-cal theory,”
Nonlinearity , vol. 24, no. 6, p. R19, 2011.[11] L. Arnold, “On the consistency of the mathematical mod-els of chemical reactions,” in
Dynamics of synergetic sys-tems . Springer, 1980, pp. 107–118.[12] J. Feng, “The hydrodynamic limit for the reaction diffu-sion equation—an approach in terms of the gpv method,”
J. Theor. Probab. , vol. 9, no. 2, pp. 285–299, 1996.[13] S. Hellander, A. Hellander, and L. Petzold, “Reactionrates for mesoscopic reaction-diffusion kinetics,”
Phys.Rev. E , vol. 91, no. 2, p. 023312, 2015.[14] S. A. Isaacson, “Relationship between the reaction–diffusion master equation and particle tracking models,”
J. Phys. A: Math. Theor. , vol. 41, no. 6, p. 065003, 2008.[15] ——, “A convergent reaction-diffusion master equation,”
The Journal of chemical physics , vol. 139, no. 5, p.054101, 2013.[16] M. Kostr`e and M. J. del Razo, “Hydrodynamic limits ofreaction-diffusion processes,” (in preparation) . [17] B. Franz, M. B. Flegg, S. J. Chapman, and R. Erban,“Multiscale reaction-diffusion algorithms: Pde-assistedbrownian dynamics,” SIAM J. on Applied Mathematics ,vol. 73, no. 3, pp. 1224–1247, 2013.[18] C. A. Smith and C. A. Yates, “The auxiliary re-gion method: a hybrid method for coupling pde-and brownian-based dynamics for reaction–diffusion sys-tems,”
R. Soc. Open Sci. , vol. 5, no. 8, p. 180920, 2018.[19] R. Erban and S. J. Chapman, “Stochastic modelling ofreaction–diffusion processes: algorithms for bimolecularreactions,”
Phys. Biol. , vol. 6, no. 4, p. 046001, 2009.[20] M. Kostr´e and M. J. del Razo, “Code and datato accompany this paper,” 2020. [Online]. Available:https://github.com/MargKos/multiscaleRD[21] M. Hoffmann, C. Fr¨ohner, and F. No´e, “ReaDDy 2: Fastand flexible software framework for interacting-particlereaction dynamics,”
PLOS Comput. Biol. , vol. 15, no. 2,p. e1006830, 2019.[22] S. S. Andrews, N. J. Addy, R. Brent, and A. P. Arkin,“Detailed simulations of cell biology with smoldyn 2.1,”
PLoS Comput. Biol. , vol. 6, no. 3, 2010.[23] I. I. Moraru, J. C. Schaff, B. M. Slepchenko, M. Bli-nov, F. Morgan, A. Lakshminarayana, F. Gao, Y. Li,and L. M. Loew, “Virtual cell modelling and simulationsoftware environment,”
IET Syst. Biol. , vol. 2, no. 5, pp.352–362, 2008.[24] J. R. Stiles, T. M. Bartol et al. , “Monte carlo methodsfor simulating realistic synaptic microphysiology usingmcell,”
Computational neuroscience: realistic modelingfor experimentalists , pp. 87–127, 2001.[25] C. Fr¨ohner and F. No´e, “Reversible interacting-particlereaction dynamics,”
J. Phys. Chem. B , vol. 122, no. 49,pp. 11 240–11 250, 2018.[26] M. Doi, “Stochastic theory of diffusion-controlled reac-tion,”
Journal of Physics A: Mathematical and General ,vol. 9, no. 9, p. 1479, 1976.[27] E. Teramoto and N. Shigesada, “Theory of bimolecu-lar reaction processes in liquids,”
Progress of TheoreticalPhysics , vol. 37, no. 1, pp. 29–51, 1967.[28] N. Agmon, “Diffusion with back reaction,”
The Journalof chemical physics , vol. 81, no. 6, pp. 2811–2817, 1984.[29] N. Agmon and A. Szabo, “Theory of reversible diffusion-influenced reactions,”
J. Chem. Phys , vol. 92, no. 9, pp.5270–5284, 1990.[30] J. C. Cavallo and M. B. Flegg, “Reversible doi and smolu-chowski kinetics for high-order reactions,”
SIAM Jour-nal on Applied Mathematics , vol. 79, no. 2, pp. 594–618,2019.[31] M. J. Del Razo, W. Pan, H. Qian, and G. Lin, “Fluores-cence correlation spectroscopy and nonlinear stochasticreaction–diffusion,”
J. Phys. Chem. B , vol. 118, no. 25,pp. 7037–7046, 2014.[32] M. J. del Razo and H. Qian, “A discrete stochastic for-mulation for reversible bimolecular reactions via diffusionencounter,”
Comm. Math. Sci. , vol. 14, no. 6, pp. 1741–1772, 2016.[33] I. V. Gopich and A. Szabo, “Kinetics of reversible dif-fusion influenced reactions: the self-consistent relaxationtime approximation,”
J. Chem. Phys , vol. 117, no. 2, pp.507–517, 2002.[34] S. S. Khokhlova and N. Agmon, “Comparison of alternateapproaches for reversible geminate recombination,”
Bull.Korean Chem. Soc , vol. 33, no. 3, p. 1021, 2012. [35] H. Kim and K. J. Shin, “Exact solution of the reversiblediffusion-influenced reaction for an isolated pair in threedimensions,”
Phys. Rev. Lett. , vol. 82, pp. 1578–1581,1999.[36] S. S. Andrews and D. Bray, “Stochastic simulation ofchemical reactions with spatial resolution and singlemolecule detail,”
Phys. Biol. , vol. 1, no. 3, p. 137, 2004.[37] D. A. Beard and H. Qian,
Chemical biophysics: quantita-tive analysis of cellular systems . Cambridge UniversityPress, 2008.[38] R. J. LeVeque,
Finite difference methods for ordinaryand partial differential equations: steady-state and time-dependent problems . Siam, 2007, vol. 98.[39] M. Kostr`e, “Hybrid models and simulation of reaction-diffusion processes,”
Master thesis , 2019.[40] S. Salsa, F. Vegni, A. Zaretti, and P. Zunino,
A primer onPDEs: models, methods, simulations . Springer Science& Business Media, 2013.[41] M. J. del Razo, H. Qian, and F. No´e, “Grand canon-ical diffusion-influenced reactions: A stochastic theorywith applications to multiscale reaction-diffusion simula-tions,”
J. Chem. Phys. , vol. 149, no. 4, p. 044102, 2018.[42] W. J. Heuett and H. Qian, “Grand canonical markovmodel: a stochastic theory for open nonequilibrium bio-chemical networks,”
J. Chem. Phys. , vol. 124, no. 4, p.044110, 2006.[43] J. Grandell,
Mixed poisson processes . CRC Press, 1997,vol. 77.[44] D. Schnoerr, R. Grima, and G. Sanguinetti, “Cox pro-cess representation and inference for stochastic reaction–diffusion processes,”
Nature communications , vol. 7, p.11729, 2016.[45] H. Wang, C. S. Peskin, and T. C. Elston, “A robustnumerical algorithm for studying biomolecular transportprocesses,”
J. Theor. Biol. , vol. 221, no. 4, pp. 491–511,2003.[46] S. MacNamara and G. Strang, “Operator splitting,” in
Splitting Methods in Communication, Imaging, Science,and Engineering . Springer, 2016, pp. 95–114.[47] D. J. Higham, “An algorithmic introduction to numeri-cal simulation of stochastic differential equations,”
SIAMreview , vol. 43, no. 3, pp. 525–546, 2001.[48] P. Deuflhard and M. Weiser,
Numerische Mathematik:Adaptive L¨osung partieller Differentialgleichungen. 3 .Walter de Gruyter, 2011, vol. 3.[49] B. S. Jovanovi´c and E. S¨uli,
Analysis of finite differ-ence schemes: for linear partial differential equationswith generalized solutions . Springer Science & BusinessMedia, 2013, vol. 46.
Appendix A: Finite difference solution of thereaction-diffusion PDE
We will show how to solve numerically the Poissonequation D ∇ u = f, with D the diffusion constant, f an arbitrary functionand homogeneous Neumann boundary conditions dudη = 06on a two dimensional finite domain (this corresponds toa reflecting boundary in a PBS) with a Finite Differencescheme [48].As the name already suggests we discretize the domainin cells of length and wide h , such that we can set up alinear system Au = f (A.1)by using the five-point stencil. This gives us the solu-tion of u for each cell. Matrix A is the so called discrete Laplace operator that depends on the boundary condi-tions. In our case, we will use homogeneous Neumannboundary conditions, so it is defined by A := − Dh B − I · · · − I . . . − I ...... − I . . . − I · · · − I B , where I is the corresponding identity matrix and B = − · · ·− − · · · − − · · · − and h is the length of the grid cell. The solution is ob-tained by inverting the matrix A in Eq. (A.1). It is possible to extend the Finite-Difference scheme tosolve the time dependent Heat-equation ∂ t u = D ∇ u. (A.2)To do so, we can combine the Finite-Difference schemewith the Euler method. Let u n denote the solution at dis-crete time t n . Then we can discretize in time the equation(A.2) for every time step n and small ∆ t to obtain u n +1 = u n + ∆ tAu n = [ I + ∆ tA ] u n . (A.3)If we use the implicit Euler scheme, we obtain u n +1 = u n + ∆ tAu n +1 ⇒ [ I − ∆ tA ] u n +1 = u n , If we combine half a time step with Euler method and halfa time step with the implicit Euler method, we obtaina very numerically robust method that is second orderaccurate in both space and time. It is called the Crank-Nicolson method [38, 49], (cid:20) I −
12 ∆ tA (cid:21) u n +1 = (cid:20) I + 12 ∆ tA (cid:21) u n ⇒ u n +1 = (cid:20) I −
12 ∆ tA (cid:21) − (cid:20) I + 12 ∆ tA (cid:21) u n ..