1. Introduction
Reaction networks are directed graphs that can be used to describe the general dynamics of population models. Their use is widespread in population ecology, chemistry, and cell biology. An example of a reaction network is
\begin{equation} \includegraphics{1.eps}\end{equation}
where the possible transformations of the species A, B, and C are depicted: in the reaction 
$A\to B+C$
, a single individual of species A can be transformed into one individual of species B and one individual of species C. In the reaction 
$B+C\to 0$
, individuals of species B and C can annihilate each other. In the opposite reaction 
$0 \to B+C$
, one individual of B and one individual of C are introduced to the system at the same time, and so on. Typically, in biochemistry the species are thought of as different chemical components that react with each other. The nodes of the graph are termed complexes. Different mathematical models can be associated with such a network in order to describe the dynamics of the system of interest. Regardless of the modeling choice, a natural question is the following: when is it possible to prove theorems that relate the structure of the reaction diagrams, which are developed by domain scientists, to general dynamical properties of the associated mathematical models? The subfield of mathematics that works on this question is often termed ‘chemical reaction network theory’, and has a history dating back to at least the early 1970s [Reference Feinberg19, Reference Horn22, Reference Horn and Jackson23]. In particular, there has been an enormous amount of work relating the possible dynamics of the associated deterministic mass-action system with the basic structure of the reaction network itself; for just a subset, see [Reference Anderson3, Reference Cappelletti and Wiuf15, Reference Craciun, Dickenstein, Shiu and Sturmfels16, Reference Craciun and Feinberg17, Reference Gopalkrishnan, Miller and Shiu20, Reference Karp, Pérez Millán, Dasgupta, Dickenstein and Gunawardena26].
There has also been work in this sense for stochastic dynamics, which are typically modeled through a continuous-time Markov chain as described in Section 3. This paper is concerned with providing structural conditions on the reaction network that guarantee that the associated continuous-time Markov chain is positive recurrent. This information is important for practical purposes; as an example, in biological models it is often the case that different chemical species evolve over different time scales. It is known that under certain assumptions, the dynamics of the slower species can be approximated by averaging out the dynamics of the faster subsystem, which is considered to be a stationary regime [Reference Ball, Kurtz, Popovic and Rempala11, Reference Kang and Kurtz25, Reference Pfaffelhuber and Popovic31]. It is therefore important to understand whether the fast subsystem admits a stationary regime in the first place. Once the existence of the stationary distribution is proven, this can be approximated by simulations, or by finite state-space projections [Reference Anderson and Ehlert9, Reference Gupta, Mikelson and Khammash21, Reference Munsky and Khammash28]. Another example concerns synthetic biology; in this field, biological circuits are designed and implemented in order to control cellular behavior. Often, it is rigorously proven that such circuits serve the desired purpose provided that a stationary regime can be reached [Reference Briat, Gupta and Khammash12, Reference Briat, Gupta and Khammash13]. Finally, from a mathematical standpoint it is inherently interesting to understand the recurrence properties of this large class of models, which are used ubiquitously in the science literature.
Despite the importance of understanding whether a stationary regime exists, checking for positive recurrence in the setting of reaction networks is usually a hard task, with the exception of a few classes of models for which positive recurrence has been shown [Reference Anderson, Cappelletti, Koyama and Kurtz5, Reference Anderson and Cotter6, Reference Anderson, Craciun and Kurtz8, Reference Anderson and Kim10, Reference Cappelletti and Wiuf14, Reference Jahnke and Huisinga24]. Probably the most important open conjecture related to stochastic models in the field is the following: if the reaction network is weakly reversible (i.e. if each connected component of the reaction graph is strongly connected; see Definition 3.2), then the associated Markov model is positive recurrent. There is a large amount of evidence that this conjecture is true. In particular, in [Reference Anderson, Craciun and Kurtz8] it was shown that if the deterministic model associated with a particular network is complex balanced (meaning that at equilibrium the total flux into each complex is equal to the total flux out), then the stochastic model associated with that same network admits a stationary distribution that is a product of Poissons. Follow-up work in [Reference Anderson, Cappelletti, Koyama and Kurtz5] showed that these models are non-explosive, implying they are positive recurrent. A key fact is that for every weakly reversible reaction network, there is a choice of rate constants that makes the model complex balanced. Hence, if there is a weakly reversible model that is not positive recurrent, the lack of recurrence must be a property of the specific choice of rate constants, and not of the network structure, something most researchers find unlikely.
In this paper we do not consider the full class of weakly reversible networks, but instead restrict ourselves to those networks that are weakly reversible and have a single connected component (termed a ‘linkage class’ in the literature). This particular set of networks was first studied in the deterministic context in [Reference Anderson2, Reference Anderson3], where it was shown that they necessarily have bounded, persistent trajectories. That is, for each trajectory 
$x(t)\in {\mathbb R}^d_{\ge 0}$
, the expression 
$\sum_{i=1}^d |\ln(x_i(t))|$
is uniformly bounded in time. This subclass of models was studied in the stochastic context in [Reference Anderson and Kim10], where positive recurrence was proven under a few additional assumptions on the network, which we will not list here due to their technical nature.
We prove in Theorem 4.1 that models that are weakly reversible and consist of a single linkage class are positive recurrent under two additional structural assumptions: (i) the system is binary (see Definition 3.3), and (ii) for each species, there is a complex that is a multiple of that species. As an example, these assumptions are met in (1).
Providing a new class of models for which positive recurrence is guaranteed is the first major contribution of this paper. The second is the method of proof, which we believe will be useful in the study of a significantly wider class of models. In particular, most previous results related to positive recurrence of stochastically modeled reaction networks that the authors are aware of utilized the Foster–Lyapunov criteria related to the continuous-time Markov chain itself. In this article, we prove positive recurrence of the continuous-time chain by studying the recurrence properties of the discrete-time Markov chain that arises by observing the continuous-time Markov chain after a fixed number of jumps.
The paper has the following outline. In Section 2 we briefly provide a catalog of notation we use throughout the paper. In Section 3, we provide the basic mathematical objects we will study. In Section 4 we state our main result, and in Section 5 we introduce the concept of tier sequences, which is fundamental to showing our main result. The connection between tiers and positive recurrence is detailed in Section 6, at the end of which the proof of the main structural result is given.
2. Notation
We denote by 
${\mathbb R}$
, 
${\mathbb R}_{>0}$
, and 
${\mathbb R}_{\ge 0}$
the real, positive, and non-negative real numbers, respectively. Similarly, we denote by 
${\mathbb Z}$
, 
${\mathbb Z}_{>0}$
, and 
${\mathbb Z}_{\ge 0}$
the integers, positive integers, and non-negative integers. Any sum of the form 
$\sum_{i = 1}^0 a_i$
is taken to be zero, as is usual.
Given two vectors 
$x,z\in{\mathbb R}^d$
, we write 
$x\geq y$
if 
$x_i\geq z_i$
for all 
$1\leq i\leq d$
. When 
$x,y \in {\mathbb R}^d_{\ge0}$
, we denote by 
$x^y$
the vector whose ith component is given by 
$x_i^{y_i}$
, with the usual convention that 
$0^0$
is taken to be 1. Moreover, we denote by 
$x\vee 1$
the vector in 
${\mathbb R}_{>0}^d$
whose ith component is 
$\max\{x_i, 1\}$
.
Given a stochastic process 
$\{X(t)\,:\,t\in{\mathbb R}_{\geq0}\}$
with state space 
$\Gamma$
, and 
$x\in\Gamma$
, we denote
\begin{equation} \mathbb{P}_x(X(t)\in A)=\mathbb{P}(X(t)\in A \mid X(0)=x) , \ \mathbb{E}_x[f(X(t))]=\mathbb{E}[f(X(t)) \mid X(0)=x]\end{equation}
for all measurable sets A, all measurable functions 
$f\colon\Gamma\to{\mathbb R}$
(with respect to the 
$\sigma$
-algebra generated by X), and all 
$t\in{\mathbb R}_{>0}$
.
We assume basic knowledge of Markov chain theory. In particular, we will not define what the generator of a Markov chain is, what it means for a Markov chain to have a stationary distribution, and what it means for a Markov chain to be positive recurrent, transient, or explosive. These concepts can be reviewed in [Reference Ethier and Kurtz18, Reference Norris29]. We will also assume that the Foster–Lyapunov function criterion for continuous-time and discrete-time Markov chains is known to the reader [Reference Meyn and Tweedie27].
3. Reaction networks and stochastic dynamics
In Section 3.1 we formally introduce reaction networks. In Section 3.2 we introduce the associated stochastic models. In Section 3.3 we introduce the embedded discrete-time Markov chain for our model.
3.1. Reaction networks
A reaction network describes the set of possible interactions among constituent chemical species, and is defined as follows.
Definition 3.1. A reaction network is given by a triple of finite sets 
$({\mathcal S},{\mathcal C},{\mathcal R})$
where:
(i) the set
${\mathcal S}=\{S_1,S_2,\ldots,S_d\}$
is a set of d symbols, called the species of the network;(ii) the set
${\mathcal C}$
is a set of linear combinations of species on
${\mathbb Z}_{\geq 0}$
, called complexes;(iii) the set
${\mathcal R}$
is a non-empty subset of
${\mathcal C}\times{\mathcal C}$
, called reactions, whose elements (y, y ′) are denoted by
$y\to y^{\prime}$
. The cardinality of
${\mathcal R}$
is denoted by r.
Given a reaction 
$y\to y^{\prime}\in{\mathcal R}$
, the complexes y and y ′ are termed the source and product complex of the reaction, respectively. It is common to assume that 
$y\to y\notin{\mathcal R}$
for all 
$y\in{\mathcal C}$
, and we will do so here. Throughout, we abuse notation slightly by letting a complex 
$y\in{\mathcal C}$
denote both a linear combination of the form
\begin{equation*}y=\sum_{i=1}^d y_iS_i\end{equation*}
as defined in Definition 3.1, and the vector whose ith component is 
$y_i$
, i.e. 
$y=(y_1,y_2,\ldots,y_d)^\top \in \mathbb{Z}^d_{\ge 0}$
. For example, if 
$y=2S_1+S_2$
then y will also be used to denote the vector 
$(2,1,0,0,\dots,0) \in {\mathbb Z}^d_{\ge 0}$
. We denote by 0 the complex y for which 
$y_i=0$
for each i.
The directed graph whose nodes are 
${\mathcal C}$
and whose directed edges are 
${\mathcal R}$
can be associated with a reaction network. Such a graph is called the reaction graph of the reaction network.
Definition 3.2. Let 
$({\mathcal S},{\mathcal C},{\mathcal R})$
be a reaction network. The connected components of the associated reaction graph are termed linkage classes. We say that a linkage class is weakly reversible if it is strongly connected, i.e. if for any two complexes 
$y_1,y_2$
in the linkage class there is a directed path from 
$y_1$
to 
$y_2$
and a directed path from 
$y_2$
to 
$y_1$
.
If all linkage classes in a reaction network are weakly reversible, then the reaction network is said to be weakly reversible.
Example 1. Consider the reaction network with the following reaction graph:
\begin{equation*}\includegraphics{2.eps}\end{equation*}
This network has three linkage classes. The right-most linkage class is weakly reversible, whereas the other two are not. Because not all linkage classes are weakly reversible, the network is not weakly reversible.
The main result of this paper pertains to binary reaction networks.
Definition 3.3. A reaction network 
$({\mathcal S},{\mathcal C},{\mathcal R})$
is called binary if 
$\sum_{i=1}^d y_i \le 2$
for all 
$y \in {\mathcal C}$
.
Most networks found in the science literature are binary.
3.2. Stochastic model
In this section we introduce a stochastic dynamical model associated with reaction networks, which is often utilized in biochemistry if few molecules of the different species are present and stochastic fluctuations cannot be ignored. First, we associate with each reaction 
$y\to y^{\prime}$
a rate function 
$\lambda_{y\to y^{\prime}}\colon {\mathbb Z}^d\to {\mathbb R}_{\geq0}$
, expressing the propensity of a reaction to occur. We further assume that a reaction 
$y\to y^{\prime}$
cannot take place if not enough reacting molecules are present, which is expressed by the condition 
$\lambda_{y\to y^{\prime}}(x)>0$
only if 
$x\geq y$
. A popular choice of rate functions is given by mass action kinetics, where it is assumed that the propensity of a reaction to occur is proportional to the number of combinations of the present molecules that can give rise to the reaction. Specifically, mass action rate functions take the form
\begin{align*} \lambda_{y\rightarrow y^{\prime}}(x)= \kappa_{y\rightarrow y^{\prime}} \lambda_y(x), \qquad\text{where}\ \ \ \lambda_y(x)=\prod_{i=1}^d \frac{x_i !}{(x_i-y_{i})!}\textbf{1}_{\{x_i \ge y_{i}\}},\end{align*}
for some positive constants 
$\kappa_{y\to y^{\prime}}$
called rate constants, and where 
$\textbf{1}_{\{x_i \ge y_{i}\}}$
is the indicator function on the event 
$\{x_i \ge y_i\}$
. A reaction network together with a choice of rate constants 
$({\mathcal S}, {\mathcal C}, {\mathcal R}, \kappa)$
is called a mass action system. We say that a mass action system is weakly reversible, binary, or has a single linkage class if the related reaction network does.
For every 
$t\geq0$
, we denote by X(t) a vector in 
${\mathbb Z}^d_{\geq0}$
whose ith entry represents the count of species 
$S_i$
at time t. We assume that X is a continuous-time Markov chain over 
${\mathbb Z}^d_{\geq0}$
, with transition rate from state 
$x\in {\mathbb Z}^d_{\geq0}$
to state 
$z \in {\mathbb Z}^d_{\geq0}$
given by
\begin{equation} q(x, z)=\sum_{\substack{y\to y^{\prime}\in {\mathcal R}\\ y^{\prime}-y=z-x}}\lambda_{y\to y^{\prime}}(x).\end{equation}
The generator 
$\mathcal{A}$
of the process X is the operator whose action on a function 
$V\,:\,{\mathbb Z}^d_{\ge 0} \to {\mathbb R}$
is given by [Reference Ethier and Kurtz18]
\begin{align*} \mathcal{A}V(x) = \sum_{y \rightarrow y^{\prime} \in {\mathcal R}} \lambda_{y\rightarrow y^{\prime}}(x)(V(x+y^{\prime}-y)-V(x)).\end{align*}
Moreover, assume 
$X(0)=x_0$
. We then define the set of reachable states as
\begin{equation*}{\mathbb{S}_{x_0}}=\{x\in{\mathbb Z}^d_{\geq0}\,:\,\mathbb{P}_{x_0}(X(t)=x)>0\ \text{for some }t\geq0\},\end{equation*}
where 
$\mathbb{P}_{x_0}$
is as in (2).
3.3. The embedded discrete-time Markov chain
Let the continuous-time Markov chain X be as defined in the previous section. In the present paper we will utilize the embedded discrete-time Markov chain of X, defined as follows [Reference Norris29].
Definition 3.4. Let 
$\tau_0=0$
and for 
$n \in {\mathbb Z}_{>0}$
let 
$\tau_n$
be the time of the nth jump of X. Then the embedded discrete-time Markov chain of X is the discrete-time Markov chain 
$\{\tilde{X}_n\}_{n=0}^\infty$
defined by 
$\tilde X_n=X(\tau_n)$
for all 
$n\in{\mathbb Z}_{\geq0}$
.
It is straightforward to show that the transition probabilities of the embedded discrete-time Markov chain 
$\tilde{X}$
are given by
\begin{equation}\mathbb{P}\big(\tilde{X}_1=z \mid \tilde{X}_0=x\big)=\frac{q(x,z)}{\overline{\lambda}(x)},\end{equation}
where q is as in (3) and
\begin{equation} \overline{\lambda}(x)=\sum_{y\to y^{\prime}\in {\mathcal R}}\lambda_{y\to y^{\prime}}(x). \end{equation}
The following result, which can be inferred from [Reference Norris29, Theorem 3.5.1], will be used.
Theorem 3.1. Suppose that 
$\inf_{x\in {\mathbb{S}_{x_0}}}\overline{\lambda}(x) > 0$
and that 
$\tilde{X}$
restricted to 
${\mathbb{S}_{x_0}}$
is positive recurrent. Then X restricted to 
${\mathbb{S}_{x_0}}$
is non-explosive and positive recurrent.
The condition 
$\inf_{x\in {\mathbb{S}_{x_0}}}\overline{\lambda}(x) > 0$
holds when the process X is associated with a weakly reversible mass action system, and the initial condition 
$X(0)=x_0$
is not an absorbing state. This is expressed by the following result, which can be derived from [Reference Paulevé, Craciun and Koeppl30, Lemma 4.6]. Note that the word ‘recurrent’ in [Reference Paulevé, Craciun and Koeppl30] is not to be intended in the usual stochastic sense, which is considered here.
Lemma 3.1. Let 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
be a weakly reversible mass action system, and assume 
$\overline\lambda(x_0)>0$
. Then 
$\inf_{x\in {\mathbb{S}_{x_0}}}\overline{\lambda}(x) > 0$
.
4. Main result
Our main result is the following.
Theorem 4.1. Consider a weakly reversible, binary mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
that has a single linkage class, and let X be the associated continuous-time Markov chain. Assume that for each 
$S\in {\mathcal S}$
, 
$\{S, 2S\} \cap {\mathcal C}\neq \emptyset$
. Then X is positive recurrent.
To prove Theorem 4.1 we will use Foster–Lyapunov functions, the notion of tiers (explained and defined in the next section), and properties of an associated n-step embedded discrete-time Markov chain. Specifically, consider the entropy-like function 
$V\colon {\mathbb R}^d_{\geq0}\to {\mathbb R}_{\geq0}$
defined by
\begin{equation}V(x) = \sum_{i=1}^d \big( x_i(\ln{x_i}-1)+1\big),\end{equation}
with the convention 
$0 \ln{0}=0$
. Notably, this function has been shown to be a Lyapunov function for the deterministic models of complex balanced mass action systems [Reference Horn22, Reference Horn and Jackson23], and, in connection with tier sequences, for those models that are weakly reversible and have a single linkage class [Reference Anderson2, Reference Anderson3]. It has been further utilized to study the stability of stochastic models of mass action systems in different works, including [Reference Agazzi, Dembo and Eckmann1, Reference Anderson, Cappelletti, Kim and Nguyen4, Reference Anderson, Craciun, Gopalkrishnan and Wiuf7, Reference Anderson and Kim10]. In this paper, we will show that V can be utilized as a Foster–Lyapunov function for the n-step discrete-time embedded chain, for a specially chosen positive integer n, even though it may not be a Foster–Lyapunov function for the continuous-time Markov chain X itself.
In the following example we show how Theorem 4.1 can be applied, and how the function V cannot be readily used as a Foster–Lyapunov function for X.
Example 2. Consider the mass action system
\begin{equation*}\includegraphics{3.eps}\end{equation*}
where each rate constant is placed next to its corresponding reaction. The system is weakly reversible and has a single linkage class, and 
$A,2B,C\in {\mathcal C}$
. Then, Theorem 4.1 applies and the associated Markov process X is positive recurrent for any choice of rate constants. To the best of our knowledge, there are no other results in the literature that allow one to quickly reach this conclusion for this model, and finding a tailored Foster–Lyapunov function may not be easy. In particular, the function V as defined in (6) is not a Foster–Lyapunov function for X. To show this, it suffices to show that there exists an initial condition 
$X(0)=x_0$
and a sequence of points 
$\{x_n\}$
in 
${\mathbb{S}_{x_0}}$
whose norm tends to infinity and for which
\begin{equation*}\limsup_{n\to\infty}\mathcal{A}V(x_n)>0.\end{equation*}
To this end, consider 
$x_0=(1,0,0)$
and the sequence defined by 
$x_n=(n,1,0)$
. Note that by sequential occurrences of reactions 1 and 5, we may conclude that 
$x_n \in {\mathbb{S}_{x_0}}$
. Moreover, we have
\begin{equation*}\mathcal{A}V(x_n)=\kappa_1n(2\ln{2}-1),\end{equation*}
which tends to infinity as n tends to infinity.
We will give here a very high-level motivation for the arguments that will follow. First, loosely speaking, the assumption that either 
$S \in {\mathcal C}$
or 
$2S \in {\mathcal C}$
ensures that, for any sequence of points 
$x_n$
going to 
$\infty$
, one of the source complexes with largest intensities is of the form S or 2S. The intensity functions for these reactions are ‘smooth’ away from the origin, in the sense that they cannot change dramatically due to a single instance of a reaction. This is unlike, say, the complex 
$S_1+S_2$
, whose intensity can change dramatically through small changes, as seen by considering the values (n, 1) and (n, 0). Second, the weak reversibility assumption is helpful as it guarantees that each product complex is also a source complex. Hence, whenever a reaction has the largest intensity, the reactions starting from the product complex are either of the same order of magnitude, or smaller. By following a sequence of such reactions (leading to the analysis of the n-step chain), we are guaranteed to eventually find a source complex whose reaction rates are of a lower order of magnitude, and we will show how this guarantees positive recurrence by the Foster–Lyapunov criterion.
5. Tiers
In this section we introduce both D-type and S-type tier sequences, which have now appeared in the literature a number of times [Reference Anderson2, Reference Anderson3, Reference Anderson, Cappelletti, Kim and Nguyen4, Reference Anderson and Kim10], and state some known relevant results. Typically, tier sequences are used in connection with the function V as described in (6). The usual strategy consists of proving that along any tier sequence 
$\{x_n\}$
we have
\begin{equation*}\limsup_{n\to\infty}\mathcal{A} V(x_n)\leq -1,\end{equation*}
implying positive recurrence by the Foster–Lyapunov criterion [Reference Meyn and Tweedie27] and from the fact that any sequence of states has a tier subsequence.
The above approach does not always work. In particular, Example 2 shows how it may fail for models in the class described in Theorem 4.1. Hence, in order to prove Theorem 4.1, in Section 5.3 we will introduce a new concept of tiers along sequences of paths. With this new concept, we will be able to prove that even if 
$\limsup_{n\to \infty} \mathcal{A}V(x_n) >0$
for a particular tier sequence 
$\{x_n\}$
, positive recurrence can still be shown by consideration of how V changes over a fixed number of transitions.
5.1. Definitions
In this section we formally introduce the concept of ‘tiers’. Tiers are partitions of the set of complexes, and they are based on the different orders of magnitude of the different rate functions along a sequence of points. We define two different types of tiers, specifically D-type and S-type tiers. We will describe the role of these objects after their formal definitions.
Definition 5.1. Let 
$( {\mathcal S},{\mathcal C},{\mathcal R})$
be a reaction network and let 
$\{x_n\}$
be a sequence in 
${\mathbb Z}^d_{\ge 0}$
, with n varying from 1 to 
$\infty$
. We say that 
${\mathcal C}$
has a D-type partition along 
$\{x_n\}$
if there exists a partition 
$\big\{T^{{\rm D},1}_{\{x_n\}},\dots, T^{{\rm D},P}_{\{x_n\}}\big\}$
of 
${\mathcal C}$
such that:
(i) if
$y,y^{\prime} \in T^{{\rm D},i}_{\{x_n\}}$
, then
exists, and
\begin{equation*}C=\lim_{n\rightarrow \infty} \frac{(x_n \vee 1)^{y}}{(x_n\vee 1)^{y^{\prime}}}\end{equation*}
$0<C<\infty$
;(ii) if
$y \in T^{{\rm D},i}_{\{x_n\}}$
and
$y^{\prime} \in T^{{\rm D},k}_{\{x_n\}}$
with
$i < k$
then
\begin{align*}\lim_{n\rightarrow \infty} \dfrac{(x_n \vee 1)^{y^{\prime}}}{(x_n\vee 1)^{y}} = 0.\end{align*}
The sets of a D-type partition are called D-type tiers along 
$\{x_n\}$
. We will say that y is in a higher D-type tier than y’ along 
$\{x_n\}$
if 
$y \in T^{{\rm D},i}_{\{x_n\}}$
and 
$y^{\prime} \in T^{{\rm D},j}_{\{x_n\}}$
with 
$i < j$
. In this case, we will write 
$y \succ_{\rm D} y^{\prime}$
. If y and y’ are in the same D-type tier, then we will write 
$y\sim_{\rm D} y^{\prime}$
.
Definition 5.2. Let 
$({\mathcal S},{\mathcal C},{\mathcal R})$
be a reaction network and let 
$\{x_n\}$
be a sequence in 
${\mathbb Z}^d_{\ge 0}$
, with n varying from 1 to 
$\infty$
. We say that 
${\mathcal C}$
has an S-type partition along 
$\{x_n\}$
if there exists a partition 
$\big\{T^{{\rm S},1}_{\{x_n\}}, \ldots , T^{{\rm S},P}_{\{x_n\}}, T^{\rm S,\infty}_{\{x_n\}}\big\}$
of 
${\mathcal C}$
such that:
(i)
$y \in T^{{\rm S},\infty}_{\{x_n\}}$
if and only if
$\lambda_{y}(x_n) = 0$
for all
$n\geq 1$
;(ii)
$y \in \bigcup_{i=1}^P T^{{\rm S},i}_{\{x_n\}}$
if and only if
$\lambda_y(x_n) \ne 0$
for all
$n\geq1$
;(iii) if
$y,y^{\prime} \in T^{{\rm S},i}_{\{x_n\}}$
, with
$i \in \{1,\dots,P\}$
, then
exists, and
\begin{align*}C=\lim_{n\rightarrow \infty} \dfrac{\lambda_{y^{\prime}}(x_n)}{\lambda_{y}(x_n)}\end{align*}
$0<C<\infty$
;(iv) if
$y \in T^{{\rm S},i}_{\{x_n\}}$
and
$y^{\prime} \in T^{{\rm S},k}_{\{x_n\}}$
with
$1\le i < k\le P$
, then
\begin{align*}\lim_{n\rightarrow \infty} \dfrac{\lambda_{y^{\prime}}(x_n)}{\lambda_{y}(x_n)} = 0.\end{align*}
The sets of an S-type partition are called S-type tiers along 
$\{x_n\}$
. We will say that y is in a higher S-type tier than y ′ along 
$\{x_n\}$
if 
$y \in T^{{\rm S},i}_{\{x_n\}}$
and 
$y^{\prime} \in T^{{\rm S},j}_{\{x_n\}}$
with 
$i < j$
.
The above tier structures make hierarchies for the complexes with respect to the sizes of 
$(x_n \vee 1)^y$
and 
$\lambda_y(x_n)$
, respectively, along a sequence 
$\{x_n\}$
. Specifically, a D-type tier structure provides an ordering for the order of magnitudes of deterministic mass action rate functions of different reactions: the deterministic mass action rate of the reaction 
$y_1\to y_1^{\prime}$
is of a higher order of magnitude with respect to the rate of 
$y_2\to y_2^{\prime}$
, along a sequence of points 
$\{x_n\}$
, if 
$y_1 \succ_{\rm D} y_2$
. Similarly, an S-type tier structure defines the order of magnitudes of the stochastic mass action rate functions of reactions with different source complexes. The reason why in Definition 5.1 the quantity 
$(x_n\vee 1)^y$
appears, rather than the mass action rate 
$x_n^y$
, is technical: this is due to the necessity of considering sequences of points 
$\{x_n\}$
with possibly null components.
Definition 5.3. Let X be the Markov process associated with a mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Let 
$\{x_n\}$
be a sequence in 
${\mathbb Z}^d_{\ge 0}$
, with n varying from 1 to 
$\infty$
. The sequence 
$\{x_n\}$
is called a proper tier sequence of X if the following holds:
(i)
$x_n \in {\mathbb{S}_{x_0}}$
for each
$n\ge 1$
;(ii)
$\lim_{n\rightarrow\infty}x_{n,i} \in [0,\infty]$
for each i, and there exists an i for which
$\lim_{n\rightarrow \infty}x_{n,i} = \infty$
;(iii)
${\mathcal C}$
has both a D-type partition and an S-type partition along
$\{x_n\}$
.
5.2. Known results
The next lemma points out that any unbounded sequence of points contains a proper tier sequence as a subsequence.
Lemma 5.1. Let X be the Markov process associated with a mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Let 
$\{x_n\}_{n=1}^\infty \subset {\mathbb{S}_{x_0}}$
be a sequence such that 
$\limsup_{n\rightarrow \infty}|x_n|$
= 
$\infty$
. Then there exists a subsequence of 
$\{x_n\}$
which is a proper tier sequence of X.
The proof of Lemma 5.1 is similar to that of Lemma 4.2 in [Reference Anderson3]. In particular, the relevant tiers can each be constructed sequentially by repeatedly taking subsequences.
The following lemma states that for each proper tier sequence, there exists at least one reaction whose source and product complex are not in the same D-type tier. This lemma was originally stated in [Reference Anderson, Cappelletti, Kim and Nguyen4] as Lemma 4.1.
Lemma 5.2. Let X be the Markov process associated with a mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. For a proper tier sequence 
$\{x_n\}$
of X, there exists a reaction 
$y \to y^{\prime} \in {\mathcal R}$
such that 
$y \in T^{{\rm D},i}_{\{x_n\}}$
and 
$y^{\prime} \in T^{{\rm D},j}_{\{x_n\}}$
for some 
$i\neq j$
.
5.3. Tiers along sequences of paths
Here, we define tier structures on a set of reactions, which will be utilized to describe the behavior of the embedded Markov chain after a fixed number of transitions. We recall here that the classical strategy of studying the behavior of the embedded chain after each step, using the function V described in (6) as a potential Lyapunov function, may not lead to a proof of positive recurrence. This is the case for Example 2, and this is the reason why we extend the method to consider mutliple transitions of the embedded chain.
Definition 5.4. Let X be the Markov process associated with a mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Let k be a finite, positive integer and let 
$R=\{y_1\to y^{\prime}_1,y_2\to y^{\prime}_2,\dots,y_k \to y^{\prime}_k \}$
be a sequence of elements of 
${\mathcal R}$
of length k. Let 
$\{x_n\}$
be a sequence such that 
$\{x_n + \sum_{j=1}^{i-1}(y^{\prime}_j -y_j)\}$
is a proper tier sequence of X for each 
$i=1,2,\dots,k$
. Then,
(i)
$R \in \mathbb{D}_{\{x_n\}}$
if
\begin{align*} y_i & \in T^{{\rm D},1}_{\big\{x_n+\sum_{j=1}^{i-1}\big(y^{\prime}_j-y_j\big)\big\}} \text{ for each } i=1,2,\dots,k \text{ and } \\
y^{\prime}_\ell & \notin T^{{\rm D},1}_{\big\{x_n+\sum_{j=1}^{\ell-1}(y^{\prime}_j-y_j)\big\}} \text{ for some } \ell \in \{1,2,\dots,k\};\end{align*}
(ii)
$R\in \mathbb{T}^{{\rm S},1}_{\{x_n\}}$
if
$y_i \in T^{{\rm S},1}_{\{x_n+\sum_{j=1}^{i-1}(y^{\prime}_j-y_j)\}}$
for each
$i=1,2,\dots,k$
.
Note that (ii) in Definition 5.4 implies that this particular sequence of reactions is one of the most likely sequences of k reactions to be observed if the process starts at state 
$x_n$
. This will be made precise in the next lemma. First, we introduce the following notation: consider a mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
and let 
$\tilde X$
be the associated embedded Markov chain, with 
$\tilde{X}_0=x$
. Denote the event that the reactions 
$y_1\to y_1^{\prime}$
, 
$y_2 \to y_2^{\prime}$
, 
$\dots$
, 
$y_k \to y_k^{\prime}$
are the first k reactions to occur in order by
\begin{equation*} x \xrightarrow{y_1\to y_1^{\prime},\dots,y_k\to y^{\prime}_k} x+\sum_{j=1}^k\big(y^{\prime}_j-y_j\big).\end{equation*}
Lemma 5.3. Let X be the Markov process associated with a weakly reversible mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Let 
$\tilde{X}$
be the embedded discrete-time Markov chain of X. Let k be a finite, positive integer and let 
$R=\{y_1\to y^{\prime}_1$
, 
$y_2\to y^{\prime}_2$
, …, 
$y_k\to y^{\prime}_k\}$
be a sequence of elements of 
${\mathcal R}$
of length k. Let 
$\{x_n\}$
be a sequence such that 
$\{x_n + \sum_{j=1}^{\ell-1} (y^{\prime}_j-y_j)\}$
is a proper tier sequence of X for each 
$\ell=1,2,\dots,k$
. Then,
(i) if
$R \in \mathbb{T}^{{\rm S},1}_{\{x_n\}}$
,
\begin{equation*}{\lim_{n\rightarrow \infty}} \mathbb{P}_{x_n}\Bigg(x_n \xrightarrow{y_1\to y_1^{\prime},\dots,y_k\to y^{\prime}_k} x_n+\sum_{j=1}^k\big(y^{\prime}_j-y_j\big) \Bigg) > 0, \quad {and}\end{equation*}
(ii) if
$R \not \in \mathbb{T}^{{\rm S},1}_{\{x_n\}}$
,
\begin{equation*}{\lim_{n\rightarrow \infty}} \mathbb{P}_{x_n}\Bigg(x_n \xrightarrow{y_1\to y_1^{\prime},\dots,y_k\to y^{\prime}_k} x_n+\sum_{j=1}^k\big(y^{\prime}_j-y_j\big)\Bigg)=0.\end{equation*}
Proof. For each 
$n\geq 1$
and 
$1\leq m\leq k$
, let 
$z_n(m)=x_n+\sum_{j=1}^{m-1}(y^{\prime}_j-y_j)$
. By (4), we have
\begin{equation}\mathbb{P}_{x_n}\Bigg (x_n \xrightarrow{y_1\to y_1^{\prime},\dots,y_k\to y^{\prime}_k} x_n+\sum_{j=1}^k\big(y^{\prime}_j-y_j\big)\Bigg )= \prod_{m=1}^{k} \frac{\lambda_{y_m\to y^{\prime}_m}(z_n(m))}{\overline{\lambda}(z_n(m))}\textbf{1}_{\{\overline{\lambda}(z_n(m))>0\}}.\end{equation}
Assume that 
$R \in \mathbb{T}^{{\rm S},1}_{\{x_n\}}$
. Then, for each 
$1\leq m \leq k$
we have 
$y_m \in T^{{\rm S},1}_{\{z_n(m)\}}$
, and, by the definition of S-type tiers,
\begin{align*}{\lim_{n\rightarrow \infty}} \dfrac{\lambda_{y_m\to y^{\prime}_m}(z_n(m))}{\overline{\lambda}(z_n(m))} = {\lim_{n\rightarrow \infty}} \dfrac{\kappa_{y_m\to y^{\prime}_m}}{{\displaystyle \sum\limits_{y\to y^{\prime}\in {\mathcal R}}}\kappa_{y\to y^{\prime}}\lambda_y(z_n(m))/\lambda_{y_m}(z_n(m))} > 0.\end{align*}
Hence, (i) holds.
Now assume that 
$R \not \in \mathbb{T}^{{\rm S},1}_{\{x_n\}}$
. Hence, by the definition of 
$\mathbb{T}^{{\rm S},1}_{\{x_n\}}$
there exists 
$1\leq \ell \le k$
with 
$y_\ell \in T^{{\rm S},j}_{\{z_n(\ell)\}}$
for 
$j > 1$
. We now consider the following single term from (7):
\begin{equation} \frac{\lambda_{y_\ell\to y^{\prime}_\ell}(z_n(\ell))}{\overline{\lambda}(z_n(\ell))} \textbf{1}_{\{\overline{\lambda}(z_n(\ell))>0\}}.\end{equation}
There are two cases. If all intensity functions are identically equal to zero at 
$z_n(\ell)$
, then the term (8) is zero. If not all of the intensity functions are zero, then, because the network is weakly reversible, there is a source complex 
$y\in T^{{\rm S},1}_{\{z_n(\ell)\}}$
, in which case
\begin{equation*}\lim_{n\to \infty} \frac{\lambda_{y_\ell\to y^{\prime}_\ell}(z_n(\ell))}{\overline{\lambda}(z_n(\ell))}\textbf{1}_{\{\overline{\lambda}(z_n(\ell))>0\}} \le \lim_{n\to \infty} \frac{\lambda_{y_\ell\to y^{\prime}_\ell}(z_n(\ell))}{\lambda_{y\to y^{\prime}}(z_n(\ell))}= 0 \end{equation*}
by the definition of S-type tiers. Hence, (ii) is shown since all the terms in (7) are less than or equal to one.
Further useful properties of tiers along sequences of paths are given in the following lemmas.
Lemma 5.4. Let X be the Markov process associated with a mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Suppose that w is a fixed vector and that both 
$\{x_n\}$
and 
$\{x_n+w\}$
are proper tier sequences for X. Suppose their D-type tier partitions are 
$\{T^{{\rm D},j}_{\{x_n\}}\}_{j=1}^{J_1}$
and 
$\{T^{{\rm D},j}_{\{x_n+w\}}\}_{j=1}^{J_2}$
, respectively. Then 
$J_1 =J_2$
and 
$T^{{\rm D},j}_{\{x_n\}} = T^{{\rm D},j}_{\{x_n+w\}}$
for each j.
Proof. To show that the D-type partitions along 
$\{x_n\}$
and 
$\{x_n+w\}$
coincide, it suffices to show that if 
$y\succ_{\rm D} y^{\prime}$
for the D-type partitions along 
$\{x_n\}$
, then 
$y\succ_{\rm D} y^{\prime}$
for the D-type partitions along 
$\{x_n+w\}$
. This holds because the limits defined in part (ii) of Definition 5.1 do not change if 
$x_n$
is replaced by 
$x_n+w$
.
Lemma 5.5. Let X be the Markov process associated with a mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Let 
$y\to y^{\prime} \in {\mathcal R}$
and assume that both 
$\{x_n\}$
and 
$\{x_n+y^{\prime}-y\}$
are proper tier sequences of X. Moreover, assume that 
$y \not \in {T^{{\rm S},\infty}_{\{x_n\}}}$
and 
$y^{\prime} \in {T^{{\rm D},1}_{\{x_n\}}}$
. Then 
$y^{\prime} \in T^{{\rm S},1}_{\{x_n+y^{\prime}-y\}}$
.
To prove Lemma 5.5 we will make use of the following result, which follows immediately from [Reference Anderson and Kim10, Corollary 7].
Lemma 5.6. Let 
$\{x_n\}$
be a proper tier sequence of a reaction network 
$({\mathcal S},{\mathcal C},{\mathcal R})$
. Suppose 
$y \notin T^{{\rm S},\infty}_{\{x_n\}}$
and 
$y \in T^{{\rm D},1}_{\{x_n\}}$
. Then 
$y \in T^{{\rm S},1}$
and 
$\lim_{n\to \infty} \lambda_y(x_n) = \infty$
.
We are now ready to prove Lemma 5.5.
Proof of Lemma 5.5: We have that 
$y \not \in {T^{{\rm S},\infty}_{\{x_n\}}}$
, which is equivalent to 
$x_n \ge y$
, which is in turn equivalent to 
$x_n + y^{\prime} - y \ge y^{\prime}$
, which implies 
$y^{\prime} \not \in T^{{\rm S},\infty}_{\{x_n+y^{\prime}-y\}}$
. Now we must just show that y ′ actually yields one of the largest intensities on the sequence 
$\{x_n+y^{\prime}-y\}$
.
Because 
$y^{\prime} \in {T^{{\rm D},1}_{\{x_n\}}}$
, Lemma 5.4 implies 
$y^{\prime} \in T^{{\rm D},1}_{\{x_n+y^{\prime}-y\}}$
. Since 
$y^{\prime} \notin T^{{\rm S},\infty}_{\{x_n+y^{\prime}-y\}}$
, Lemma 5.6 implies 
$y^{\prime}\in T^{{\rm S},1}_{\{x_n+y^{\prime}-y\}}$
, and the result is shown.
6. Tiers and positive recurrence
The goal of this section is to prove the following theorem.
Theorem 6.1. Consider a weakly reversible mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
that has a single linkage class, and let X be the associated continuous-time Markov chain, with 
$X(0)=x_0$
. Suppose
\begin{equation}T^{{\rm S},1}_{\{z_n\}} \subset T^{{\rm D},1}_{\{z_n\}}\ for\ any\ proper\ tier\ sequence\ \{z_n\}\ \text{of}\ X.\end{equation}
Then X is positive recurrent.
To prove this result, we require a number of preliminary lemmas.
Lemma 6.1. Let 
$\{x_n\}$
be a sequence of points in 
${\mathbb Z}^d_{\ge 0}$
such that 
$\lim_{n\to \infty} |x_n|=\infty$
. Let 
$I = \{ i \mid x_{n,i} \rightarrow \infty \ \textrm{as} \ n \rightarrow \infty \}$
, and assume that 
$I \ne \emptyset$
and that 
$I^{\rm c} = \{i \mid \sup_{n} x_{n,i} < \infty\}$
. Let 
$\eta \in {\mathbb R}^d$
be such that 
$x_n + \eta \in {\mathbb R}^d_{\ge 0}$
for each 
$n \ge 1$
, and let V be defined as in (6). Then there exists a positive constant C such that, for all n,
\begin{align*} V(x_n+\eta)-V(x_n) \le \ln{( (x_n\vee 1)^\eta )}+C.\end{align*}
Proof. Since the terms associated with 
$I^{\rm c}$
are uniformly bounded, there exists a 
$C_1\ge 0$
such that
\begin{align}V(x_n+\eta)-V(x_n) &\le {\displaystyle \sum\limits_{i\in I}}[ (x_{n,i}+\eta_i)(\ln{(x_{n,i}+\eta_i)}-1) - x_{n,i}(\ln{(x_{n,i})}-1)] + C_1\notag\\&= {\displaystyle \sum\limits_{i\in I}}\bigg[ x_{n,i}\ln \left(1+\frac{\eta_i}{x_{n,i}}\right)-\eta_i + \eta_i\ln{(x_{n,i}+\eta_i)}\bigg] + C_1\notag.\end{align}
Using the fact that 
$\lim_{t \to \infty} \big(1 + \tfrac{\alpha}{t}\big)^t = {\rm e}^\alpha$
for any 
$\alpha$
, we have that
\begin{equation*}x_{n,i}\ln \left(1+\frac{\eta_i}{x_{n,i}}\right)-\eta_i \rightarrow 0 \ \textrm{as} \ n \rightarrow \infty \end{equation*}
for each 
$i\in I$
. Hence, there is a 
$C_2>0$
for which
\begin{equation*}V(x_n+ \eta)-V(x_n) \le {\displaystyle \sum\limits_{i \in I}} \eta_i\ln (x_{n,i}+\eta_i) + C_2 \end{equation*}
for each n. The limit
\begin{equation*} \ln \left[ (x_{n,i}+\eta_i)^{\eta_i}\right] - \ln\big(x_{n,i}^{\eta_i}\big) \to 0 \end{equation*}
as 
$n\to \infty$
then implies the existence of a 
$C_3>0$
for which
\begin{align*}V(x_n+\eta)&-V(x_n) \le \sum_{i \in I}\ln \left( x_{n,i}^{\eta_i}\right) + C_3 = \ln \left( \prod_{i\in I} x_{n,i}^{\eta_i}\right) + C_3.\end{align*}
Finally, we may add back in the bounded terms to conclude that there is a 
$C>0$
for which
\begin{align*}V(x_n+\eta)&-V(x_n) \le \ln{\left({\displaystyle \prod\limits_{i=1}^d} (x_{n,i} \vee 1)^{\eta_i}\right)} + C = \ln{(x_n\vee 1)^\eta} +C,\end{align*}
and the proof is complete.
Lemma 6.2. Let X be the Markov process associated with a weakly reversible mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Suppose (9) holds, and let 
$\tilde{X}$
be the embedded discrete-time Markov chain of X. Let k be a finite, positive integer, and let 
$R=\{y_1\to y^{\prime}_1$
, 
$y_2\to y^{\prime}_2$
, … , 
$y_k\to y^{\prime}_k\}$
be a sequence of elements of 
${\mathcal R}$
of length k. Let 
$\{x_n\}$
be a fixed sequence such that 
$\{x_n + \sum_{j=1}^{i-1}(y^{\prime}_j-y_j)\}$
is a proper tier sequence of X for each 
$i=1,2,\dots,k$
. Then, for the function V defined as (6):
(i) there is a constant
$K>0$
such that
for all
\begin{align*} \mathbb{P}_{x_n}\Bigg (x_n \xrightarrow{y_1\to y_1^{\prime},\dots,y_k\to y^{\prime}_k} x_n+\sum_{j=1}^k\big(y^{\prime}_j-y_j\big)\Bigg )\Bigg(V\Bigg(x_n+\sum_{j=1}^k \big(y^{\prime}_j-y_j\big)\Bigg)-V(x_n)\Bigg) \le K \end{align*}
$n \ge 1$
, and(ii) if
$R \in \mathbb{T}^{{\rm S},1}_{\{x_n\}}\cap \mathbb{D}_{\{x_n\}}$
, then
\begin{multline*} {\lim_{n\rightarrow \infty}} \mathbb{P}_{x_n}\Bigg (x_n \xrightarrow{y_1\to y_1^{\prime},\dots,y_k\to y^{\prime}_k} x_n+\sum_{j=1}^k\big(y^{\prime}_j-y_j\big)\Bigg ) \times \\
\Bigg(V\Bigg(x_n+\sum_{j=1}^k \big(y^{\prime}_j-y_j\big)\Bigg)-V(x_n)\Bigg) = -\infty.\end{multline*}
Proof. For each 
$n\geq 1$
and 
$1\leq m\leq k$
, let 
$z_n(m)=x_n+\sum_{j=1}^{m-1}(y^{\prime}_j-y_j)$
. Note that each 
$\{z_n(m)\}$
is a proper tier sequence, which implies 
$z_n(m) \in {\mathbb R}^d_{\ge 0}$
for each n and m. Hence, Lemma 6.1 implies that for each proper tier sequence 
$\{x_n\}$
, there is a constant 
$C>0$
such that, for all n,
\begin{align}V\Bigg(x_n+\sum_{m=1}^k \big(y^{\prime}_j-y_j\big)\Bigg)-V(x_n) &= \sum_{m=1}^k\Big ( V(z_n(m)+y^{\prime}_m-y_m)-V(z_n(m)) \Big) \notag\\[5pt]
&\le \ln{\left ( \prod_{m=1}^k \frac{(z_n(m) \vee 1)^{y^{\prime}_m}}{(z_n(m) \vee 1)^{y_m}} \right ) }+ C.\end{align}
First, we suppose 
$y_\ell \in T^{{\rm S},\infty}_{\{z_n(\ell)\}}$
for some 
$\ell$
; then, by (7),
\begin{align*}\mathbb{P}_{x_n}\Bigg (x_n \xrightarrow{y_1\to y_1^{\prime},\dots,y_l\to y^{\prime}_k} x_n+\sum_{j=1}^k\big(y^{\prime}_j-y_j\big)\Bigg )=0,\end{align*}
and (i) is shown.
Now assume 
$y_m \not \in T^{{\rm S},\infty}_{\{z_n(m)\}}$
for all 
$1\le m\le k$
. Then, there is at least one complex 
$\hat y_m\in T^{{\rm S},1}_{\{z_n(m)\}}$
for each 
$1\le m\le k$
. By (9) and [Reference Anderson and Kim10, Lemma 6], for any 
${y}\in {\mathcal C}$
and any 
$1\leq m\leq k$
,
\begin{equation}{\lim_{n\rightarrow \infty}} \frac{(z_n(m)\vee 1)^{{y}}}{\lambda_{\hat y_m}(z_n(m))}= {\lim_{n\rightarrow \infty}} \frac{(z_n(m)\vee 1)^{{y}}}{(z_n(m)\vee 1)^{\hat y_m}}\frac{(z_n(m)\vee 1)^{\hat y_m}}{\lambda_{\hat y_m}(z_n(m))}
< \infty. \end{equation}
\begin{align}\mathbb{P}_{x_n}\Bigg (x_n& \xrightarrow{y_1\to y_1^{\prime},\dots,y_l\to y^{\prime}_k} x_n+\sum_{j=1}^k(y^{\prime}_j-y_j)\Bigg )\Bigg(V\Bigg(x_n+\sum_{j=1}^k (y^{\prime}_j-y_j)\Bigg)-V(x_n)\Bigg) \notag\\[4pt]
\le &\left(\prod_{m=1}^{k}\frac{\kappa_{y_m\to y^{\prime}_m} \lambda_{y_m}(z_n(m))}{\overline{\lambda}(z_n(m))}\right)\left ( \ln{\left ( \prod_{m=1}^k \frac{(z_n(m) \vee 1)^{y^{\prime}_m}}{(z_n(m) \vee 1)^{y_m}} \right ) }+ C \right )\notag\\[4pt]
= &\left(\prod_{m=1}^k\frac{\kappa_{y_m\to y^{\prime}_m} \lambda_{y_m}(z_n(m))}{(z_n(m)\vee 1)^{y_m}}\right) \psi_n\left(\ln{\left(\frac{1}{\psi_n}\right)}+\ln{(\phi_n)}+C\right),\end{align}
where
\begin{equation*} \psi_n = \prod_{m=1}^k \frac{(z_n(m) \vee 1)^{y_m}}{\overline{\lambda}(z_n(m))} \quad \text{and} \quad \phi_n =\prod_{m=1}^k \frac{(z_n(m) \vee 1)^{y^{\prime}_m}}{\overline{\lambda}(z_n(m))},\end{equation*}
and we recall the definition of 
$\overline \lambda$
from (5). Since the network is weakly reversible, any complex is a source complex, which in turn implies that 
$\lambda_{y}(x)\leq \overline{\lambda}(x)$
for any 
$y\in{\mathcal C}$
and 
$x\in{\mathbb Z}^d$
. By combining this with (11), it follows that there is a constant 
$C^{\prime\prime}>0$
such that
\begin{align} 0<\psi_n \le C^{\prime\prime} \quad \text{and} \quad 0<\phi_n \le C^{\prime\prime} \qquad \text{for all} \ n. \end{align}
Hence, (i) is proved by combining (12), (13), the fact that 
$\psi\ln{(1/\psi)}$
is bounded from above if 
$\psi$
is bounded from above, and finally the fact that 
$\lambda_y(x)\le x^y \le (x \vee 1)^y$
for any 
$y\in{\mathcal C}$
and 
$x\in{\mathbb Z}^d$
.
In order to prove (ii), suppose 
$R \in \mathbb{T}^{{\rm S},1}_{\{x_n\}}\cap \mathbb{D}_{\{x_n\}}$
. Since 
$R \in \mathbb{T}^{{\rm S},1}_{\{x_n\}}$
, Lemma 5.3 implies that
\begin{equation*} {\lim_{n\rightarrow \infty}} \mathbb{P}_{x_n}\Bigg (x_n \xrightarrow{y_1\to y_1^{\prime},\dots,y_k\to y^{\prime}_k} x_n+\sum_{j=1}^k(y^{\prime}_j-y_j)\Bigg )>0. \end{equation*}
Hence, (ii) is shown if we prove that the right-hand side of (10) goes to 
$-\infty$
as 
$n \to \infty$
. Because 
$R \in \mathbb{D}_{\{x_n\}}$
, which implies each 
$y_m \in T^{{\rm D},1}_{\{z_n(m)\}}$
, we may conclude that each term of the form
\begin{equation*} \frac{(z_n(m) \vee 1)^{y^{\prime}_m}}{(z_n(m) \vee 1)^{y_m}}\end{equation*}
is uniformly bounded from above in n. Moreover, and again because 
$R \in \mathbb{D}_{\{x_n\}}$
, there exists 
$1\leq \ell\leq k$
such that 
$y_\ell \in T^{{\rm D},1}_{z_n(m)}$
, but 
$y^{\prime}_\ell \not \in T^{{\rm D},1}_{\{z_n(m)\}}$
. Hence,
\begin{align*}{\lim_{n\rightarrow \infty}} \ln{\left ( \prod_{m=1}^k \frac{(z_n(m) \vee 1)^{y^{\prime}_m}}{(z_n(m) \vee 1)^{y_m}} \right ) } = -\infty,\end{align*}
and (ii) is shown.
Lemma 6.3. Let 
$k \in {\mathbb Z}_{>0}$
be fixed and let X be the Markov process associated with a weakly reversible mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
, with 
$X(0)=x_0$
. Suppose (9) holds and that for any proper tier sequence 
$\{x_n\}$
of X there is a sequence of k reactions 
$R=\{y_1\to y^{\prime}_1,\dots,y_k \to y^{\prime}_k\}$
such that, by potentially considering a subsequence of 
$\{x_n\}$
, the following holds:
(i)
$\big\{x_n+\sum_{j=1}^{m-1}(y^{\prime}_j-y_j)\big\}$
is a proper tier sequence for each
$m=1,2,\dots,k$
, and(ii)
$R \in \mathbb{T}^{{\rm S},1}_{\{x_n\}} \cap \mathbb{D}_{\{x_n\}}$
.
Then X is positive recurrent.
Proof. Let 
$\tilde X_n$
, 
$n\ge 0$
, be the embedded chain for X. For 
$n \ge 0$
, define 
$Z^k_n = \tilde X_{nk}$
. Note that positive recurrence of 
$Z^k_n$
implies positive recurrence of 
$\tilde X_n$
, which in turn implies positive recurrence of X via Theorem 3.1.
Suppose, in order to find a contradiction, that 
$Z^k_n$
, 
$n \ge 0$
, is not positive recurrent for some choice of 
$Z^k_0 = \tilde X_0 = x_0$
. Then, by the usual Foster–Lyapunov criteria [Reference Meyn and Tweedie27], there exists a sequence of points 
$x_n\in{\mathbb{S}_{x_0}}$
with 
$\lim_{n\to\infty}|x_n| = \infty$
and for which
\begin{equation} \mathbb{E}_{x_n}\big[V\big(Z^k_{1}\big)\big] - V(x_n) > -1,\end{equation}
where V is defined in (6). By Lemma 5.1, we may choose a subsequence of 
$\{x_n\}$
that is a proper tier sequence of X. By assumption, there is then a further subsequence, which we will also denote by 
$\{x_n\}$
, and an ordered sequence of reactions 
$\{y_1 \to y_1^{\prime}, \dots, y_k \to y_k^{\prime}\}$
for which both (i) and (ii) in the statement of the lemma hold. Denote by 
${\mathcal R}_k$
all possible sequences of reactions of length k chosen from 
${\mathcal R}$
. We will denote elements of 
${\mathcal R}_k$
by 
$\bar{y}_1\to\bar{y}^{\prime}_1, \dots, \bar{y}_k\to\bar{y}^{\prime}_k$
. Perhaps after consideration of another subsequence, by applying Lemma 5.1 we may further assume that 
$\big\{x_n+\sum_{j=1}^{m-1}(\bar y^{\prime}_j-\bar y_j)\big\}$
is a proper tier sequence for each 
$m=1,2,\dots,k$
and for each element of 
${\mathcal R}_k$
, which are finitely many.
We may now observe that
\begin{multline*}\mathbb{E}_{x_n}\big[V\big(Z^k_{1}\big)\big] - V(x_n)= \sum_{{\mathcal R}_k} \mathbb{P}_{x_n}\Bigg(x_n \xrightarrow{\bar{y}_1\to\bar{y}^{\prime}_1, \dots, \bar{y}_k\to\bar{y}^{\prime}_k} x_n+\sum_{j=1}^k(\bar{y}^{\prime}_j-\bar{y}_j) \Bigg) \times \\
\Bigg( V\bigg(x_n+\sum_{j=1}^k(\bar{y}^{\prime}_j-\bar{y}_j)\bigg)-V(x_n) \Bigg),\end{multline*}
where 
$\mathbb{P}_{x}$
is a probability measure of 
$\tilde X$
with an initial point x. Note that by Lemma 6.2 there exists a constant 
$\overline K$
such that
\begin{align*}&\sup_n \mathbb{P}_{x_n}\Bigg(x_n \xrightarrow{\bar{y}_1\to\bar{y}^{\prime}_1, \dots, \bar{y}_k\to\bar{y}^{\prime}_k} x_n+\!\!\sum_{j=1}^k(\bar{y}^{\prime}_j-\bar{y}_j) \Bigg)\Bigg(V\Bigg(x_n+\!\!\sum_{j=1}^k(\bar{y}^{\prime}_j-\bar{y}_j)\Bigg)-V(x_n)\Bigg) \le \overline K \end{align*}
for each element of 
${\mathcal R}_k$
. Note that the particular K from the conclusion of Lemma 6.2 depends upon the particular element of 
${\mathcal R}_k$
under consideration. Therefore, 
$\overline K$
is the maximum of all such K taken over these elements, which are finitely many. Next, note that for the particular sequence 
$y_1\to y_1^{\prime}, \dots, y_k \to y_k^{\prime}$
satisfying (i) and (ii), Lemma 6.2 implies that
\begin{equation*}\lim_{n\to \infty} \mathbb{P}_{x_n}\Bigg(x_n \xrightarrow{{y}_1\to{y}^{\prime}_1, \dots, {y}_k\to{y}^{\prime}_k} x_n+\sum_{j=1}^k({y}^{\prime}_j-{y}_j) \Bigg)\Bigg(V\Bigg(x_n+\sum_{j=1}^k({y}^{\prime}_j-{y}_j)\Bigg)-V(x_n)\Bigg) =-\infty.\end{equation*}
We immediately see that the above three equations contradict (14), and so the result is shown.
The next lemma states that assumption (9) actually ensures that conditions (i) and (ii) in Lemma 6.3 hold.
Lemma 6.4. Consider a weakly reversible mass action system 
$({\mathcal S},{\mathcal C},{\mathcal R}, \kappa)$
that has a single linkage class, and let X be the associated continuous-time Markov chain. Let 
$X(0)=x_0$
and assume that 
$x_0$
is not an absorbing state. Let 
$r=|{\mathcal R}|$
, and suppose (9) holds. Let 
$\{x_n\}$
be a proper tier sequence of X. Then, by potentially considering a subsequence of 
$\{x_n\}$
, there exists a sequence of reactions 
$R=\{y_1\to y^{\prime}_1,y_2\to y^{\prime}_2,\dots,y_r \to y^{\prime}_r \}$
such that
(i)
$\big\{x_n+\sum_{j=1}^{m-1}(y^{\prime}_j-y_j)\big\}$
is a proper tier sequence for each
$m=1,2,\dots,r$
, and(ii)
$R \in \mathbb{T}^{{\rm S},1}_{\{x_n\}} \cap \mathbb{D}_{\{x_n\}}$
.
Proof. Since 
$x_0$
is not absorbing and the network is weakly reversible, it follows from Lemma 3.1 that for any 
$x\in {\mathbb{S}_{x_0}}$
there exists a reaction 
$y\to y^{\prime}$
with 
$\lambda_{y\to y^{\prime}}(x)>0$
, hence 
$T^{{\rm S},1}_{\{x_n\}}\neq \emptyset$
. Let 
$y^\star\in T^{{\rm S},1}_{\{x_n\}}$
, which by assumption is contained in 
${T^{{\rm D},1}_{\{x_n\}}}$
. Moreover, by Lemma 5.2 there exists 
$y^{\star\star}\in {\mathcal C}\setminus {T^{{\rm D},1}_{\{x_n\}}}$
. Since the network 
$({\mathcal S}, {\mathcal C}, {\mathcal R})$
is weakly reversible and 
$y^\star$
and 
$y^{\star\star}$
are in the same linkage class, there exists a sequence of reactions 
$\{y_1\to y_1^{\prime},\dots,y_H\to y^{\prime}_H\}$
with 
$y_1=y^\star\in T^{{\rm S},1}_{\{x_n\}}$
, 
$y^{\prime}_j=y_{j+1}$
for all 
$1\leq j<H$
, and 
$y^{\prime}_H=y^{\star\star}\notin{T^{{\rm D},1}_{\{x_n\}}}$
. It follows from Lemma 5.1 that by potentially considering a subsequence of 
$\{x_n\}$
, we may assume that 
$\big\{x_n+\sum_{j=1}^{m}(y^{\prime}_j-y_j)\big\}$
is a proper tier sequence of X for each 
$1\leq m\leq H$
. It follows from Lemma 5.4 that 
$y^{\prime}_H\notin T^{{\rm D},1}_{\{x_n+\sum_{j=1}^{H-1} (y^{\prime}_j-y_j)\}}$
. Let h be the smallest index such that 
$y^{\prime}_h\notin T^{{\rm D},1}_{\{x_n+\sum_{j=1}^{h-1} (y^{\prime}_j-y_j)\}}$
. Note that 
$1\leq h\leq H\leq r$
. Define
\begin{equation*}\tilde R=\{y_1\to y_1^{\prime},\dots,y_h\to y_h^{\prime}\}.\end{equation*}
By Lemma 5.5, for each 
$2\leq m\leq h$
,
\begin{equation*}y_m=y^{\prime}_{m-1}\in T^{{\rm S},1}_{\{x_n+\sum_{j=1}^{m-1} (y^{\prime}_j-y_j)\}},\end{equation*}
which implies that 
$\tilde R\in\mathbb{T}^{{\rm S},1}_{\{x_n\}}$
. Moreover, 
$\tilde R\in\mathbb{D}_{\{x_n\}}$
by (9) and by the definition of h.
If 
$h=r$
, then the the proof is complete by choosing 
$R=\tilde R$
.
If 
$h<r$
, then we may recursively add reactions to 
$\tilde R$
as follows. Since there are finitely many reactions and there is no absorbing state in 
${\mathbb{S}_{x_0}}$
, by potentially considering a subsequence of 
$\{x_n\}$
we may assume that there exists a reaction 
$\hat y_{h+1}\to \hat y^{\prime}_{h+1}$
such that
\begin{equation*}\lambda_{\hat y_{h+1}\to \hat y^{\prime}_{h+1}}\Bigg(x_n+\sum_{j=1}^{h} (y^{\prime}_j-y_j)\Bigg)=\max_{y\to y^{\prime}\in{\mathcal R}}\lambda_{y\to y^{\prime}}\Bigg(x_n+\sum_{j=1}^{h} (y^{\prime}_j-y_j)\Bigg)>0.\end{equation*}
It follows that 
$\hat y_{h+1}\in T^{{\rm S},1}_{\{x_n+\sum_{j=1}^{h} (y^{\prime}_j-y_j)\}}$
, and hence, by (9),
\begin{equation*}\{y_1\to y_1^{\prime},\dots,y_h\to y_h^{\prime}, \hat y_{h+1}\to \hat y^{\prime}_{h+1}\}\in \mathbb{T}^{{\rm S},1}_{\{x_n\}} \cap \mathbb{D}_{\{x_n\}}.\end{equation*}
Moreover, by Lemma 5.1, there exists a subsequence of 
$\{x_n\}$
such that
\begin{equation*}\Bigg\{x_n+\sum_{j=1}^{h} (y^{\prime}_j-y_j)+ \hat y^{\prime}_{h+1}-\hat y_{h+1}\Bigg\}\end{equation*}
is a proper tier sequence. If 
$h+1=r$
, the proof is complete. Otherwise, the argument can be iterated with the new set of reactions, until the cardinality of the sequence of interest is r. This concludes the proof.
We are now ready to prove Theorem 6.1.
Proof of Theorem 6.1. If 
$x_0$
is an absorbing state, then X is positive recurrent and the proof is complete. If 
$x_0$
is not absorbing, then Lemma 6.4 can be applied. Therefore, the proof is concluded by Lemma 6.3, by choosing 
$k=r$
.
6.1. Structural network conditions
In this section, we prove Theorem 4.1. We do so by showing that under the structural assumptions in the statement of Theorem 4.1, Theorem 6.1 can be applied.
Proof of Theorem 4.1. If 
$x_0$
is an absorbing state, then X is positive recurrent and the proof is concluded. Now assume that 
$x_0$
is not absorbing. In order to apply Theorem 6.1 and conclude that X is positive recurrent, we only need to show that under the assumptions of Theorem 4.1, (9) holds. By Lemma 3.1, for any 
$x\in {\mathbb{S}_{x_0}}$
there exists a reaction 
$y\to y^{\prime}$
with 
$\lambda_{y\to y^{\prime}}(x)>0$
, hence 
$T^{{\rm S},1}_{\{x_n\}}\neq \emptyset$
for any proper tier sequence 
$\{x_n\}$
. Fix any 
$y^\star\in {T^{{\rm S},1}_{\{x_n\}}}$
, 
$y^{\star\star}\in{T^{{\rm D},1}_{\{x_n\}}}$
. We will conclude the proof by showing that necessarily
\begin{equation} y^\star \sim_{\rm D} y^{\star\star}, \end{equation}
which in turn implies 
$y^\star\in {T^{{\rm D},1}_{\{x_n\}}}$
and hence (9). We consider two cases. First, assume that 
$y^{\star\star}\notin {T^{{\rm S},\infty}_{\{x_n\}}}$
. We have
\begin{equation} \frac{(x_n\vee 1)^{y^{\star\star}}}{(x_n \vee 1)^{y^\star}} = \frac{(x_n\vee 1)^{y^{\star\star}}}{\lambda_{y^{\star\star}}(x_n)}\cdot \frac{\lambda_{y^{\star\star}}(x_n)}{\lambda_{y^\star}(x_n)}\cdot \frac{\lambda_{y^\star}(x_n)}{(x_n\vee 1)^{y^\star}}.\end{equation}
By [Reference Anderson and Kim10, Corollary 7], we have 
$y^{\star\star}\in {T^{{\rm S},1}_{\{x_n\}}}$
. Hence, the second term on the right-hand side of (16) tends to a positive constant as 
$n\to\infty$
, and by [Reference Anderson and Kim10, Lemma 6] the first and the third terms on the right-hand side of (16) tend to 1 as 
$n\to\infty$
. Hence, the quantity (16) tends to a positive constant as 
$n\to\infty$
, which implies (15).
Now assume that 
$y^{\star\star}\in {T^{{\rm S},\infty}_{\{x_n\}}}$
, meaning that
\begin{equation} \lambda_{y^{\star\star}}(x_n)=0\qquad\text{for all }n\geq1.\end{equation}
Since 
$y^{\star\star}\in {T^{{\rm D},1}_{\{x_n\}}}$
, by [Reference Anderson and Kim10, Lemma 5]
\begin{equation} \lim_{n\to\infty} (x_n\vee 1)^{y^{\star\star}}=\infty.\end{equation}
Since by assumption the network is binary, the only possibility for both (17) and (18) to hold is that there exist 
$1\leq u,v\leq d$
with 
$u\neq v$
such that: (a) 
$y^{\star\star}_u=y^{\star\star}_v=1$
and 
$y^{\star\star}_i=0$
for all 
$1\leq i\leq d$
with 
$i\notin\{ u,v\}$
, (b) 
$x_{nu}=0$
for all 
$n\geq 1$
, and (c) 
$\lim_{n\to\infty} x_{nv}=\infty$
. By assumption, either 
$\tilde y=S_v \in {\mathcal C}$
or 
$\hat y=2S_v\in {\mathcal C}$
, or both inclusions hold. If 
$\hat y$
were in 
${\mathcal C}$
, then we would have
\begin{equation*}\lim_{n\to \infty}\frac{(x_n\vee 1)^{y^{\star\star}}}{(x_n\vee 1)^{\hat y}}=\lim_{n\to \infty}\frac{x_{nv}}{x_{nv}^2}=0,\end{equation*}
which would imply 
$\hat y\succ_{\rm D} y^{\star\star}$
. This is in contradiction with the assumption that 
$y^{\star\star}\in{T^{{\rm D},1}_{\{x_n\}}}$
, hence 
$\hat y\notin {\mathcal C}$
and 
$\tilde y\in{\mathcal C}$
. Note that in this case 
$(x_n\vee 1)^{y^{\star\star}}=x_{nv}=\lambda_{\tilde y}(x_n)$
for all 
$n\geq1$
. Hence,
\begin{equation*}\lim_{n\to\infty} \frac{(x_n\vee 1)^{y^{\star\star}}}{(x_n\vee 1)^{y^\star}}=\lim_{n\to\infty}\frac{\lambda_{\tilde y}(x_n)}{\lambda_{y^\star}(x_n)}\cdot \frac{\lambda_{y^\star}(x_n)}{(x_n\vee 1)^{y^\star}}<\infty,\end{equation*}
because 
$y^\star\in{T^{{\rm S},1}_{\{x_n\}}}$
and the second term on the right-hand side limit tends to 1 by [Reference Anderson and Kim10, Lemma 6]. Then, either 
$y^\star\succ_{\rm D} y^{\star\star}$
or (15) holds. Since 
$y^{\star\star}\in{T^{{\rm D},1}_{\{x_n\}}}$
, the former cannot hold and the proof is concluded.
Acknowledgements
We gratefully acknowledge the American Institute of Mathematics (AIM) for hosting our SQuaRE entitled ‘Dynamical Properties of Deterministic and Stochastic Models of Reaction Networks’ where this work was initiated. The first author was supported by Army Research Office grant W911NF-18-1-0324. The second author received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme grant agreement no. 743269 (CyberGenetics project). The third author was supported by NSF grant DMS-1616233 (PI: German Enciso).
