
The past few decades have seen robust research on questions regarding the existence, form, and properties of stationary distributions of stochastically modeled reaction networks. When a stochastic model admits a stationary distribution an important practical question is: what is the rate of convergence of the distribution of the process to the stationary distribution? With the exception of [
Citation: David F. Anderson, Jinsu Kim. Mixing times for two classes of stochastically modeled reaction networks[J]. Mathematical Biosciences and Engineering, 2023, 20(3): 4690-4713. doi: 10.3934/mbe.2023217
[1] | Linard Hoessly, Carsten Wiuf . Fast reactions with non-interacting species in stochastic reaction networks. Mathematical Biosciences and Engineering, 2022, 19(3): 2720-2749. doi: 10.3934/mbe.2022124 |
[2] | Giorgos Minas, David A Rand . Parameter sensitivity analysis for biochemical reaction networks. Mathematical Biosciences and Engineering, 2019, 16(5): 3965-3987. doi: 10.3934/mbe.2019196 |
[3] | Daniele Cappelletti, Badal Joshi . Transition graph decomposition for complex balanced reaction networks with non-mass-action kinetics. Mathematical Biosciences and Engineering, 2022, 19(8): 7649-7668. doi: 10.3934/mbe.2022359 |
[4] | Ting Kang, Yanyan Du, Ming Ye, Qimin Zhang . Approximation of invariant measure for a stochastic population model with Markov chain and diffusion in a polluted environment. Mathematical Biosciences and Engineering, 2020, 17(6): 6702-6719. doi: 10.3934/mbe.2020349 |
[5] | Luyao Li, Licheng Fang, Huan Liang, Tengda Wei . Observer-based event-triggered impulsive control of delayed reaction-diffusion neural networks. Mathematical Biosciences and Engineering, 2025, 22(7): 1634-1652. doi: 10.3934/mbe.2025060 |
[6] | H. J. Alsakaji, F. A. Rihan, K. Udhayakumar, F. El Ktaibi . Stochastic tumor-immune interaction model with external treatments and time delays: An optimal control problem. Mathematical Biosciences and Engineering, 2023, 20(11): 19270-19299. doi: 10.3934/mbe.2023852 |
[7] | Dengxia Zhou, Meng Liu, Ke Qi, Zhijun Liu . Long-time behaviors of two stochastic mussel-algae models. Mathematical Biosciences and Engineering, 2021, 18(6): 8392-8414. doi: 10.3934/mbe.2021416 |
[8] | Xuehui Ji, Sanling Yuan, Tonghua Zhang, Huaiping Zhu . Stochastic modeling of algal bloom dynamics with delayed nutrient recycling. Mathematical Biosciences and Engineering, 2019, 16(1): 1-24. doi: 10.3934/mbe.2019001 |
[9] | David F. Anderson, Tung D. Nguyen . Results on stochastic reaction networks with non-mass action kinetics. Mathematical Biosciences and Engineering, 2019, 16(4): 2118-2140. doi: 10.3934/mbe.2019103 |
[10] | Linda J. S. Allen, P. van den Driessche . Stochastic epidemic models with a backward bifurcation. Mathematical Biosciences and Engineering, 2006, 3(3): 445-458. doi: 10.3934/mbe.2006.3.445 |
The past few decades have seen robust research on questions regarding the existence, form, and properties of stationary distributions of stochastically modeled reaction networks. When a stochastic model admits a stationary distribution an important practical question is: what is the rate of convergence of the distribution of the process to the stationary distribution? With the exception of [
Stochastic models of reaction networks, which will be formally introduced in Section 2, are used ubiquitously in the biosciences to track the (integer valued) counts of different interacting species. This class of models is utilized at multiple different scales, from molecular processes to the level of populations. Numerous previous works have considered questions related to the stationary behavior of these models: Lyapunov function approaches for positive recurrence [2,3,4], algebraic approaches for closed forms of stationary distributions [5,6,7], stability for first-order reaction networks [1,8], and so on. In this paper, we study the natural followup question: if π is the stationary distribution of the model, how fast does Pt(x,⋅) converge to π(⋅), where Pt(x,A)=P(X(t)∈A|X(0)=x) is the time-dependent distribution of the process, denoted by X, given an initial condition of X(0)=x?
A natural route forward is to consider the mixing time of the process. Let ε∈(0,12). For a Markov process X with a stationary distribution π, the mixing time for the process with initial condition x is
τεx=inft≥0{‖Pt(x,⋅)−π(⋅)‖TV≤ε}, | (1.1) |
where the total variation distance between two probability measures on a measurable space (Ω,F) is defined as ‖μ−ν‖TV=supA∈F|μ(A)−ν(A)|. In our case of a discrete state space, ‖μ−ν‖TV=12∑x|μ(x)−ν(x)|, and so the total variation norm is an L1 norm. See [9] for a nice introduction to mixing times.
In Theorems 3.1 and 3.2, we will provide results on the mixing times for two distinct classes of reaction networks, with each class characterized by different graph topological conditions on the associated reaction network. The network types considered in this paper were first considered in [2], where positive recurrence (and hence existence of a stationary distribution) was established for each. Thus, this paper can naturally be viewed as a followup to [2]. For the sake of consistency, throughout this manuscript we used the same format and wordings of the key definitions and theorems as [2].
For now, we simply write Cond1 and Cond2 for the sets of conditions considered in Theorems 3.1 and 3.2, respectively (the specific conditions will be detailed in Section 3 when we have the proper terminology available). Each of our theorems is then of the following form: for a stochastically modeled reaction network whose associated graph satisfies the conditions Condi and whose stationary distribution is πi, there is an ηi>0 and a function Bi:Zd≥0→R≥0 for which
‖Pti(x,⋅)−πi(⋅)‖TV≤Bi(x)e−ηit,for all t≥0. | (1.2) |
Hence, τεx≤1η[ln(Bi(x)+ln(ε−1)]=O(ln(Bi(x)). The exponential decay rate (1.2) with ηi>0 means that the Markov process X is exponentially ergodic. Moreover, we characterize Bi(x) in each case with B1(x)=C1(|x|+1)ln(|x|+2) and B2(x)=C2, for some C1,C2>0. Hence, the convergence for the second class is shown to be uniform over the initial condition.
Our proofs rely upon the use of Foster-Lyapunov functions. Specifically, for Theorem 3.1 we will provide an appropriate function, V:Zd≥0→R≥0, with V(x)→∞ as |x|→∞, and constants a>0 and b>0 so that AV(x)≤−aV(x)+b for all x outside some compact set, where A is the generator of the process. For Theorem 3.2 we will provide an appropriate function V, and positive constants a,b, and δ so that AV(x)≤−aV(x)1+δ+b for x outside a compact. The (known) results of Section 4 related to Foster-Lyapunov functions are then utilized to complete the argument. The new results presented here should be compared with the results of [2] in which positive recurrence was proven by demonstrating the existence of an appropriate function V with AV(x)≤−1 for all x outside some compact set. The sharper estimates provided here are a result of a more careful analysis of the generator.
The remainder of the paper is outlined as follows. In Section 2, we introduce stochastic models of reaction networks. In Section 3, we precisely state our main results. In Section 4, we provide some necessary results related to tiers (in the context of reaction networks) and Foster-Lyapunov functions. In Section 5, we provide proofs for our main results using appropriate Foster-Lyapunov functions. In Section 6, we provide some generalizations to our main results.
Stochastic reaction networks are a class of continuous-time Markov chains on Zd≥0 that are used to model biochemical systems and other population processes [10,11,12,13]. Of particular interest to mathematical researchers is how the topological conditions of the underlying reaction graph relates to the qualitative behavior of the associated dynamical system (such as the mixing time). We proceed with the basic definitions from the field.
A reaction network is a graphical construct that describes a set of possible interactions among some constituent "species."
Definition 2.1. A reaction network is given by a triple of finite sets (S,C,R) where:
(i) The species set S={S1,S2,⋯,Sd}, with d<∞, contains the species of the reaction network.
(ii) The reaction set R={R1,R2,⋯,Rr}, with r<∞, consists of ordered pairs (y,y′)∈R where
y=d∑i=1yiSiandy′=d∑i=1y′iSi, | (2.1) |
and where the values yi,y′i∈Z≥0 are the stoichiometric coefficients. We will often write reactions (y,y′) as y→y′.
(iii) The complex set C consists of the linear combinations of the species in (2.1). Specifically, C={y | y→y′∈R}∪{y′ | y→y′∈R}. For the reaction y→y′, the complexes y and y′ are termed the source and product complex of the reaction, respectively.
The state space of our eventual Markov model, which will be fully described after some discussion about reaction networks, will be described by the copy numbers of the species. The transitions of the model will be determined via the reactions.
We note that in particular examples, as in most of ours below, we do not enumerate the species as S1,…,Sd and instead utilized more suggestive notation (e.g., E for Enzyme, P for a protein, etc.), or alphabetical labeling (e.g., A, B, C, etc.). The clear enumeration of {S1,S2,…,Sd} will be used primarily for theory and proofs.
The linear combinations in (2.1) are, for example, of the form 2S1+S2 or 2S1+S2+4S3. However, we will also associate each complex y with the vector whose j-th component is yj, i.e. y=(y1,y2,⋯,yd)T∈Zd≥0.
For example, when S={S1,S2,…,Sd}, the complex 2S1+S2 is associated with the vector (2,1,0,0,…,0)T∈Zd≥0. We further note that it is reasonable to consider a complex y with yi=0 for each i. This complex is denoted by ∅.
We sometimes enumerate complexes and/or reactions. Hence, we will sometimes write, for example, y1 or yr to denote particular complexes, as opposed to the 1st or rth element of the vector y. Similarly, we may refer to the rth reaction as yr→y′r∈R. In such cases, the jth coordinate of, for example, y1 and yr will be denoted y1,j and yr,j, respectively. Context will always make it clear if yr is referring to the rth element of a complex y or if yr denotes a particular complex (such as the source complex of the reaction yr→y′r).
It is most common to present a reaction network with a directed graph, termed the reaction graph of the network, in which the vertices are the complexes, each complex is written exactly one time (even if it appears in more than one reaction), and the directed edges are given by the reactions. We present an example to solidify notation.
Example 2.1. The following reaction graph is associated with a reaction network modeling a usual substrate-enzyme kinetics
S+E⇄SE→E+P. |
For this reaction network, S={S,E,SE,P}, C={S+E,SE,E+P} and R={S+E→SE,SE→S+E,SE→E+P}.
It is conditions on the reaction graph that will be important to us. It is unfortunately the case that some of the terminology used in the study of reaction networks is different than that used more commonly. Hence, we discuss here some of the differences.
A subgraph, which is connected, of the reaction graph is termed a linkage class [14]. If a linkage class is strongly connected then we say that linkage class is weakly reversible. Finally, we say that the entire network/graph is weakly reversible if each connected component (linkage class) is weakly reversible.
Example 2.2. If we consider a reaction network described with the reaction graph,
![]() |
then one can see that there are four species, eight complexes (vertices), seven reactions (directed edges), and three linkage classes (connected components). The right-most connected component is strongly connected, and hence weakly reversible, whereas the other two are not. Hence, the network as a whole is not said to be weakly reversible.
The following notation will be utilized.
1. For u,v∈Rd≥0, we define uv=∏di=1uvii, with 00 taken to be equal to 1.
2. For x∈Zd≥0, we let x!=∏di=1xi!.
We now show the most common way to stochastically model a reaction network. Given a reaction network (S,C,R), a (stochastic) kinetics is an assignment of a function λy→y′:Zd≥0→R≥0 to each reaction y→y′∈R. The time evolution of the copy numbers of species is then modeled by means of a continuous-time Markov chain with state space Zd≥0, whose transition rate from state x to state z is given by
qx,z=∑y→y′∈Ry′−y=z−xλy→y′(x), |
where the sum is over those reactions whose occurrence causes a net change that is precisely z−x. If infinitely many transitions occur within a finite-time T∞ (in other words, an explosion occurs at T∞), we let X(t)=Δ for any t≥T∞, where Δ is a cemetery state not contained in Zd≥0 [15]. The infinitesimal generator A of the associated Markov process acts on functions via the operation
Af(x)=∑y→y′∈Rλy→y′(x)(f(x+y′−y)−f(x)), | (2.2) |
for f:Zd≥0→R [16].
One of the most widely used choices of stochastic kinetics, and the one which the present paper focuses on exclusively, is given by (stochastic) mass-action kinetics, where for any reaction y→y′∈R
λy→y′(x)=κy→y′x!(x−y)!1{x≥y}=κy→y′d∏i=1xi!(xi−yi)!1{xi≥yi}, | (2.3) |
for some positive constant κy→y′, called a reaction rate constant (or rate constant for short). We denote K={κy→y′}. We may then denote the stochastic mass-action system by (S,C,R,K). To include the rate constants into the reaction graph, we place them right next to the arrow of the corresponding reaction as in yκy→y′→y′. Typical stochastic simulation algorithms for generating sample trajectories of this model are the Gillespie algorithm [12,17] or the next reaction method [18,19], or are approximated via tau-leaping [20,21]. Expectations and probabilities can be simulated by more advanced techniques [22].
In [2], Anderson and Kim studied two classes of reaction networks and proved that any stochastic mass-action system associated with a network from one of these classes is necessarily positive recurrent. Moreover, they showed that positive recurrence holds regardless of the choice of reaction rate constants. Hence, all such models admit a stationary distribution. The different sets of conditions characterizing the classes are properties of the reaction graph alone, and are easily checked. In this paper, we prove that these models are also exponentially ergodic. Specifically, Theorem 3.1 below states that networks satisfying the conditions detailed in Theorem 1 of [2] are exponentially ergodic. Theorem 3.2 below states that networks satisfying the conditions detailed in Theorem 2 of [2] are exponentially ergodic and, moreover, that the convergence is uniform over the choice of initial condition. After providing the results, with examples, we will close this section with two remarks that we hope shed light on the theorems. We also provide generalizations to our results in Section 6.
Before stating the theorems, a few definitions related to possible network structures are needed. These terms are introduced in [2], but they are also universally used in this field.
Definition 3.1. A reaction network (S,C,R) is called binary if ‖y‖ℓ1=∑di=1yi≤2 for all y∈C.
Definition 3.2. Let (S,C,R) be a reaction network with S={S1,S2,⋯,Sd}. The complex ∅ is termed the zero complex. Unary complexes and binary complexes are the complexes of the form Si and Si+Sj, respectively, for i,j∈{1,…,d}, where we could have i=j. Binary complexes of the form 2Si are called double complexes. If 2Si∈C for each i=1,2,…,d, then the reaction network (S,C,R) is said to be double-full.
Definition 3.3. We will say that (S,C,R) is open if both ∅→S∈R and S→∅∈R for each S∈S.
That is, the system is open if there are inflows and outflows for each species. Note that in the case of open systems, all states communicate with each other and so Zd≥0 is irreducible.
In the two papers [23,24], Anderson showed that deterministically modeled weakly reversible reaction networks with a single linkage class are necessarily persistent and have bounded trajectories. That is, they are "well-behaved." Our first result pertains to a stochastic analog of those models. That is, we start with a single linkage class that is weakly reversible and binary, and then add all inflow and outflow reactions (making the network "open"). Theorem 1 in [2] states that a stochastically modeled system associated with such a network is positive recurrent regardless of the choice of rate constants. The result below gives us exponential ergodicity.
Theorem 3.1. Let (S,C,R) be a weakly reversible, binary reaction network that has a single linkage class and S={S1,…,Sd}. Let ˜C=C∪{∅}∪{Si | Si∈S} and ˜R=R∪Si∈S{∅→Si,Si→∅}. Then, for any choice of reaction rate constants K={κy→y′} there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)ln(|x|+2)e−ηt,for all t≥0 and all x∈Zd≥0, | (3.1) |
where π is the unique stationary distribution of the model on Zd≥0. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|→∞.
Moreover, if there is a w∈Rd>0 for which w⋅(y′−y)=0 for all y→y′∈R, that is if the network (S,C,R) is conservative, then there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)e−ηt,for all t≥0 and all x∈Zd≥0. | (3.2) |
The mixing time of the model still satisfies τεx=O(ln(|x|)), as |x|→∞.
We provide an example.
Example 3.1. We take S={A,B,C} and the binary reaction network with a single linkage class, (S,C,R), to be
A →B. ↖↙2C |
Adding inflows and outflows yields the reaction network (S,˜C,˜R), which has reaction graph
![]() |
(3.3) |
Hence, by Theorem 3.1, for any choice of rate constants the Markov model associated to the network (S,˜C,˜R) is exponentially ergodic with mixing time bounded above by the logarithm of the initial counts.
We provide the results of a numerical example in Figure 1A. Specifically, we consider the stochastic model associated with the network (3.3) in which all rate constants are selected to be equal to one. In this case, the stationary distribution is a product of Poissons [6]. We then parameterized the model via the initial condition xm=(m,m,m)T. On the left hand side of Figure 1A we provide plots of ‖Pt(x,⋅)−π(⋅)‖TV, as a function of time, for each of m=1,10,100. On the right hand side of Figure 1A we provide a plot of τεxm, for xm=(m,m,m)T and ε=0.1, as a function of m. Note that this plot qualitatively agrees with Theorem 3.1 in that it appears logarithmic in m. The plots were made via approximating Pt(x,⋅) via Monte Carlo methods, and estimating the total variation distance via the approximation
‖Pt(x,⋅)−π(⋅)‖TV=12∑z|Pt(x,z)−π(z)|≈12∑z∈[0,200]3|Pt(x,z)−π(z)|. |
The next theorem, which allows us to conclude that the convergence is also uniform over initial states, drops both the weak reversibility and the single linkage class assumption. However, it adds the assumption that the network is double-full and an assumption related to paths through the reaction graph that start with double complexes. Another difference with the previous theorem, though not the corollay, is that now Zd≥0 may no longer be irreducible. This is also an extension of Theorem 2 in [2] that only shows ergodicity of the proposed reaction network class.
Theorem 3.2. Let (S,C,R) be a binary reaction network satisfying the following two conditions:
(i) the reaction network is double-full, and
(ii) for each double complex (of the form 2Si) there is a directed path within the reaction graph beginning with the double complex itself and ending with either a unary complex (of the form Sj) or the zero complex, ∅.
Let K={κy→y′} be a choice of reaction rate constants. Let X(0)=x∈S⊂Zd≥0 be in a closed communication class, S, and let π be the stationary distribution of the model on that class. Then, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,C,R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤Ce−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. Thus, the mixing time satisfies τεx=O(1), as |x|→∞.
We provide an example.
Example 3.2. Consider the following network
2A⇄A+B⇄B,A⇄2C⇄B+C.2B⇄∅,C⇄A+C | (3.4) |
This network is double full since 2A,2B,2C∈C. Moreover, for each double complex, there is a sequence of directed edges starting with the double complex and arriving at a complex that is either ∅ or a unary complex. Specifically, we have the following such sequences:
2A→A+B→B2C→A2B→∅ |
Hence, by Theorem 3.2, for any choice of rate constants the associated Markov model is exponentially ergodic. Moreover, the convergence time is uniform over the initial state.
As in Example 3.1, we provide the results of a numerical example in Figure 1B. Specifically, we consider the stochastic model associated with the network (3.4) in which all rate constants are selected to be equal to one. In this case, the stationary distribution is a product of Poissons [6]. We then parameterized the model via the initial condition xm=(m,m,m)T. On the left hand side of Figure 1B we provide plots of ‖Pt(x,⋅)−π(⋅)‖TV, as a function of time, for each of m=1,10,100. On the right hand side of Figure 1B we provide a plot of τεxm, for xm=(m,m,m)T and ε=0.1, as a function of m. Note that this plot qualitatively agrees with Theorem 3.2 in that it appears to be bounded above in m. The plots were made via approximating Pt(x,⋅) via Monte Carlo methods, and estimating the total variation distance via the approximation
‖Pt(x,⋅)−π(⋅)‖TV=12∑z|Pt(x,z)−π(z)|≈12∑z∈[0,200]3|Pt(x,z)−π(z)|. |
We close this section with two remarks.
Remark 3.1. We note that the exponential convergence bounds in Theorem 3.1 (given in (3.1) and (3.2)) are not uniform over the initial state, x. This can be understood by the fact that the "strength" of the convergence is governed by the linear outflow reactions Si→∅. Such linear terms lead to exponential decay and, hence, logarithmic mixing times. Conversely, for Theorem 3.2 the decay has quadratic rate, which should be compared to the ODE ˙x(t)=−x(t)2, which implodes from infinity (i.e., the time needed to hit a given compact from a particular state is bounded from above).
Remark 3.2. To prove Theorems 3.1 and 3.2 we require a number of results that we collect in the next section. However, we will provide a bit of intuition here. The main point is that we must show that the "average drift" of the model (which will be clarified in the next section through the use of Foster-Lyapunov functions) points inward toward the origin from each x∈Zd≥0 outside of a compact set. There are then two main cases to consider: (ⅰ) those x that are "away" from the boundary in that there are sufficient counts for each reaction to have positive intensity, and (ⅱ) those x near the boundary in that at least one of the reactions has zero intensity (for example, perhaps xi=0 for some i). The case of x "away" from the boundary can be handled in a manner that is similar to the deterministic case [23,24], and relies on the overall structure of the networks. The existence of the outflow reactions in Theorem 3.1 and of the double complex condition in Theorem 3.2 ensure that the drift towards the origin also takes place for the x near the boundary.
This brief section is dedicated to introducing two key analytical tools we will use in the proofs of Theorems 3.1 and 3.2: tiers and Foster-Lyapunov functions.
"Tiers" were introduced to allow for the partition of the complexes along sequences in a manner that tells us which reactions/complexes are "dominating" the dynamics of the mathematical model in different regions of the state space. Tiers were first introduced by Anderson in [23,24] in the context of deterministically modeled reaction networks, but have recently been successfully used to analyze stochastic models as well [2,25].
There are multiple ways to partition the set of complexes along a particular sequence of points. Below, and for reasons that will become clear throughout our analysis, we choose to do so via the function x↦(x∨1), which is the vector whose jth component is
(x∨1)j=xj∨1=max{xj,1}. |
This particular choice of function (and hence partition) was first used in [2] where it was called a "D-type partition."
Definition 4.1. Let (S,C,R) be a reaction network with S={S1,…,Sd}. A sequence {xn} of points in Rd≥0 is called a D-type tier-sequence if
(i) for each i∈{1,…,d} the limit limn→∞xn,i exists (it could be infinity) and limn→∞xn,i=∞ for at least one i, and
(ii) for any pair of complexes y,y′∈C, the limit
limn→∞(xn∨1)y(xn∨1)y′, |
exists (it could be infinity).
Remark 4.1. Note that, given a sequence {xn} of vectors in Rd≥0 with lim supn→∞‖xn‖∞=∞, it is always possible to find a subsequence that is a D-type tier-sequence. This is simply because we assume that only finitely many complexes are present. See Lemma 4.2 of [23].
Definition 4.2. Let (S,C,R) be a reaction network and let {xn} be a D-type tier-sequence in Rd≥0. We define a partition of the complexes C=∪Ki=1TD,i{xn} along the tier-sequence {xn} in the following recursive manner:
(i) we say that a complex y is in tier 1 (and write y∈TD,1{xn}) if for all complexes y′∈C
limn→∞(xn∨1)y(xn∨1)y′>0; |
(ii) we say that a complex y is in tier i (and write y∈TD,i{xn}) if there exists y′∈TD,i−1{xn} with
limn→∞(xn∨1)y(xn∨1)y′=0 |
and for all complexes y′∉⋃i−1j=1TD,j{xn} we have
limn→∞(xn∨1)y(xn∨1)y′>0. |
The mutually disjoint subsets TD,i{xn} are called D-type tiers along {xn}. If y∈TD,i{xn} and y′∈TD,j{xn} with i<j we will write y≻Dy′. If y,y′ are in the same D-type tier, then we will write y∼Dy′.
Thus, those complexes in tier 1 maximize (xn∨1)y along the sequence, with those in tier 2 being second largest, etc. Note that as a consequence of the above definition, if y,y′∈Ti{xn} then there exists C∈(0,∞) such that
limn→∞(xn∨1)y(xn∨1)y′=C, |
explaining the notation y∼Dy′, and if y∈Ti{xn} and y′∈Tj{xn} for some i<j, then
limn→∞(xn∨1)y′(xn∨1)y=0, |
explaining the notation y≻Dy′.
The following example demonstrates the concept.
Example 4.1. Consider the reaction graph below.
2A ⟶ B⇌∅⇌A. ↖↙A+B |
Let C be enumerated with y1=2A,y2=B,y3=∅,y4=A,y5=A+B. For xn=(n,0)⊤, (xn∨1)yi is equal to n2,1,1,n and n for i=1,2,3,4 and 5, respectively. On the other hand, if xn=(n,n+1)⊤, then (xn∨1)yi is equal to n2,n+1,1,n and n(n+1) for i=1,2,3,4 and 5, respectively. Note that both sequences are tier-sequences. The tier structures (i.e., the partition) for each tier-sequence is shown in the following table.
![]() |
We close this section with a useful lemma.
Lemma 4.1. Let {xn}, with xn∈Zd≥0, be a D-type tier-sequence of a reaction network (S,C,R) and let y→y′∈R. If limn→∞λy→y′(xn)>0, where λy→y′ is given in (2.3), then there is a C>0 with
limn→∞λy→y′(xn)(xn∨1)y=C. |
Moreover, if the source complex y satisfies |yi|≤1 for each i∈{1,…,d}, then C=κy→y′.
Proof. Since limn→∞λy→y′(xn)>0, we have that (xn∨1)y=xyn for all n large enough. The result then follows from the basic polynomial structure of λy→y′.
The following known result is key to our analysis. See, for example, [26,27].
Theorem 4.2. Let X be a continuous-time Markov chain on an irreducible, countably infinite state space S⊂Zd with generator A. Suppose there exists a positive function V on S satisfying the following.
(i) V(x)→∞, as |x|→∞, and
(ii) there are positive constants a and b such that
AV(x)≤−aV(x)+bfor all x∈S. | (4.1) |
Then there exists positive constants η and C such that
‖Pt(x,⋅)−π(⋅)‖TV≤C(V(x)+1)e−ηt,∀x∈S. |
Note that if a continuous-time Markov chain X satisfies the conditions in Theorem 4.2 with the function V defined from (5.1), then the mixing time satisfies τεx≤1ηln(C(V(x)+1)ε)=O(ln|x|).
The next result utilizes a 'super Lyapunov function.' Such functions will allow us to show that, for certain models, the convergence is uniform over the initial state. This theorem is a version of Theorem 3.2 in [28].
Theorem 4.3. Let X be a continuous-time Markov chain on an irreducible, countable state space S⊂Rd with generator A. Suppose there exists a positive function V on S satisfying the following.
(i.) V(x)→∞, as |x|→∞, and
(ii.) there are positive constants a, b and δ such that
AV(x)≤−aV(x)1+δ+bfor all x∈S. | (4.2) |
Then there exist positive constants η and C such that
‖Pt(x,⋅)−π(⋅)‖TV≤Ce−ηt. |
Note that in this case we have that the mixing time satisfies τεx=O(1).
We let V:Zd≥0→R≥0 be the following function
V(x)=d∑i=1[xi(ln(xi)−1)+1], | (5.1) |
with 0(ln(0)−1) taken to be zero. This function has a long history in the field of chemical reaction network theory: it has been used to show stability of deterministically modeled reaction networks [24,29,30,31], to show positive recurrence of stochastic reaction networks [2,4,25], and also plays a key role in the study of "non-equilibrium" systems in the physics literature (see [32] and followup references). The basic idea behind the proofs of Theorems 3.1 and 3.2 is to use the infinitesimal behavior of the Markov process, which can be described via the generator A defined in (2.2), and the Foster-Lyapunov function V defined above in (5.1), to show that the models under consideration admit an "inward" drift at every state outside some compact set, in the sense of Theorems 4.2 and 4.3. To do so we make use of the tier structures of the model.
We begin by first providing a series of lemmas.
Lemma 5.1. Let A be the generator (2.2) of the continuous-time Markov chain associated to a reaction network (S,C,R) with mass-action kinetics (2.3) and rate constants {κy→y′}. Let V be the function defined in (5.1). For each D-type tier-sequence {xn} there is a constant C1>0 for which
AV(xn)≤∑y→y′∈Rλy→y′(xn)(ln((xn∨1)y′(xn∨1)y)+C1). | (5.2) |
Furthermore, suppose that there exists a reaction y1→y′1∈R such that
(i) y1∈TD,1{xn},
(ii) y1≻y′1, and
(iii) limn→∞λy1→y′1(xn)(xn∨1)y1=C2, for some C2>0.
Then there exists a positive constant C3 such that for large enough n,
AV(xn)≤C3∑y→y′∈Ry≻Dy′λy→y′(xn)ln((xn∨1)y′(xn∨1)y), | (5.3) |
where the sum is over those reactions, y→y′∈R, with y≻Dy′.
The proof of (5.2) can be found in [2,Lemma 8]. The proof of (5.3) can be found in the proof of [2,Theorem 9].
The next lemma gives a condition for when (5.3) holds.
Lemma 5.2. Let (S,C,R), with S={S1,…,Sd}, be a weakly reversible, binary reaction network that has a single connected component. Suppose that C is not a subset of {Si+Sj:i,j∈{1,…,d}}. Let ˜C=C∪{∅}∪{Si | Si∈S} and ˜R=R∪Si∈S{∅→Si,Si→∅}. Then, for any D-type tier-sequence {xn} for (S,˜C,˜R), there exists a reaction y1→y′1∈˜R such that
(i) y1∈TD,1{xn},
(ii) y1≻Dy′1 and
(iii) limn→∞λy1→y′1(xn)(xn∨1)y1=C, for some C>0.
Proof. First note that ∅∉TD,1{xn}. This follows because there is at least one i∈{1,…,d} with xn,i→∞, since {xn} is a D-type tier sequence, and so Si≻D∅.
Next, we observe that if Si∈TD,1{xn} for some i, then we are done because in this case the reaction Si→∅ satisfies all of (i),(ii), and (iii).
Now suppose that there is no Si∈TD,1{xn}. Then it must be the case that Si+Sj∈TD,1{xn} for some i,j (we could have i=j). However, because C does not consist solely of binary complexes, we must have at least one complex, ˜y∈C, with |˜y|∈{0,1} (i.e., it is either unary or the zero complex). In particular, we know that ˜y∉TD,1{xn}. Hence, by the weak reversibility of (S,C,R), there is a sequence of reactions starting with Si+Sj and ending with a complex not contained in TD,1{xn}. This implies that there is at least one reaction, y→y′∈R, with y∈TD,1{xn} and y≻Dy′. That is, conditions (i) and (ii) are satisfied for these reactions.
We now look closer at y, the source complex of this reaction. If y=2Sk for some k, then we have that this reaction also satisfies (ⅲ) with xn,k→∞, as n→∞, and we are done. Similarly, if y=Sk+Sℓ for some k≠ℓ, and both limn→∞xn,k>0 and limn→∞xn,ℓ>0, then Lemma 4.1 may be used to imply that (ⅲ) holds as well, and we are done. The final case to consider is when y=Sk+Sℓ with limn→∞xn,k>0 and limn→∞xn,ℓ=0. However, this would imply that Sk∈TD,1{xn} as well, which is impossible by assumption. Hence, all cases are covered and the proof is complete.
While the inflow reactions are not crucial, the outflows in ˜R is necessary for Lemma 5.2. The next example demonstrates.
Example 5.1. Let (S,C,R) be a reaction network formed with the reactions in the 'triangle' within the reaction network shown in Example 4.1. That is,
2A ⟶ B ↖↙A+B. |
Let (S,˜C,˜R) be the full network in Example 4.1, where all the inflows and outflows are included. For the D-type tier-sequence xn=(n,0)⊤, the reaction 2A→B satisfies all the conditions (ⅰ), (ⅱ) and (ⅲ) in Lemma 5.2. For the D-type tier sequence xn=(0,n)⊤, the reaction A+B→2A satisfies (ⅰ) and (ⅱ) in Lemma 5.2. However (ⅲ) does not hold because the rate function λA+B→2A is zero at xn for each n. In this case, the out-flow reaction B→∅ satisfies (ⅰ), (ⅱ), and (ⅲ) in Lemma 5.2.
Now we begin the proof of Theorem 3.1.
Proof of Theorem 3.1. We first suppose the existence of a w∈Rd>0 with w⋅(y′−y)=0 for all y→y′∈R. Let W(x)=w⋅x=∑di=1wixi. Note that lim|x|→∞|W(x)|=∞ since w∈Rd>0. For each i∈{1,…,d}, let ei∈Zd≥0 be the vector with a 1 in the ith coordinate and zeros elsewhere. Since w conserves the reactions in R, we have
AW(x)=∑y→y′∈Rλy→y′(x)(W(x+y′−y)−W(x))+d∑i=1λSi→∅(x)(W(x−ei)−W(x))+d∑i=1λ∅→Si(x)(W(x+ei)−W(x))=−d∑i=1κSi→∅⋅wixi+d∑i=1κ∅→Siwi≤−aW(x)+b, |
where a=min{κSi→∅} and b=∑di=1κ∅→Siwi. Hence, (4.1) holds with V replaced with W, giving us the second statement of the theorem, by Theorem 4.2.
Now we assume that no such w exists. Note that this implies that neither C⊂{2Si:Si∈S} nor C⊂{Si:Si∈S}. We claim that there exists a compact set K⊂S such that for some positive constant a, the function V defined in (5.1) satisfies
AV(x)≤−aV(x)for all x∈S∖K. | (5.4) |
If this claim holds, then the result follows from Theorem 4.2 after setting b=(a+1)maxx∈KAV(x).
We prove the claim by contradiction. We therefore suppose that there is not a compact K, together with an a>0, so that (5.4) holds. In this case, there exists a sequence {xn}⊂S such that
limn→∞||xn||∞=∞andAV(xn)≥−1nV(xn), for all n. | (5.5) |
By Remark 4.1, we can find a subsequence of {xn} that is a D-type tier-sequence. For ease of notation, we also denote this subsequence by {xn}. By considering a further subsequence if needed, we also can assume that there exists an i such that
limn→∞xn,ixn,j>0,for any j, | (5.6) |
where we denote the ith coordinate of xn by xn,i. Because of (5.6), we will say that xn,i is a "maximal coordinate, " since no other element is asymptotically larger. Without loss of generality, we assume xn,1 is such a maximal coordinate of xn. That is, there exists a positive constant C′ such that for each i
xn,i≤C′xn,1,for all n large enough. |
By Lemmas 5.1 and 5.2 there is a C1>0 for which (5.3) holds for all n large enough. That is,
AV(xn)≤C1∑y→y′∈Ry≻Dy′λy→y′(xn)ln((xn∨1)y′(xn∨1)y). | (5.7) |
Note that because we are assuming that xn,1 is a maximal coordinate, it holds that S1≻∅, and so the reaction S1→∅ is included in the sum above. Since each term in the summation of the right-hand side of (5.7) is negative for large n, the inequality stands after dropping other terms except for the term related to S1→∅. Hence we have
AV(xn)≤C1λS1→∅(xn)ln((xn∨1)(0,0,…,0)(xn∨1)(1,0,…,0))=−C1λS1→∅(xn)ln(xn,1∨1). |
However, since we assumed xn,1 is a maximal coordinate of xn, and limn→∞xn,1=∞ by construction, we can find some positive constant C2 such that V(xn)≤C2xn,1ln(xn,1) for all n large enough. Combining the above, there is a C3>0 so that for all n large enough we have
AV(xn)≤−C1λS1→∅(xn)ln(xn,1∨1).≤−C3V(xn). |
This contradicts (5.5), and the result is shown.
Remark 5.1. Since (S,˜C,˜R) is weakly reversible in Theorem 3.1, the associated Markov process X with X(0)=x is irreducible. Moreover, because ˜R is open, the irreducible state space is all of Zd≥0.
We now turn to the proof of Theorem 3.2. We begin with a lemma that characterizes tiers of double-full binary reaction networks.
Lemma 5.3. Let (S,C,R) be a binary reaction network satisfying the same conditions given in Theorem 3.2. Let {xn} be a D-type tier-sequence. Then the following hold.
(i) For y→y′∈R with y∈TD,1{xn}, we have that λy→y′(xn)>0 for all n large enough,
(ii) 2Si∈TD,1{xn} for some i, and
(iii) There exists y→y′∈R such that y∈TD,1{xn} and y≻y′.
Proof. There are no complexes with ‖y‖ℓ1>2, and our assumptions imply that each double complex is the source complex for some reaction. Hence, (ⅰ) and (ⅱ) follow.
To see that (ⅲ) holds, you simply have to note that no unary complex (or the zero complex) can be in TD,1{xn}. (ⅲ) then follows from the existence of a directed path from any double complex to a unary complex or the zero complex.
We now prove Theorem 3.2.
Proof of Theorem 3.2. We will show that V defined in (5.1) satisfies the conditions of Theorem 4.3 with δ=1/2 (and value δ∈(0,1) would work). As in the proof of Theorem 3.1, we proceed via a proof by contradiction. Therefore, in order to find a contradiction, we assume that there is no compact set K⊂S, together with an a>0, for which
AV(x)≤−aV(x)3/2for all x∈S∖K. |
Then there exists a sequence {xn}, with xn∈S, so that
limn→∞‖xn‖∞=∞andAV(xn)≥−1nV(xn)3/2, | (5.8) |
for all n. By Remark 4.1, we may extract a subsequence that is a D-type tier sequence and that still satisfies (5.8). As before, for ease of notation we denote this subsequence by {xn}.
By result (ⅱ) in Lemma 5.3, there is a species Si for which 2Si∈TD,1{xn}. We will now show that there must be a C>0 with
AV(xn)≤−Cx2n,i, |
which would contradict (5.8).
By Lemma 5.3, there is a y1→y′1∈R for which y1∈TD,1{xn} and y1≻y′1. By Lemma 4.1 and Lemma 5.1 (ⅰ) we also have that the limit in (ⅲ) of Lemma 5.1 holds. Hence, we may utilize (5.3), with a bound similar to (5.7) in the proof of Theorem 3.1, to conclude that there is a C2>0 so that for n large enough we have
AV(xn)≤C2λy1→y′1(xn)ln((xn∨1)y′1(xn∨1)y1)≤−C2λy1→y′1(xn), | (5.9) |
where the second inequality holds since ln((xn∨1)y′1(xn∨1)y1)→−∞, as n→∞. Since we know that 2Si∈TD,1{xn}, we can then conclude that there is a C>0 so that
AV(xn)≤−C2λy1→y′1(xn)≤−Cx2n,i. |
This was our goal, so we are done.
In this section, we generalize Theorem 3.1 to provide four more classes of reaction networks for which the associated Markov chains with intensity functions (2.3) are exponentially ergodic. Note that we no longer have previous results guaranteeing the models under consideration in this section are positive recurrent. Hence, existence of π is now part of each result. Note, however, that this fact follows immediately since the existence of positive a and b for which AV(x)≤−aV(x)+b outside a compact immediately implies that AV(x)≤−1 outside that same compact (since V is positive and bounded from below).
The first generalization arises by noting that the proof of Theorem 3.1 only utilizes the presence of the outflow reactions. Thus, the theorem can be generalized to only allow for a portion of (or even none of) the inflow reactions. However, in this case the system is not open and so Zd≥0 may not be a closed communication class, and so the generalized statement must restrict to a particular closed class.
Corollary 6.1. Let (S,C,R) be a weakly reversible, binary reaction network that has a single linkage class and S={S1,…,Sd}. Let ˜C=C∪{∅}∪{Si | Si∈S} and ˜R=R∪Si∈S{Si→∅}∪Si∈I{∅→Si}, where I⊂S is any subset of the species (including the empty set). Let K={κy→y′} be a choice of reaction rate constants. Let X(0)=x∈Zd≥0 be in a closed communication class, S. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)ln(|x|+2)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|→∞.
Suppose now that there is a w∈Rd>0 for which w⋅(y′−y)=0 for all y→y′∈R. Let X(0)=x∈Zd≥0 be in a closed communication class, S. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|→∞.
Example 6.1. Consider the system with reaction graph in (3.3) except C→∅ is removed, i.e. the network with reaction graph
![]() |
For this model Zd≥0 is irreducible, and there is a unique stationary distribution on Zd≥0. Moreover, because of Corollary 6.1, it is also exponentially ergodic with a mixing time that is bounded above by the logarithm of the initial counts.
As noted in Section 5, the results of Lemma 5.2 are crucial in proving Theorem 3.1. Hence it is possible to find other structural conditions for the reaction network that guarantee exponential ergodicity by ensuring Lemma 5.2 holds. Essentially, in Lemma 5.2 the existence of a directed path of reactions starting with a complex y∈TD,1{xn} and ending with a complex y′∉TD,1{xn} was crucial. Hence we can generalize Theorem 3.1 with alternative conditions that guarantee the existence of such a path of reactions.
Using the observation above, the corollary below allows us to drop the single linkage class assumption. However, π could be a point mass (for example, when the network is 2S1→S1→∅ and no inflow of S1 is added).
Corollary 6.2. Let (S,C,R) be a binary reaction network with S={S1,…,Sd}. Suppose that for each binary complex y∈C, there is a directed path within the reaction graph beginning with y and ending with either a unary complex (of the form Sj) or the zero complex, ∅. Let ˜C=C∪{∅}∪{Si | Si∈S} and ˜R=R∪Si∈S{Si→∅}∪Si∈I{∅→Si}, where I⊂S is any subset of the species. Let K={κy→y′} be a choice of reaction rate constants. Let X(0)=x∈Zd≥0 be in a closed communication class, S. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)ln(|x|+2)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|→∞.
Suppose now that there is a w∈Rd>0 for which w⋅(y′−y)=0 for all y→y′∈R. Let X(0)=x∈Zd≥0 be in a closed communication class, S. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|→∞.
Proof. It suffices to show that the results (ⅰ), (ⅱ), and (ⅲ) in Lemma 5.2 hold. However, in the proof of Lemma 5.2, the weak reversibility and single linkage class assumptions of (S,C,R) were only used to show the existence of a directed path of reactions beginning with some complex y∈TD,1{xn} and ending with either a unary complex or the zero complex when TD,1{xn} only contains binary complexes. Hence by directly assuming the existence of such a sequence, the result follows.
Example 6.2. Let a reaction network (S,C,R) be described with the following reaction graph.
2B→A+B→C,2C→A←∅→B. |
Then each binary complex in C has a directed path of reactions beginning with itself and ending with a unary complex. Then adding all the out-flows, we have the new reaction network (S,˜C,˜R), which has reaction graph
2B→A+B→ C↓2C→A⇌ ∅⇌B. |
By Collorary 6.2, for any choice of rate constants the Markov model associated to the network (S,˜C,˜R) is exponentially ergodic with mixing time bounded above by the logarithm of the initial counts.
We may also generalize by not requiring that all the species have an out-flow reaction (Si→∅). However, in this case we add the assumption that there is a "chain" of first order reactions (of the form Si→Sj) connecting a species with no out-flow reaction to one with an out-flow reaction. We first consider the case where only one species, S1, say, has an outflow.
Corollary 6.3. Let (S,C,R) be a binary reaction network with S={S1,…,Sd} and {S1,S2,…,Sd}⊂C. Suppose that
1. for each binary complex y∈C, there is a directed path within the reaction graph beginning with y and ending with either a unary complex (of the form Sj) or the zero complex, ∅.
2. for each unary complex Si∈C such that i≠1, there is a directed path within the reaction graph such that it begins with Si, ends with S1, and consists only of unary complexes (i.e. the path is of the form Si→Si1→Si2→⋯→S1).
Let ˜C=C∪{∅} and ˜R=R∪{S1→∅}∪Si∈I{∅→Si}, where I⊂S is any subset of the species. Let K={κy→y′} be a choice of reaction rate constants. Let X(0)=x∈Zd≥0 be in a closed communication class, S. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)ln(|x|+2)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. Therefore the mixing time of the model satisfies τεx=O(ln(|x|)), as |x|→∞.
Suppose now that there is a w∈Rd>0 for which w⋅(y′−y)=0 for all y→y′∈R. Let X(0)=x∈Zd≥0 be in a closed communication class, S. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. The mixing time in this case still satisfies τεx=O(ln(|x|)), as |x|→∞.
Proof. We first show that the results in Lemma 5.2 hold for (S,˜C,˜R). Let {xn} be an arbitrary tier-sequence. We only need to consider the case that at least one unary complex is in TD,1{xn}. The proof of the other cases is the same as the proof of Lemma 5.2. Let y=Sℓ∈TD,1{xn}. Then on a direct path Sℓ→Si1→⋯→S1→∅, there must exist a reaction of the form either Sj→Sk for some j≠k such that Sj∈TD,1{xn} and Sk∉TD,1{xn} or S1→∅ with S1∈TD,1{xn} and ∅∉TD,1{xn} (as always). Then letting this reaction be y1→y′1, Lemma 5.2 follows.
To complete the proof, it suffices to start at (5.7) because the previous parts in the proof of Theorem 3.1 stand without further modification. Let {xn} be an arbitrary tier-sequence. Note that i) if there exists a reaction y→y′ such that y=Sr for some r and y≻Dy′ and ii) if xn,r is a maximal coordinate of {xn}, then the proof is done in the same way as the proof of Theorem 3.1. Let xn,ℓ be the maximal coordinate of {xn}. Then over a directed path Sℓ=Si0→Si1→Si2→⋯→S1→∅, we have that either i) there exists Sr→Sm∈R for some r,m such that xn,r is a maximal coordinate and Sr≻DSm or ii) S1→∅ satisfies that xn,1 is a maximal coordinate and S1≻D∅. This completes the proof.
Example 6.3 Let (S,C,R) be a reaction network that is described with the following reaction graph
![]() |
Then (S,C,R) satisfies conditions (ⅰ) and (ⅱ) of Corollary 6.3 where S1 can be taken to be A, B, or C. Hence, (S,˜C,˜R) can be made with an additional one of the out-flows {A→∅,B→∅,C→∅,D→∅}. Then, by Corollary 6.3, for any choice of rate constants the Markov model associated to the network (S,˜C,˜R) is exponentially ergodic with mixing time bounded above by the logarithm of the initial counts.
Finally, Corollary 6.3 can be generalized further by adding additional out-flows to (S,C,R) when some unary complex does not satisfy the directed path condition. In the corollary below, Sp plays the role of {S2,…,Sd} and Se plays the role of {S1} in Corollary 6.3.
Corollary 6.4. Let (S,C,R) be a binary reaction network with S={S1,…,Sd}. We assume the following conditions.
1. For each binary complex y∈C, there is a directed path within the reaction graph beginning with y and ending with either a unary complex (of the form Sj) or the zero complex, ∅.
2. There exist Sp⊆S and a non-empty subset Se⊆S such that (i) S=Sp∪Se, and (ii) for each Si∈Sp there is a directed path within the reaction graph such that it begins with Si and ends with some Sk∈Se. Furthermore the directed path consists only of unary complexes (i.e. the path is of the form Si→Si1→Si2→⋯→Sk).
Let ˜C=C∪{Si | Si∈Se}∪{∅} and ˜R=R∪S∈Se{S→∅}∪Si∈I{∅→Si}, where I⊂S is any subset of the species. Let K={κy→y′} be a choice of reaction rate constants. Let X(0)=x∈Zd≥0 be in a closed communication class. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)ln(|x|+2)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. Therefore the mixing time of the model satisfies τεx=O(ln(|x|)), as |x|→∞.
Suppose now that there is a w∈Rd>0 for which w⋅(y′−y)=0 for all y→y′∈R. Let X(0)=x∈Zd≥0 be in a closed communication class, S. Then there is a stationary distribution π on that class. Moreover, there exist C,η>0 so that the Markov process X with intensity functions (2.3) associated to the reaction network (S,˜C,˜R) satisfies
‖Pt(x,⋅)−π(⋅)‖TV≤C(|x|+1)e−ηt,for all t≥0, |
where neither C nor η depend upon x∈S. The mixing time in this case still satisfies τεx=O(ln(|x|)), as |x|→∞.
Proof. If Se≠ and Sp≠∅, the proof is exactly the same as the proof of Corollary 6.3 after replacement of S1 by Sk∈Se. If Sp=∅, then Scp=S. In this case the result still follows since (S,˜C,˜R) has the same condition of Corollary 6.2.
Example 6.4. We give an example for Corollary 6.4 with a version of enzyme-substrate kinetics. Consider the reaction network described by the reaction graph below:
S⇌P←∅→E↓↑S+E⇌ SE⇌E+P. | (6.1) |
We put Sp={S} and Se={P,E,SE}. Then (S,˜C,˜R) is made by adding reactions E→∅ and P→∅ to this reaction network (note that SE→∅ is already in the network and so does not need to be added). Similar to Figure 1, in Figure 2 we provide the results of a numerical simulation for the continuous-time Markov chain X associated with (S,˜C,˜R) in which all rate constants are selected to be equal to one. In this case, the stationary distribution of X is a product form of Poissons [6].
DFA is supported from NSF-DMS-2051498 and Army Research Office grant W911NF-18-1-0324. JK is supported from the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT)(No. 2022R1C1C1008491) and the Basic Science Research Institute Fund, whose NRF grant number is 2021R1A6A1A10042944. We also thank the three reviewers of our paper, whose comments improved the final product appreciably.
The authors declare there is no conflict of interest.
[1] | C. Xu, M. C. Hansen, C. Wiuf, Full classification of dynamics for one-dimensional continuous time markov chains with polynomial transition rates, 2021, arXiv: 2006.10548. https://doi.org/10.48550/arXiv.2006.10548 |
[2] |
D. F. Anderson, J. Kim, Some network conditions for positive recurrence of stochastically modeled reaction networks, SIAM J. Appl. Math., 78 (2018), 2692–2713. https://doi.org/10.1137/17M1161427 doi: 10.1137/17M1161427
![]() |
[3] |
A. Agazzi, J. C. Mattingly, Seemingly stable chemical kinetics can be stable, marginally stable, or unstable, Commun. Math. Sci., 18 (2020), 1605–1642. https://dx.doi.org/10.4310/CMS.2020.v18.n6.a5. doi: 10.4310/CMS.2020.v18.n6.a5
![]() |
[4] |
D. F. Anderson, D. Cappelletti, J. Kim, Stochastically modeled weakly reversible reaction networks with a single linkage class, J. Appl. Probab., 57 (2020), 792–810. https://doi.org/10.1017/jpr.2020.28 doi: 10.1017/jpr.2020.28
![]() |
[5] |
D. F. Anderson, S. L. Cotter, Product-form stationary distributions for deficiency zero networks with non-mass action kinetics, Bull. Math. Biol., 78 (2016), 2390–2407. https://doi.org/10.1007/s11538-016-0220-y doi: 10.1007/s11538-016-0220-y
![]() |
[6] |
D. F. Anderson, G. Craciun, T. G. Kurtz, Product-form stationary distributions for deficiency zero chemical reaction networks, Bull. Math. Biol., 72 (2010), 1947–1970. https://doi.org/10.1007/s11538-010-9517-4 doi: 10.1007/s11538-010-9517-4
![]() |
[7] |
B. Pascual-Escudero, L. Hoessly, An algebraic approach to product-form stationary distributions for some reaction networks, SIAM J. Appl. Dynam. Syst., 21 (2022), 588–615. https://doi.org/10.1137/21M1401498 doi: 10.1137/21M1401498
![]() |
[8] |
C. Gadgil, C. H. Lee, H. G. Othmer, A stochastic analysis of first-order reaction networks, Bull. Math. Biol., 67 (2005), 901–946. ttps://doi.org/10.1016/j.bulm.2004.09.009 doi: 10.1016/j.bulm.2004.09.009
![]() |
[9] | D. A. Leven, Y. Peres, Markov Chains and Mixing Times, American Mathematical Society, 107 2017. https://doi.org/10.1007/s00283-018-9839-x |
[10] | D. F. Anderson, T. G. Kurtz, Continuous time Markov chain models for chemical reaction networks, in Design and Analysis of Biomolecular Circuits: Engineering Approaches to Systems and Synthetic Biology (ed. H. K. Et al.), Springer, 2011, 3–42. ttps://doi.org/10.1007/978-1-4419-6766-4_1 |
[11] | D. F. Anderson, T. G. Kurtz, Stochastic analysis of biochemical systems, vol. 1.2 of Stochastics in Biological Systems, 1st edition, Springer International Publishing, Switzerland, 2015. https://doi.org/10.1007/978-3-319-16895-1 |
[12] |
D. T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, J. Comput. Phys., 22 (1976), 403–434. https://doi.org/10.1016/0021-9991(76)90041-3 doi: 10.1016/0021-9991(76)90041-3
![]() |
[13] | D. J. Wilkinson, Stochastic Modelling for Systems Biology, Chapman and Hall/CRC Press, 2006. ttps://doi.org/10.1201/9781351000918 |
[14] | J. Clark, D. A. Holton, A first look at graph theory, World Scientific, 1991. https://doi.org/10.1142/1280 |
[15] | J. Norris, Markov Chains, Cambridge University Press, 1997. https://doi.org/10.1017/CBO9780511810633 |
[16] | S. N. Ethier, T. G. Kurtz, Markov processes: Characterization and convergence, vol. 282, John Wiley & Sons, 2009. https://doi.org/10.1002/9780470316658 |
[17] |
D. T. Gillespie, Exact Stochastic Simulation of Coupled Chemical Reactions, J. Phys. Chem., 81 (1977), 2340–2361. https://doi.org/10.1021/j100540a008 doi: 10.1021/j100540a008
![]() |
[18] |
D. F. Anderson, A modified Next Reaction Method for simulating chemical systems with time dependent propensities and delays, J. Chem. Phys., 127 (2007), 214107. https://doi.org/10.1063/1.2799998 doi: 10.1063/1.2799998
![]() |
[19] |
M. A. Gibson, J. Bruck, Efficient exact stochastic simulation of chemical systems with many species and many channels, J. Phys. Chem. A, 105 (2000), 1876–1889. https://doi.org/10.1021/jp993732q doi: 10.1021/jp993732q
![]() |
[20] |
D. F. Anderson, Incorporating postleap checks in tau-leaping, J. Chem. Phys., 128 (2008), 54103. https://doi.org/10.1063/1.2819665 doi: 10.1063/1.2819665
![]() |
[21] |
D. T. Gillespie, Approximate accelerated simulation of chemically reaction systems, J. Chem. Phys., 115 (2001), 1716–1733. https://doi.org/10.1063/1.1378322 doi: 10.1063/1.1378322
![]() |
[22] |
D. F. Anderson, D. J. Higham, Multi-level Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics, SIAM Mult. Model. Simul., 10 (2012), 146–179. https://doi.org/10.1137/110840546 doi: 10.1137/110840546
![]() |
[23] |
D. F. Anderson, A proof of the global attractor conjecture in the single linkage class case, SIAM J. Appl. Math, 71 (2011), 1487–1508. https://doi.org/10.1137/11082631X doi: 10.1137/11082631X
![]() |
[24] |
D. F. Anderson, Boundedness of trajectories for weakly reversible, single linkage class reaction systems, J. Math. Chem., 49 (2011), 2275–2290. https://doi.org/10.1007/s10910-011-9886-4 doi: 10.1007/s10910-011-9886-4
![]() |
[25] |
D. F. Anderson, D. Cappelletti, J. Kim, T. D. Nguyen, Tier structure of strongly endotactic reaction networks, Stochast. Process. Appl., 130 (2020), 7218–7259. https://doi.org/10.1016/j.spa.2020.07.012 doi: 10.1016/j.spa.2020.07.012
![]() |
[26] | S. P. Meyn, R. L. Tweedie, Stability of Markovian Processes Ⅲ: Foster-Lyapunov Criteria for Continuous-Time Processes, Adv. Appl. Probab., 25 (1993), 518–548, http://www.jstor.org/stable/10.2307/1427522. |
[27] |
R. L. Tweedie, Criteria for ergodicity, exponential ergodicity and strong ergodicity of Markov processes, J. Appl. Prob., 18 (1981), 122–130. https://doi.org/10.2307/3213172 doi: 10.2307/3213172
![]() |
[28] |
A. Athreya, T. Kolba, J. C. Mattingly, Propagating Lyapunov functions to prove noise-induced stabilization, Electron. J. Probab., 17 (2012), 1–38. https://doi.org/10.1214/EJP.v17-2410 doi: 10.1214/EJP.v17-2410
![]() |
[29] |
M. Feinberg, Complex balancing in general kinetic systems, Arch. Rational Mech. Anal., 49 (1972), 187–194. https://doi.org/10.1007/BF00255665 doi: 10.1007/BF00255665
![]() |
[30] |
M. Feinberg, Chemical reaction network structure and the stability of complex isothermal reactors - Ⅰ. The Deficiency Zero and Deficiency One theorems, Review Article 25, Chem. Eng. Sci., 42 (1987), 2229–2268. https://doi.org/10.1016/0009-2509(87)80099-4 doi: 10.1016/0009-2509(87)80099-4
![]() |
[31] |
F. J. M. Horn, R. Jackson, General Mass Action Kinetics, Arch. Rat. Mech. Anal., 47 (1972), 81–116. https://doi.org/10.1007/BF00251225 doi: 10.1007/BF00251225
![]() |
[32] |
D. F. Anderson, G. Craciun, M. Gopalkrishnan, C. Wiuf, Lyapunov functions, stationary distributions, and non-equilibrium potential for reaction networks, Bull. Math. Biol., 77 (2015), 1744–1767. https://doi.org/10.1007/s11538-015-0102-8 doi: 10.1007/s11538-015-0102-8
![]() |
1. | David F. Anderson, Daniele Cappelletti, Wai-Tong Louis Fan, Jinsu Kim, A New Path Method for Exponential Ergodicity of Markov Processes on \({{\mathbb Z^{d}}}\), with Applications to Stochastic Reaction Networks, 2025, 24, 1536-0040, 1668, 10.1137/24M1665933 |