Macroscopic limit for stochastic chemical reactions involving diffusion and spatial heterogeneity

To model bio-chemical reaction systems with diffusion one can either use stochastic, microscopic reaction-diffusion master equations or deterministic, macroscopic reaction-diffusion system. The connection between these two models is not only theoretically important but also plays an essential role in applications. This paper considers the macroscopic limits of the chemical reaction-diffusion master equation for first-order chemical reaction systems in highly heterogeneous environments. More precisely, the diffusion coefficients as well as the reaction rates are spatially inhomogeneous and the reaction rates may also be discontinuous. By carefully discretizing these heterogeneities within a reaction-diffusion master equation model, we show that in the limit we recover the macroscopic reaction-diffusion system with inhomogeneous diffusion and reaction rates.


Introduction
Modeling chemical reactions and diffusion have been investigated extensively in the literature due to the frequent appearance of these processes in nature.There are different levels of these models.At the microscopic scale, one can use the stochastic reaction-diffusion master equations (RDME), while at the macroscopic, the reaction-diffusion systems (RDS) are more common.The RDME is useful for simulation as well as when molecular stochasticity due to small quantities is taken into account, while the latter is more convenient to investigate qualitative behaviour of the system.The connection between these models have been investigated since the seventies, with the pioneering works of Kurtz [Kur70,Kur71] where this issue was studied for the homogeneous (or well-stirred) case, i.e.only chemical reactions are taken into consideration.In these two works, the Markov jump process, which describes chemical reactions in the microscopic level, is shown to converge in probability to the corresponding differential system, which models the same reactions at the macroscopic level.When a diffusion process is also present, attempts to connect microscopic RDME and macroscopic RDS have been made in [BVdB80,NP77], and a first rigorous proof was provided in [AT80].In [AT80], the authors studied the reaction and diffusion processes of a single chemical, and proved the convergence in the sense of probability of the corresponding Markov jump process to the reaction-diffusion equation.This result was later improved by [Kot86a,Kot86b,Kot88,Blo91,Blo93] where, among other things, central limit theorems were established for local diffusion, and recently by [WT22] where the diffusion is assumed to be non-local.
We highlight that most of existing works have dealt with scalar reaction-diffusion equations.An exception is the recent work [DP22] where the authors showed exponential convergence to equilibrium for systems with with degenerate reaction rataes; i.e. reactions only occur in a subdomain of positive measure.Our current work extends the literature by considering general first-order, or unimolecular, reaction systems with high levels of heterogeneity, for which we show convergence in probability of microscopic RDME to corresponding RDS.More precisely, we consider chemical reaction networks with an arbitrary number of chemicals, the diffusion and reaction coefficients are spatially inhomogeneous, and in particular, the reaction coefficients can be discontinuous.
Another motivation of this paper stems from applications in molecular communication, an emerging field spanning telecommunication, bio-and nano-technologies.The core idea of molecular communication is to send and receive messages using chemical signals mimicking how cells communicate, which was first proposed in [NSM + 05, NEH13].Understanding dynamics of bio-chemical systems is essential for designing molecular communication schemes.In our recent work [AET20], we proposed a novel detection scheme called "equilibrium signaling" which is robust to the geometry and heterogeneity of considered system.In that work, the microscopic RDME was used for computational purpose while the macroscopic PDE provided quantitative dynamics of the bio-chemical system, and the connection between these models has been assumed therein.Our current work provides rigorous convergence from RDME to PDE for general first-order reaction systems, of which the problem in [AET20] is a special case.
Finally, we remark that the study of RDME for higher order reactions is sparse in the literature.One particular reason for this is that, in the rescaling limit of RDME the biomolecular reactions are lost in two or more dimensions, see [Isa09].Attempts to circumvent this problem includes the convergent RDME, see [Isa13].A general approach has been proposed in [MSW23] where the gradient flow structure of detailed balanced reaction systems has been investigated in a formal setting.Let us also mention related works considering the limit from system of interacting particles to reaction-diffusion equations [LLN20] or integral-differential equations [IMS22].
Our paper is organized as follows: in the next section, in order to clearly present ideas and methods, we start with a reversible system of two chemicals in heterogeneous medium in one dimension.For this system, we first consider in Subsection 2.1 its macroscopic PDE description and basic properties such as global existence, boundedness and large time behaviour of solutions.The corresponding microscopic RDMS equation is written down in Subsection 2.2.The main result concerning the limit from RDME to PDE model is stated in Subsection 2.3 together with its proof.These results are then extended to general first-order reaction networks in Section 3 in arbitrary dimension.
Notation: In the this paper, we will use the following notation: • For 1 ≤ p ≤ ∞, L p (Ω) denotes the usual Lebesgue space with the associated norm • The Sobolev space H 1 0 (Ω) is defined as where the gradient is understood in the distributional sense and the value on the boundary is understood in the trace sense.
• For a function u : Ω → R, Ω ⊂ R n , the gradient of is defined by For a vector field F = (F 1 , . . ., F n ) : Ω → R n , the divergence of F is defined by • Consider the random variable X defined on the measurable space (X , B) with underlying probability space (Ω X , F , P).The probability distribution of X is denoted by P and the probability of an event {X ∈ B}, B ∈ B is denoted by P(X ∈ B).
2 The case of two chemicals

Macroscopic PDE Model
Consider the following reversible reaction of two chemical species S 1 and S 2 S 1 where λ 1 and λ 2 are reaction rates.The reactions take place in a bounded vessel Ω ⊂ R n , with Lipschitz boundary ∂Ω.In general, the reaction rates λ 1 and λ 2 are spatially inhomogeneous; that is, the reaction rate depends on the spatial position in the medium with λ i : Ω → [0, ∞).
Denote by u i (x, t) the concentration of S i at time t > 0 and position x ∈ Ω.The chemicals species diffuse through the vessel with the diffusive fluxes By applying the law of mass action, with the diffusive flux given in (2), we obtain the following macroscopic reaction-diffusion system x ∈ Ω, i = 1, 2. (3) Here, ∇ • F is the divergence of the vector field F .
Lemma 2.1.Assume that Ω ⊂ R n be a bounded domain with Lipschitz boundary ∂Ω.
Then there exists a constant C P > 0 such that for all u ∈ H 1 0 (Ω).Theorem 2.1 (Global bounded weak solution of (3)).Assume the following for each i ∈ {1, 2} and λ i ≥ 0 a.e. in Ω.Then, for any non-negative initial data u 0 = (u i,0 ) ∈ L ∞ + (Ω) 2 , there exists a unique global non-negative, weak solution to (3) (in the sense of (9)) which satisfies the mass dissipation Moreover, the solution is bounded uniformly in time, i.e. there is ρ > 0 depending on initial data, In addition, if , the solution to (3) can be represented by where and the semigroup {T (t) = e L t } t≥0 is generated by the operator Furthermore, if λ 1 (x) = δλ 2 (x), x ∈ Ω, for some positive constant δ, then the solution to (3) decays exponentially to zero, i.e. for any p ∈ [2, ∞) there are α p , C p > 0 such that Proof.A weak solution to (3) on (0, T ), T > 0, is a pair of functions It is noted that the sum of the right-hand side of the equations in (3) equals to zero.Therefore, the global existence and uniform-in-time boundedness of a unique non-negative, weak solution to (3) follows immediately from [FMTY21, Theorem 1.1].
It remains to show the decay of the solution in case Thanks to the homogeneous Dirichlet boundary conditions and the Poincaré inequality in Lemma 2.1, it follows that there is β > 0 satisfying and consequently Now for any 2 < p < ∞, we have The proof for u 2 follows the same way.

When the diffusion coefficients
that the operator L defined in (7) generates an analytic contraction semigroup {T (t) = e L t } t≥0 , using which the solution to (3) can be written as (6).Note that this solution is also a weak solution and therefore enjoys the uniform-in-time bounds (5) as well as the exponential decay (8).
Since the system is linear, the global existence is not surprising, and one can get a bound in L ∞ (Ω)-norm which might depend on time (typically exponential).The uniform-intime bound (5), which is a consequent of the dissipation of mass (4), plays an important role in defining the stopping process in the sequel (see ( 21)).

Reaction-diffusion master equations
Let Ω = (0, 1), and consider a volume in one dimension, i.e.V = (0, L) for some L > 0. The case of a cube Ω = (0, 1) n and V = (0, L) n in R n will be discussed in Section 3 for general systems.At the microscopic level, the system (1) with diffusion can be modeled via a continuous-time Markov jump process [AT80,Kot86a].Let N ∈ N and consider a uniform partition of V into N cells of equal size w = L/N and define the j-th cell by (x j−1 , x j ] = ((j − 1)w, jw], j = 1, . . ., N. It is remarked that this partition is not a discretization of the domain Ω in the macroscopic model.In fact, later on, it will be assumed that the volume V , as well as each cell, will "blow up", i.e. become unbounded (see ( 25)).To account for the boundary conditions detailed in the sequel, the volume is extended to an interval of length L + 2w consisting of N + 2 cells of length w.The number of molecules of each species l ∈ {1, 2} in cell j = 0, . . ., N + 1 at time t is denoted by X l j (t).We assume that the number of molecules for the cells j = 0 and j = N + 1 are always zero, which corresponds to the homogeneous Dirichlet boundary condition in (3).The total state of the system is denoted by X(t), which forms a continuous-time Markov jump process on N 2(N +2) .Suppose that at time t, X l j (t ) and denote the j-th unit vector in R N +2 by e j .The transition rates for the process X(t) are given by and q k,k+m = 0 otherwise.Here, D l j and λ l j denote the diffusion coefficient and reaction rate constant in the j-th cell for species l, which are calculated from the functions D l and λ l via It is remarked that we only need integrability of D l and λ l for these voxel diffusion and reaction rates to be meaningful.
The evolution of X(t) is described by the reaction-diffusion master equation (RDME); i.e., the Kolmogorov forward equation for the Markov jump process described by the transition rates in (10).Let ) and e l j = e (l−1)N +j ∈ N 2(N +2) .The RDME is then given by Recall that, for all t > 0, the boundary layers satisfy Xl j (t) = 0, j = 0, N + 1 to enforce the zero boundary condition.
We now define the concentration process C(t) via where the last two conditions account for the zero boundary condition.The process C(t) forms a continuous-time Markov jump process, a property inherited from the process X(t).Observe that where f : R It follows that C(t) is a density dependent continuous Markov process [Kur71] with waiting time parameter and the jump distribution function given by We remark that due to the definition (16), the sum over all m ∈ N 2(N +2) in ( 17) and ( 18) are indeed finite sums.

The macroscopic limit
To allow for a comparison between the microscopic and macroscopic models, we interpret C(t) as a function on [0, 1].In particular, let The function u(t) lies in a subspace X N , which is a subspace of Define τ S as the first exit time of u(t) from with arbitrary ρ > ρ, where ρ is the uniform-in-time bound (5) of solutions in Theorem 2.1.The stopped process ũ(t) is then defined as Let Then, ũ(t) is a jump Markov process with The main result of the current work is the following theorem.
The rest of this section is devoted to prove Theorem 2.2.We start off with the behaviour of a martingale related to the space-time Markov jump process ũ.Then discretisation of the inhomogeneous diffusion and its corresponding semigroup are considered, which in turn will be used in comparing ũ and u in the last part.

An accompanying martingale
Consider the operator on X N which is viewed as a stochastic approximation of the average infinitesimal rate of change of u.We now establish some properties of Lemma 2.2.The process (28) is a martingale with respect to P u (the distribution induced by an initial value u).
Proof.With τ defined in (24) we observe that and sup τ (u) By [Kur71, Proposition 2.1], the process (28) is then a martingale with respect to P u .
Observe that for an arbitrary initial distribution, Proof.Since f is the infinitesimal generator of the process ũ, the result follows by definition.

Discretization of diffusion and reaction terms
Due to the inhomogeneity of diffusion coefficients and reaction rate constants, we approximate them by the following piecewise constants functions: for l = 1, 2, D l (r) = D l j and λ l (x) = λ l j , for x ∈ I j = [x j−1 , x j ), where D l j and λ l j are defined as in (12).Now to study the corresponding discrete evolution we define the discrete reactions and the discrete diffusion operator where the difference quotients are defined as Using the definition of D l , we can rewrite for x ∈ I j , 1 ≤ j ≤ N − 1.Note that from the definition of f in (16), F in (27), we have Remark 2.1.The form of L l N (u)(x) in (35) generalizes the analogous expression for spatially homogeneous diffusion, as considered in, e.g., [AT80,Kot86a].It is not the only potential generalization; however, as we will see, is the correct choice to establish the desired convergence result.
Proof.For simplicity, we omit the superscript l in L l N .We show that 1 0 For u, v ∈ H 1 0 (Ω), we have where we used v| [−h,0] = u| [1,1+h] = 0 at the last step.Now we apply this "integration by parts" to obtain which is the self-adjointness of the operator L N .
Remark 2.2.We note that, unlike the case of spatially homogeneous diffusion coefficients as considered in [AT80, Kot86a], the operator L N is not self-adjoint under Neumann boundary conditions (i.e., a reflective boundary).
Proof.We will apply the Hille-Yosida Theorem.Due to the definition, we have D(−L N ) = {u ∈ C 0 (Ω) : u| ∂Ω = 0} which is dense in L 2 (Ω).Let ω > 0. We claim that for each v ∈ L 2 (Ω), there exists a unique solution to the equation Indeed, by defining the bilinear form a a(u, z) = ω u, z + −L N u, z we have Therefore, thanks to Lax-Milgram theorem, (39) has a unique solution u ∈ L 2 (Ω), which means that (ωI − L N ) −1 is one-to-one and onto.It remains to show that By multiplying (39) by u in L 2 (Ω) and using −L N u, u ≥ 0 we have which shows (40).Therefore, we can apply Hille-Yosida theorem to obtain that T N (t) = exp(−L N t) is an analytic contraction semigroup in L 2 (Ω).The properties (i) and (ii) therefore follow immediately from contraction semigroup.

Proof of Theorem 2.2
Using (6), (28), and (37), we can split the difference ũ − u as Let y(t) be the solution to which is equivalent to Therefore, we have Inserting this into (42) and taking the norm of both sides give We estimate the terms (I), (II) and (III) separately.

Estimate of (II).
Lemma 2.6.Let {T N (t)} t≥0 and {T (t)} t≥0 be the semigroups generated by L N and L respectively.We have where for each t > 0.
Proof.By triangle inequality, it holds where we used that {T N (t)} t≥0 is a contraction semigroup at the last step.It remains to show ǫ 1 (t) := (T N (t) − T (t))u(0) → 0 as N → ∞ for each t > 0. This is nothing else but the L 2 (Ω)-convergence of finite difference scheme for the diffusion equation ∂ t u − L u = 0 where coefficients of L are smooth, in this case, C 1 , and therefore the convergence follows from e.g.[Jov89].
Estimate of (I).We show that Indeed, we compute using ( 28) and ( 43) It follows that which is (47) due to w(0) = 0. Further calculations lead to where we used Lemma 2.5 (ii) at the last estimate.
Estimate of (III).Now consider Using the contraction property of the semigrioup T N (t) and boundedness of the reaction rate constants λ 1 (•) and λ 2 (•), we have where the constant c 0 depends on λ 1 L ∞ (Ω) and λ 2 L ∞ (Ω) .For (III 2 ), we use the contraction of {T N (t)} again to have with the property lim for each t > 0, thanks to the definition of R N in (34).Finally, by using similar arguments to the proof of Lemma 2.6, we have for each t > 0. These estimates yield

General first order chemical reaction networks
In this section, we show that the techniques presented above can be applied to any first order chemical reaction network with heterogeneities.Let S 1 , . . ., S K be K ≥ 1 chemicals which react through the following reactions Here λ ji : Ω → R + denotes the reaction constant from S i to S j which depends on the position x ∈ Ω.
Macroscopic PDE model: Let u i (x, t) be the concentration of the chemical S i at position x ∈ Ω and time t > 0. Assume moreover that D i (x) is the diffusion coefficient of S i .The macroscopic reaction-diffusion system reads as where subject to the homogeneous Dirichlet boundary conditions and non-negative initial data The diffusion coefficients, which could be possibly discontinuous, are merely assumed to be bounded from above and below, i.e. there D * , D * > 0 such that and the reaction rate functions are bounded, non-negative, i.e.
• If we consider the reaction network (54) as a directed graph where all the species are nodes and reactions are directed edges, then the weak reversibility in Definition 3.1 is equivalent to the strongly connectedness of the corresponding directed graph.
• The weak reversibility for general chemical reaction networks is in fact more general, especially for higher order chemical reactions, in which the corresponding directed graph might have several disjoint strongly connected components.However, since the network in (54) is of first order, each strongly connected component can be treated separately, and therefore we adapt the definition in Definition 3.1 for convenience.
Proof.Denote by f i (u) the right hand side of (55).Thanks to (56), it follows that The global existence of bounded solutions and uniform-in-time boundedness then follow from [FMTY21] immediately.It remains to show the exponential decay in the case of weakly reversible networks and a ij (x) = γ ij ϕ(x).We rewrite the equation of u i as Due to the weak reversibility of the network, there exists a unique positive state u Thanks to [FPT17, Lemma 2.3], we have for some δ > 0, which implies E(u)(t) ≤ e −δt E(u 0 ) for all t > 0, and therefore, for all i = 1, . . ., K, The convergence in L p (Ω)-norm follows from (59) and L 2 (Ω)-convergence u i (t) p L p (Ω) ≤ u i (t) p−2 L ∞ (Ω) u i (t) 2 L 2 (Ω) , thus u i (t) L p (Ω) ≤ Ce −(2/p)t , ∀t > 0.
On the boundary of V , additional cells are added.Molecules in cells (j 1 , . . ., j i , 1, j i+1 , . . ., j n ) and (j 1 , . . ., j i , N, j i+1 , . . ., j n ) may diffuse into these boundary cells in which they are absorbed.That is, the number of molecules in these cells are always zero.The total number of cells in the system is then (N + 2) n .The number of molecules of each species l ∈ {1, . . ., K} in cell j = (j 1 , . . ., j n ) at time t is denoted by Xl j (t).The number of molecules in cells (j 1 , . . ., j i , 0, j i+1 , . . ., j n ) and (j 1 , . . ., j i , N + 1, j i+1 , . . ., j n ) are always zero.The total state of the system is denoted by X(t), which is modeled by a continuous-time Markov jump process on N K(N +2) n .
We remark that the main difference from the case of K = 2 chemical species arises in the analog of Lemma 2.3.In the case of general K, for the martingale process z(t) we have where λ * = max i,j λ ij L ∞ (Ω) .Hence, if we seek E[ z(t) 2 ] → 0 as N, w → 0, we require N 2 w n → 0, which corresponds to assumption (H1').