Scaling limits of spatial compartment models for chemical reaction networks
aa r X i v : . [ m a t h . P R ] O c t The Annals of Applied Probability (cid:13)
Institute of Mathematical Statistics, 2015
SCALING LIMITS OF SPATIAL COMPARTMENT MODELS FORCHEMICAL REACTION NETWORKS
By Peter Pfaffelhuber and Lea Popovic Albert-Ludwigs-Universit¨at Freiburg and Concordia University
We study the effects of fast spatial movement of molecules on thedynamics of chemical species in a spatially heterogeneous chemicalreaction network using a compartment model. The reaction networkswe consider are either single- or multi-scale. When reaction dynam-ics is on a single-scale, fast spatial movement has a simple effect ofaveraging reactions over the distribution of all the species. When re-action dynamics is on multiple scales, we show that spatial movementof molecules has different effects depending on whether the movementof each type of species is faster or slower than the effective reactiondynamics on this molecular type. We obtain results for both whenthe system is without and with conserved quantities, which are linearcombinations of species evolving only on the slower time scale.
1. Introduction.
When chemical species react, they are present in some(open or closed) system with a spatial dimension. Most models of chemi-cal reaction systems describe the evolution of the concentration of chemicalspecies and ignore both stochastic and spatial effects inherent in the system.This can be justified by law of large number results when both: the num-ber of species across all molecular types is large, and when the movementof molecules within the system is much faster than the chemical reactionsthemselves. In applications where these assumptions hold, the system isspatially homogeneous, and the use of the deterministic law of mass actionkinetics is approximately appropriate [Kurtz (1970)].However, in biological cells, low numbers of certain key chemical speciesinvolved in the reaction systems result in appreciable noise in gene expression
Received February 2013; revised July 2014. Supported by the University Faculty award and Discovery grant of NSERC (NaturalSciences and Engineering Council of Canada).
AMS 2000 subject classifications.
Primary 60J27, 60J25; secondary 60F17, 92C45,80A30.
Key words and phrases.
Chemical reaction network, multiple time scales, stochastic av-eraging, scaling limits, quasi-steady state assumption.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Applied Probability ,2015, Vol. 25, No. 6, 3162–3208. This reprint differs from the original inpagination and typographic detail. 1
P. PFAFFELHUBER AND L. POPOVIC and many regulatory functions of the cell, and lead to cell–cell variabilityand different cell fate decisions [McAdams and Arkin (1997), Elowitz et al.(2002)]. In order to understand the key effects of intrinsic noise in chem-ical reaction networks on the overall dynamics, one needs to derive newapproximations of stochastic models describing the evolution of molecularcounts of chemical species. Furthermore, the cell is not a spatially homoge-neous environment, and it has been repeatedly demonstrated that spatialconcentration of certain molecules plays an important role in many cellu-lar processes [Howard and Rutenberg (2003), Takahashi, Tanase-Nicola andTen Wolde (2010)]. Deterministic spatial models (reaction–diffusion PDEs)are insufficient for this purpose, since local fluctuations of small molecularcounts can have propagative effects, even when the overall total number ofmolecules in the cell of the relevant species is high.Numerical simulations of biochemical reactions are indispensable sincethe stochastic spatial dynamics of most systems of interest is analyticallyintractable. A number of different simulation methods have been developedfor this purpose, ranging from exact methods (accounting for each stochasticevent) to approximate methods (replacing exact stochastics for some aspectsof the system with approximate statistical distributions); see, for example,Fange et al. (2010), Drawert et al. (2010), Jeschke, Ewald and Uhrmacher(2011). For effective computations, a mesoscopic form of the full stochasticspatial reaction model is necessary. These are compartment models in whicha heterogeneous system is divided into homogeneous subsystems in eachof which a set of chemical reactions are performed. Molecular species aredistributed across compartments and their diffusion is modeled by movesbetween neighboring compartments [Burrage et al. (2011)].An important mathematical feature of models of biochemical reactionslies in the essential multi-scale nature of the reaction processes. In somecases, all chemical species are present in a comparable amount and changetheir concentrations on the same temporal scale. We call this a single-scale reaction network. However, if at least one chemical species changes its abun-dance (relative to its abundance) on a much faster time scale, we call this a multi-scale reaction network. The fast change of the concentrations of somechemical species then has an impact on the dynamics of the slow species.When one adds spatial movement of species into the system, then the specieswith fast movement can have an averaging effect on the dynamics as well.The overall dynamics depends on how these two averaging factors interact,and subsequently determines the evolution of the slow species on the finaltime scale of interest.In this paper, we analyze the effects of movement of molecules on thedynamics of the molecular counts of chemical species. We consider a fi-nite number of compartments in which different reactions can happen, withspecies moving between compartments, and derive results for the evolution
CALING LIMITS OF SPATIAL COMPARTMENT MODELS of the sum total of molecules in all compartments. We consider the followingfor chemical reactions within compartments: (i) single-scale chemical reac-tions, (ii) multi-scale chemical reactions without conserved quantities and(iii) multi-scale chemical reactions with conserved quantities. We focus onthe derivation of simplified models obtained as limits of rescaled versions ofthe original model [in the spirit of Kang, Kurtz and Popovic (2014), Ballet al. (2006), Franz, Liebscher and Zeiser (2012)]. We stress that our resultsare for mesoscopic models of spatial systems, as opposed to models in whichthe number of compartments increases and the size of compartments shrinks[Blount (1994), Kouritzin and Long (2002), Kotelenez (1988)].Our goal is to find reduced models that capture all relevant stochasticfeatures of the original model, while focusing only on the quantities that areeasy to measure (sum total of molecular numbers for each species) in thesystem. Our results for this reduced dynamics make stochastic simulationsalmost trivial, while at the same time differentiating (being able to detect)between different cases of system heterogeneity and between different rela-tive time scales of species movement.1.1. Outline of results.
After introducing a model for chemical reac-tions in Section 2.1, we provide results in nonspatial systems which weneed later for spatial results: asymptotics for a single-scale reaction net-work (Lemma 2.4), asymptotics for two-scale systems without (Lemma 2.7)and with conserved quantities on the fast time scale (Lemma 2.12), as wellas extensions to three-scale systems (Lemma 2.10). Examples of each typeof system are provided as well (Examples 2.5, 2.8, 2.13). In Section 3, wegive results for spatial compartment models: for single-scale spatial systems(Theorem 3.7), for two-scale spatial systems (Theorem 3.13) in the absenceof conserved quantities, and in the presence of conserved quantities (The-orem 3.23). We also place the earlier examples in a spatial setting (Exam-ples 3.10, 3.19, 3.27). We conclude the paper with a discussion of possibleimplications and extensions.
Remark 1.1 (Notation). For some Polish space E , we denote the setof continuous (bounded and continuous, continuous with compact support),real-valued functions by C ( E ) [ C b ( E ), C c ( E )]. In general, we write x := ( x k ) k for vectors and x := ( x ik ) i,k for matrices. In addition, x i · is the i th lineand x · k is the k th column of x . We denote by D ( I ; E ) the space of c`adl`agfunctions I ⊆ R → E , which is equipped, as usual, with the Skorohod topol-ogy, metrized by the Skorohod metric, d Sk . For sets F, F ′ ⊆ E , we write F − F ′ := { f ∈ F : f / ∈ F ′ } . P. PFAFFELHUBER AND L. POPOVIC
2. Chemical reactions in a single compartment.
Before we present ourmain results on spatial systems in the next section, we provide basic resultson reaction networks in a single compartment, which are special cases oftheorems given in Kang and Kurtz (2013). Our results for the spatial caseboth rely on them and are proved using similar techniques.We consider a set I of different chemical species, involved in K differentreactions of the form ( ν ik ) i ∈I → ( ν ′ ik ) i ∈I (2.1)with ν = ( ν ik ) i ∈I ,k ∈K , ν ′ = ( ν ′ ik ) i ∈I ,k ∈K ∈ Z I×K + and ν ik = l if l molecules ofthe chemical species i take part in reaction k and ν ′ ik = l if reaction k pro-duces l molecules of species i . In the chemical reaction literature, ζ = ν ′ − ν is called the stoichiometric matrix of the system, and P i ∈I ν ik the order ofreaction k . In addition, we set ζ · k := ( ζ ik ) i ∈I .2.1. The Markov chain model and the rescaled system.
Denoting by X i ( t )the number of molecules of species i at time t , we assume that ( X ( t )) t ≥ with X ( t ) = ( X i ( t )) i ∈I is solution of X i ( t ) = X i (0) + X k ∈K ζ ik Y k (cid:18)Z t Λ CR k ( X ( u )) du (cid:19) , (2.2)where the Y k ’s are independent (rate 1) Poisson processes and Λ CR k ( X ( u )) isthe reaction rate of reaction k at time u , k ∈ K . We will assume throughoutthe following. Assumption 2.1 (Dynamics of unscaled single compartment reaction net-work). The reaction network dynamics satisfies the following conditions:(i) The reaction rate x Λ CR k ( x ) is a nonnegative locally Lipshitz, lo-cally bounded function and Λ CR k = 0, k ∈ K .(ii) Given ( Y k ) k ∈K , the time-change equation (2.2) has a unique solution.The most important chemical reaction kinetics is given by mass action,that is, Λ CR k ( x ) = κ ′ k Y i ∈I ν ik ! (cid:18) x i ν ik (cid:19) (2.3)for constants κ k . In other words, the rate of reaction k is proportional tothe number of possible combinations of reacting molecules. Solutions to (2.2)can be guaranteed by using, for example, Ethier and Kurtz (1986), Theo-rem 6.2.8; see also their Remark 6.2.9(b). Note, however, that Assump-tion 2.1(i) does not suffice to guarantee a global solution to (2.2), since it CALING LIMITS OF SPATIAL COMPARTMENT MODELS has to be certain—usually by imposing some growth condition—that thesolution does not become infinite in finite time.Chemical reaction networks in many applications involve chemical specieswith vastly differing numbers of molecules and reactions with rate constantsthat also vary over several orders of magnitude Ball et al. (2006), Examples.This wide variation in number and rate yield phenomena that evolve onvery different time scales. Recognizing that the variation in time scales isdue both to variation in species number and to variation in rate constants,we normalize species numbers and rate constants by powers of a param-eter N which we assume to be large, and consider a sequence of models,parametrized by N ∈ N . Rescaled versions of the original model, under cer-tain assumptions, have a limit as N → ∞ . We will use stochastic equationsof the form (2.2) driven by independent Poisson processes to show conver-gence, exploiting the law of large numbers and martingale properties of thePoisson processes. We rely heavily on the stochastic averaging methods thatdate back to Khasminskii, for which we follow the formalism in terms ofmartingale problems from Kurtz (1992).We rescale the system as follows: consider the solution ( X N ( t )) t ≥ of (2.2)with the chemical reaction rates Λ CR k replaced by Λ CR ,Nk . For real-valued α = ( α i ) i ∈I , β = ( β k ) k ∈K , γ (2.4)with α i ≥ , i ∈ I , we denote the ( α, β, γ )-rescaled system by V Ni ( t ) := N − α i X Ni ( N γ t ) , λ CR ,Nk ( v ) := N − β k Λ CR ,Nk (( N α i v i ) i ∈I ) , (2.5)where ( α, β, γ ) is chosen so that: V Ni = O (1) , i ∈ I , for all time (a.s. doesnot go infinity in finite time, but also does not have a.s. zero limit for alltime), and λ CR ,Nk = O (1) , k ∈ K (for all values of v when it is not equal tozero). We will restrict to the case γ = 0 which can always be achieved whenconsidering β ′ k = β k + γ, k ∈ K . From (2.2), we see that the ( α, β, γ )-rescaledsystem V N := ( V N ( t )) t ≥ is a solution to the system of stochastic equations V Ni ( t ) = V Ni (0) + X k ∈K N − α i ζ ik Y k (cid:18) N β k + γ Z t λ CR ,Nk ( V N ( u )) du (cid:19) , (2.6)which has a unique solution thanks to Assumption 2.1. In vector notation,we use the diagonal matrix N − α with i th diagonal entry N − α i and write V N ( t ) = V N (0) + X k ∈K N − α ζ · k Y k (cid:18) N β k + γ Z t λ CR ,Nk ( V N ( u )) du (cid:19) . (2.7)The reaction rates satisfy the following. P. PFAFFELHUBER AND L. POPOVIC
Assumption 2.2 (Dynamics of scaled single compartment reaction net-work). There exist locally Lipshitz functions λ CR k : R I + → R + , k ∈ K with N − β k Λ CR ,Nk (( N α i v i ) i ∈I ) N →∞ −→ λ CR k ( v )(2.8)uniformly on compacts. [Without loss of generality, we will assume thatconvergence in (2.8) is actually an identity; our results easily generalize bythe assumed uniform convergence on compacts.]In the special case of mass action kinetics (2.3), if α i = 1 for all i ∈ I and κ k = κ ′ k N − ( P i ν ik )+1 with β k = 1 and some κ ′ k > k ∈ K , then N − β k Λ CR ,Nk (( N α i v i ) i ∈I ) N →∞ −→ κ ′ k Y i ∈I v ν ik i . The polynomial on the right-hand side is known in the literature for deter-ministic chemical reaction systems as the mass action kinetic rate.2.2.
Single scale systems.
For i ∈ I , the set of reactions which changethe number of species i is K i := { k ∈ K : ζ ik = 0 } (a reaction of the form A + B → A + C does not change the number ofspecies A ). A chemical reaction network is a single scale system if ( α, β, γ )from (2.4) satisfy max k ∈K i β k + γ = α i , i ∈ I . (2.9)For i ∈ I , let K ∗ i ⊆ K i be the set of reactions such that β k + γ = α i , and let K ∗ = S i ∈I K ∗ i . Define ζ ∗ by ζ ∗ ik = lim N →∞ N − α i N β k + γ ζ ik . Then ζ ∗ is the matrix whose i ∈ I , k ∈ K ∗ i entries are ζ ik and its i ∈ I , k ∈K i − K ∗ i entries are zero. Let I ◦ be the subset of species with α i = 0, calledthe discrete species , and let K ∗◦ = S i ∈I ◦ K ∗ i , called the slow reactions . Let I • be the subset of species with α i >
0, called the continuous species , and let K ∗• = S i ∈I • K ∗ i , called the fast reactions . Then K ∗ = K ∗◦ ∪ K ∗• . Note that bydefinition I ◦ and I • are disjoint, and by definition of K ∗ i (and as reactionrates come with a single scaling N β k + γ ), K ∗◦ and K ∗• are also disjoint. In thelimit of the rescaled system, the species indexed by I ◦ are Z + -valued (hencethe name discrete species), while the species indexed by I • are R + -valued(continuous species). See Table 1 for an overview of these definitions. Wenext assume the following. CALING LIMITS OF SPATIAL COMPARTMENT MODELS Table 1
An overview of different sets and possibilities in the case γ = 0 . The set I is split intodiscrete ( I ◦ ) and continuous ( I • ) chemical species, while the set K ∗ is split into slow( K ∗◦ ) and fast ( K ∗• ) reactions. The gray boxes give the reactions which still appear in thelimit dynamics. A special feature of single-scale systems is that discrete species areexactly changed through slow reactions, and continuous species are changed by fastreactions. In particular, discrete species are not changed by fast reactions. This isdifferent in multi-scale networks; see Table 2 Slow reactions Fast reactions k ∈ K ∗◦ , β k = 0 k ∈ K ∗• , β k > Discrete species, α i = 0 , i ∈ I ◦ ζ ∗ ik (cid:26) = 0 , k ∈ K ∗ i ,= 0 , else ζ ik = ζ ∗ ik = 0Continuous species, α i > , i ∈ I • ζ ik = 0 or β k + γ − α i < ⇒ ζ ∗ ik = 0 ζ ∗ ik (cid:26) = 0 , k ∈ K ∗ i ,= 0 , else Assumption 2.3 (Dynamics of the reaction network). For Poisson pro-cesses ( Y k ) k ∈K ∗◦ , the time-change equation V ( t ) = V (0) + X k ∈K ∗◦ ζ ∗· k Y k (cid:18)Z t λ CR k ( V ( u )) du (cid:19) (2.10) + X k ∈K ∗• ζ ∗· k Z t λ CR k ( V ( u )) du has a unique solution V := ( V ( t )) t ≥ .Actually, the last display is shorthand notation for V i ( t ) = V i (0) + X k ∈K ∗◦ ζ ∗ ik Y k (cid:18)Z t λ CR k ( V ( u )) du (cid:19) , i ∈ I ◦ ,V i ( t ) = V i (0) + X k ∈K ∗• ζ ∗ ik Z t λ CR k ( V ( u )) du, i ∈ I • . Lemma 2.4 (Convergence of single-scale reaction networks).
Let V N be the vector process of rescaled species amounts for the reaction networkwhich is the unique solution to (2.6). Assume ( α, β, γ ) from (2.5) satisfy thesingle scale system assumptions (2.9), and Assumptions 2.2 and 2.3 for therescaled reaction network are satisfied. Then, if V N (0) = ⇒ V (0) , the processof rescaled amounts V N converges weakly to the solution V of (2.10) in theSkorohod topology. P. PFAFFELHUBER AND L. POPOVIC
Table 2
As in the single-scale case, the set I is split into discrete and continuous chemicalspecies. In addition, discrete and continuous species are either changed on the fast orslow time scale. (This means that I f ◦ , I f • , I s ◦ , I s • are disjoint sets.) The set of reactions issplit into several categories, which can overlap. Here, k ∈ K f ◦ is a reaction which changesa discrete species on the fast time scale, etc. Note that such a reaction can as well changea continuous species on the slow time scale. The separation of fast and slow time scalesis determined by (2.11) with ε = 1 . As in Table 1, we mark the cells which finallydetermine the dynamics of the limiting object Fast time scale
N dt
Reactions ofdiscrete species onfast scale k ∈ K f ◦ , β k = 1 Reactions ofcontinuous specieson fast scale k ∈ K f • , β k > Other reactions k ∈ K − K f ◦ − K f • , β k ≥ Discrete species, α i = 0 , i = I f ◦ ζ fik (cid:26) = 0 , k ∈ K fi , = 0 , else ζ ik = ζ fik = 0 ζ fik = 0Continuous species, α i > , i = I f • ζ ik = 0 or β k − α i − < ⇒ ζ fik = 0 ζ fik (cid:26) = 0 , k ∈ K fi , = 0 , else ζ ik = 0 or β k − α i − > ⇒ ζ fik = 0 Slow time scale dt Reactions of discrete species onslow scale k ∈ K s ◦ , β k = 0 Reactions of continuous species onslow scale k ∈ K s • , β k > Discrete species, α i = 0 , i = I s ◦ ζ sik (cid:26) = 0 , k ∈ K si = 0 , else ζ ik = ζ sik = 0Continuous species, α i > , i = I s • ζ ik = 0 or β k − α i < ⇒ ζ sik = 0 ζ sik (cid:26) = 0 , k ∈ K si , = 0 , else The proof of Lemma 2.4 is an extension of the classical theorem for con-vergence of Markov chains to solutions of ODEs; see Kurtz (1970, 1981), orEthier and Kurtz (1986). It is essentially shown in Kang and Kurtz (2013),Theorem 4.1. For other recent related results, see Franz, Liebscher and Zeiser(2012).
Example 2.5 (Self-regulating gene). We give a simple example of asingle-scale reaction network which leads to a piecewise deterministic solu-tion [similar to Franz, Liebscher and Zeiser (2012), Section 5]. Consider aself-regulating gene modeled by the set of reactions : G + P κ ′ −→ G ′ + P, : G ′ κ ′ −→ G, CALING LIMITS OF SPATIAL COMPARTMENT MODELS : G ′ κ ′ −→ G ′ + P, : P κ ′ −→ ∅ , where G is the inactivated gene, G ′ is the activated gene (hence G, G ′ sumsto 1 and is conserved by the reactions), and P is the protein expressed by thegene. Here, describes activation of the gene by the protein, is spontaneousdeactivation of the gene, is production of the protein by the activated geneand is degradation of the protein. Let x = ( x G , x G ′ , x P ) = (1 − x G ′ , x G ′ , x P )and let the reaction rates beΛ CR ( x ) = κ ′ x G x P = κ ′ (1 − x G ′ ) x P , Λ CR ( x ) = κ ′ x G ′ , Λ CR ( x ) = κ ′ x G ′ , Λ CR ( x ) = κ ′ x P with scaling α G = α G ′ = 0 , α P = 1, that is, I ◦ = { G, G ′ } and I • = { P } , aswell as β = 0 , β = 0 , β = 1 , β = 1 ,κ ′ = N − κ , κ ′ = κ , κ ′ = N κ , κ ′ = κ . Then v G = x G = 1 − v G ′ , v G ′ = x G ′ , v P = N − x P , and see (2.8), λ CR ( v ) = κ v G v P = κ (1 − v G ′ ) v P , λ CR ( v ) = κ v G ′ ,λ CR ( v ) = κ v G ′ , λ CR ( v ) = κ v P . Here, K ∗◦ = K G = K G ′ = { , } and K ∗• = K P = { , } . In this example, thematrices ζ and ζ ∗ are given by ζ = ζ ∗ = GG ′ P − − − . Moreover, according to Lemma 2.4, the limit ( V ( t )) t ≥ of ( V N ( t )) t ≥ solves V G ′ ( t ) = V G ′ (0) + Y (cid:18) κ Z t (1 − V G ′ ( u )) V P ( u ) du (cid:19) − Y (cid:18) κ Z t V G ′ ( u ) du (cid:19) ,V P ( t ) = V P (0) + κ Z t V G ′ ( u ) du − κ Z t V P ( u ) du. P. PFAFFELHUBER AND L. POPOVIC
Multi-scale systems.Two-scale systems.
We say that the chemical network (2.1) is a twoscale system if ( α, β, γ ) from (2.4) are such that: there is a partition of I into (disjoint) I f , called the fast species , and I s , called the slow species ,such that, for some ε > k ∈K i β k + γ = α i + ε, i ∈ I f , (2.11) max k ∈K i β k + γ = α i , i ∈ I s . Without loss of generality, we assume that γ = 0, and that our choice of N is such that ε = 1 in (2.11), so the relative change of fast species happens atrate O ( N ) and the relative change of slow species happens at rate O (1).We first consider what happens on the faster time scale N dt . For each i ∈ I f , let K fi ⊆ K i be the set of reactions with β k = α i + 1. Define K f = { k ∈ K : ∃ i ∈ I f : k ∈ K i , β k = α i + 1 } , (2.12)and a matrix ζ f with |I f | rows and |K f | columns defined by ζ fik = lim N →∞ N − ( α i +1) N β k ζ ik , i ∈ I f , k ∈ K fi . (2.13)This matrix identifies a subnetwork of reactions and their effective changeon the faster time scale N dt . Let I f ◦ ⊆ I f be the subset of fast species forwhich α i = 0, and let K f ◦ = S i ∈I f ◦ K fi be the subset of reactions changingthese fast discrete species on this time scale. Let I f • ⊆ I f be the subset offast species for which α i >
0, which are continuous species on the fast timescale, and let K f • = S i ∈I f • K fi be the subset of reactions changing continuousspecies on this time scale. Since I f ◦ and I f • are disjoint in I f = I f ◦ ∪ I f • , andsince β k is unique for each reaction in K f = K f ◦ ∪ K f • , it follows that K f ◦ and K f • are disjoint as well.We next consider what happens on the slower time scale dt . For each i ∈ I s , let K s = { k ∈ K : ∃ i ∈ I s , β k = α i } , (2.14)be the set of reactions such that β k = α i , and ζ s = ( ζ sik ) i ∈I s ,k ∈K s defined by ζ sik = lim N →∞ N − α i N β k ζ ik , i ∈ I s , k ∈ K si . (2.15)This matrix identifies the subnetwork of reactions and their effective changeon the slower time scale dt . Let I s ◦ ⊆ I s be the subset of (discrete) slowspecies for which α i = 0, and K s ◦ = S i ∈I s ◦ K si the subset of reactions changing CALING LIMITS OF SPATIAL COMPARTMENT MODELS discrete species on this time scale. Let I s • ⊆ I be the subset of (continuous)slow species for which α i >
0, and K s • = S i ∈I s • K si the subset of reactionschanging continuous species on this time scale. As before, I s ◦ and I s • beingdisjoint in I s = I s ◦ ∪ I s • implies that K s ◦ and K s • are disjoint in K s = K s ◦ ∪ K s • as well.Note, however, that there is no reason for K s to be disjoint from K f . Infact, there may be reactions in K s with parameter β k that make an effectivechange on the time scale dt in a slow species of high enough abundance α i = β k , that also effectively change some fast species on the time scale N dt , that is, for some j ∈ I f with β k = α j + 1. The important factor forlimiting results is that we identify contributions from reactions on each ofthe two scales independently, and make assumptions on their stability.Our initial division of species into fast and slow may include some con-served quantities , that is, linear combinations of fast species that remainunchanged on the faster time scale N dt . Let N (( ζ f ) t ) be the null space of( ζ f ) t . If its dimension is > N dt . In spatial systems, the fast speciesare counts of species in a single compartment (which evolves due to both,movement and chemical reactions) while conserved quantities are the sumtotal of the coordinates in all compartments (which evolves only accordingto chemical reactions). For now—unless stated otherwise—we assume thatthe basis for the species is such that dim( N (( ζ f ) t )) = 0.Define the fast process V Nf := ( V Nf ( t )) t ≥ as V Nf ( t ) := ( V Ni ( t )) i ∈I f and the slow process V Ns := ( V Ns ( t )) t ≥ as V Ns ( t ) := ( V Ni ( t )) i ∈I s . We give necessaryassumptions on the dynamics of V Nf on the time scale N dt conditionalon V Ns ( t ) ≡ v s being constant, on the dynamics of V Ns , and on the overallbehavior of V N in order to obtain a proper limiting dynamics of slow species, V Ns . Assumption 2.6 (Dynamics of a two-scale reaction network). Recall λ CR k from (2.8). The two-scale reaction network (2.11) with effective change ζ f as in (2.13) on time scale N dt and ζ s as in (2.15) on time scale dt satisfiesthe following conditions:(i) For each v s ∈ R |I s | + there exists a well-defined process V f | v s , giving thedynamics of the fast species given the vector of slow species, that is,the solution of V f | v s ( t ) = V f | v s (0) + X k ∈K f ◦ ζ f · k Y k (cid:18)Z t λ CR k ( V f | v s ( u ) , v s ) du (cid:19) (2.16) P. PFAFFELHUBER AND L. POPOVIC + X k ∈K f • ζ f · k Z t λ CR k ( V f | v s ( u ) , v s ) du with a unique stationary probability measure µ v s ( dz ) on R |I f | + , suchthat ¯ λ CR k ( v s ) = Z R |I f | + λ CR k ( z, v s ) µ v s ( dz ) < ∞ , k ∈ K s . (2.17)(ii) There exists a well-defined process V s that is the solution of V s ( t ) = V s (0) + X k ∈K s ◦ ζ s · k Y k (cid:18)Z t ¯ λ CR k ( V s ( u )) du (cid:19) (2.18) + X k ∈K s • ζ s · k Z t ¯ λ CR k ( V s ( u )) du with ¯ λ CR k given by (2.17).(iii) There exists a locally bounded function ψ : R |I| + → R , ψ ≥ ψ ( x ) → ∞ as x → ∞ , and(iii-a) for each t > N E (cid:20)Z t ψ ( V N ( u )) du (cid:21) < ∞ ;(iii-b) for all k ∈ K lim K →∞ sup | x | >K λ CR k ( x ) ψ ( x ) = 0 . Lemma 2.7 (Convergence of two-scale reaction networks).
Let V N bethe vector process of rescaled species amounts for the reaction network whichis the unique solution to (2.6) [or (2.7)]. Assume that ( α, β, γ = 0) satisfythe two-scale system assumptions (2.11) for some I f , I s and ε = 1 , and theAssumptions 2.6 are satisfied. Then, if V Ns (0) N →∞ = ⇒ V s (0) , the process ofrescaled amounts of the slow species V Ns ( · ) converges weakly to the solution V s ( · ) of (2.18) with rates given by (2.17) in the Skorokhod topology. The proof of Lemma 2.7 is given in Kang and Kurtz (2013), Theorem 5.1.
Example 2.8 (Production from a fluctuating source). We present herean example of a reaction network with two time scales with no conservedspecies on the fast time scale. In our example, two species A and B react CALING LIMITS OF SPATIAL COMPARTMENT MODELS and produce species C . Source B fluctuates as it is quickly transported intothe system and degrades very fast. We have the set of reactions : A + B κ ′ −→ C, : ∅ κ ′ −→ B, : B κ ′ −→ ∅ . Here, the sum of the numbers of molecules A and C is constant (but bothwill turn out to be slow species), so we only need to consider the dynamics ofthe A molecules. We denote molecules numbers by x A and x B , respectively,set x = ( x A , x B ) and consider the reaction rates as given by mass actionkinetics, Λ CR ( x ) = κ ′ x A x B , Λ CR ( x ) = κ ′ , Λ CR ( x ) = κ ′ x B . (2.19)For the scaled system, we use α A = α C = 1, α B = 0. So, setting the rescaledspecies counts v A = N − x A , v B = x B and β = 0 , β = 1 , β = 1 , (2.20) κ = κ ′ , κ = N − κ ′ , κ = N − κ ′ , (2.21)we write λ CR ( v ) = κ v A v B , λ CR ( v ) = κ , λ CR ( v ) = κ v B . (2.22)Now, the process V N = ( V NA , V NB ) is given by (2.6) as V NA ( t ) = V NA (0) − N − Y (cid:18) N Z t κ V NA ( u ) V NB ( u ) du (cid:19) ,V NB ( t ) = V NB (0) − Y (cid:18) N Z t κ V NA ( u ) V NB ( u ) du (cid:19) + Y ( N κ t )(2.23) − Y (cid:18) N Z t κ V NB ( u ) du (cid:19) . From this representation, it should be clear that V B is fast while V A is aslow species. For γ = 0, ε = 1, we have K s = { } , K f = { , , } (in particular K s ∩ K f = ∅ ), K A = { } , K B = { , , } and I f = I ◦ f = { B } , I s = I • s = { A } .The matrices describing the reaction dynamics on both scales are ζ = AB (cid:18) − − − (cid:19) , ζ f = B ( − − , ζ s = A ( − , where the three columns in ζ and ζ f give reactions , and , and ζ s isa 1 × A is involved. Notethat N (( ζ f ) t ) = { } , indicating that there are no conserved quantities onthe fast time scale. In order to study the dynamics of the slow species, P. PFAFFELHUBER AND L. POPOVIC V s := V A , we apply Lemma 2.7 and have to check Assumption 2.6. Here, forPoisson processes Y , Y and Y , and fixed V s = V A = v A , from (2.16), V B | v A ( t ) − V B | v A (0) = − Y (cid:18)Z t κ v A V B | v A ( u ) du (cid:19) + Y ( κ t ) − Y (cid:18)Z t κ V B | v A ( u ) du (cid:19) d = Y ( κ t ) − Y + (cid:18)Z t ( κ v A + κ ) V B | v A ( u ) du (cid:19) for some Poisson process Y + which is independent of Y . Note that V B | v A ( · )is a birth–death process with constant birth rate κ and linear death rates,proportional to κ v A + κ . It is well known that in equilibrium, V B | v A d = X with X ∼ Poi (cid:18) κ κ + κ v A (cid:19) , which gives the desired µ v A ( dv B ). Hence, (2.17) gives¯ λ CR ( v A ) = E [ κ v A X ] = κ κ v A κ + κ v A . Finally, Lemma 2.7 implies that in the limit N → ∞ , we obtain the dynamics V A ( t ) = V A (0) − Z t ¯ λ CR ( V s ( u )) du = V A (0) − Z t κ κ V A ( u ) κ + κ V A ( u ) du. Three time scales.
Chemical reaction networks with more than two timescales also appear in the literature; see E, Liu and Vanden-Eijnden (2007)for a simulation algorithm for such systems. One example is the heat shockresponse in
Escherichia coli , introduced by Srivastava, Peterson and Bentley(2001) and studied in detail by Kang (2012). Here, we state an extensionof Lemma 2.7 to reaction networks with more than two time scales [seeKang, Kurtz and Popovic (2014)]. Namely, suppose that for some γ ∈ R theparameters α, β in (2.5) and (2.8) are such that: there is a partition of I into disjoint sets I f , I m , I s such that, for some ε > ε > k ∈K i β k + γ = α i + ε , i ∈ I f , max k ∈K i β k + γ = α i + ε , i ∈ I m , (2.24) max k ∈K i β k + γ = α i , i ∈ I s . We will assume, as before, that γ = 0, and that our choice of N is such that ε = 1 in (2.24), so the relative change of fastest species I f happens at rate CALING LIMITS OF SPATIAL COMPARTMENT MODELS O ( N ), the relative change of the middle species I m happens at rate O ( N ε ),0 < ε <
1, and the relative change of slow species I s happens at rate O (1).Again, we need to consider what happens on each single time scale sepa-rately. In addition to earlier definitions, for each i ∈ I m we let K mi ⊆ K bethe set of reactions with β k = α i + ε , K m = S i ∈I m K mi , and a matrix ζ m with |I m | rows and |K m | columns defined by ζ mik = lim N →∞ N − ( α i + ε ) N β k ζ ik , i ∈ I m , k ∈ K mi , (2.25)which identifies a subnetwork of reactions and their effective change on themiddle time scale N ε dt , and we let I m ◦ ⊂ I m be the subset of middle speciesfor which α i = 0, I m • ⊆ I m be the subset of fast species for which α i >
0, andfinally K m ◦ = S i ∈I m ◦ K mi , and K m • = S i ∈I m • K mi . We now need an additionalset of assumptions on the dynamics of V Nf on the time scale N dt conditionalon ( V Nm ( t ) , V Ns ( t )) = ( v m , v s ) being constant, and on the dynamics of V Nm on the time scale N ε dt conditional on V Ns ( t ) = v s being constant. Assumption 2.9 (Dynamics of a three-scale reaction network). Thethree time scale reaction network (2.24) with effective change (2.13) on timescale
N dt , (2.25) on time scale N ε dt and (2.15) on time scale dt satisfiesthe following conditions:(i-a) For each ( v m , v s ) ∈ R |I m | + |I s | + there exists a well-defined process thatis the solution of V f | ( v m ,v s ) ( t ) = V f | ( v m ,v s ) (0) + X k ∈K f ◦ ζ f · k Y k (cid:18)Z t λ CR k ( V f | ( v m ,v s ) ( u ) , v m , v s ) du (cid:19) + X k ∈K f • ζ f · k Z t λ CR k ( V f | ( v m ,v s ) ( u ) , v m , v s ) du with a unique stationary probability measure µ ( v m ,v s ) ( dz ) on R |I f | + , such that˜ λ CR k ( v m , v s ) = Z R |I f | + λ CR k ( z, v m , v s ) µ ( v m ,v s ) ( dz ) < ∞ , k ∈ K m . (2.26)(i-b) For each v s ∈ R |I s | + there exists a well-defined process that is thesolution of V m | v s ( t ) = V m | v s (0) + X k ∈K m ◦ ζ m · k Y k (cid:18)Z t ˜ λ CR k ( V m | v s ( u ) , v s ) du (cid:19) + X k ∈K m • ζ m · k Z t ˜ λ CR k ( V m | v s ( u ) , v s ) du, P. PFAFFELHUBER AND L. POPOVIC which has a unique stationary probability measure µ v s ( dv m ) on R |I m | , suchthat ¯ λ CR k ( v s ) = Z R |I m | + ˜ λ CR k ( v m , v s ) µ v s ( dv m ) < ∞ , k ∈ K s . (2.27)(ii) There exists a well-defined process that is the solution of (2.18) with¯ λ CR k given by (2.27).(iii) see Assumption 2.6(iii).The extension of Lemma 2.7 then becomes the following. Lemma 2.10 (Convergence of three-scale reaction networks).
Let V N bethe vector process of rescaled species amounts for the reaction network whichis the unique solution to (2.6) [or (2.7)]. Assume that ( α, β, γ = 0) satisfy thethree time scale system assumptions (2.24) for some I f , I m , I s and < ε <ε = 1 , and the Assumptions 2.9 are satisfied. Then, if V N (0) N →∞ = ⇒ V (0) , theprocess of rescaled amounts of the slow species V Ns converges weakly to thesolution V s of (2.18) with rates given by (2.27) in the Skorokhod topology. The proof of Lemma 2.10 is along the same lines as the proof of Lemma 2.7,this time applying the stochastic averaging twice; see Kang (2012) for thesame approach.
Conserved quantities.
We turn now to the problem of conserved quanti-ties . Suppose we have a two-scale reaction network with dim( N (( ζ f ) t )) =: n f >
0. Then there exist linearly independent R -valued vectors θ c i = ( θ c i , . . . ,θ c i |I f | ), i = 1 , . . . , n f such that t
7→ h θ c i , V f | v s ( t ) i where V f | v s from (2.16) isconstant. In other words, the change of h θ c i , V Nf ( t ) i on the time scale N dt goes to 0. We set Θ f := ( θ c i ) i =1 ,...,n f , that is, N (( ζ f ) t ) = span(Θ f )(2.28)and note that the construction implies that θ c i has a unique parameter α i associated with it, which we denote by α c i , i = 1 , . . . , | Θ f | (all the species inthe support of θ c i must have the same scaling parameter α c i , as any specieswith a greater value of the scaling parameter does not effectively contributeto the conservation law in the limit).We assume that the effective changes for these combinations are on thetime scale dt , that is, sup k ∈K : h θ ci ,ζ · ,k i6 =0 β k ≤ α c i . In other words, we excludethe possibility that they create a new time scale, or that they effectivelyremain constant as then we do not need to worry about their dynamics.This will be all we need for our main results on the compartment model of CALING LIMITS OF SPATIAL COMPARTMENT MODELS multi-scale reaction networks. If they change on the time scale dt , we needto consider their behavior together with that of the slow species.We let V Nc = ( V Nc i ) i =1 ,..., | Θ f | be the vector of rescaled conserved quantities.For each i = 1 , . . . , | Θ f | , let K θ ci be the set of slow reactions such that β k = α c i and h θ c i , ζ · k i 6 = 0, and let K c := S | Θ f | i =1 K θ ci . Note that K c ∩ K f = ∅ byconstruction. Let ζ c be the matrix with | Θ f | rows and |K c | columns definedby ζ cθ ci k = lim N →∞ N − α ci N β k h θ c i , ζ · k i , i = 1 , . . . , | Θ f | , k ∈ K ci . (2.29)Let Θ f ◦ ⊆ Θ f be the subset of conserved quantities for which α c i = 0, K c ◦ = S θ ci ∈ Θ f ◦ K θ ci . Let Θ f • ⊆ Θ f be the subset of conserved quantities for which α c i > K c • = S θ ci ∈ Θ f • K θ ci . As before, K c ◦ and K c • are disjoint.We extend our results, under obvious modifications of our earlier assump-tions below. Note that the dynamics of conserved quantities depends on thatof the fast species in the same way as the dynamics of the slow species does. Assumption 2.11 (Dynamics of a two-scale reaction network with con-served quantities). The two-scale reaction network (2.11) with effectivechange (2.13) on time scale
N dt and (2.15) and (2.29) on time scale dt satisfies the following conditions:(i) For each ( v s ; v c ) ∈ R |I s | + | Θ f | + , v c := ( v c i ) i =1 ,..., | Θ f | , there exists a well-defined process that is the solution of V f | ( v s ; v c ) ( t ) = V f | ( v s ; v c ) (0)+ X k ∈K f ◦ ζ f · k Y k (cid:18)Z t λ CR k ( V f | ( v s ; v c ) ( u ) , v s ) du (cid:19) + X k ∈K f • ζ f · k Z t λ CR k ( V f | ( v s ; v c ) ( u ) , v s ) du, which satisfies the constraints h θ c i , V f | ( v s ; v c ) i = v c i , θ c i ∈ Θ f , (2.30)and which has a unique stationary probability measure µ ( v s ; v c ) ( dz ) on R |I f | + [concentrated on the linear subspace such that (2.30) is satisfied], such that¯ λ CR k ( v s , v c ) = Z R |I f | + λ CR k ( z, v s , v c ) µ ( v s ; v c ) ( dz ) < ∞ , k ∈ K s . (2.31) P. PFAFFELHUBER AND L. POPOVIC (ii) In addition to a well-defined solution of (2.18), there exists a well-defined process that is the solution of V c ( t ) = V c (0) + X k ∈K c ◦ ζ c · k Y k (cid:18)Z t ¯ λ CR k ( V s ( u ) , V c ( u )) du (cid:19) (2.32) + X k ∈K c • ζ c · k Z t ¯ λ CR k ( V s ( u ) , V c ( u )) du, that is, V c i ( t ) = V c i (0) + X k ∈K c ◦ ζ cθ ci k Y k (cid:18)Z t ¯ λ CR k ( V s ( u ) , V c ( u )) du (cid:19) (2.33) + X k ∈K c • ζ cθ ci k Z t ¯ λ CR k ( V s ( u ) , V c ( u )) du, where the rates in both (2.18) and (2.32), (2.33) are given by (2.31).(iii) Same as in Assumptions 2.6.The following Lemma 2.12 is again in Theorem 5.1 of Kang and Kurtz(2013). Lemma 2.12 (Convergence of two-scale reaction networks with conservedfast quantities).
Let V N be the process of rescaled species amounts (2.6)for a two-scale reaction network, with α, β satisfying (2.11), γ = 0 , ε = 1 ,and with conserved quantities Θ f = ( θ c i ) i =1 ,..., | Θ f | [which is a basis of thenull space of (( ζ ik ) i ∈I f ,k ∈K f ) t ] whose effective change is on time scale dt ,with Assumptions 2.11 satisfied. Then, if V N (0) N →∞ = ⇒ V (0) , we have jointconvergence of the process of rescaled amounts of the slow and conservedquantities ( V Ns ( · ) , V Nc ( · )) N →∞ = ⇒ ( V s ( · ) , V c ( · )) in the Skorohod topology, with V s the solution of (2.18) and V c the solution of (2.32) with rates given by(2.31). It is clear that the result on the limiting dynamics of the conserved quan-tities which change on the time scale dt holds even if we do not have anyslow species on this time scale. We then only have the dynamics of conservedquantities following (2.33) with the rates ¯ λ CR k obtained using the station-ary probability measure for the fast species µ v c ( · ) which depends on theconserved quantities only. Analogously, it is possible that the dynamics ofconserved species on time scale dt is trivial in which case we have the dynam-ics of slow quantities following (2.18) with V c ( u ) = v c , u >
0. Furthermore,
CALING LIMITS OF SPATIAL COMPARTMENT MODELS if we have a reaction network on three scales, it is obvious how to writethe analogous result for the conserved quantities on whichever slower timescale their dynamics is. Both of these situations appear in the dynamics ofcompartment reaction network models, and the above lemmas provide allthe tools we need for our results on models with movement between com-partments. Example 2.13 (Michaelis–Menten kinetics). One of the simplest multi-scale reaction systems with conserved quantities on the fast time scale leadsto the well-known Michaelis–Menten kinetics. A substrate S is transformedinto a product P with the help of an enzyme E via a complex ES formedby enzyme and substrate. The set of reactions is : E + S κ ′ −→ ES, − : ES κ ′− −→ E + S, (2.34) : ES κ ′ −→ E + P. The sum of numbers of free and bound enzymes
E, ES is a constant, whichwill be denoted m (also the sum of numbers S, ES, P molecules is a constantbut
S, P will be more abundant and they will not effectively contribute toa conserved quantity on the fast time scale). We denote molecules numbersby x S , x E , x ES and x P , let x = ( x S , x E , x ES , x P ) and let the reaction ratesbe given by mass action kinetics,Λ CR ( x ) = κ ′ x S x E , Λ CR − ( x ) = κ ′− x ES , Λ CR ( x ) = κ ′ x ES . (2.35)For the scaled system, we use α S = α P = 1 , α E = α ES = 0. Setting therescaled species counts v S = N − x S , v E = x E , v ES = x ES , v P = N − x P , and β = 0 , β − = 1 , β = 1 , (2.36) κ = κ ′ , κ − = N − κ ′− , κ = N − κ ′ , we write λ CR ( v ) = κ v S v E , λ CR − ( v ) = κ − v ES , λ CR ( v ) = κ v ES . Note that the rescaling gives that V NS + V NP = O (1), such that we only needto describe V NS . The process V N = ( V NS , V NE , V NES ) is given by V NS ( t ) = V NS (0) − N − Y (cid:18) N Z t κ V NS ( u ) V NE ( u ) du (cid:19) + N − Y − (cid:18) N Z t κ − V NES ( u ) du (cid:19) , P. PFAFFELHUBER AND L. POPOVIC V NE ( t ) = V NE (0) − Y (cid:18) N Z t κ V NS ( u ) V NE ( u ) du (cid:19) + Y − (cid:18) N Z t κ − V NES ( u ) du (cid:19) + Y (cid:18) N Z t κ V NES ( u ) du (cid:19) ,V NES ( t ) = V NES (0) + Y (cid:18) N Z t κ V NS ( u ) V NE ( u ) du (cid:19) − Y − (cid:18) N Z t κ − V NES ( u ) du (cid:19) − Y (cid:18) N Z t κ V NES ( u ) du (cid:19) . From this, we see V E , V ES are fast while V S , V P are slow species. For γ = 0, ε = 1, we have K f = K s = { , − , } , K S = { , − } , K E = K ES = { , − , } , K P = { } and I f = I ◦ f = { E, ES } , I s = I • s = { S, P } . The matrices describingthe reaction dynamics on both scales are ζ = SEESP − − − −
10 0 1 , ζ f = EES (cid:18) − − − (cid:19) , (2.37) ζ s = SP (cid:18) − (cid:19) , where the three columns give reactions , − and . Note that N (( ζ f ) t ) =span((1 , V c := V E + V ES is constant (at least on thefast time scale). In order to study the dynamics of the slow species, V s :=( V S , V P ), we have to check Assumption 2.11 and apply Lemma 2.12. Here,for Poisson processes Y and Y − , and fixed V s = ( V S , V P ) = ( v S , v P ) = v s and V c = V E + V ES =: m , from (2.16), V E | ( v s ,v c ) ( t ) − V E | ( v s ,v c ) (0)= − Y (cid:18)Z t κ v S V E ( u ) du (cid:19) + Y − (cid:18)Z t κ − V ES ( u ) du (cid:19) + Y (cid:18)Z t κ V ES ( u ) du (cid:19) d = − Y (cid:18)Z t κ v S V E ( u ) du (cid:19) + Y − + (cid:18)Z t ( κ − + κ ) V ES ( u ) du (cid:19) ,V ES | ( v s ,v c ) ( t ) − V ES | ( v s ,v c ) (0)= Y (cid:18)Z t κ v S V E ( u ) du (cid:19) + Y − + (cid:18)Z t ( κ − + κ ) V ES ( u ) du (cid:19) CALING LIMITS OF SPATIAL COMPARTMENT MODELS for some Poisson process Y − + which is independent of Y . Note that( V E | ( v s ,v c ) ( · ) , V ES | ( v s ,v c ) ( · )) behaves like an Ehrenfest urn with two compart-ments, where each E turns to ES at rate κ v S , and each ES turns to E atrate ( κ − + κ ). It is well known that ( V E | ( v s ,v c ) , V ES | ( v s ,v c ) ) d = ( X, m − X )has an equilibrium X ∼ Binom (cid:18) m, κ − + κ κ − + κ + κ v S (cid:19) , which gives the desired µ ( v s ,m ) ( dv E , dv ES ) concentrated on V c = V E + V ES = m . The conserved species V c do not change even on the time scale dt (thiswill no longer be the case in a heterogeneous compartment model in thenext section). Hence, (2.31) gives¯ λ ( v s ) = E [ κ v S X ] = κ v S m ( κ − + κ ) κ − + κ + κ v S , ¯ λ − ( v s ) = E [ κ − ( m − X )] = κ − mκ v S κ − + κ + κ v S , ¯ λ ( v s ) = E [ κ ( m − X )] = κ mκ v S κ − + κ + κ v S . Lemma 2.12 implies that in the limit N → ∞ , we obtain the dynamics V S ( t ) = V S (0) − Z t ¯ λ ( V s ( u )) du + Z t ¯ λ − ( V s ( u )) du (2.38) = V S (0) − Z t mκ κ V S ( u ) κ − + κ + κ V S ( u ) du, for V S , which is the classical Michaelis–Menten kinetics.
3. Chemical reactions in multiple compartments.
We now assume thatthe chemical system is separated into a set of D compartments , and chemicalspecies can migrate within these compartments. For species i ∈ I , movement happens from compartment d ′ to d ′′ at rate Λ M i,d ′ ,d ′′ .3.1. The Markov chain model.
Denoting by X id ( t ) the number of mol-ecules of species i in compartment d at time t , we assume that ( X ( t )) t ≥ with X ( t ) = ( X id ( t )) i ∈I ,d ∈D is solution of X id ( t ) = X id (0) + X k ∈K ζ ik Y kd (cid:18)Z t Λ CR kd ( X · d ( u )) du (cid:19) (3.1) + X d ′ ,d ′′ ∈D ( δ d ′′ ( d ) − δ d ′ ( d )) Y i,d ′ ,d ′′ (cid:18)Z t Λ M i,d ′ ,d ′′ X id ′ ( u ) du (cid:19) , P. PFAFFELHUBER AND L. POPOVIC where δ d ( · ) is a Dirac delta function, X · d = ( X id ) i ∈I and all the Y · ’s areindependent (rate 1) Poisson processes. We assume the following. Assumption 3.1 (Dynamics of un-scaled multi-compartment reactionnetwork). The reaction network dynamics satisfies the following conditions:(i) Same as Assumption 2.1(i) in each compartment and for all k thereis at least one d with Λ CR kd = 0.(ii) Given ( Y kd ) k ∈K ,d ∈D , and ( Y i,d ′ ,d ′′ ) i ∈I ,d ′ ,d ′′ ∈D , the time change equa-tion (3.1) has a unique solution.For X i ( t ) := P d ∈D X id ( t ), X i ( t ) d = X i (0) + X k ∈K ζ ik Y k (cid:18)Z t X d ∈D Λ CR kd ( X · d ( u )) du (cid:19) , for some independent (rate 1) Poisson processes ( Y k ) k ∈K . However, sincethe rate P d ∈D Λ CR kd ( X · d ( s )) depends on all entries in X ( s ), the process(( X i ( t )) i ∈I ) t ≥ is not in general Markov.3.2. The rescaled system.
Consider the solution of (3.1) with the chem-ical reaction rates Λ CR kd and movement rates Λ M i,d ′ ,d ′′ replaced by Λ CR ,Nkd andΛ M ,Ni,d ′ ,d ′′ , respectively. For real-valued α = ( α i ) i ∈I , β = ( β k ) k ∈K , γ, η = ( η i ) i ∈I , with α i ≥ , i ∈ I , we denote the ( α, β, γ, η )-rescaled system by V Nid ( t ) := N − α i X Nid ( N γ t ) , i ∈ I , d ∈ D ,λ CR ,Nk ( v ) := N − β k Λ CR ,Nk (( N α i v i ) i ∈I ) , k ∈ K ,λ M ,Ni,d ′ ,d ′′ ( v ) := N − η i Λ M ,Ni,d ′ ,d ′′ (( N α i v i ) i ∈I ) , i ∈ I , d ′ , d ′′ ∈ D , where α, β, γ, η is chosen so that V Nid = O (1) , λ CR ,Nk = O (1) , λ M ,Ni,d ′ ,d ′′ = O (1)(reactions of the same type and species of the same type are scaled by thesame parameters in each compartment). Again, we will restrict to the case γ = 0. Assumption 3.2 (Dynamics of scaled multiple compartment reaction net-work). In addition to Assumption 2.2 within each compartment, there exist λ M i,d ′ ,d ′′ , i ∈ I , d ′ , d ′′ ∈ D with N − η i Λ M ,Ni,d ′ ,d ′′ N →∞ −→ λ M i,d ′ ,d ′′ . (3.2)Again, we will assume that this convergence is actually an identity. CALING LIMITS OF SPATIAL COMPARTMENT MODELS The ( α, β, γ, η )-rescaled system V N ( t ) = N − α X ( N γ t ) is the unique solu-tion to the system of stochastic equations V Nid ( t ) = V Nid (0) + X k ∈K N − α i ζ ik Y kd (cid:18) N β k + γ Z t λ CR kd ( V N · d ( u )) du (cid:19) (3.3) + X d ′ ,d ′′ ∈D N − α i ( δ d ′′ ( d ) − δ d ′ ( d )) Y i,d ′ ,d ′′ (cid:18) N α i + η i + γ Z t λ M i,d ′ ,d ′′ V Nid ′ ( s ) (cid:19) . In addition, define S N = ( S Ni ) i ∈I with S Ni := X d ∈D V Nid , (3.4)then S Ni solves S Ni ( t ) = S Ni (0) + X k ∈K N − α i ζ ik X d ∈D Y kd (cid:18) N β k + γ Z t λ CR kd ( V N · d ( u )) du (cid:19) . (3.5) Remark 3.3 (Heterogeneity of the reaction network). Our set-up doesnot preclude the option that different compartments may have different re-action networks all together. We contain all possible reactions in the stoi-chiometric matrix ζ , then setting individual compartment rates λ CR kd to zeroin desired compartments can achieve this.3.3. Spatial single-scale systems.
We can now examine the effect of het-erogeneity on the chemical reaction systems via compartmental models. Weassume that (2.9) holds within every compartment. The sets I , I ◦ , I • and K , K ∗◦ , K ∗• and ζ ∗ are used as in Section 2. We assume that movement ofspecies is fast, η i > , i ∈ I and it has a unique equilibrium. Assumption 3.4 (Equilibrium for movement). For each species i ∈ I ,the movement Markov chain, given through the jump rates λ M i,d ′ ,d ′′ from d ′ to d ′′ , has a unique stationary probability distribution denoted by ( π i ( d )) d ∈D . Lemma 3.5 (Movement equilibrium).
Let Assumption 3.4 hold. (1)
Let i ∈ I be such that α i = 0 (i.e., i ∈ I ◦ ). Consider the Markov chainof only the movement of molecules of species i , that is, the solution of V id ( t ) = V id (0) + X d ′ ,d ′′ ∈D ( δ d ′′ ( d ) − δ d ′ ( d )) Y i,d ′ ,d ′′ (cid:18)Z t λ M i,d ′ ,d ′′ V id ′ ( u ) du (cid:19) started with P d ∈D V id (0) = s i . Then, the unique equilibrium probability dis-tribution of this Markov chain is given as the multinomial distribution withparameters ( s i ; ( π i ( d )) d ∈D ) . P. PFAFFELHUBER AND L. POPOVIC (2)
Let i ∈ I be such that α i > (i.e., i ∈ I • ). Consider the limitingdeterministic process of only the movement of molecules of species i , that is,the solution of V id ( t ) = V id (0) + X d ′ ∈D Z t ( λ M i,d ′ ,d V id ′ ( u ) − λ M i,d,d ′ V id ( u )) du started with P d ∈D V id (0) = s i . Then the unique equilibrium of this processis given by ( s i π i ( d )) d ∈D .We denote the equilibrium probability distribution of movement of all species,started in ( s i ) i ∈I by P s and by E s the corresponding expectation operator.From the above, P s is a product of multinomial and point mass distributions. Proof.
In (1), we have an Ehrehfest urn model with |D| boxes; due itsreversibility it is easy to check we have the correct equilibrium. In (2), wehave a deterministic system of |D| equations whose equilibrium is equallyeasy to obtain. (cid:3)
We start with the simplest results for chemical reaction networks whichare on a single scale, and describe the effect of mixing on the heterogeneouschemical reaction system.
Assumption 3.6 (Dynamics of the spatial single-scale reaction network).The spatial single-scale reaction network on time scale dt , where Assump-tion 3.4 holds, satisfies the following conditions:(i) Given ( Y k ) k ∈K ∗◦ , the time change equation S ( t ) = S (0) + X k ∈K ∗◦ ζ ∗· k Y k (cid:18)Z t ¯ λ CR k ( S ( u )) du (cid:19) (3.6) + X k ∈K ∗• ζ ∗· k Z t ¯ λ CR k ( S ( u )) du has a unique solution S := ( S ( t )) t ≥ , where¯ λ CR k ( s ) := E s (cid:20)X d ∈D λ CR kd ( V · d ) (cid:21) . (3.7)(ii) Same as (iii) in Assumption 2.6 for all d ∈ D . Theorem 3.7 (Heterogeneous single-scale system).
Let V N be the vec-tor process of rescaled species amounts for the reaction network which isthe unique solution to (3.3). Assume that ( α, β, γ = 0) satisfy single scale CALING LIMITS OF SPATIAL COMPARTMENT MODELS assumption (2.9) within compartments and η i = η > , i ∈ I . Let S N ( t ) =( S Ni ( t )) i ∈I be S Ni ( t ) := P d ∈D V Nid ( t ) , and suppose Assumptions 3.4 and 3.6for the rescaled network hold. If S N (0) N →∞ = ⇒ S (0) , then the process of rescaledsums S N ( · ) converges weakly to the unique solution S ( · ) of (3.6) in the Sko-rohod topology. Proof.
In the heterogeneous reaction network, we have |I|× |D| species;one for each type and each compartment, with rescaled amounts V Nid . Move-ment between compartments can be viewed as (at most) |D| × |D| first-orderreactions involving only species of the same type i ∈ I in different compart-ments, with net change in compartment d of ( δ d ′′ ( d ) − δ d ′ ( d )) ( d ′ ,d ′′ ) ∈D×D atrate { Λ M i,d ′ ,d ′′ , d ′ , d ′′ ∈ D} . This set of reactions together with the original re-actions within each compartment give an overall network in which all thespecies V Nid with ( i, d ) ∈ I × D are fast, whose conserved quantities is a vec-tor of sums over all the compartments for each species type, which are givenby S Ni := P d ∈D V Nid , i ∈ I . Since η i > , i ∈ I the movement reactions changeall the species amounts on the time scale N η dt , and its effective changes onthis time scale are still ( δ d ′′ ( d ′ ) − δ d ′ ( d ′′ )) ( d ′ ,d ′′ ) ∈D×D while the original withincompartment reactions effectively change only the conserved sum quantitieson the time scale dt and its effective changes on this time scale are given by ζ ∗ .In order to apply Lemma 2.12, set ε := η and we need to check Assump-tions 2.11. In this special case, there are no slow species only fast speciesand conserved quantities. Condition (i) is simply the requirement that—inthe limit N → ∞ —for fixed given vector of sums of species movement leadsto a well-defined process on the species amounts in different compartments,which for each value of the vector of sums s has a unique stationary prob-ability measure P s , which is concentrated on P d ∈D v id = s i . This is exactlyimplied by Lemma 3.5 under Assumption 3.4. Conditions (ii) and (iii) inAssumptions 2.11 is assumed in the statement of the theorem.Let us consider the dynamics of the conserved quantities. Here, θ c i =(1 j = i ) j ∈I ,d ∈D is the i th conserved quantity. On the time scale dt , the reactiondynamics of these conserved sums is a Markov chain whose effective change isgiven by the matrix ζ c = ζ ∗ with overall rate equal to a sum of the individualcompartment rates.Since the equilibrium for the movement dynamics P s is given by P s ( dv ) = Y i ∈I ◦ (cid:18) s i v i · · · v i |D| (cid:19) π i (1) v i · · · π i ( |D| ) v i |D| (3.8) × Y i ∈I • δ π i (1) s i ( dv i ) · · · δ π i ( D ) s i ( dv i |D| ) , P. PFAFFELHUBER AND L. POPOVIC the averaged rates for reaction dynamics in each compartment under theequilibrium probability measure as considered in (3.7) are exactly of theform (2.31),¯ λ CR k ( s ) = X v · : P d v d = s · · · X v |I|· : P d v |I| d = s |I| X d ∈D λ CR kd ( v · d ) × Y i ∈I ◦ (cid:18) s i v i · · · v i |D| (cid:19) π i (1) v i · · · π i ( |D| ) v i |D| (3.9) × Y i ∈I • δ π i (1) s i ( v i ) · · · δ π i ( |D| ) s i ( v i |D| )= X d ∈D E s [ λ CR kd ( V · d )] . (cid:3) Corollary 3.8 (Mass-action kinetics).
Let α, β, γ, η be as in Theo-rem 3.7. If the reaction rates are given by mass action kinetics for some κ kd , k ∈ K , d ∈ D λ CR kd ( v · d ) = κ kd Y i ∈I ◦ ν ik ! (cid:18) v id ν ik (cid:19) · Y i ∈I • v ν ik id , (3.10) then the limit of S N ( · ) in the Skorohod topology is the solution of (3.6) withrates given by ¯ λ CR k ( s ) = X d ∈D κ kd Y i ∈I ◦ ν ik ! (cid:18) s i ν ik (cid:19) π i ( d ) ν ik · Y i ∈I • ( π i ( d ) s i ) ν ik . If α i = 0 for all i ∈ I , the limit process for the sums is a Markov chain modelfor reaction networks with mass action kinetics (2.3) whose rate parametersare ¯ κ k = X d ∈D κ kd Y i ∈I π i ( d ) ν ik . (3.11) If α i > for all i ∈ I , the limit process for the sums is the deterministicsolution to an ordinary differential equation dS ( t ) = X k ∈K ζ ∗· k ¯ λ CR k ( S ( t )) dt with mass action kinetics (2.3) whose rate parameters are (3.11). Remark 3.9 (Different time scales for the movement). From the point ofview of the limit on time scale dt , the parameters for time scale of movementof different species types do not have to all be equal η i = η ; as long as η i > i ∈ I , it is easy to show that the limit dynamics of S N ( · ) is as above. CALING LIMITS OF SPATIAL COMPARTMENT MODELS Proof of Corollary 3.8.
We plug (3.10) into (3.9). This gives¯ λ CR k ( s ) = X x · : P d x d = s if α =0 · · · X x |I|· : P d x |I| d = s |I| if α |I| =0 X d ∈D κ kd × Y i ∈I ◦ ν ik ! (cid:18) x id ν ik (cid:19) (cid:18) s i x i · · · x i |D| (cid:19) π i (1) x i · · · π i ( |D| ) x i |D| × Y i ∈I • ( π i ( d ) s i ) ν ik = X d ∈D κ kd X x d =0 ,...,s if α =0 · · · X x |I| d =0 ,...,s |I| if α |I| =0 × Y i ∈I ◦ ν ik ! (cid:18) s i x id (cid:19) (cid:18) x id ν ik (cid:19) π i ( d ) x id (1 − π i ( d )) s i − x id · Y i ∈I • ( π i ( d ) s i ) ν ik = X d ∈D κ kd Y i ∈I ◦ ν ik ! (cid:18) s i ν ik (cid:19) π i ( d ) ν ik X x d =0 ,...,s if α =0 · · · X x |I| d =0 ,...,s |I| if α |I| =0 × Y i ∈I ◦ (cid:18) s i − ν ik x id − ν ik (cid:19) π i ( d ) x id − ν ik (1 − π i ( d )) s i − x id · Y i ∈I • ( π i ( d ) s i ) ν ik = X d ∈D κ kd Y i ∈I ◦ ν ik ! (cid:18) s i ν ik (cid:19) π i ( d ) ν ik · Y i ∈I • ( π i ( d ) s i ) ν ik . When α i = 0 for all i ∈ I only the first sum in (3.6) exists, whereas when α i > i ∈ I only the second sum in (3.6) exists. (cid:3) Example 3.10 (Self-regulating gene in multiple compartments). Weplace the reaction kinetics from Example 2.5 in a spatial multi-compartmentsetting. Let the dynamics initiate with S G (0) + S G ′ (0) = 1 active and inac-tive genes and S P (0) = s P proteins in the whole space. If η G , η G ′ , η P >
0, themovement is faster than any effective reaction dynamics, and the limitingdynamics of the rescaled sums in the whole system solves S G ′ ( t ) = S G ′ (0) + Y (cid:18) ¯ κ Z t (1 − S G ′ ( u )) S P ( u ) du (cid:19) − Y (cid:18) ¯ κ Z t S G ′ ( u ) du (cid:19) ,S P ( t ) = S P (0) + ¯ κ Z t S G ′ ( u ) du − ¯ κ Z t S P ( u ) du, P. PFAFFELHUBER AND L. POPOVIC where ¯ κ k are given by (3.11) and π G , π G ′ , π P are the equilibrium distribu-tions of the movement of G, G ′ and P , respectively. Given the values ofsystem sums S G ′ ( t ) , S P ( t ) the molecules of G ′ , P will then be distributed incompartments according to V G ′ d ( t ) ∼ Multinom( S G ′ ( t ) , π G ′ ( d )) , V P d ( t ) ∼ δ S P ( t ) π P ( d ) . Spatial multi-scale systems.
We next consider heterogeneous reac-tion networks on multiple time scales, with interplay between time scales onwhich the reaction network dynamics evolves and time scales on which thespecies move between compartments. We give results for chemical reactionson two time scales, extensions to more are obvious.We stick to our notation from Section 2.3. In particular, we assume the re-action dynamics (within each compartment) has a separation of time scales(2.11) with ε = 1 , γ = 0. We set K f and K s as in (2.12) and (2.14), respec-tively, and I f and I s for the sets of fast and slow species, if only chemicalreactions within compartments are considered. The scaling parameters formovement of all fast species is η i = η f for i ∈ I f while for all slow species is η i = η s , i ∈ I s . We assume both η f , η s >
0. In order to assess the interplay ofdynamics on different time scales, we need to consider all possible orderingsof ε = 1 , η f and η s . In the sequel, we assume that η f , η s = 1 for simplic-ity. Moreover, the cases η s ≤ η f < η f < η s <
1, as well as 1 < η s ≤ η f and 1 < η f < η s lead to the same limiting behavior, because the movementprocesses occurring on the time scale N η f dt and N η s dt are independent (amovement of one species type on a time scale that is different from that ofreactions depends only on its own molecular counts and is independent ofother species types). Therefore, we are left with the four cases(1) 1 < η s , η f ; (2) η s < < η f ;(3.12) (3) η f < < η s ; (4) η f , η s < . As in the nonspatial situation, we also need to distinguish the cases when(i) there are no conserved quantities on the time scale of fast species [mean-ing that N (( ζ f ) t ) is the null space], and (ii) when some quantities areconserved [i.e., N (( ζ f ) t ) = span(Θ f ), where Θ f = ( θ c i ) i =1 ,...,n f is a linearlyindependent family of R |I f | -valued vectors]. In the latter case, the quantities h θ c i , ( V Nid ( · )) i ∈I f i also change on the time scale N η f dt for d ∈ D by move-ment of the fast species, but h θ c i , ( S Ni ( · )) i ∈I f i is constant on the time scale N η f dt . We start with the case of N (( ζ f ) t ) = null space. No conserved quantities on the fast time scale.
We need to consider dif-ferent processes of possible effective reaction dynamics for fast species and
CALING LIMITS OF SPATIAL COMPARTMENT MODELS their sums, conditional on knowing the values of the slow species. In each ofthe four cases above we need to consider different intermediate processes andassumptions on them. We write here, distinguishing fast and slow species, v = ( v f , v s ) with v f = ( v id ) i ∈I f ,d ∈D , v s = ( v id ) i ∈I s ,d ∈D , as well as s = ( s f , s s ), s f = ( s i ) i ∈I f , s s = ( s i ) i ∈I s . Assumption 3.11 (Dynamics of the spatial multi-scale reaction network).In each case (1)–(4), the spatial two-scale reaction network on time scale
N dt , where Assumption 3.4 holds, satisfies the following conditions:(i) (1) Given ( Y k ) k ∈K f ◦ , the time-change equation of the dynamics of S f given the value of S s = s s S f | s s ( t ) = S f | s s (0) + X k ∈K f ◦ ζ f · k Y k (cid:18)Z t ˜ λ CR(1) k ( S f | s s ( u ) , s s ) du (cid:19) (3.13) + X k ∈K f • ζ f · k Z t ˜ λ CR(1) k ( S f | s s ( u ) , s s ) du has a unique solution, where for all s f , s s ˜ λ CR(1) k ( s f , s s ) = Z R |I f |×|D| + × R |I s |×|D| + X d ∈D λ CR kd ( v · d,f , v · d,s ) P ( s f ,s s ) ( dv f , dv s )(3.14) < ∞ , where v · d,f = ( v id ) i ∈I f , v · d,s = ( v id ) i ∈I s , for P ( s f ,s s ) a product of multinomialand point mass probability distributions for both v f and v s defined in (3.8).In addition, S f | s s ( · ) has a unique stationary probability measure µ s s ( ds f )on R |I f | + .(2) Given ( Y k ) k ∈K f ◦ , the time-change equation of the dynamics of S f giventhe value of V s = v s S f | v s ( t ) = S f | v s (0) + X k ∈K f ◦ ζ f · k Y k (cid:18)Z t ˜ λ CR(2) k ( S f | v s ( u ) , v s ) du (cid:19) (3.15) + X k ∈K f • ζ f · k Z t ˜ λ CR(2) k ( S f | v s ( u ) , v s ) du has a unique solution, where for all s f , v s ˜ λ CR(2) k ( s f , v s ) = Z R |I f |×|D| + X d ∈D λ CR kd ( v · d,f , v · d,s ) P s f ( dv f ) < ∞ (3.16) P. PFAFFELHUBER AND L. POPOVIC for P s f a product of multinomial and point mass probability distributionsas in (3.8), where I is replaced by I f , s by s f and v by v f . In addition, S f | v s ( · ) has a unique stationary probability measure µ v s ( ds f ) on R |I f | + .(3) Given ( Y kd ) k ∈K f ◦ ,d ∈D , the time-change equation of the dynamics of V f given the value of S s = s s V · d,f | s s ( t ) = V · d,f | s s (0) + X k ∈K f ◦ ζ f · k Y kd (cid:18)Z t ˜ λ CR(3) kd ( V · d,f | s s ( u ) , s s ) du (cid:19) (3.17) + X k ∈K f • ζ fik Z t ˜ λ CR(3) kd ( V · d,f | s s ( u ) , s s ) du has a unique solution, where for all v f , s s ˜ λ CR(3) kd ( v · d,f , s s ) = Z R |I s |×|D| + λ CR kd ( v · d,f , v · d,s ) P s s ( dv s ) < ∞ (3.18)for P s s a product of multinomial and point mass probability distributionsas in (3.8), where I is replaced by I s , s by s s and v by v s . In addition, V f | s s ( · ) = ( V id,f | s s ( · )) i ∈I f ,d ∈D has a unique stationary probability measure µ s s ( dv f ) on R |I f |×|D| + .(4) Given ( Y kd ) k ∈K f ◦ ,d ∈D , the time-change equation of the dynamics of V f given the value of V s = v s V · d,f | v s ( t ) = V · d,f | v s (0)+ X k ∈K f ◦ ζ f · k Y kd (cid:18)Z t ˜ λ CR(4) kd ( V · d,f | v s ( u ) , v · d,s ) du (cid:19) (3.19) + X k ∈K f • ζ fik Z t ˜ λ CR(4) kd ( V · d,f | v s ( u ) , v · d,s ) du has a unique solution with unique stationary probability measure µ v s ( dv f )on R |I f |×|D| + . Here, we set ˜ λ CR(4) kd := λ CR kd . (3.20)(ii) There exists a well-defined process S s ( · ) that is the unique solutionof S s ( t ) = S s (0) + X k ∈K s ◦ ζ s · k Y k (cid:18)Z t ¯ λ CR( ℓ ) k ( S s ( u )) du (cid:19) CALING LIMITS OF SPATIAL COMPARTMENT MODELS (3.21) + X k ∈K s • ζ s · k Z t ¯ λ CR( ℓ ) k ( S s ( u )) du, where rates (¯ λ CR( ℓ ) k ) ℓ =1 , , , are given from (˜ λ CR( ℓ ) k ) ℓ =1 , , , in each case as¯ λ CR(1) k ( s s ) = Z R |I f | + ˜ λ CR(1) k ( s f , s s ) µ s s ( ds f )(3.22) = X d ∈D Z Z λ CR kd ( v · d,f , v · d,s ) P ( s f ,s s ) ( dv f , dv s ) µ s s ( ds f ) < ∞ ;¯ λ CR(2) k ( s s ) = Z R |I s |×|D| + Z R |I f | + ˜ λ CR(2) k ( s f , v s ) µ v s ( ds f ) P s s ( dv s )(3.23) = X d ∈D Z Z Z λ CR kd ( v · d,f , v · d,s ) P s f ( dv f ) µ v s ( ds f ) P s s ( dv s ) < ∞ ;¯ λ CR(3) k ( s s ) = X d ∈D Z R |I f |×|D| + ˜ λ CR(3) kd ( v f , s s ) µ s s ( dv f )(3.24) = X d ∈D Z Z λ CR kd ( v · d,f , v · d,s ) P s s ( dv s ) µ s s ( dv f ) < ∞ ;¯ λ CR(4) k ( s s ) = X d ∈D Z R |I s |×|D| + Z R |I f |×|D| + ˜ λ CR(4) k ( v f , v s ) µ v s ( dv f ) P s s ( dv s )(3.25) = X d ∈D Z Z λ CR kd ( v · d,f , v · d,s ) µ v s ( dv f ) P s s ( dv s ) < ∞ . (iii) Same as (iii) in Assumption 2.6 in each compartment. Remark 3.12 (Equivalent formulation). For the dynamics under theabove assumption, the following is immediate: In each case (1)–(4), the spa-tial two-scale reaction network on time scale dt , where Assumption 3.11holds, satisfies the following condition: given ( Y k ) k ∈K s ◦ , the time-changeequation (3.21) has a unique solution, where for all s s ¯ λ CR( ℓ ) k ( s s ) := E s s (cid:20)X d ∈D λ CR kd ( V · d ) (cid:21) < ∞ (3.26)and the distribution of ( V id ) i ∈I ,d ∈D in (3.26) depends on the parameters η s , η f as follows:( ℓ ) = (1) P s s ( dv f , dv s ) = Z R |I f | + µ s s ( ds f ) P ( s f ,s s ) ( dv f , dv s ) , (3.27) P. PFAFFELHUBER AND L. POPOVIC ( ℓ ) = (2) P s s ( dv f , dv s ) = P s s ( dv s ) Z R |I f | + µ v s ( ds f ) P s f ( dv f ) , (3.28) ( ℓ ) = (3) P s s ( dv f , dv s ) = P s s ( dv s ) µ s s ( dv f ) , (3.29) ( ℓ ) = (4) P s s ( dv f , dv s ) = P s s ( dv s ) µ v s ( dv f ) . (3.30)We can now state our results for the limiting behavior of S Nf := ( S Ni , i ∈I f ) and S Ns := ( S Ni , i ∈ I s ) on the time scales N dt and dt . Theorem 3.13 (Two-scale system without conserved fast quantities).
Let V N be the vector process of rescaled species amounts for the reactionnetwork which is the unique solution to (3.3). Assume that ( α, β, γ = 0) satisfy two-scale system assumption (2.11) for some I f , I s with ε = 1 and N (( ζ f ) t ) = 0 [with ζ f from (2.13)] within compartments without conservedquantities on the fast time scale. In addition, η i = η f > for all i ∈ I f and η i = η s > for all i ∈ I s , one of the cases (1)–(4) holds, and Assumption 3.11holds. Then, if S N (0) N →∞ = ⇒ S (0) , the rescaled sums of slow species S Ns ( · ) from(3.4) converges weakly to the unique solution S s ( · ) of (3.21) in the Skorohodtopology. Remark 3.14 (Interpretation). The rates in Theorem 3.13 have an in-tuitive interpretation. In order to compute E s [ λ CR kd ( V · d )], we have to knowthe distribution of V given S s . Consider case (1) as an example: since move-ment of particles are the fastest reactions in the system, given the valueof S s = s s , (i) V s are distributed according to P s s ( dv s ) from (3.8), and(ii) S f is distributed according to the probability measure µ s s ( ds f ) fromAssumption 3.11(i)(1); then, given the value of S f = s f , the values of V f are distributed according to P s f ( dv f ) from (3.8). In case (3), fast reactionswithin compartments intertwine with movement between them: given thevalue of S s = s s , (i) V s are again distributed according to P s s from (3.8),but (ii) V f are distributed according to µ s s ( dv f ). Proof of Theorem 3.13.
The proof relies on use of Lemmas 2.10 and2.12. Let us first consider the case (1): η f , η s >
1. In this case on the twofastest time scales N η f dt and N η s dt , we have movement of fast and slowspecies, respectively, whose sums are unchanged on any time scale fasterthan N dt . We view these two fastest time scales as one, since dynamics ofmovement of fast and slow species are in this case independent of each otherand can be combined. Regarding all of the movement as a set of first-orderreactions as in proof of Theorem 3.7, we have a three time scale dynamics:
CALING LIMITS OF SPATIAL COMPARTMENT MODELS movement of all species is the fast process on time scales N η f dt, N η s dt ,effective change of fast species is the medium process on time scale N dt , andeffective change of slow species is the slow process on the time scale dt . Thefast process of movement of all species has a stationary probability measurethat is a product of multinomial and point mass probability distributions P ( s f ,s s ) from (3.8). Arguments from Lemma 2.12 imply that on the timescale N dt all rates for the reaction network dynamics ˜ λ CR(1) k are sums overcompartments of rates averaged with respect to P ( s f ,s s ) as in (3.14), and themedium process of the sums of fast species S f has effective change given by ζ f . Condition (i)(1) of Assumption 3.11 ensures that on time scale N dt themedium process S f ( · ) is well defined and has a unique stationary probabilitydistribution µ s s ( ds f ). Condition (ii) of Assumption 3.11 then ensures that, inaddition to conditions (i-a) and (i-b), also condition (ii) in the assumptionsfor Lemma 2.10 is met, and consequently the limiting dynamics of the slowprocess S s ( · ) with effective change given by ζ f is well defined and given bythe solution of (3.21) with rates as in (3.22).Next, consider the case (2): η f > , η s <
1. In this case we have a four timescale dynamics: movement of fast species is the fast process on time scale N η f dt , effective change of fast species is the medium-fast process on timescale N dt , movement of slow species is the medium-slow process on timescale N η s dt , and finally effective change of slow species is the slow processon time scale dt . The fast process on time scale N η f dt of movement of allspecies has a stationary probability measure that is P s f (3.8) (over i ∈ I f only). Lemma 2.12 implies that on the next time scale N dt rates ˜ λ CR(2) k are averaged with respect to P s f as in (3.16), and the medium-fast process S f has an effective change given by ζ f . We now have that this process iswell defined and has a unique stationary probability distribution µ v s ( ds f ).Furthermore, on the next time scale N η s dt we only have the movement ofslow species, which has a stationary probability measure that is P s s from(3.8) (over i ∈ I s only). Finally, the limiting dynamics of the slow process S s ( · ) on time scale dt , by an extension of Lemma 2.10 to four time scales,is well defined and given by the solution of (3.21) with rates as in (3.23).Let us next consider the case (3): η f < , η s >
1. We again have a four timescale dynamics: movement of slow species is the fast process on time scale N η s dt , effective change of fast species is the medium-fast process on timescale N dt , movement of fast species is the medium-slow process on timescale N η f dt , and finally effective change of slow species is the slow processon time scale dt . Lemma 2.12 implies that on the medium-fast time scale N dt rates ˜ λ CR(3) kd are averaged with respect to P s s (3.8) (over i ∈ I s only)as in (3.18), and the medium-fast process V f has an effective change given P. PFAFFELHUBER AND L. POPOVIC by ζ f in each compartment d ∈ D . This process is well defined and has aunique stationary probability distribution µ s s ( dv f ). On the next time scale N η s dt we only have the movement of slow species, which has a stationaryprobability measure that is P s f (3.8) (over i ∈ I f only). Finally, on time scale dt , Lemma 2.10 extended to four time-scales implies the limiting dynamicsof the slow process S s ( · ) is well defined and given by the solution of (3.21)with rates as in (3.24).Finally, we consider case (4): η f < , η s <
1. On the fast time scale
N dt in each compartment d ∈ D , independently we have reaction dynamics ofthe fast species, with a unique equilibrium µ v s ( v f ) that must be a productdistribution over the different compartments. Similarly, to case (1) when allthe movement is fastest, now all the movement is on the medium time scale,and the movement of all molecules (fast and slow) is independent and canbe viewed as combined on one time scale with unique stationary probabilitydistribution P ( s f ,s s ) (3.8) (over all i ∈ I ). This implies that on the slow timescale dt , the effective change of S s is due to reaction dynamics with rates˜ λ CR , kd ( v Nf , v Ns ) that have been averaged over P ( s f ,s s ) , and is given by ζ s .Lemma 2.10 implies that S s ( · ) is well defined and given by the solution of(3.21) with rates as in (3.25). (cid:3) If the movement of fast species is slower than fast reactions [i.e., we con-sider cases (3) or (4)], the equilibria for reactions is always attained beforemovement of fast species can change this equilibrium, as stated in the nextcorollary.
Corollary 3.15 (Irrelevance of movement of fast species).
In cases (3)or (4) of Theorem 3.13, the limiting dynamics of S s is independent of P s f . Proof.
The assertion can be seen directly from (3.29) and (3.30), sincethe right-hand sides do not depend on P s f . (cid:3) If all slow species are continuous, the limiting dynamics for cases (1), (2)and (3), (4) are equal. The key to this observation is the following lemma.
Lemma 3.16.
In the situation of Theorem 3.13, assume all slow speciesare continuous, I s ◦ = ∅ , and let πs s := ( π i ( d ) s i ) i ∈I s ,d ∈D . (i) For stationary probability measures µ s s ( ds f ) of S f | s s from Assump-tion 3.11 (i)(1) and µ v s ( ds f ) of S f | v s from (i)(2) , we have µ s s ( ds f ) = µ πs s ( ds f ) . CALING LIMITS OF SPATIAL COMPARTMENT MODELS (ii) Likewise, for stationary probability measures µ s s ( dv f ) of V f | s s from (i)(3) and µ v s ( dv f ) of V f | v s from (i)(4) , we have µ s s ( dv f ) = µ πs s ( dv f ) . Proof. (i) It suffices to show that µ πs s is a stationary probability mea-sure for the process S f | s s ( · ) from (3.13), since we assumed that this processhas a unique stationary probability distribution. Note that, by independenceof the movement of fast and slow species, for k ∈ K f ,˜ λ CR(1) k ( s f , s s ) = Z X d ∈D λ CR kd ( v · d,f , v · d,s ) P ( s f ,s s ) ( dv f , dv s )= Z X d ∈D λ CR kd ( v · d,f , v · d,s ) P s f ( dv f ) P s s ( dv s )= Z ˜ λ CR(2) k ( s f , v s ) P s s ( dv s )= ˜ λ CR(2) k ( s f , πs s ) . Since these rates are equal, the corresponding equilibrium distributions mustbe equal as well. In other words, the equilibrium distribution µ s s ( ds f ) fromAssumption 3.11(i)(1) must equal the equilibrium distribution µ πs s ( ds f )from Assumption 3.11(i)(2). (ii) follows along similar lines. (cid:3) Corollary 3.17.
Suppose in Theorem 3.13 all slow species are contin-uous, I s ◦ = ∅ . Then, the dynamics of (3.21) is the same among the first twocases (1) and (2), and among the last two cases (3) and (4). Proof.
Since P s s ( dv s ) is the delta-measure on πs s , all assertions canbe read directly from (3.27)–(3.30) together with Lemma 3.16. (cid:3) We note that the case (1) (where all species move faster than the fastreactions occur) plays a special role under mass action kinetics.
Corollary 3.18 (Homogeneous mass action kinetics).
Suppose that inTheorem 3.13 the reaction rates are given by mass action kinetics with con-stants satisfying the homogeneity condition κ k := |D| κ kd Y i ∈I π i ( d ) ν ik . Then the dynamics of S s in case (1) is the same as for the system regardedas a single compartment. P. PFAFFELHUBER AND L. POPOVIC
Proof.
For (1) from (3.14), we only need to calculate the average withrespect to equilibrium of the movement dynamics for both slow and fastspecies. The same calculation as for mass action kinetics in Corollary 3.8 weget the first equality in˜ λ CR(1) k ( s f , s s ) = X d ∈D Z λ CR kd ( v · d,f , v · d,s ) P ( s f ,s s ) ( dv f , dv s )= X d ∈D κ kd Y i ∈I ◦ ν ik ! (cid:18) s i ν ik (cid:19) π i ( d ) ν ik · Y i ∈I • ( π i ( d ) s i ) ν ik = κ k Y i ∈I ◦ ν ik ! (cid:18) s i ν ik (cid:19) Y i ∈I • s ν ik i . Since the right-hand side gives the reaction rates for mass action kinet-ics within a single compartment, as given through V f | v s from (2.16), theequilibrium µ v s ( dz ) from Assumption 2.6(i) and µ s s ( ds f ) from Assump-tion 3.11(i)(1) must be the same and the assertion follows. (cid:3) Example 3.19 (Production from fluctuating source in multiple compart-ments). We consider reaction kinetics from Example 2.8 and extend it toa spatial multi-compartment setting. Recall that the chemical reaction net-work is given (within compartments) by the set of reactions : A + B κ ′ d −→ C, : ∅ κ ′ d −→ B, : B κ ′ d −→ ∅ . We consider Λ CR k ( x ) as in (2.19) with κ ′ k replaced by κ ′ kd , k ∈ { , , } . Wehave x = ( x · d ) d ∈D , x · d = ( x Ad , x Bd ), and the dynamics are given byΛ CR d ( x · d ) = κ ′ d x Ad x Bd , Λ CR d ( x · d ) = κ ′ d , Λ CR d ( x · d ) = κ ′ d x Bd . Movement of species is given as in (3.3). Scaling in each compartment is asin the nonspatial setting (2.20), (2.21) and (2.22), so rescaled species countsare v Ad = N − x Ad , v Bd = x Bd , and rates are λ CR d ( v · d ) = κ d v Ad v Bd , λ CR d ( v · d ) = κ d , λ CR ( v · d ) = κ d v Bd . The process V N = ( V NAd , V
NBd ) is given as in (2.23) and additional move-ment terms. We set η s = η A for movement of slow species and η f = η B formovement of fast species. We assume (as in Assumption 3.4) that move-ment of species A, B have stationary probability distributions ( π A ( d )) d ∈D and ( π B ( d )) d ∈D . We derive the dynamics of S A as S A ( t ) = S A (0) − Z t ¯ λ CR ( S A ( u )) du CALING LIMITS OF SPATIAL COMPARTMENT MODELS for appropriate λ . Since the slow species A are continuous, we are in theregime of Corollary 3.17 and we distinguish the following two cases: Dynamics in the cases (1) + (2). We have S B | s A ( t ) − S B | s A (0)= − Y (cid:18)Z t Z X d ∈D κ d v Ad v Bd P ( s A ,S B | sA ( u )) ( dv A , dv B ) du (cid:19) + Y (cid:18)X d ∈D κ d t (cid:19) − Y (cid:18)Z t Z X d ∈D κ d v Bd P ( s A ,S B | sA ( u )) ( dv A , dv B ) du (cid:19) d = − Y + (cid:18)Z t (cid:18)X d ∈D κ d π A ( d ) π B ( d ) | {z } =:¯ κ s A + X d ∈D κ d π B ( d ) | {z } =:¯ κ (cid:19) S B | s A ( u ) du (cid:19) + Y (cid:18)X d ∈D κ d | {z } =:¯ κ t (cid:19) . Hence, the equilibrium of the above process is as in Example 2.8 given by X ∼ µ s A ( ds B ) = Pois (cid:18) ¯ κ ¯ κ + ¯ κ s A (cid:19) . We can now compute ¯ λ CR(1)+(2) from (3.22) as¯ λ CR(1)+(2) ( s A ) = − Z X d ∈D κ d π A ( d ) s A π B ( d ) s B µ s A ( ds B )(3.31) = − X d ∈D κ d π A ( d ) s A π B ( d ) ¯ κ ¯ κ + ¯ κ s A = − ¯ κ ¯ κ s A ¯ κ + ¯ κ s A . Dynamics in the cases (3) + (4). We have in each compartment d ∈ D V B d | s A ( t ) − V B d | s A (0)= − Y (cid:18)Z t Z κ d v Ad V B d | s A ( u ) P s A ( dv A d ) du (cid:19) + Y ( κ d t ) − Y (cid:18)Z t Z κ d V B d | s A ( u ) P s A ( dv A d ) du (cid:19) d = − Y + (cid:18)Z t ( κ d π A ( d ) s A + κ d ) V B d | s A ( u ) du (cid:19) + Y ( κ d t ) . P. PFAFFELHUBER AND L. POPOVIC
Hence, the equilibrium of the above process is X ∼ µ s A ( dv B d ) = Pois (cid:18) κ d κ d + κ d π A ( d ) s A (cid:19) and for ¯ λ CR(3)+(4) from (3.24) we have¯ λ CR(3)+(4) ( s A ) = − X d ∈D Z κ d v Ad v Bd µ s A ( dv B d ) P s A ( dv A d )(3.32) = − X d ∈D κ d κ d π A ( d ) s A κ d + κ d π A ( d ) s A . Note that we are in the regime of Corollary 3.15, which shows that ¯ λ CR(3)+(4) is independent of π B . Comparison of dynamics in cases (1) + (2) and (3) + (4). Let us comparethe case (1) + (2), when the turnover rate of A is given by (3.31), and(3) + (4), when the rate is given by (3.32). First note that even when thenetwork is spatially homogeneous (the chemical constants satisfy assumptionin Corollary 3.18) there is a marked difference between the dynamics of cases(1) + (2) (as in single compartment case) and cases (3) + (4) depending onthe movement equilibria π A and π B . However, if we additionally suppose theslow species A are equidistributed π A ( d ) = 1 / |D| then all four cases have thesame dynamics. Conserved quantities on the fast time scale.
Now, we include conservedquantities in our two-scale system in multiple compartments, that is, wehave a two-scale reaction network with dim( N (( ζ f ) t )) =: n f >
0. We willuse the same notation as in Section 2.3. In particular, Θ f := ( θ c j ) j =1 ,...,n f are linearly independent vectors which span the null space of ( ζ f ) t . Ev-ery θ c j has a unique parameter α i associated with it, j = 1 , . . . , | Θ f | = n f .Here, Θ f ◦ is the subset of conserved quantities for which α c j = 0, and Θ f • is the subset of conserved quantities for which α c j >
0. Conservation meansthat t
7→ h θ c j , S f | s s ( t ) i with S f | s s from (3.13) is constant, j = 1 , . . . , | Θ f | . Welet S Nc j = h θ c j , S Nf i and S Nc = ( S Nc j ) i =1 ,..., | Θ f | be the vector of rescaled con-served quantities. Again, K θ cj is the set of reactions such that β k = α c j and h θ c j , ζ · k i 6 = 0, and let K c := S | Θ f | j =1 K θ ci , K c ◦ := S | Θ f ◦ | j =1 K θ cj and K c • := S | Θ f • | j =1 K θ cj . We still let ζ c be the matrix defined by (2.29).Again, we consider the four cases as given in (3.12). In addition, we assumethat h θ c j , S N i changes on the time scale dt . We write here, distinguishingfast species, conserved quantities and slow species, v = ( v f , v c , v s ) with v f = CALING LIMITS OF SPATIAL COMPARTMENT MODELS ( v id ) i ∈I f ,d ∈D , v c = ( h θ c j , v · d,f i ) j =1 ,..., | Θ f | ,d ∈D , v s = ( v id ) i ∈I s ,d ∈D , as well as s = ( s f , s c , s s ), s f = ( s i ) i ∈I f , s c = ( h θ c j , s f i ) j =1 ,..., | Θ f | , and s s = ( s i ) i ∈I s . Remark 3.20 (Conserved quantities as new species). In light of theprevious results, one would guess that conserved quantities on the fast timescale can be handled as if they are new chemical species, evolving on theslow time scale. However, an important distinction between slow species andconserved quantities does exist: movement of conserved quantities occurs onthe time scale N η f dt rather than on the time scale N η s dt , on which it oc-curs for the slow species. This implies an important distinction between slowspecies and conserved quantities occurs in cases (2) and (3); the averagingmeasures over the intermediate time scales treat conserved and slow speciesdifferently.Although what follows resembles our previous results, in order to be ableto use them, we do have to state the assumptions and results for systems withconserved species explicitly. We omit all the proofs as they follow analogoussteps to those for systems without conserved species. Assumption 3.21 (Dynamics of the spatial multi-scale reaction networkwith conserved quantities). In each case (1)–(4), the spatial two-scale re-action network on time scale
N dt , where Assumption 3.4 holds, satisfies thefollowing conditions:(i) (1) Given ( Y k ) k ∈K f ◦ , the time-change equation of the dynamics of S f given the values of S s = s s and S c = s c , denoted ( S f | ( s s ,s c ) ( t )) t ≥ , givenby (3.13) with S f | s s replaced by S f | ( s s ,s c ) , has a unique solution, where˜ λ CR(1) k ( s f , s s ) is given by (3.14). In addition, S f | ( s s ,s c ) ( · ) has a unique station-ary probability measure µ ( s s ,s c ) ( ds f ) on R |I f | + with h θ c j , s f i = s c j , µ ( s s ,s c ) -almost surely, j = 1 , . . . , | Θ f | .(2) Given ( Y k ) k ∈K f ◦ , the time-change equation of the dynamics of S f giventhe value of V s = v s and S c = s c , denoted ( S f | ( v s ,s c ) ( t )) t ≥ , given by (3.15)with S f | v s replaced by S f | ( v s ,s c ) , has a unique solution, where ˜ λ CR(2) k ( s f , v s )is given by (3.16). In addition, S f | ( v s ,s c ) ( · ) has a unique stationary probabil-ity measure µ ( v s ,s c ) ( ds f ) on R |I f | + with h θ c j , s f i = s c j , µ ( v s ,s c ) -almost surely, j = 1 , . . . , | Θ f | .(3) Given ( Y kd ) k ∈K f ◦ ,d ∈D , the time-change equation of the dynamics of V f given the values of S s = s s and V c = v c , denoted by ( V f | ( s s ,v c ) ( t )) t ≥ ,given by (3.17) with V · d,f | s s replaced by V · d,f | ( s s ,v c ) , has a unique solution, P. PFAFFELHUBER AND L. POPOVIC where ˜ λ CR(3) kd ( v f , s s ) is given by (3.18). In addition, V f | ( s s ,v c ) ( · ) has a uniquestationary probability measure µ ( s s ,v c ) ( dv f ) on R |I f |×|D| + with h θ c j , v f i = v c j , µ ( s s ,v c ) -almost surely, j = 1 , . . . , | Θ f | .Moreover, given s s and s c , the movement dynamics of V c | ( s s ,s c ) is a uniquesolution V c | ( s s ,s c ) ( · ) = ( V · d,c | ( s s ,s c ) ) d ∈D of the time-change equations h θ c j , V · d,f i | ( s s ,s c ) ( t ) − h θ c j , V · d,f i | ( s s ,s c ) (0)= X i ∈I f ◦ θ c j i X d ′ ,d ′′ ∈D ( δ d ′′ ( d ) − δ d ′ ( d )) × Y i,d ′ ,d ′′ (cid:18)Z t λ M i,d ′ ,d ′′ Z v id ′ µ ( s s ,V c | ( ss,sc ) ( u )) ( dv f ) du (cid:19) ,θ c j ∈ Θ f ◦ , (3.33) h θ c j , V · d,f i | ( s s ,s c ) ( t ) − h θ c j , V · d,f i | ( s s ,s c ) (0)= X i ∈I f • θ c j i X d ′ ∈D Z t Z ( λ M i,d ′ ,d v id ′ − λ M i,d,d ′ v id ) × µ ( s s ,V c | ( ss,sc ) ( u )) ( dv f ) du, θ c j ∈ Θ f • , with an equilibrium probability distribution of movement P ( s s ,s c ) ( dv c ) with P d ∈D v · d,c = s c , P ( s s ,s c ) ( dv c )-almost surely.(4) Given ( Y kd ) k ∈K f ◦ ,d ∈D , the time-change equation of the dynamics of V f given the values of V s = v s and V c = v c , denoted by ( V f | ( v s ,v c ) ( t )) t ≥ , givenby (3.19) with V · d,f | v s replaced by V · d,f | ( v s ,v c ) , has a unique solution, where˜ λ CR(4) kd ( v f , v s ) is given by (3.20). In addition, V c | ( v s ,v c ) has a unique sta-tionary probability measure µ ( v s ,v c ) ( dv f ) with h θ c j , v f i = v c j , µ ( v s ,v c ) -almostsurely, j = 1 , . . . , | Θ f | .Moreover, given v s and s c , the movement dynamics of V c | ( v s ,s c ) is aunique solution V c | ( v s ,s c ) = ( V · d,c | ( v s ,s c ) ) d ∈D of the time-change equations(3.33) with V c | ( s s ,s c ) replaced by V c | ( v s ,s c ) with an equilibrium probabilitydistribution of movement P ( v s ,s c ) ( dv c ) with P d ∈D v · d,c = s c , P ( v s ,s c ) -almostsurely.(ii) From { ˜ λ CR( ℓ ) k } ℓ =1 , , , , we set in each case¯ λ CR(1) k ( s s , s c ) = Z R |I f | + ˜ λ CR(1) k ( s f , s s ) µ ( s s ,s c ) ( ds f );(3.34) CALING LIMITS OF SPATIAL COMPARTMENT MODELS ¯ λ CR(2) k ( s s , s c ) = Z R |I s |×|D| + Z R |I f | + ˜ λ CR(2) k ( s f , v s ) µ ( v s ,s c ) ( ds f ) P s s ( dv s );(3.35)¯ λ CR(3) k ( s s , s c ) = Z R | Θ f |×|D| + Z R |I f |×|D| + ˜ λ CR(3) k ( v f , s s ) µ ( s s ,v c ) ( dv f )(3.36) × P ( s s ,s c ) ( dv c );¯ λ CR(4) k ( s s , s c ) = Z R |I s |×|D| + Z R | Θ f |×|D| + Z R |I f |×|D| + ˜ λ CR(4) k ( v f , v s ) µ ( v s ,v c ) ( dv f )(3.37) × P ( v s ,s c ) ( dv c ) P s s ( dv s ) . For j = 1 , , ,
4, there exists a well-defined process ( S s ( · ) , S c ( · )) that is theunique solution of S s ( t ) = S s (0) + X k ∈K s ◦ ζ s · k Y k (cid:18)Z t ¯ λ CR( ℓ ) k ( S s ( u ) , S c ( u )) du (cid:19) (3.38) + X k ∈K s • ζ s · k Z t ¯ λ CR( ℓ ) k ( S s ( u ) , S c ( · )) du, and for j = 1 , . . . , | Θ f | , S c j ( t ) = S c j (0)+ X k ∈K c ◦ X i ∈I f θ ic j ζ sik Y k (cid:18)Z t ¯ λ CR( ℓ ) k ( S s ( u ) , S c ( u )) du (cid:19) (3.39) + X k ∈K c • X i ∈I f θ ic j ζ sik Z t ¯ λ CR( ℓ ) k ( S s ( u ) , S c ( u )) du. (iii) Same as (iii) in Assumption 2.6 in each compartment. Remark 3.22 (Equivalent formulation). For the dynamics under theabove assumption, the following is immediate: In each case (1)–(4), the spa-tial two-scale reaction network on time scale dt , where Assumption 3.21holds, satisfies the following condition: given ( Y k ) k ∈K s ◦ ∪K c ◦ , the time changeequations (3.38) and (3.39) have a unique solution, with¯ λ CR( ℓ ) k ( s s ) := E ( s s ,s c ) (cid:20)X d ∈D λ CR kd ( V · d ) (cid:21) < ∞ . (3.40)The distribution of ( V id ) i ∈I ,d ∈D in (3.40) depends on the parameters η s , η f as follows:(1) P ( s s ,s c ) ( dv f , dv s ) = P s s ( dv s ) Z R |I f | + P s f ( dv f ) µ ( s s ,s c ) ( ds f ) , P. PFAFFELHUBER AND L. POPOVIC (2) P ( s s ,s c ) ( dv f , dv s ) = P s s ( dv s ) Z R |I f | + P s f ( dv f ) µ ( v s ,s c ) ( ds f ) , (3) P ( s s ,s c ) ( dv f , dv s ) = P s s ( dv s ) Z R | Θ f |×|D| + µ ( s s ,v c ) ( dv f ) P ( s s ,s c ) ( dv c ) , (4) P ( s s ,s c ) ( dv f , dv s ) = P s s ( dv s ) Z R | Θ f |×|D| + µ ( v s ,v c ) ( dv f ) P ( v s ,s c ) ( dv c ) . Theorem 3.23 (Heterogeneous two-scale system with conserved fast quan-tities).
Let V N be the vector process of rescaled species amounts for the re-action network which is the unique solution to (3.3). Assume that ( α, β, γ =0) satisfy two-scale system assumptions (2.11) for some I f , I s with ε = 1 and N (( ζ f ) t ) = span(Θ f ) [with ζ f from (2.13) and Θ f from (2.28)] withincompartments with conserved quantities ( θ c j ) i =1 ,..., | Θ f | on the fast time scale.In addition, η i = η f > , i ∈ I f , η i = η s > , i ∈ I s , one of the cases (1)–(4)holds and Assumption 3.21 holds. Then, if ( S Ns (0) , S Nc (0)) N →∞ = ⇒ ( S s (0) , S c (0)) ,we have joint convergence of the process of rescaled amounts of slow andconserved quantities ( S Ns ( · ) , S Nc ( · )) to ( S s ( · ) , S c ( · )) in the Skorohod topol-ogy, with S s the solution of (3.38) and S c the solution of (3.39) with ratesgiven by (3.34)–(3.37). Results analogous to Corollary 3.15 stating the irrelevance of the move-ment of fast species in cases (3) and (4) does not carry over to the case withconserved quantities, since on the time scale N η f dt conserved quantities arestill preserved and their movement equilibria affects the end result. Lemma 3.24.
In the situation of Theorem 3.23, assume I s ◦ = ∅ , that is,all slow species are continuous, and let πs s := ( π i ( d ) s i ) i ∈I s ,d ∈D . (i) For stationary probability measures µ ( s s ,s c ) ( ds f ) of S f | ( s s ,s c ) from As-sumption 3.21 (i) (1) and µ ( v s ,s c ) ( ds f ) of S f | ( v s ,s c ) from (i)(2) , we have µ ( s s ,s c ) ( ds f ) = µ ( πs s ,s c ) ( ds f ) . (ii) Likewise, for stationary probability measures µ ( s s ,v c ) ( dv f ) of V f | ( s s ,v c ) from (i)(3) and µ ( v s ,v c ) ( dv f ) of V f | ( v s ,v c ) from (i)(4) , we have µ ( s s ,v c ) ( dv f ) = µ ( πs s ,v c ) ( dv f ) . Corollary 3.25.
Suppose in Theorem 3.23 all slow species are contin-uous, I s ◦ = ∅ . Then dynamics of (3.38) is the same in cases (1), (2) andalso in cases (3), (4). CALING LIMITS OF SPATIAL COMPARTMENT MODELS Corollary 3.26 (Homogeneous mass action kinetics).
Corollary 3.18carries over to the same situation as in Theorem 3.23.
Example 3.27 (Michaelis–Menten kinetics in multiple compartments).We place Michaelis–Menten reaction kinetics from Example 2.13 in a spatialmulti-compartment setting. The chemical reaction network is given (withincompartments) by the set of reactions from (2.34), with κ ′ k replaced by κ ′ kd incompartment d . We have x = ( x · d ) d ∈D , x · d = ( x Sd , x Ed , x ESd , x
P d ), and thedynamics in each compartment d is given by rates (2.35) with κ ′ k replacedby κ ′ kd . Movement of species is given as in (3.3). Again, we set α S = α P =1 , α E = α ES = 0, and κ d = κ ′ d , κ − d = N − κ ′− d and κ d = N − κ ′ d as in(2.37) so, setting the rescaled species counts v Sd = N − x Sd , v Ed = x Ed , v ESd = x ESd , v
P d = N − x P d , and β = 1 , β − = 1 , β = 1 as in (2.36). We write λ CR d ( v · d ) = κ d v Sd v Ed , λ CR − d ( v · d ) = κ − d v ESd , λ CR d ( v · d ) = κ d v ESd . The process V N = ( V NSd , V
NEd , V
NESd , V
NP d ) is given as in Example (2.13) plusadditional movement terms. We set η s = η S = η P for movement of slowspecies and η f = η E = η ES for movement of fast species. We assume (as inAssumption 3.4) that movement of species i has a stationary probability dis-tribution ( π i ( d )) d ∈D . We have I f = I f ◦ = { E, ES } and I s = I s • = { S, P } and K f = K s = { , − , } , K S = { , − } , K E = K ES = K and ζ, ζ f , ζ s as in (2.37).For conserved quantities within compartments, we set V Cd := V Ed + V ESd and note that while movement changes the values of V Cd the overall sum S C = P d ∈D V Cd := m is a conserved quantity for all times, and thus, thedynamics of S C is trivial.We derive the dynamics of S S (as in Example 2.13, S S + S P is a conservedquantity on the slow time scale). We have from (3.21) that S S ( t ) = S S (0) − Z t ¯ λ CR ( S S ( u )) du for appropriate ¯ λ . Since all slow species are continuous, we are in the regimeof Corollary 3.25 and we only need to distinguish the following two cases: Dynamics in cases (1) + (2). From (3.13) with S f | s s replaced by S f | ( s s ,s c ) ,we have S E | s S ( t ) − S E | s S (0)= − Y (cid:18)Z t Z X d ∈D κ d v Sd v Ed P ( s S ,S E | sS ( u )) ( dv S , dv E ) du (cid:19) P. PFAFFELHUBER AND L. POPOVIC + Y − + (cid:18)Z t Z X d ∈D ( κ − d + κ d ) v ESd P ( m − S E | sS ( u )) ( dv ES ) du (cid:19) = − Y (cid:18)Z t (cid:18)X d ∈D κ d π S ( d ) π E ( d ) | {z } =:¯ κ (cid:19) s S S E | s S ( u ) du (cid:19) + Y − + (cid:18)Z t (cid:18)X d ∈D ( κ − d + κ d ) π ES ( d ) | {z } =:¯ κ − +¯ κ (cid:19) ( m − S E | s S ( u )) du (cid:19) . Hence, the equilibrium of the above process is as in Example 2.13 given by X ∼ µ ( s S ,m ) ( ds E ) where X ∼ Binom (cid:18) m, ¯ κ − + ¯ κ ¯ κ − + ¯ κ + ¯ κ s S (cid:19) and S ES | s S has equilibrium m − X . We next compute ¯ λ CR(1)+(2) from (3.34)as ¯ λ CR(1)+(2) ( s S ) = ¯ λ CR(1)+(2) ( s S ) − ¯ λ CR(1)+(2) − ( s S )= X d ∈D Z ( κ d π S ( d ) π E ( d ) s S s E − κ − d π ES ( d )( m − s E )) µ ( s S ,m ) ( ds E )= X d ∈D κ d π S ( d ) π E ( d ) s S m ¯ κ − + ¯ κ ¯ κ − + ¯ κ + ¯ κ s S − κ − d π ES ( d ) m ¯ κ s S ¯ κ − + ¯ κ + ¯ κ s S = m ¯ κ ¯ κ s S ¯ κ − + ¯ κ + ¯ κ s S . Comparing this with (2.38), we see that in cases (1) + (2) Michaelis–Mentenkinetics in multiple compartments equals the same kinetics in a single com-partment, when κ i is exchanged by ¯ κ i , i = − , , ; compare also with Corol-lary 3.25. Dynamics in cases (3) + (4). For simplicity, we assume that λ M d,d ′ := λ M E,d,d ′ = λ M ES,d,d ′ , that is, movement of E and ES is the same, and hence π E ( d ) = π ES ( d ) , d ∈ D [we will use this property for deriving P ( s S ,m ) ( dv C )and P ( v S ,m ) ( dv C ) below]. We will treat the cases (3) and (4) separately andshow the result of Corollary 3.25 which states that these two cases lead tothe same limiting dynamics. CALING LIMITS OF SPATIAL COMPARTMENT MODELS (3) From Assumption 3.21(i)(3), for v C with P d ∈D v Cd = m , we have V Ed | ( s S ,v C ) ( t ) − V Ed | ( s S ,v C ) (0)= − Y d (cid:18)Z t Z κ d V Ed | ( s S ,v C ) ( u ) v Sd P s S ( dv Sd ) du (cid:19) + Y ( − + ) d (cid:18)Z t ( κ − d + κ d )( v Cd − V Ed | ( s S ,v C ) ( u )) du (cid:19) = − Y d (cid:18)Z t κ d V Ed | ( s S ,v C ) ( u ) s S π S ( d ) du (cid:19) + Y ( − + ) d (cid:18)Z t ( κ − d + κ d )( v Cd − V Ed | ( s S ,v C ) ( u )) du (cid:19) . Hence, the equilibrium of V Ed | ( s S ,v C ) is as in Example 2.13 given by X d ∼ µ ( v Sd ,v Cd ) ( ds Ed ) = Binom (cid:18) v Cd , κ − d + κ d κ − d + κ d + κ d s S π S ( d ) (cid:19) ( s Ed )and V ESd | ( s S ,v C ) has equilibrium v Cd − X d . We compute ¯ λ CR(3) from (3.36)as¯ λ CR(3) ( s S ) = ¯ λ CR(3) ( s S ) − ¯ λ CR(3) − ( s S )= X d ∈D Z κ d v Ed s S π S ( d ) − κ − d ( v Cd − v Ed ) µ ( v Sd ,v Cd ) ( dv Ed ) P ( s S ,m ) ( dv Cd )= X d ∈D Z (cid:18) κ d v Cd κ − d + κ d κ − d + κ d + κ d s S π S ( d ) s S π S ( d ) − κ − d v Cd κ d s S π S ( d ) κ − d + κ d + κ d s S π S ( d ) (cid:19) P ( s S ,m ) ( dv Cd )= X d ∈D Z v Cd κ d κ d s S π S ( d ) κ − d + κ d + κ d s S π S ( d ) P ( s S ,m ) ( dv Cd ) . Consider the equilibrium P ( s S ,m ) ( dv Cd ) of movement dynamics for conservedspecies V Cd = V Ed + V ESd . Since we assume the same migration dynamics for E and ES , the equilibrium P ( s S ,m ) ( dv C ) is given by a multinomial distribu-tion with parameters m, ( π E ( d )) d ∈D and R v Cd P ( s S ,m ) ( dv Cd ) = mπ E ( d ) , d ∈D .For (4), the overall rate ¯ λ CR(4) from (3.37) has the same form as ¯ λ CR(3) except that P ( s S ,m ) ( dv Cd ) is replaced by P ( v Sd ,m ) ( dv Cd ). We first derive the P. PFAFFELHUBER AND L. POPOVIC equilibrium probability distribution µ ( v Sd ,v Cd ) ( dv Ed ) as above. Here, we findthat s S π S ( d ) is replaced by v Sd , leading to X d ∼ µ ( v Sd ,v Cd ) ( dv Ed ) = Binom (cid:18) v Cd , κ − d + κ d κ − d + κ d + κ d v Sd (cid:19) ( v Ed ) . The conserved quantities V Cd ( · ) follow the same dynamics as in case (3) ex-cept that s S π S ( d ) is replaced by v Sd and, therefore, P ( v S ,m ) ( dv C ) is a multi-nomial distribution with parameters m, ( π E ( d )) d ∈D as in case (3). Hence,in the equation for the rates we have R v Cd P ( v Sd ,m ) ( dv Cd ) = mπ E ( d ) , d ∈ D ,and the limiting dynamics in cases (3) and (4) is given by¯ λ CR(3)+(4) ( s S ) = X d ∈D mπ E ( d ) κ d κ d s S π S ( d ) κ − d + κ d + κ d s S π S ( d ) . Although we have seen that the dynamics for cases (1) + (2), as well as for(3) + (4) is the same, in general they are quite different from each other,unless some very special relationships between the chemical constants andmovement equilibria in different compartments are assumed. See the resultsin Pfaffelhuber and Popovic (2014) on notions of dynamical homogeneitythat allow one to make some interesting conclusions.
4. Discussion.
Specific features and extensions of spatial chemical reac-tion models. (a)
Heterogeneous reaction and migration rates.
The reaction rates Λ CR kd in general depend on the compartment d . For the same reason, the outflowof species i from compartment d ′ , P d ′′ ∈D Λ M i,d ′ ,d ′′ might depend on i and d ′ . Moreover, it is possible that Λ CR kd ( x · d ) is zero for some compartments,that is, our model is flexible enough to restrict some reactions to a subsetof compartments. Analogously, movement of certain species types can berestricted to only a subset of compartments, that is, Λ M i,d,d ′ can also be set tozero for some i, d, d ′ . The only thing which is required is that every reaction k happens within at least one compartment.(b) Geometry of space.
The geometry of the spatial system has not beenexplicitly relevant for our results. The reason is that movement dynamics isassumed to happen at a different time scale (either faster or slower) than theeffective reaction dynamics of either the slow or fast species. This impliesthat only the equilibrium of the movement is relevant for any dynamicsoccurring on the respectively slower scale.(c)
Chemical conformations.
Our model can be extended in order to modeldifferent chemical conformations of chemical species instead of spatial com-partments. For this, let D i be the set of possible conformations of species i . Then any molecule of species i performs a Markov chain on D i due tochanges in conformation. Moreover, in this case for each type of reaction k CALING LIMITS OF SPATIAL COMPARTMENT MODELS its reaction rate Λ CR k,d,d ′ might then depend on all conformations of reactingand produced molecules d = ( d i ) i ∈I and d ′ = ( d ′ i ) i ∈I , respectively. For exam-ple, our results can be applied to Michaelis–Menten kinetics with multipleconformations of the enzyme and of the enzyme-substrate complex [see Kou(2008)].(d) Other density dependent processes.
The model can also be appliedto other density dependent Markov chain models, such as epidemic or eco-logical models. Analogous results can also be made for density dependentstochastic differential models of stochastic population growth in spatiallyheterogeneous environments [see Evans et al. (2013)].
Conclusions.
The main conclusion of our paper is the following algorithmfor determining the dynamics of a spatial chemical reaction network: assumewe are given a network of the form (2.1) in a spatial context, that is, (3.1)holds with reaction rates as in Assumption 3.1; introduce a (large) scalingconstant N and rewrite the dynamics of all species in the form (3.3) (forsome α i ’s, η i ’s and β k ’s) assuming (2.8) and (3.2) hold (admittedly, thechoice of N , α i ’s and β k ’s is rather an art than a science—for simplicity,we are assuming here that this step has been done already); in addition,suppose every species moves between compartments as in Assumption 3.4;the goal is to understand the dynamics of overall normalized sums of speciesover compartments as given in (3.5).There are two cases: either the system is on a single-scale, that is, (2.9)holds, or the system is two-scale, that is, (2.11) holds. (We do not treathigher order scales in this paper.)(i) In the single-scale case Theorem 3.7 applies. Essentially, one has toaverage all reaction rates of reactions affecting slow species over the equilib-rium distribution of movement of all species. If reaction rates are given bymass action kinetics, Corollary 3.8 applies.(ii) The two-scale case is considerably more complicated. Here, everyspecies is either fast or slow and we have to consider all orders of the timescale of fast reactions and movement of fast and slow species. We call S f theoverall sum of normalized fast species and S s the overall sum of normalizedslow species. Consider the submatrices of slow and fast reactions, ζ f and ζ s from (2.13) and (2.15), respectively. A conserved quantity for the fastreaction subnetwork is a nontrivial element of the null-space of ( ζ f ) t .(ii-a) If there is no conserved quantity, we can use Theorem 3.13. Here,there are up to four time scales to consider, movement of fast and slowspecies, the time scale of the fast reactions and the time scale of the slowspecies. In all cases, in order to determine the effective rate on S s on a slowertime scale, one has to average over the equilibrium of all higher time scales.Interestingly, if all slow species are continuous (i.e., have a deterministic P. PFAFFELHUBER AND L. POPOVIC process as a limit), it only matters if the fast species move faster or slowerthan fast reactions. The speed of the movement of slow species does notmatter (see Corollary 3.17).(ii-b) If there are conserved quantities for the fast reaction subnetwork,these conserved quantities can still change on a slower time scale. Here, weare assuming that this time scale is the same as the time scale of the slowspecies. The main difference from the case without conserved quantities isthat on the fast time scale, the equilibria we need to consider for averagingare concentrated on a fixed conserved quantity. Then, basically, the con-served quantity can be treated as new species with its own dynamics (whichchanges on the timescale of slow species by assumption). Again, there arefour cases to consider; see Theorem 3.23. Also, if all slow quantities are con-tinuous, it only matters if the fast species move faster or slower than thefast reactions; see Corollary 3.25.REFERENCES
Ball, K. , Kurtz, T. G. , Popovic, L. and
Rempala, G. (2006). Asymptotic analysisof multiscale approximations to reaction networks.
Ann. Appl. Probab. Blount, D. (1994). Density-dependent limits for a nonlinear reaction–diffusion model.
Ann. Probab. Burrage, K. , Burrage, P. M. , Leier, A. , Marquez-Lago, T. and
Nicolau, D. V. Jr. (2011). Stochastic simulation for spatial modelling of dynamic processes in a living cell.In
Design and Analysis of Biomolecular Circuits: Engineering Approaches to Systemsand Synthetic Biology
Drawert, B. , Lawson, M. J. , Petzold, L. and
Khammash, M. (2010). The diffu-sive finite state projection algorithm for efficient simulation of the stochastic reaction–diffusion master equation.
Journal of Chemical Physics
E, W. , Liu, D. and
Vanden-Eijnden, E. (2007). Nested stochastic simulation algorithmsfor chemical kinetic systems with multiple time scales.
J. Comput. Phys.
Elowitz, M. B. , Levine, A. J. , Siggia, E. D. and
Swain, P. S. (2002). Stochastic geneexpression in a single cell.
Science Signalling
Ethier, S. N. and
Kurtz, T. G. (1986).
Markov Processes: Characterization and Con-vergence . Wiley, New York. MR0838085
Evans, S. N. , Ralph, P. L. , Schreiber, S. J. and
Sen, A. (2013). Stochastic populationgrowth in spatially heterogeneous environments.
J. Math. Biol. Fange, D. , Berg, O. G. , Sjoberg, P. and
Elf, J. (2010). Stochastic reaction–diffusionkinetics in the microscopic limit.
Proc. Natl. Acad. Sci. USA
Franz, U. , Liebscher, V. and
Zeiser, S. (2012). Piecewise-deterministic Markovprocesses as limits of Markov jump processes.
Adv. in Appl. Probab. Howard, M. and
Rutenberg, A. D. (2003). Pattern formation inside bacteria: Fluctu-ations due to the low copy number of proteins.
Phys. Rev. Lett. Jeschke, M. , Ewald, R. and
Uhrmacher, A. M. (2011). Exploring the performance ofspatial stochastic simulation algorithms.
J. Comput. Phys. Kang, H.-W. (2012). A multiscale approximation in a heat shock response model of E.coli.
BMC Systems Biology Kang, H.-W. and
Kurtz, T. G. (2013). Separation of time-scales and model reductionfor stochastic reaction networks.
Ann. Appl. Probab. Kang, H.-W. , Kurtz, T. G. and
Popovic, L. (2014). Central limit theorems and diffu-sion approximations for multiscale Markov chain models.
Ann. Appl. Probab. Kotelenez, P. (1988). High density limit theorems for nonlinear chemical reactions withdiffusion.
Probab. Theory Related Fields Kou, S. C. (2008). Stochastic modeling in nanoscale biophysics: Subdiffusion within pro-teins.
Ann. Appl. Stat. Kouritzin, M. A. and
Long, H. (2002). Convergence of Markov chain approximations tostochastic reaction–diffusion equations.
Ann. Appl. Probab. Kurtz, T. G. (1970). Solutions of ordinary differential equations as limits of pure jumpMarkov processes.
J. Appl. Probab. Kurtz, T. G. (1981).
Approximation of Population Processes . CBMS-NSF Regional Con-ference Series in Applied Mathematics . SIAM, Philadelphia, PA. MR0610982 Kurtz, T. G. (1992). Averaging for martingale problems and stochastic approximation.In
Applied Stochastic Analysis (New Brunswick, NJ, 1991) . Lecture Notes in Controland Inform. Sci.
McAdams, H. H. and
Arkin, A. (1997). Stochastic mechanisms in gene expression.
Proc.Natl. Acad. Sci. USA Pfaffelhuber, P. and
Popovic, L. (2014). How spatial heterogeneity shapes multi-scale biochemical reaction network dynamics.
Journal of the Royal Society Interface Srivastava, R. , Peterson, M. and
Bentley, W. (2001). Stochastic kinetic analysis ofthe escherichia coli stress circuit using sigma(32)-targeted antisense.
Biotechnology andBioengineering Takahashi, K. , Tanase-Nicola, S. and
Ten Wolde, P. R. (2010). Spatio-temporalcorrelations can drastically change the response of a mapk pathway.
Proc. Natl. Acad.Sci. USA
Abteilung f¨ur Mathematische StochastikAlbert-Ludwigs-Universit¨at FreiburgEckerstr. 179104 FreiburgGermanyE-mail: [email protected]