Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Mixing times for two classes of stochastically modeled reaction networks

  • 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 [1] pertaining to models whose state space is restricted to the non-negative integers, there has been a notable lack of results related to this rate of convergence in the reaction network literature. This paper begins the process of filling that hole in our understanding. In this paper, we characterize this rate of convergence, via the mixing times of the processes, for two classes of stochastically modeled reaction networks. Specifically, by applying a Foster-Lyapunov criteria we establish exponential ergodicity for two classes of reaction networks introduced in [2]. Moreover, we show that for one of the classes the convergence is uniform over the initial state.

    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

    Related Papers:

    [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 [1] pertaining to models whose state space is restricted to the non-negative integers, there has been a notable lack of results related to this rate of convergence in the reaction network literature. This paper begins the process of filling that hole in our understanding. In this paper, we characterize this rate of convergence, via the mixing times of the processes, for two classes of stochastically modeled reaction networks. Specifically, by applying a Foster-Lyapunov criteria we establish exponential ergodicity for two classes of reaction networks introduced in [2]. Moreover, we show that for one of the classes the convergence is uniform over the initial state.



    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=inft0{Pt(x,)π()TVε}, (1.1)

    where the total variation distance between two probability measures on a measurable space (Ω,F) is defined as μνTV=supAF|μ(A)ν(A)|. In our case of a discrete state space, μνTV=12x|μ(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:Zd0R0 for which

    Pti(x,)πi()TVBi(x)eηit,for all t0. (1.2)

    Hence, τεx1η[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:Zd0R0, 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 Zd0 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=di=1yiSiandy=di=1yiSi, (2.1)

    and where the values yi,yiZ0 are the stoichiometric coefficients. We will often write reactions (y,y) as yy.

    (iii) The complex set C consists of the linear combinations of the species in (2.1). Specifically, C={y | yyR}{y | yyR}. For the reaction yy, 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)TZd0.

    For example, when S={S1,S2,,Sd}, the complex 2S1+S2 is associated with the vector (2,1,0,0,,0)TZd0. 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 yryrR. 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 yryr).

    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+ESEE+P.

    For this reaction network, S={S,E,SE,P}, C={S+E,SE,E+P} and R={S+ESE,SES+E,SEE+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,vRd0, we define uv=di=1uvii, with 00 taken to be equal to 1.

    2. For xZd0, 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 λyy:Zd0R0 to each reaction yyR. The time evolution of the copy numbers of species is then modeled by means of a continuous-time Markov chain with state space Zd0, whose transition rate from state x to state z is given by

    qx,z=yyRyy=zxλyy(x),

    where the sum is over those reactions whose occurrence causes a net change that is precisely zx. If infinitely many transitions occur within a finite-time T (in other words, an explosion occurs at T), we let X(t)=Δ for any tT, where Δ is a cemetery state not contained in Zd0 [15]. The infinitesimal generator A of the associated Markov process acts on functions via the operation

    Af(x)=yyRλyy(x)(f(x+yy)f(x)), (2.2)

    for f:Zd0R [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 yyR

    λyy(x)=κyyx!(xy)!1{xy}=κyydi=1xi!(xiyi)!1{xiyi}, (2.3)

    for some positive constant κyy, called a reaction rate constant (or rate constant for short). We denote K={κyy}. 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κyyy. 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 y1=di=1yi2 for all yC.

    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 2SiC 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 SR and SR for each SS.

    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 Zd0 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 | SiS} and ˜R=RSiS{Si,Si}. Then, for any choice of reaction rate constants K={κyy} 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,)π()TVC(|x|+1)ln(|x|+2)eηt,for all t0 and all xZd0, (3.1)

    where π is the unique stationary distribution of the model on Zd0. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|.

    Moreover, if there is a wRd>0 for which w(yy)=0 for all yyR, 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,)π()TVC(|x|+1)eηt,for all t0 and all xZd0. (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=12z|Pt(x,z)π(z)|12z[0,200]3|Pt(x,z)π(z)|.
    Figure 1.  A and B (Left). For each initial condition X(0)=(m,m,m), the images on the left provide plots of Pt(x,)π()TV over time for the reaction networks (3.3) (top) and (3.4) (bottom), where π represents their stationary distributions. A and B (Right). The images on the right provide plots of the mixing times τεxm with ε=0.1 and xm=(m,m,m), as functions of m. The top image is for the reaction network (3.3) whereas the bottom image is for the the reaction network (3.4). Note that the top image appears to show logarithmic growth whereas the bottom appears to show an asymptotic upper bound. Both observations agree with our theory.

    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 Zd0 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={κyy} be a choice of reaction rate constants. Let X(0)=xSZd0 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,)π()TVCeηt,for all t0,

    where neither C nor η depend upon xS. Thus, the mixing time satisfies τεx=O(1), as |x|.

    We provide an example.

    Example 3.2. Consider the following network

    2AA+BB,A2CB+C.2B,CA+C (3.4)

    This network is double full since 2A,2B,2CC. 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:

    2AA+BB2CA2B

    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=12z|Pt(x,z)π(z)|12z[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 xZd0 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(x1), which is the vector whose jth component is

    (x1)j=xj1=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 Rd0 is called a D-type tier-sequence if

    (i) for each i{1,,d} the limit limnxn,i exists (it could be infinity) and limnxn,i= for at least one i, and

    (ii) for any pair of complexes y,yC, the limit

    limn(xn1)y(xn1)y,

    exists (it could be infinity).

    Remark 4.1. Note that, given a sequence {xn} of vectors in Rd0 with lim supnxn=, 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 Rd0. 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 yTD,1{xn}) if for all complexes yC

    limn(xn1)y(xn1)y>0;

    (ii) we say that a complex y is in tier i (and write yTD,i{xn}) if there exists yTD,i1{xn} with

    limn(xn1)y(xn1)y=0

    and for all complexes yi1j=1TD,j{xn} we have

    limn(xn1)y(xn1)y>0.

    The mutually disjoint subsets TD,i{xn} are called D-type tiers along {xn}. If yTD,i{xn} and yTD,j{xn} with i<j we will write yDy. If y,y are in the same D-type tier, then we will write yDy.

    Thus, those complexes in tier 1 maximize (xn1)y along the sequence, with those in tier 2 being second largest, etc. Note that as a consequence of the above definition, if y,yTi{xn} then there exists C(0,) such that

    limn(xn1)y(xn1)y=C,

    explaining the notation yDy, and if yTi{xn} and yTj{xn} for some i<j, then

    limn(xn1)y(xn1)y=0,

    explaining the notation yDy.

    The following example demonstrates the concept.

    Example 4.1. Consider the reaction graph below.

    2A    BA. A+B

    Let C be enumerated with y1=2A,y2=B,y3=,y4=A,y5=A+B. For xn=(n,0), (xn1)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 (xn1)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 xnZd0, be a D-type tier-sequence of a reaction network (S,C,R) and let yyR. If limnλyy(xn)>0, where λyy is given in (2.3), then there is a C>0 with

    limnλyy(xn)(xn1)y=C.

    Moreover, if the source complex y satisfies |yi|1 for each i{1,,d}, then C=κyy.

    Proof. Since limnλyy(xn)>0, we have that (xn1)y=xyn for all n large enough. The result then follows from the basic polynomial structure of λyy.

    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 SZd 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  xS. (4.1)

    Then there exists positive constants η and C such that

    Pt(x,)π()TVC(V(x)+1)eηt,xS.

    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 τεx1η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 SRd 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  xS. (4.2)

    Then there exist positive constants η and C such that

    Pt(x,)π()TVCeηt.

    Note that in this case we have that the mixing time satisfies τεx=O(1).

    We let V:Zd0R0 be the following function

    V(x)=di=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 {κyy}. 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)yyRλyy(xn)(ln((xn1)y(xn1)y)+C1). (5.2)

    Furthermore, suppose that there exists a reaction y1y1R such that

    (i) y1TD,1{xn},

    (ii) y1y1, and

    (iii) limnλy1y1(xn)(xn1)y1=C2, for some C2>0.

    Then there exists a positive constant C3 such that for large enough n,

    AV(xn)C3yyRyDyλyy(xn)ln((xn1)y(xn1)y), (5.3)

    where the sum is over those reactions, yyR, with yDy.

    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 | SiS} and ˜R=RSiS{Si,Si}. Then, for any D-type tier-sequence {xn} for (S,˜C,˜R), there exists a reaction y1y1˜R such that

    (i) y1TD,1{xn},

    (ii) y1Dy1 and

    (iii) limnλy1y1(xn)(xn1)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 SiD.

    Next, we observe that if SiTD,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 SiTD,1{xn}. Then it must be the case that Si+SjTD,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, ˜yC, with |˜y|{0,1} (i.e., it is either unary or the zero complex). In particular, we know that ˜yTD,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, yyR, with yTD,1{xn} and yDy. 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 limnxn,k>0 and limnxn,>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 limnxn,k>0 and limnxn,=0. However, this would imply that SkTD,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 2AB satisfies all the conditions (ⅰ), (ⅱ) and (ⅲ) in Lemma 5.2. For the D-type tier sequence xn=(0,n), the reaction A+B2A satisfies (ⅰ) and (ⅱ) in Lemma 5.2. However (ⅲ) does not hold because the rate function λA+B2A 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 wRd>0 with w(yy)=0 for all yyR. Let W(x)=wx=di=1wixi. Note that lim|x||W(x)|= since wRd>0. For each i{1,,d}, let eiZd0 be the vector with a 1 in the ith coordinate and zeros elsewhere. Since w conserves the reactions in R, we have

    AW(x)=yyRλyy(x)(W(x+yy)W(x))+di=1λSi(x)(W(xei)W(x))+di=1λSi(x)(W(x+ei)W(x))=di=1κSiwixi+di=1κSiwiaW(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:SiS} nor C{Si:SiS}. We claim that there exists a compact set KS such that for some positive constant a, the function V defined in (5.1) satisfies

    AV(x)aV(x)for all  xSK. (5.4)

    If this claim holds, then the result follows from Theorem 4.2 after setting b=(a+1)maxxKAV(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

    limnxn,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,iCxn,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)C1yyRyDyλyy(xn)ln((xn1)y(xn1)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((xn1)(0,0,,0)(xn1)(1,0,,0))=C1λS1(xn)ln(xn,11).

    However, since we assumed xn,1 is a maximal coordinate of xn, and limnxn,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,11).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 Zd0.

    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 yyR with yTD,1{xn}, we have that λyy(xn)>0 for all n large enough,

    (ii) 2SiTD,1{xn} for some i, and

    (iii) There exists yyR such that yTD,1{xn} and yy.

    Proof. There are no complexes with y1>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 KS, together with an a>0, for which

    AV(x)aV(x)3/2for all  xSK.

    Then there exists a sequence {xn}, with xnS, so that

    limnxn=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 2SiTD,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 y1y1R for which y1TD,1{xn} and y1y1. 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λy1y1(xn)ln((xn1)y1(xn1)y1)C2λy1y1(xn), (5.9)

    where the second inequality holds since ln((xn1)y1(xn1)y1), as n. Since we know that 2SiTD,1{xn}, we can then conclude that there is a C>0 so that

    AV(xn)C2λy1y1(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 Zd0 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 | SiS} and ˜R=RSiS{Si}SiI{Si}, where IS is any subset of the species (including the empty set). Let K={κyy} be a choice of reaction rate constants. Let X(0)=xZd0 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,)π()TVC(|x|+1)ln(|x|+2)eηt,for all t0,

    where neither C nor η depend upon xS. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|.

    Suppose now that there is a wRd>0 for which w(yy)=0 for all yyR. Let X(0)=xZd0 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,)π()TVC(|x|+1)eηt,for all t0,

    where neither C nor η depend upon xS. 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 Zd0 is irreducible, and there is a unique stationary distribution on Zd0. 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 yTD,1{xn} and ending with a complex yTD,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 2S1S1 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 yC, 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 | SiS} and ˜R=RSiS{Si}SiI{Si}, where IS is any subset of the species. Let K={κyy} be a choice of reaction rate constants. Let X(0)=xZd0 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,)π()TVC(|x|+1)ln(|x|+2)eηt,for all t0,

    where neither C nor η depend upon xS. Thus, the mixing time satisfies τεx=O(ln(|x|)), as |x|.

    Suppose now that there is a wRd>0 for which w(yy)=0 for all yyR. Let X(0)=xZd0 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,)π()TVC(|x|+1)eηt,for all t0,

    where neither C nor η depend upon xS. 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 yTD,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.

    2BA+BC,2CAB.

    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

    2BA+B  C2CA 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 SiSj) 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 yC, 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 SiC such that i1, 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 SiSi1Si2S1).

    Let ˜C=C{} and ˜R=R{S1}SiI{Si}, where IS is any subset of the species. Let K={κyy} be a choice of reaction rate constants. Let X(0)=xZd0 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,)π()TVC(|x|+1)ln(|x|+2)eηt,for all t0,

    where neither C nor η depend upon xS. Therefore the mixing time of the model satisfies τεx=O(ln(|x|)), as |x|.

    Suppose now that there is a wRd>0 for which w(yy)=0 for all yyR. Let X(0)=xZd0 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,)π()TVC(|x|+1)eηt,for all t0,

    where neither C nor η depend upon xS. 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=STD,1{xn}. Then on a direct path SSi1S1, there must exist a reaction of the form either SjSk for some jk such that SjTD,1{xn} and SkTD,1{xn} or S1 with S1TD,1{xn} and TD,1{xn} (as always). Then letting this reaction be y1y1, 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 yy such that y=Sr for some r and yDy 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=Si0Si1Si2S1, we have that either i) there exists SrSmR for some r,m such that xn,r is a maximal coordinate and SrDSm or ii) S1 satisfies that xn,1 is a maximal coordinate and S1D. 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 yC, 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 SpS and a non-empty subset SeS such that (i) S=SpSe, and (ii) for each SiSp there is a directed path within the reaction graph such that it begins with Si and ends with some SkSe. Furthermore the directed path consists only of unary complexes (i.e. the path is of the form SiSi1Si2Sk).

    Let ˜C=C{Si | SiSe}{} and ˜R=RSSe{S}SiI{Si}, where IS is any subset of the species. Let K={κyy} be a choice of reaction rate constants. Let X(0)=xZd0 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,)π()TVC(|x|+1)ln(|x|+2)eηt,for all t0,

    where neither C nor η depend upon xS. Therefore the mixing time of the model satisfies τεx=O(ln(|x|)), as |x|.

    Suppose now that there is a wRd>0 for which w(yy)=0 for all yyR. Let X(0)=xZd0 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,)π()TVC(|x|+1)eηt,for all t0,

    where neither C nor η depend upon xS. 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 SkSe. 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:

    SPE↓↑S+E SEE+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].

    Figure 2.  Plots for the model provided in Example 6.4 in (6.1). For each initial condition X(0)=(m,m,m), the image on the left provides a plot of Pt(x,)π()TV over time for the reaction network (6.1), where π represents the stationary distribution. The image on the right provides a plot of the mixing times τεxm with ε=0.1 and xm=(m,m,m), as a function of m. Note that the growth on the right appears to be logarithmic, which agrees with our theory.

    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
  • This article has been cited by:

    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
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2035) PDF downloads(94) Cited by(1)

Figures and Tables

Figures(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog