1. Introduction
Continuous-time Markov chains (CTMCs) on a countable state space are widely used in applications, for example, in genetics [Reference Ewens20], epidemiology [Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani37], ecology [Reference Gardiner24], biochemistry and systems biology [Reference Wilkinson45], sociophysics [Reference Weidlich and Haag44], and queueing theory [Reference Gross and Harris26]. For a CTMC on a countable state space, criteria for dynamical properties (explosivity, recurrence, certain absorption, positive recurrence, etc.) are among the fundamental topics and areas of interest.
A primary source of inspiration for our work comes from stochastic reaction network (SRN) theory, where examples are abundant. In the present context, SRNs are CTMC models of (chemical) reaction networks with polynomial transition rates [Reference Anderson and Kurtz4] (see Section 4.1 for a precise definition). In particular, we are interested in one-species reaction networks, where the reactions take the form $n\text{S}\ce{->[\kappa]}m\text{S}$ for two nonnegative integers, n, m, and $\kappa>0$ , a positive reaction rate constant. Here S represents a (chemical) species common to all reactions in the network, and the reaction represents the conversion of n molecules of the species S into m molecules of the same species. Each reaction has a transition rate, a propensity to ‘fire’. The transition rate of $n\text{S}\ce{->[\kappa]}m\text{S}$ is $\eta(x)=\kappa x(x-1)\ldots(x-n+1)$ , $x\in\mathbb{N}_0$ . Whenever the reaction fires, the corresponding Markov chain on the state space $\mathbb{N}_0$ jumps from the current state x to the state $x+m-n$ , the number of S molecules after the firing of the reaction. Different reactions may contribute to the same transition in the state space.
While the chemical terminology may suggest that the usage of such models is restricted, this is by far not so. In fact, SRNs have widespread use in the sciences with species interpreted as agents, individuals, and similar entities, and reactions as interactions among these [Reference Anderson and Kurtz4]. One might emphasize the susceptible–infectious–recovered (SIR) model in epidemiology as a particular example [Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani37].
Consider the following two examples of one-species SRNs from the recent literature, consisting of seven and five reactions, respectively [Reference Anderson, Cappelletti, Koyama and Kurtz1]:
A key issue is to understand whether the graphical representation of the reaction networks determines the dynamics of the corresponding CTMCs, irrespective of their initial values. The first network is explosive (except if the initial state is 0, which forms a singleton communicating class), while the second is positive recurrent on the positive integers (again 0 forms a singleton class) [Reference Anderson, Cappelletti, Koyama and Kurtz1], which might be inferred from known birth–death process (BDP) criteria [Reference Anderson5]. However, these criteria are not computationally simple and blind to the graphical structure of the networks. A simple explanation for the drastic difference in the dynamics of these two random walks on $\mathbb{N}_0$ is desirable but remains unknown [Reference Anderson, Cappelletti, Koyama and Kurtz1].
Motivated by the above concern, we provide criteria for dynamical properties of CTMCs on $\mathbb{N}_0$ with polynomial-like transition rates and without the possibility of arbitrary large negative jumps, as in the examples above. These CTMCs are ubiquitous in applications [Reference Barbour7]. Specifically, we provide simple threshold criteria for the existence and nonexistence of moments of hitting times, positive recurrence and null recurrence, and exponential ergodicity of stationary distributions and quasi-stationary distributions (QSDs) in terms of four easily computable parameters, derived from the transition rates. Additionally, we provide necessary and sufficient conditions for explosivity, recurrence versus transience, certain absorption, and implosivity. These conditions provide simple explanations for the dynamical discrepancies between the two SRNs in (1).
Our approach is to apply the classical semimartingale approach used in Lamperti’s problem [Reference Lamperti29], as well as Lyapunov–Foster theory [Reference Champagnat and Villemonais12, Reference Menshikov and Petritis32, Reference Meyn and Tweedie33] with delicately constructed Lyapunov functions (in particular, we make use of the techniques in [Reference Menshikov and Petritis32]). The problem of finding neat and desirable necessary and sufficient conditions for dynamical properties of CTMCs has existed for a long time [Reference Meyn and Tweedie33, Reference Meyn and Tweedie34, Reference Champagnat and Villemonais12]. However, the fact that this has not been accomplished yet indicates that it might be a nontrivial task. A main contribution of this paper is to identify a large class of CTMCs (without a built-in detailed balanced structure) for which computationally simple sufficient and necessary criteria can be established for dynamical properties of interest. Our criteria save the effort of constructing Lyapunov functions and applying Lyapunov–Foster theory case by case. Also, a case-by-case approach is ignorant of the underlying graphical structure of the Markov chain.
The simple necessary and sufficient conditions for the dynamical properties are determined by calculating up to four parameters, R, $\alpha$ , $\beta$ , and $\gamma$ , that are expressed in terms of the coefficients of the first two terms of the polynomial-like transition rate functions (the specific assumptions are given in Section 2). For illustration, let $\Omega$ be the set of jump sizes, and let
be the transition rate functions, where $d_\omega$ is the degree of $\lambda_{\omega}$ and O is Landau’s symbol. Define $R=\max_{\omega\in\Omega}d_{\omega}$ and
Based on these four parameters, a full classification of the dynamical properties can be achieved (see Theorems 1, 3, 7, and 9) and is summarized in Table 1 below. The parameters $\alpha,\beta,\gamma$ depend only on the coefficients of the monomials of degree R and $R-1$ . Furthermore, the parameter $\alpha$ might be interpreted as a sum over the jump sizes, weighted by the coefficients of the monomials of degree R. Similarly, $\gamma$ might be interpreted as a sum over the jump sizes, weighted by the coefficients of the monomials of degree $R-1$ .
To see the power of our results, consider the following SRN, which is not a BDP:
where m is a positive integer, and $\kappa_i$ , $i=1,\ldots$ , is a positive rate constant. Then $R=m+1$ , $\alpha=2\kappa_5-\kappa_{4}$ ,
and
The criteria established in Section 3 (and collected in Table 1) imply that the SRN is (in the sense of the underlying irreducible CTMC on $\mathbb{N}_0$ )
-
(a) explosive almost surely (a.s.) if and only if (i) $\alpha>0$ or (ii) $\alpha=0$ , $\beta>0$ , $R>2$ , and non-explosive if either (i) or (ii) fails;
-
(b) recurrent if and only if (iii) $\alpha<0$ or (iv) $\alpha=0$ , $\beta\le0$ , and transient if and only if both (iii) and (iv) fail;
-
(c) positive recurrent if and only if (iii), (v) $\alpha=0$ , $\beta<0$ , or (vi) $\alpha=0$ , $\beta=0$ , $R>2$ holds, and null recurrent if and only if (vii) $\alpha=0$ , $\beta=0$ , $R=2$ ;
-
(d) implosive if and only if (iii) or (viii) $\alpha=0$ , $\beta\le0$ , $R>2$ holds, and non-implosive if and only if both (ii) and (viii) fail.
Here implosive means positive recurrent with uniformly bounded expected first return time. (See Subsection 3.5 for the precise definition.) The above example shows the applicability and simplicity of our results; in fact, the computations could easily be implemented in a software program that takes a reaction network as input and outputs the network dynamical properties. Furthermore, the example illustrates the richness of the dynamical properties that might reside within a single example by varying the parameters of the model. All possibilities for $\alpha, \beta$ and R are covered (the parameter $\gamma$ is irrelevant for SRNs [Reference Wiuf and Xu46]). The stability of the chain depends only on $\alpha$ (which is independent of m), unless $\alpha=0$ , in which case the sign of $\beta$ determines the stability. If so, then $\beta=\kappa_3-m\kappa_2-3\kappa_5$ , which depends on m. Thus, if $\kappa_3-\kappa_2-3\kappa_5>0$ , then if $m>1$ is chosen large enough, the stability of the chain flips. The parameter $\alpha$ plays a role similar to that of the largest Lyapunov exponent for $\alpha\not=0$ . Analogously, when $\alpha=0$ , the parameter $\beta$ determines the stochastic stability and hence plays a role similar to that of the second-largest Lyapunov exponent in the critical case.
Brief description of our approach
Although our approach is essentially based on Lyapunov–Foster-type results, the sharp criteria for diverse dynamical properties of CTMCs are established by combining a mixture of results [Reference Aspandiiarov and Iasnogorodski6, Reference Champagnat and Villemonais12, Reference Chen16, Reference Menshikov and Petritis32, Reference Meyn and Tweedie33]; in particular [Reference Menshikov and Petritis32] provides useful criteria.
The most prominent difficulties in deriving necessary and sufficient conditions for dynamical properties of general CTMCs, with multiple jump sizes, lie in the non-calculability of stationary distributions/measures, as well as the nonexistence of orthogonal polynomials [Reference Karlin and McGregor28]. This also explains why, in general, a partial result in terms of a sufficient but not a necessary condition, by construction of a Lyapunov function, is likely. Here we discover that Lyapunov–Foster theory and the semimartingale approach are indeed enough to derive necessary and sufficient conditions. To obtain conditions that are not only sufficient but also necessary, we check whether the negation of a condition is also sufficient for the reverse dynamical property. Moreover, to show null recurrence, we also rely on estimates for stationary measures in [Reference Aspandiiarov and Iasnogorodski6]. Finally, we would like to point out that some of the Lyapunov functions we use appear to be rarely used in the literature.
Comparison with results in the literature
Complete classification of dynamical properties seems quite rare in the literature. Here, we summarize relevant results together with the methods applied.
Reuter provided necessary and sufficient conditions for explosivity of CTMCs (known as Reuter’s criterion) [Reference Reuter38], but these conditions are difficult to check except in special cases, e.g. for BDPs [Reference Karlin and McGregor28] and competition processes [Reference Reuter39]. This is due to the fact that the conditions involve infinitely many algebraic equations.
Karlin and McGregor established threshold results for explosivity and recurrence, as well as certain absorption of BDPs with and without absorbing states, by means of the so-called Karlin–McGregor integral representation formula [Reference Karlin and McGregor28]. The existence of such a formula is essentially due to the tridiagonal structure of the Q-matrix. For the same reason, it is delicate to extend such an approach to generalized BDPs: pure birth processes [Reference Chen14], one-sided skip-free CTMCs [Reference Chen, Pollet, Zhang and Cairns13, Reference Chen15, Reference Chen16], and recently (higher-dimensional) quasi-birth–death processes (QBDPs) with tridiagonal block structure of the Q-matrix [Reference Fernández and de la Iglesia22]. There are no restrictions on the types of transition rates for the results in [Reference Karlin and McGregor28] to hold, except the BDP assumption. In contrast, for our results to hold, we require the polynomial transition rates to have uniformly bounded degrees, but we do not require the BDP assumption.
In the context of QSDs, there are few threshold results for certain absorption, existence and uniqueness, and quasi-ergocidity of QSDs. Van Doorn [Reference Van Doorn42, Reference Van Doorn43] obtained ergodicity, existence and nonexistence, and uniqueness of QSDs for absorbed BDPs, also building on the Karlin–McGregor integral representation formula. Later, Ferrari et al. [Reference Fernández, Kesten, Martinez and Picco23] generalized the results in [Reference Van Doorn43]. They derived a necessary and sufficient condition for the existence of a QSD on the positive integers for which zero is an absorbing state, using the so-called renewal dynamical approach, assuming the CTMC is non-explosive, and that the absorption time is finite and unbounded with probability one. Then the existence of a QSD is equivalent to finiteness of the exponential moment of the absorption time, for one initial transient state (and hence all of them). But such a moment condition is again not straightforward to verify, pending the assumptions.
To sum up, general checkable threshold criteria for dynamical properties of CTMCs (absorbed and non-absorbed), other than generalized BDPs, are few. We identify a class of CTMCs with polynomial-like transition rates and without arbitrary large backward jumps for which simple, checkable criteria for absorbed and non-absorbed CTMCs are found, based on the coefficients of the polynomials. The price for this is to impose some further mild regularity conditions, in addition to the two requirements mentioned above.
Impact of our work and further extensions
The following are some possible extensions of our work from the theoretical perspective:
-
The sufficient condition for the existence of ergodic stationary distributions and QSDs allows us to further investigate the tail asymptotics of these distributions [Reference Xu, Hansen and Wiuf47] and the computation of these distributions (in a forthcoming paper).
-
The novel combination of the approaches presented here can further be extended to establish criteria for the dynamics of one-dimensional CTMCs with asymptotic polynomial transition rates, and higher-dimensional CTMCs with Q-matrix of a certain block structure, in analogy with QBDPs and BDPs.
-
A deeper understanding of the threshold parameters may provide insight into the dynamics of higher-dimensional CTMCs on lattices.
From the perspective of applications, we have the following:
-
The criteria can be applied to completely classify the dynamics of one-dimensional mass-action SRNs, and in particular, we can prove the so-called positive recurrence conjecture [Reference Anderson and Kim2] for weakly reversible reaction networks in one dimension [Reference Wiuf and Xu46].
-
The criteria can be used to establish bifurcations of one-dimensional SRNs (in a forthcoming paper).
Outline
In Section 2, the notation and standing assumptions are introduced. Section 3 develops threshold criteria for dynamical properties of CTMCs. Applications to SRNs, a class of branching processes, a general bursty single-cell stochastic gene expression model, and population processes of non-BDP type are provided in Section 4. Proofs of the main results are provided in Section 5. Additional tools used in the proofs as well as proofs of some elementary propositions are in the appendix.
2. Preliminaries and assumptions
Let $\mathbb{R}$ , $\mathbb{R}_{\ge 0}$ , $\mathbb{R}_{>0}$ be the sets of real, nonnegative real, and positive real numbers, respectively. Let $\mathbb{Z}$ be the set of integers, $\mathbb{N}=\mathbb{Z}\cap \mathbb{R}_{>0}$ , and $\mathbb{N}_0=\mathbb{N}\cup\{0\}$ . For $x,y\in \mathbb{N}$ , let $x^{\underline{y}}=x(x-1)\cdots(x-y+1)$ be the descending factorial of x.
Let $(Y_t\,:\, t\ge 0)$ (or $Y_t$ for short) be a CTMC on a closed, infinite state space $\mathcal{Y}\subseteq\mathbb{N}_0$ with conservative transition rate matrix $Q=(q_{x,y})_{x,y\in\mathcal{Y}}$ ; that is, every row sums to zero. A set $A\subseteq\mathcal{Y}$ is closed if $q_{x,y}=0$ for all $x\in A$ and $y\in\mathcal{Y}\setminus A$ [Reference Norris36]. Assume the absorbing set $\partial\subsetneq\mathcal{Y}$ is finite (potentially empty) and closed. Hence, $\mathcal{Y}\setminus\partial$ is unbounded.
Let $\Omega=\{y-x\,:\, q_{x,y}>0,\ \text{for some}\ x,y\in\mathcal{Y}\}$ be the set of jump sizes. For $\omega\in\Omega$ , define the transition rate function by
Let $\Omega_{\pm}=\{\omega\in\Omega\,:\, {\text{sgn}}(\omega)=\pm1\}$ be the sets of forward and backward jump sizes, respectively. Throughout, we assume the following regularity conditions:
-
(A1) $\Omega_+\neq\varnothing$ , $\Omega_-\neq\varnothing$ .
-
(A2) $\#\Omega_-<\infty$ .
-
(A3) $\sum_{\omega\in\Omega}\lambda_{\omega}(x)|\omega|<\infty,$ for all $x\in\mathcal{Y}$ .
-
(A4) There exist $u, M\in\mathbb{N}$ such that $\lambda_{\omega}$ is a strictly positive polynomial of degree $\le M$ on the set $\mathcal{Y}\setminus\{0,\ldots,u-1\}$ , for all $\omega\in\Omega$ .
-
(A5) $\mathcal{Y}\setminus\partial$ is irreducible.
If either $\Omega_+=\varnothing$ or $\Omega_-=\varnothing$ , then $Y_t$ is a pure birth or death process (possibly with multiple jump sizes). The classification of states and the dynamics of such processes are simpler than under ( $\textbf{A1}$ ). Indeed, one can derive parallel results from the corresponding results under ( $\textbf{A1}$ ). Assumption ( $\textbf{A2}$ ) implies $Y_t$ cannot make arbitrary large negative jumps. Assumption ( $\textbf{A3}$ ) is a regularity condition that ensures that functions like x, $\log x$ , and $\log\log x$ are in the domain of the infinitesimal generator of the CTMC, in order to serve as Lyapunov functions. If $\Omega$ is finite, then ( $\textbf{A2}$ ) and ( $\textbf{A3}$ ) are automatically fulfilled. In that case, the sums above are trivially polynomials for large x.
Assumption ( $\textbf{A4}$ ) implies that the Markov chain can make all jumps in $\Omega$ with positive probability from any ‘large’ state $x\in\mathcal{Y}$ , $x\ge u$ . Moreover, ( $\textbf{A3}$ ) and ( $\textbf{A4}$ ) together imply that $\sum_{\omega\in\Omega}\lambda_{\omega}(x)$ and $\sum_{\omega\in\Omega}\lambda_{\omega}(x)\omega$ are polynomials of degree $\le M$ for $x\in\mathcal{Y}\setminus\{0,\ldots,u-1\}$ (Proposition 1). If ( $\textbf{A4}$ ) fails, simple examples show that $\sum_{\omega\in\Omega}\lambda_{\omega}(x)$ and $\sum_{\omega\in\Omega}\lambda_{\omega}(x)\omega$ may not be polynomials. That these sums are polynomials for large states is an essential property that we rely on in proofs.
Assumption ( $\textbf{A4}$ ) is common in applications, especially in the context of chemical reaction networks and population processes [Reference Anderson, Kurtz, Koeppl and Setti3, Reference Ethier and Kurtz19]. Assumption ( $\textbf{A5}$ ) is also standard and is generally satisfied in applications [Reference Champagnat and Villemonais12, Reference Menshikov and Petritis32], potentially by restricting the state space. One can show that ( $\textbf{A4}$ ) and ( $\textbf{A5}$ ) together imply that $\mathcal{Y}\setminus\partial$ is infinite; thus the assumptions are not compatible with a finite state space.
With the above assumptions, the following three parameters are well-defined and finite:
If $R=0$ , then trivially $\gamma=0$ . In particular, if $\Omega$ is finite, the following additional parameter is also well-defined and finite:
with $\beta<\gamma$ . The parameter $\alpha$ encodes the sign of the average jump size of the chain. It is straightforward to verify that (3) is a consequence of the above parameter definitions, owing to the asymptotic expansions (2) of the transition rate functions.
Example 1. Recall the example (4) in the introduction,
Then $\Omega=\{1,2,-1,m,-m\}$ , $m\ge 1$ , and
Hence, $R=m+1$ , and
by (3). As mentioned in the introduction, the sign of $\alpha$ determines the stochastic stability of the CTMC [Reference Meyn and Tweedie33]. Hence, $\alpha$ plays a similar role as the largest Lyapunov exponent. Analogously, when $\alpha=0$ , the parameter $\beta$ determines the stochastic stability and hence plays a role similar to that of the second-largest Lyapunov exponent in the critical case.
3. Criteria for dynamical properties
In this section, we provide threshold criteria for various dynamical properties in terms of $R, \alpha, \beta, \gamma$ . Proofs are relegated to Section 5. For ease of comparison, we collect all parameter conditions used in the main theorems below. These are listed in the order in which they appear in the main theorems; see Table 2. Figure 1 shows implications among the 21 conditions.
3.1. Explosivity and non-explosivity
The sequence $J=(J_n)_{n\in\mathbb{N}_0}$ of jump times of a CTMC $Y_t$ are defined by $J_0=0$ and $J_n=\inf\{t\ge J_{n-1}\,:\, Y_t\neq Y_{J_{n-1}}\}$ , $n\ge1$ , where $\inf\varnothing=\infty$ by convention. The lifetime is denoted by $\zeta=\sup_{n}J_n$ . The process $Y_t$ is said to explode (with positive probability) at $y\in\mathcal{Y}$ if $\mathbb{P}_y(\{\zeta<\infty\})>0$ . In particular, $Y_t$ explodes a.s. at $y\in\mathcal{Y}$ if $\mathbb{P}_y(\{\zeta<\infty\})=1$ , and does not explode at $y\in\mathcal{Y}$ if $\mathbb{P}_y(\{\zeta<\infty\})=0$ [Reference Meyn and Tweedie33]. Hence, $Y_t$ does not explode if $Y_0\in\partial$ (since $\partial$ is closed and finite), and $\mathbb{E}_y(\zeta)<\infty$ implies that $Y_t$ explodes at y a.s. Recall that non-explosivity and explosivity are class properties. They hold for either all or no states in $\mathcal{Y}\setminus\partial$ . Hence, we simply say $Y_t$ is explosive (respectively, explosive a.s.) if it explodes with positive probability (respectively, explodes a.s.) at some state in $\mathcal{Y}\setminus\partial$ , and $Y_t$ is non-explosive if it does not explode at some state in $\mathcal{Y}\setminus\partial$ .
We present necessary and sufficient conditions for explosivity and non-explosivity.
Theorem 1. Assume $(\textbf{A1})$ – $(\textbf{A5})$ , and that $\Omega$ is finite. Then $Y_t$ is explosive with positive probability if and only if either (C1) or (C2) holds. Moreover, $Y_t$ is explosive a.s. whenever it is explosive, provided $\partial=\varnothing$ .
Theorem 2. Assume $(\textbf{A1})$ – $(\textbf{A5})$ and that $\Omega$ is infinite. Then $Y_t$ is explosive if (C1) holds, and it is non-explosive if one of the three conditions (C3), (C4), (C5) holds. Moreover, $Y_t$ is explosive a.s. whenever it is explosive, provided $\partial=\varnothing$ .
Explosion might occur with probability less than one for CTMCs with non-polynomial transition rates and $\partial=\varnothing$ [Reference Menshikov and Petritis32]. Reuter’s criterion and generalizations of it provide necessary and sufficient conditions for explosivity (with positive probability) for general CTMCs in terms of convergence or divergence of a series [Reference Chen, Pollet, Zhang and Cairns13, Reference Li and Li30, Reference Reuter38]. However, these conditions are not easy to check. In comparison, for CTMCs with polynomial transition rates, Theorem 1 provides an explicit and checkable necessary and sufficient condition.
3.2. Recurrence versus transience, and certain absorption
For a nonempty subset $A\subseteq\mathcal{Y}$ , let $\tau_A=\inf\{t\ge0\,:\, Y_t\in A\}$ be the hitting time of A, with the convention that $\inf\varnothing=\infty$ . If $Y_0\in A$ , then $\tau_A=0$ . Let $\tau^+_A=\inf\{t\ge J_1\,:\, Y_t\in A\}$ be the first return time to A. Obviously, $\tau_A=\tau^+_A$ if and only if $Y_0\notin A$ . The process $Y_t$ has certain absorption if the hitting time of $\partial$ is finite a.s. for all $Y_0\in\mathcal{Y}$ .
Theorem 3. Assume $(\textbf{A1})$ – $(\textbf{A5})$ and that $\Omega$ is finite.
(i) Assume $\partial=\varnothing$ . Then $Y_t$ is recurrent if either (C3) or (C6) holds, while it is transient if neither of them holds.
(ii) Assume $\partial\neq\varnothing$ . Then $Y_t$ has certain absorption if and only if either (C3) or (C6) holds.
The results show that CTMCs with polynomial transition rates cannot have an infinite series of critical transitions from recurrence to transience, for varying parameter values. This stands in contrast to the case of CTMCs with non-polynomial transition rates, as discovered in [Reference Menshikov and Petritis32]. One might hope that this phenomenon carries over to CTMCs with polynomial transition rates in dimensions higher than one.
Theorem 4. Assume $(\textbf{A1})$ – $(\textbf{A5})$ and that $\Omega$ is infinite.
(i) Assume $\partial=\varnothing$ . Then $Y_t$ is recurrent if (C3) holds, and it is transient if (C7) holds.
(ii) Assume $\partial\neq\varnothing$ . Then $Y_t$ has certain absorption if (C3) holds, while it does not have certain absorption if (C7) holds.
3.3. Moments of hitting times
Below we present threshold results on the existence of moments of hitting times for recurrent states only, as transient states have infinite return time. Therefore, in light of Theorem 3, we investigate the existence and nonexistence of moments of hitting times only for $\alpha<0$ and for $\alpha=0,\ \beta\le0$ . Moreover, limited by the tools we apply, we do not discuss existence and nonexistence of moments of absorption times for $\partial\neq\varnothing$ . Hence, we assume $Y_t$ is irreducible on $\mathcal{Y}$ (equivalently, $\partial=\varnothing$ ) and provide existence and nonexistence of moments of hitting times for states in $\mathcal{Y}$ .
Theorem 5. Assume $(\textbf{A1})$ – $(\textbf{A5})$ , $\partial=\varnothing$ , and that $\Omega$ is finite. Then the following hold:
(i) There exists a finite nonempty subset $B\subseteq\mathcal{Y}$ such that
for $\delta>0$ , provided one of the conditions (C3), (C9), (C10) holds; for $0<\delta<1/2$ , provided (C13) holds; for $0<\delta<\tfrac\beta{\beta-\gamma}$ , provided (C14) holds; and for $0<\delta<1$ , provided (C16) holds. In particular, $\mathbb{E}_x(\tau_B)<+\infty$ , provided one of the conditions (C3), (C9), (C10), (C11) holds.
(ii) There exists a finite nonempty subset $B\subseteq\mathcal{Y}$ such that
for $\delta>1$ , provided (C13) holds; for $\delta>\tfrac\beta{\beta-\gamma}$ , provided (C15) holds; and for $\delta>1$ , provided (C16) holds. In particular, $\mathbb{E}_x(\tau_B)=+\infty$ provided (C12) holds.
Theorem 6. Assume $(\textbf{A1})$ – $(\textbf{A5})$ , $\partial=\varnothing$ , and that $\Omega$ is infinite. If (C3) is fulfilled and $0<\delta\le1$ , then (5) holds.
3.4. Positive recurrence and null recurrence
We provide sharp criteria for positive and null recurrence, as well as exponential ergodicity of stationary distributions and QSDs.
If $\partial=\varnothing$ , then $\tau_{\partial}=\infty$ a.s., and the conditional process $(Y_t\,:\, \tau_{\partial}>t)$ reduces to $(Y_t\,:\, t\ge 0)$ . If $\partial\neq\varnothing$ , and $\tau_{\partial}<\infty$ a.s. (that is, $Y_t$ has certain absorption), then the process conditioned to never be absorbed, $(Y_t\,:\, \tau_{\partial}>t)$ , is referred to as the Q-process [Reference Champagnat and Villemonais11, Reference Collet, Martínez and San Martín18].
The process $(Y_t\,:\, \tau_{\partial}>t)$ on $\mathcal{Y}\setminus\partial$ is said to be exponentially ergodic if there exist a probability measure $\mu_*$ and $0<\delta<1$ such that for all probability measures $\mu$ on $\mathcal{Y}\setminus\partial$ , there exists a constant $C_{\mu}>0$ such that
(see [Reference Gibbs and Su25]). The measure $\mu^*$ is also said to be exponentially ergodic. In particular, if $C_{\mu}$ can be chosen independently of $\mu$ , then $(Y_t\,:\, \tau_{\partial}>t)$ and $\mu^*$ is said to be uniformly exponentially ergodic. Moreover, if $\partial=\varnothing$ , then $\mu_*$ is the unique ergodic stationary distribution; if $\partial\neq\varnothing$ , then $\mu_*$ is a quasi-limiting distribution (QLD) [Reference Collet, Martínez and San Martín18].
If $\partial\neq\varnothing$ , a probability measure $\nu$ on $\mathcal{Y}\setminus\partial$ is a QSD for $Y_t$ if for all $t\ge0$ and all sets $B\subseteq \mathcal{Y}\setminus\partial$ ,
Any QLD is a QSD [Reference Collet, Martínez and San Martín18]. The existence of a QSD implies certain absorption, and exponential ergodicity of the Q-process implies existence of a unique QSD [Reference Collet, Martínez and San Martín18]. A probability measure $\nu$ on $\mathcal{Y}\setminus\partial$ is a quasi-ergodic distribution if, for any $x\in\mathcal{Y}\setminus\partial$ and any bounded function f on $\mathcal{Y}\setminus\partial$ [Reference Breyer and Roberts10, Reference He, Zhang and Zhu27], the following limit holds:
A quasi-ergodic distribution is in general different from a QSD [Reference He, Zhang and Zhu27].
Theorem 7. Assume $(\textbf{A1})$ – $(\textbf{A5})$ and that $\Omega$ is finite.
(i) Assume $\partial=\varnothing$ and that $Y_t$ is recurrent. Then $Y_t$ is positive recurrent and there exists a unique stationary distribution $\pi$ on $\mathcal{Y}$ , if and only if one of the conditions (C3), (C9), (C10), (C11) holds, while $Y_t$ is null recurrent if and only if none of the conditions (C3), (C9), (C10), (C11) hold. Moreover, $Y_t$ is exponentially ergodic if either (C17) or (C18) holds.
(ii) Assume $\partial\neq\varnothing$ and that $Y_t$ has certain absorption. Then there exist no QSDs if none of the conditions (C3), (C9), (C10), (C11) hold. In contrast, there exists a unique uniformly exponentially ergodic QLD, supported on $\mathcal{Y}\setminus\partial$ , if either (C18) or (C19) holds. Morevover, it is also a unique quasi-ergodic distribution and the unique stationary distribution of the Q-process.
Theorem 8. Assume $(\textbf{A1})$ – $(\textbf{A5})$ and that $\Omega$ is infinite.
-
(i) Assume $\partial=\varnothing$ . Then $Y_t$ is positive recurrent and there exists a unique stationary distribution $\pi$ on $\mathcal{Y}$ if (C3) holds. Moreover, $\pi$ is exponentially ergodic if (C17) holds.
-
(ii) Assume $\partial\neq\varnothing$ . Then there exist no QSDs if (C7) holds, while there exists a unique uniformly exponentially ergodic QLD supported on $\mathcal{Y}\setminus\partial$ if (C19) holds.
We provide some perspectives:
-
The convergence (or ergodicity) in Theorem 8(ii) is uniform with respect to the initial distribution, while in contrast, the convergence in Theorem 8(i) is not uniform. Indeed, for the subcritical linear BDP, the stationary distribution is exponentially ergodic but not uniformly so [Reference Anderson5].
-
Indeed, one can obtain uniform exponential ergodicity in Theorem 7(i) with (C18) or (C19) by choosing a non-reachable absorption set (potentially empty), hence imposing that the time to extinction is infinite. In this case, the QSD is in fact a stationary distribution [Reference Champagnat and Villemonais12].
-
The subtle difference between the conditions for positive recurrence and for exponential ergodicity of QSDs lies in the fact that we have no a priori estimate of the decay parameter
\begin{equation*}\psi_0=\inf\bigl\{\psi>0\,:\, \liminf_{t\to\infty}e^{\psi t}\mathbb{P}_x(X_t=x)>0\bigr\} \end{equation*}(which is independent of x) [Reference Champagnat and Villemonais12]. We cannot compare $\psi_0$ with $-\alpha$ when $R>1$ , or with $-\beta$ when $R>2$ and $\alpha=0$ . Refer to the constructive proofs (using Lyapunov functions) in Appendix A for details. Hence, one may believe that the condition we provide for quasi-ergodicity generically is stronger than that for ergodicity.The only gap cases that remain for QSDs are (C11), (C20), (C21), where neither existence of a QSD nor exponential ergodicity of the Q-process are known to occur, provided one of the three conditions hold.
3.5. Implosivity
Assume $\partial=\varnothing$ . Then $Y_t$ is irreducible on $\mathcal{Y}$ . Let $B\subsetneq\mathcal{Y}$ be a nonempty proper subset. Then $Y_t$ implodes towards B [Reference Menshikov and Petritis32] if there exists $t_*>0$ such that
Implosion towards a single state $x\in\mathcal{Y}$ implies finite expected first return time to the state, and thus positive recurrence of x. Indeed,
where $\tau_x^+=\tau_{\{x\}}^+$ , $\tau_x=\tau_{\{x\}}$ , and $J_1$ has finite expectation since x is not absorbing. Hence, $Y_t$ does not implode towards any transient or null recurrent state.
The process $Y_t$ is implosive if $Y_t$ implodes towards any state of $\mathcal{Y}$ , and otherwise, $Y_t$ is non-implosive. Hence, implosivity implies positive recurrence. If $Y_t$ implodes towards a finite nonempty subset of $\mathcal{Y}$ , then $Y_t$ is implosive (see Proposition 17).
Theorem 9. Assume $(\textbf{A1})$ – $(\textbf{A5})$ , $\partial=\varnothing$ , and that $\Omega$ is finite. Then $Y_t$ is implosive, and there exists $\epsilon>0$ such that for every nonempty finite subset $B\subseteq\mathcal{Y}$ and every $x\in \mathcal{Y}\setminus B$ ,
if either (C18) or (C19) holds, while $Y_t$ is non-implosive otherwise.
Theorem 10. Assume $(\textbf{A1})$ – $(\textbf{A5})$ , $\partial=\varnothing$ , and that $\Omega$ is infinite. Then $Y_t$ is implosive if (C19) holds.
3.6. Relations between absorbed CTMCs and non-absorbed CTMCs
It is worth comparing the properties of absorbed CTMCs to those of non-absorbed CTMCs, as first discussed by Karlin and McGregor [Reference Karlin and McGregor28]. Indeed, in the case of a finite set of absorbing states $\partial\neq\varnothing$ and an irreducible state space $\mathcal{Y}\setminus \partial$ , one can add positive transition rates to the transition matrix from the states in $\partial$ to a finite subset of states in $\mathcal{Y}\setminus\partial$ , such that $\mathcal{Y}$ is irreducible for the new chain. Conversely, if $\partial=\varnothing$ , then one can prescribe a finite set $\partial^{\prime}\subseteq\mathcal{Y}$ , and delete all transitions from $\partial^{\prime}$ to $\mathcal{Y}^{\prime}=\mathcal{Y}\setminus\partial$ , so that $\mathcal{Y}^{\prime}$ is irreducible and $\partial^{\prime}$ an absorbing set for the new chain. These operations can be viewed as simple extensions to those of [Reference Karlin and McGregor28], proposed in the context of BDPs. As the dynamical properties we discuss generically are determined by transitions among states of large values, the operations provide a way to link the dynamics of an absorbed CTMC with that of a corresponding non-absorbed CTMC, and vice versa.
Figure 2 shows the implications among the properties, in agreement with the parameter conditions derived in the main theorems. In Examples 2 and 3 below, counterexamples are given. It remains unknown whether exponential ergodicity of the Q-process implies implosivity, and whether ergodicity implies existence of a QSD; see Figure 2.
Example 2. (i) Consider the sublinear BDP on $\mathbb{N}_0$ with birth rates $\lambda_j=a$ and death rates $\mu_j=b$ for $j\in\mathbb{N}$ . We have $\partial=\{0\}$ , $R=0$ , and $\alpha=a-b$ . Hence, the process is non-explosive for any initial state by Theorem 1. By [Reference Van Doorn43], the process has certain absorption with decay parameter $\psi_0=\big(\sqrt{a}-\sqrt{b}\big)^2$ when $a\le b$ , and it admits a continuum family of QSDs when $\alpha<0$ . This shows that existence of a QSD does not imply exponential ergodicity of the Q-process in general.
(ii) Consider the linear BDP on $\mathbb{N}_0$ with birth rates $\lambda_j=aj$ and death rates $\mu_j=bj$ for $j\in\mathbb{N}$ . Assume $a\le b$ . We have $\partial=\{0\}$ , $R=1$ , and $\alpha=a-b$ . Hence, the process is non-explosive for any initial state by Theorem 1. By [Reference Van Doorn43], the process has certain absorption with decay parameter $\psi_0=\big(\sqrt{a}-\sqrt{b}\big)^2$ when $a\le b$ , and it admits no QSDs for $\alpha=0$ (hence, $\beta=-a<0=\gamma$ ), while it admits a continuum family of QSDs for $\alpha<0$ . This shows that the process has certain absorption, but no QSDs for $\alpha=0$ , which is also justified by Theorem 7(ii). Moreover, it also shows that ergodicity of the non-absorbed process does not imply uniqueness of a QSD or exponential ergodicity of the Q-process.
(iii) Consider the superlinear BDP on $\mathbb{N}_0$ with birth rates $\lambda_j=j^2$ and death rates $\mu_j=j^2$ for $j\in\mathbb{N}$ . We have $\partial=\{0\}$ , $R=2$ , $\alpha=0$ , and $\beta=-1<0$ . Hence, the process is non-explosive and has certain absorption by Theorem 1 and Theorem 3(ii). By [Reference Van Doorn43], the process admits either no QSDs or a continuum family of QSDs. This shows that certain absorption does not imply exponential ergodicity of the Q-process.
Implosivity is indeed a stronger property than positive recurrence (e.g., when $R\le1$ , $\alpha<0$ ), as shown in the following example (see also Table 1).
Example 3. Let $Y_t$ be an irreducible BDP on $\mathbb{N}_0$ with $\Omega=\{1,-1\}$ and
In this case, $R=1$ and $\alpha=-1$ . By Theorems 7 and 9, $Y_t$ is positive recurrent and admits an ergodic stationary distribution, but $Y_t$ is non-implosive.
4. Applications
4.1. Stochastic reaction networks
A reaction network $(\mathcal{C},\mathcal{R})$ on a finite set $\mathcal{S}=\{S_1,\ldots,S_n\}$ (with elements called species) is an edge-labeled finite digraph with node set $\mathcal{C}$ (with elements called complexes) and edge set $\mathcal{R}$ (with elements called reactions), such that the elements of $\mathcal{C}$ are nonnegative linear combinations of species, $y=\sum_{i=1}^ny^i S_i$ , identified with vectors $y=\big(y^1,\ldots,y^n\big)$ in $\mathbb{N}_0^n$ . Reactions are directed edges between complexes, written as $y\to y^{\prime}$ . We assume that every species has a positive coefficient in some complex, and that every complex is in some reaction. Hence, the reaction network can be deduced from the reactions alone, and it is customary simply to list (or draw) the reactions. If $n=1$ , the reaction network is a one-species reaction network.
A stochastic reaction network (SRN) is a reaction network together with a CTMC X(t), $t\ge 0$ , on $\mathbb{N}_0^n$ , modeling the number of molecules of each species over time. A reaction $ y\to y^{\prime}$ fires with transition rate $\eta_{y\to y^{\prime}}(x)$ , in which case the chain jumps from $X(t)=x$ to $x+y^{\prime}-y$ [Reference Anderson and Kurtz4]. The Markov process with transition rates $\eta_{y\to y^{\prime}}\,:\,\mathbb{N}^n_0\to \mathbb{R}_{\geq 0}$ , ${y\to y^{\prime}}\in\mathcal{R}$ , has Q-matrix
Hence, the transition rate from the state x to $x+\omega$ is
For (stochastic) mass-action kinetics, the transition rate for $y\to y^{\prime}$ is
(in accordance with the transition rate introduced in the introduction for one-species reaction networks), where $x!\,:\!=\,\prod_{i=1}^nx_i!$ , and $\kappa_{y\to y^{\prime}}$ is a positive reaction rate constant [Reference Anderson, Kurtz, Koeppl and Setti3, Reference Anderson and Kurtz4]. Generally, we number the reactions and write $\kappa_1,\kappa_2,\ldots$ for convenience.
In this section, we apply the results developed in Section 3 to some examples of SRNs.
Example 4. Consider the following two reaction networks:
with $\Omega=\{-1,1\}$ in both cases and with transition rates
respectively. By Theorem 7, the first is positive recurrent and admits an exponentially ergodic stationary distribution on $\mathbb{N}_0$ , since $\alpha=-\kappa^A_2$ and $R=1$ , while by Theorem 1, the second reaction network is explosive for any initial state, since $\alpha=\kappa^B_3>0$ and $R=2$ . Indeed, these two reaction networks are structurally equivalent in the sense that there is only one irreducible component $\mathbb{N}_0$ [Reference Xu, Hansen and Wiuf48].
Example 5. Consider the following pair of SRNs from the introduction:
For the first reaction network, $R=4$ , $\alpha=0$ , and $\beta=1$ , and for the second, $R=3$ , $\alpha=0$ , and $\beta=0$ . By Theorem 1, the first is explosive for any initial state and the second does not explode for any initial state.
Example 6. (i) Consider a strongly connected reaction network:
For the underlying CTMC $Y_t$ , $\Omega=\{1,-2\}$ , and
Hence $Y_t$ is irreducible on $\mathbb{N}$ with 0 a neutral state. Moreover, $R=3$ and $\alpha=-2\kappa_3$ . By Theorem 7, there exists a unique exponentially ergodic stationary distribution on $\mathbb{N}$ .
(ii) Consider a similar reaction network including direct degradation of S:
The threshold parameters are the same as in (i). Let $\partial=\{0\}$ with $\mathbb{N}_0\setminus\partial=\mathbb{N}$ an irreducible component. By Theorem 7, the network has a uniformly exponentially ergodic QSD on $\mathbb{N}$ .
4.2. An extended class of branching processes
Consider an extended class of branching processes [Reference Chen16] with transition rate matrix $Q=(q_{x,y})_{x,y\in\mathbb{N}_0}$ :
where $\mu$ is a probability measure on $\mathbb{N}_0$ , $q_0=\sum_{y\in\mathbb{N}}q_{0,y}$ , and r(x) is a positive finite function on $\mathbb{N}_0$ . Assume the following:
-
(H1) $\mu(0)>0$ , $\mu(0)+\mu(1)<1$ ;
-
(H2) $\sum_{y\in\mathbb{N}}q_{0,y}y<\infty$ , ${E}=\sum_{k\in\mathbb{N}_0}k\mu(k)<\infty$ ;
-
(H3) r(x) is a polynomial of degree $R\ge1$ for large x.
The next theorem follows from the results in Section 3. We would like to mention that the results below provide conditions for different dynamical regimes in terms of only R and M. In contrast, the condition for positive recurrence in [Reference Chen16] also depends on the integrability of a definite integral as well as summability of a series, which nonetheless never appear.
Theorem 11. Assume ${(\textbf{H1})}$ – ${(\textbf{H3})}$ . Let $Y_t$ be a process generated by the Q-matrix given in (7) and $Y_0\neq0$ . Then $Y_t$ is non-explosive if one of the following conditions holds: (1) $R\le1$ , (2) ${E}<1$ , (3) $R=2$ , ${E}=1$ , while it is explosive with positive probability if (4) ${E}>1$ , $R>1$ . Furthermore, the following hold:
-
(i) If $q_0>0$ , then $Y_t$ is irreducible on $\mathbb{N}_0$ and is
-
(i-1) recurrent if ${E}<1$ , and transient if ${E}>1$ ;
-
(i-2) positive recurrent and exponentially ergodic if ${E}<1$ ;
-
(i-3) implosive if $R>1$ and ${E}<1$ .
-
-
(ii) If $q_0=0$ , then $\partial=\{0\}$ , and $Y_t$ has certain absorption if ${E}<1$ , while it does not if ${E}>1$ . Moreover, the process admits no QSDs if ${E}>1$ , while it admits a uniformly exponentially ergodic QSD on $\mathbb{N}$ if $R>1$ and ${E}<1$ .
Proof. For all $k\in\mathbb{N}\cup\{-1\}$ , let
By ${(\textbf{H1})}$ , $\mu(k)>0$ for some $k\in\mathbb{N}$ . Note that $\mathcal{Y}\setminus\partial$ is irreducible, with $\partial=\varnothing$ if $q_0>0$ and $\partial=\{0\}$ if $q_0=0$ . Hence, regardless of $q_0$ , by positivity of r, ${(\textbf{A1})}$ – ${(\textbf{A2})}$ are satisfied with $\Omega_-=\{-1\}$ and $\Omega_+=\{y\,:\, q_{0y}>0\}\cup(\textrm{supp}\ \mu\setminus\{0,1\}-1)$ . Moreover, ${(\textbf{H2})}$ – ${(\textbf{H3})}$ imply ${(\textbf{A3})}$ – ${(\textbf{A4})}$ . Let $r(x)=ax^R+bx^{R-1}+\textrm{O}\big(x^{R-2}\big)$ with $a>0$ . Since $R\ge1$ , it coincides with $\max\{\deg\!(\lambda_{\omega})\,:\, \omega\in\Omega\}$ . It is straightforward to verify that
where ${E^{\prime}}=\sum_{k\in\mathbb{N}}k(k-1)\mu(k)>0$ . Hence $\alpha$ has the same sign as ${E}-1$ , and $\beta<0$ whenever ${E}=1$ (or equivalently, $\alpha=0$ ). Furthermore, $\alpha=0$ implies $\gamma=0$ . In addition, in light of the fact that $R\ge1$ , the condition ${E}\ge 1$ decomposes into three possibilities:
Then the conclusions follow directly from Theorems 1, 3, 7 and 9.
Corollary 1. Assume ${(\textbf{H1})}$ – ${(\textbf{H3})}$ , that $\mu$ has finite support, and that $\{y\,:\, q_{0y}>0\}$ is finite. Then $Y_t$ is non-explosive if and only if either $R=1$ or ${E}\le1$ . Furthermore, the following hold:
-
(i) If $q_0>0$ , then $Y_t$ is irreducible and is
-
(i-1) recurrent if ${E}\le1$ , and transient otherwise;
-
(i-2) positive recurrent if and only if ${E}<1$ , or ${E}=1$ , $R>1$ , while it is null recurrent if and only if ${E}=R=1$ ; furthermore, $Y_t$ is exponentially ergodic if ${E}<1$ , or if ${E}=1$ , $R>2$ ;
-
(i-3) implosive if and only if either of the two conditions $R>1$ , ${E}<1$ , or $R>2$ , ${E}=1$ holds.
-
-
(ii) If $q_0=0$ , then $\partial=\{0\}$ , and $Y_t$ has certain absorption if and only if ${E}\le1$ . Moreover, the process admits no QSDs if ${E}>1$ , or ${E}=R=1$ , while it admits a uniformly exponentially ergodic QSD on $\mathbb{N}$ if either $R>1$ , ${E}<1$ , or $R>2$ , ${E}=1$ .
The extended branching process under more general assumptions (allowing more general forms of r) is addressed in [Reference Chen16]. In that reference, the conditions given for the dynamic behavior of the process seem more involved than here and even become void in some situations (e.g., in [Reference Chen16, Corollary 1.5(iii)], where the definite integral indeed is always infinite under ( $\textbf{H1}$ )–( $\textbf{H3}$ ).)
4.3. A general single-cell stochastic gene expression model
To model single-cell stochastic gene expression with bursty production, we propose the following one-species generalized reaction network (consisting potentially of infinitely many reactions) with mass-action kinetics:
where $c_m\ge0$ for $m=0,\ldots,J_1$ , $r_m\ge0$ for $m=1,\ldots,J_2$ , $J_1\in\mathbb{N}_0,\ J_2\in\mathbb{N}$ , and $\mu_m$ , for $m=0,\ldots,J_1$ , are probability distributions on $\mathbb{N}$ . Assume the following:
-
(H4) $J_1\le J_2$ , $c_0>0$ , $c_{J_1}>0$ , $r_1>0$ , and $r_{J_2}>0$ .
-
(H5) ${E}_m=\sum_{k=1}^{\infty}k\mu_m(k)<\infty$ , for $m=0,\ldots,J_1$ .
This network encompasses several single-cell stochastic gene expression models in the presence of bursting; see e.g. [Reference Be’er and Assaf8, Reference Chen and Jia17, Reference Mackey, Tyran-Kamińska and Yvinec31, Reference Shahrezaei and Swain41]. Ergodicity and an exact formula for the ergodic stationary distribution (when it exists) are the main concerns of these references. The first set of $J_1$ reactions accounts for bursty production of mRNA copies with transcription rate $c_m$ and burst size distribution $\mu_m$ . The second set of $J_2$ reactions accounts for degradation of mRNA with degradation rate $r_m$ [Reference Chen and Jia17, Reference Shahrezaei and Swain41].
The network (8) reduces to the specific model studied in the following references:
-
[Reference Chen and Jia17, Section 4] (see also [Reference Mackey, Tyran-Kamińska and Yvinec31, Section 3.2]), when $J_1=0$ , $J_2=1$ , and $\mu_0$ is a geometric distribution.
-
[Reference Schwabe, Rybakova and Bruggeman40], when $J_1=0$ , $J_2=1$ , and $\mu_0$ is a negative binomial distribution.
-
[Reference Mackey, Tyran-Kamińska and Yvinec31, Example 3.6], when $J_1=J_2=1$ , and $\mu_0=\mu_1$ are geometric distributions.
-
[Reference Falk, Mendler and Drossel21] when $J_1=2$ , $J_2=3$ , $\mu_0=\delta_1$ , $\mu_2=\delta_k$ for some $k\in\mathbb{N}$ , and $c_1=r_2=0$ . Here $\delta_i$ is the Dirac delta measure at i.
Theorem 12. Assume ${(\textbf{H4})}$ – ${(\textbf{H5})}$ , and that $\mu_m$ has finite support whenever $c_m>0$ for $m=0,\ldots,J_1$ . Then the process $Y_t$ associated with the network (8) is irreducible on $\mathbb{N}_0$ , and it is positive recurrent and there exists an ergodic stationary distribution on $\mathbb{N}_0$ if and only if one of the following conditions holds:
-
(i) $J_1<J_2$ ,
-
(ii) $J_1=J_2$ and $c_{J_2}{E}_{J_2}<r_{J_2}$ ,
-
(iii) $J_1=J_2>2$ , $c_{J_2}{E}_{J_2}=r_{J_2}$ , and $c_{J_2-1}{E}_{J_2-1}\le r_{J_2-1}+\frac{1}{2}c_{J_2}\big({E}_{J_2}+{E}^{\prime}_{J_2}\big)$ ,
-
(iv) $J_1=J_2=2$ , $c_{J_2}{E}_{J_2}=r_{J_2}$ , and $c_{J_2-1}{E}_{J_2-1}<r_{J_2-1}+\frac{1}{2}c_{J_2}\big({E}_{J_2}+{E}^{\prime}_{J_2}\big)$ ,
-
(v) $J_1=J_2=1$ , $c_{J_2}{E}_{J_2}=r_{J_2}$ , and $c_{J_2-1}{E}_{J_2-1}<r_{J_2-1}$ ,
where ${E}^{\prime}_m=\sum_{k=1}^{\infty}k^2\mu_m(k)$ . Moreover, the stationary distribution is exponentially ergodic if one of (i), (ii), and (iii) holds. Furthermore, the process $Y_t$ is implosive if and only if (iii) or $J_2>1$ with (i) or (ii).
Proof. We have $\Omega=\{-1\}\cup\Big(\cup_{j=0}^{J_1}\textrm{supp}\ \mu_j\Big)$ , and
for $k\in\mathbb{N}$ and $x\in\mathbb{N}_0$ , where $x^{\,\underline{j}}=\prod_{i=0}^{j-1}(x-i)$ is the descending factorial. By ${(\textbf{H4})}$ , ${(\textbf{A1})}$ – ${(\textbf{A2})}$ are satisfied; moreover, the irreducibility of $Y_t$ also follows from [Reference Xu, Hansen and Wiuf48]. Under ${(\textbf{H5})}$ , the mass-action kinetics yield ${(\textbf{A3})}$ – ${(\textbf{A4})}$ . Since $J_1\le J_2$ by ${(\textbf{H4})}$ , we have $R=J_2\ge1$ . Since
we have $\alpha=c_{J_1}{E}_{J_1}\delta_{J_1,J_2}-r_{J_2}$ , where $\delta_{i,j}$ is the Kronecker delta. When $\alpha=0$ , we have
The combination of the conditions (i) and (ii) is equivalent to $\alpha<0$ ; the condition (iii) is equivalent to $R>2$ , $\alpha=0$ , $\beta\le0$ ; the condition (iv) is equivalent to $R=2$ , $\alpha=0$ , $\beta<0$ ; the condition (v) is equivalent to $R=1$ , $\alpha=0$ , $\gamma<0$ . The conclusions then follow from Theorems 7 and 9.
Corollary 2. Assume ${(\textbf{H4})}$ – ${(\textbf{H5})}$ , and that $\mu_m$ has infinite support and $c_m>0$ for some $m=0,\ldots,J_1$ . Then $Y_t$ is irreducible on $\mathbb{N}_0$ , and is positive recurrent with an ergodic stationary distribution if one of Theorem 12(i), Theorem 12(ii), and Theorem 12(v) holds. Moreover, the stationary distribution is exponentially ergodic if either Theorem 12(i) or Theorem 12(ii) holds. In addition, the process $Y_t$ is implosive if $J_2>1$ .
Proof. Based on the proof of Theorem 12, the conclusions follow directly from Corollaries 8 and 10.
4.4. Stochastic populations under bursty reproduction
Two stochastic population models with bursty reproduction are investigated in [Reference Be’er and Assaf8].
The first model is a Verhulst logistic population process with bursty reproduction. The process $Y_t$ is a CTMC on $\mathbb{N}_0$ with transition rate matrix $Q=(q_{x,y})_{x,y\in\mathbb{N}_0}$ satisfying
where $c>0$ is the reproduction rate, $K\in\mathbb{N}$ is the typical population size in the long-lived metastable state prior to extinction [Reference Be’er and Assaf8], and $\mu$ is the burst size distribution. Assume the following:
(H6) ${E_b}=\sum_{k=1}^{\infty}k\mu(k)<\infty$ .
Approximations of the mean time to extinction and QSD are discussed in [Reference Be’er and Assaf8] against various burst size distributions of finite mean (e.g., Dirac measure, Poisson distribution, geometric distribution, negative-binomial distribution). Nevertheless, the existence of a QSD is not proved there. Here we prove the certain absorption and ergodicity of the QSD for this population model.
Theorem 13. Assume ${(\textbf{H6})}$ . The Verhulst logistic model $Y_t$ with bursty reproduction has certain absorption. Moreover, there exists a uniformly exponentially ergodic QSD on $\mathbb{N}$ trapped to zero.
Proof. We have $\Omega=\textrm{supp}\ \mu\cup\{-1\}$ , $\lambda_{-1}(x)=\frac{c}{K}x^2+x$ , $\lambda_k(x)=c\mu(k)x$ , for $k\in\mathbb{N}$ and $x\in\mathbb{N}$ . Let $\partial=\{0\}$ ; then $\mathbb{N}_0\setminus\partial=\mathbb{N}$ is irreducible [Reference Xu, Hansen and Wiuf48]. Hence ${(\textbf{A1})}$ – ${(\textbf{A5})}$ are satisfied. Moreover, $R=2$ , $\alpha=-\frac{c}{K}<0$ , and thus the conclusions follow from Theorems 3 and 7, together with Corollaries 4 and 8 for finite $\textrm{supp}\ \mu$ and infinite $\textrm{supp}\ \mu$ , respectively.
The second model is a runaway model of a stochastic population including bursty pair reproduction [Reference Be’er and Assaf8]. This model can be described as a generalized reaction network, where c, K, and $\mu$ are defined as in the first model. The survival probability of this population model is addressed in [Reference Be’er and Assaf8]. Nevertheless, it turns out that this model is explosive for any initial state.
Theorem 14. Assume ${(\textbf{H6})}$ . The runaway model is explosive.
Proof. We have $\Omega=\textrm{supp}\ \mu\cup\{-1\}$ , $\lambda_k(x)=\frac{c}{K}\mu(k)x(x-1)$ , $\lambda_{-1}(x)=x$ , for $k\in\mathbb{N}$ and $x\in\mathbb{N}$ . Let $\partial=\{0,1\}$ . Then $\mathbb{N}_0\setminus\partial=\mathbb{N}\setminus\{1\}$ is irreducible [Reference Xu, Hansen and Wiuf48]. Hence ${(\textbf{A5})}$ is valid. Moreover, it is easy to verify that ${(\textbf{A1})}$ – ${(\textbf{A4})}$ are also satisfied. In addition, $R=2$ , $\alpha=\frac{c}{K}{E_b}>0$ , and thus the conclusions follow from Theorem 1 and Corollary 2 for finite $\textrm{supp}\ \mu$ and infinite $\textrm{supp}\ \mu$ , respectively.
5. Proofs
5.1. Proof of Theorem 1
Hereafter, we use the notation $[m,n]_1$ ( $[m,n[_1$ , etc.) for the set of consecutive integers from m to n, with $m,n\in\mathbb{N}_0\cup\{+\infty\}$ . The notation is adopted from [Reference Xu, Hansen and Wiuf48].
Moreover, throughout the proofs, we assume without loss of generality that $\mathcal{Y}=\mathbb{N}_0$ , and $\partial\subseteq\{0\}$ for ease of exposition. Indeed, the dynamical properties of the CTMCs discussed in this paper depend only on the transition structure of the states $x\in\mathcal{Y}$ with large values of x. By Assumptions ( $\textbf{A4}$ )–( $\textbf{A5}$ ), all jumps in $\Omega$ are possible. When $\partial\neq\varnothing$ , it is standard to ‘glue’ all states in $\partial$ to be a single state 0, since $\partial$ is finite and the set of states in $\mathcal{Y}\setminus\partial$ one jump away from $\partial$ is also finite, by ( $\textbf{A2}$ ).
We prove the conclusions case by case.
(a) Assume $\partial=\varnothing$ . Then $Y_t$ is irreducible on $\mathbb{N}_0$ , and one can directly apply Propositions 2 and 3 with appropriate Lyapunov functions to be determined.
(b) Assume $\partial\neq\varnothing$ . Let $Z_t$ be the irreducible CTMC on the state space $\mathcal{Y}\setminus\partial$ with $Z_0=Y_0$ and transition operator $\widetilde{Q}$ being Q restricted to $\mathcal{Y}\setminus\partial$ :
In the following, we show that $Z_t$ is explosive if and only if $Y_t$ is explosive, and hence the case (b) reduces to the case (a). This equivalence is not quite trivial. There is a positive probability that, starting from any non-absorbing state, the chain will jump to an absorbing state in a finite number of steps. So we need to show that this will not happen with probability one. Otherwise, explosivity is not possible. Assume first that $Z_t$ is explosive. Then $\widetilde{Q}v=v$ for some bounded nonnegative nonzero v. Let $u_x=v_x \unicode{x1D7D9}_{\mathcal{Y}\setminus\partial}(x),\ \forall x\in\mathcal{Y}$ . It is straightforward to verify that $Qu=u$ . By Proposition 4, $Y_t$ is also explosive. Conversely, assume that $Y_t$ is explosive; then $Qu=u$ for some bounded nonnegative nonzero u. Let $w=u|_{\partial}$ , i.e., $w_x=u_x$ for all $x\in\partial$ . Since $Q|_{\partial}$ is a lower-triangular matrix with nonpositive diagonal entries, and $w\ge0$ , it is readily deduced that $w=0$ by Gaussian elimination. This implies from $Qu=u$ that $\widetilde{Q}v=v$ with $v=u|_{\mathcal{Y}\setminus\partial}$ . Hence $Z_t$ is explosive. To sum up, $Z_t$ is explosive if and only if $Y_t$ is explosive.
Based on the above analysis, it remains to prove the conclusions for the case (a) using Propositions 2 and 3. We first prove the conclusions assuming $\Omega$ is finite.
(i) We prove explosivity by Proposition 2. Let the lattice interval $A=[0,x_0-\min\!\Omega_-[_1$ for some $x_0>1$ to be determined. Since $\#\Omega_-<\infty$ , $A\subseteq\mathbb{N}_0$ is finite. Let f be decreasing and bounded such that $f(x)=\unicode{x1D7D9}_{[0,x_0[_1}(x)+x^{-\delta}\unicode{x1D7D9}_{[x_0,\infty[_1}(x)$ for all $x\ge x_0$ , with $\delta>0$ to be determined. Obviously, Proposition 2(i) is satisfied for the set A. Next we verify the conditions in Proposition 2(ii). It is easy to verify by straightforward calculation that
where $\epsilon=\delta\alpha/2$ provided (C1) holds with $\delta<R-1$ , or $\epsilon=\delta\!\left(\beta-\delta\vartheta\right)/2$ provided (C2) holds with $\delta<\min\{\beta/\vartheta,R-2\}$ , and $x_0$ is chosen large enough. Since $\delta>0$ can be arbitrarily small, in either case, there exist $\delta$ and $\epsilon$ such that the conditions in Proposition 2 are fulfilled, and thus $\mathbb{E}_x\zeta<+\infty$ for all $x\in\mathcal{Y}$ . In particular, $Y_t$ is explosive a.s.
(ii) Now we prove non-explosivity using Proposition 3. Let $f(x)=\log\log\!(x+1)$ and $g(x)=(|\alpha|+|\beta|+1)(x+M)$ for all $x\in\mathbb{N}_0$ , with some $M>0$ to be determined. One can show that all the conditions in Proposition 3 are satisfied with some large constant $M>0$ , provided neither (C1) nor (C2) holds. Hence $Y_t$ is non-explosive.
5.2. Proof of Theorem 2
We first prove explosivity under the condition (C1). Let f be as in the proof of Theorem 1(i), and let
Then $\alpha=\alpha_++\alpha_-$ . Since $\alpha>0$ , we have $R_+=R$ and there exists $\epsilon_0\in]0,1[$ such that $\alpha_-+(1-\epsilon_0)\alpha_+>0$ . By Proposition 1, there exist $N_0,\ u^{\prime}\in\mathbb{N}$ such that
By ( $\textbf{A3}$ ), $\Omega\setminus]N_0,\infty[_1$ is finite. Hence, choosing $x_0\ge u^{\prime}$ large, we have for all $x\in\mathcal{Y}\setminus A$
where
and $\delta<R-1$ . The rest of the argument is the same as that of the proof of Theorem 1(i).
Next, we prove non-explosivity under (C3), (C4), or (C5). Let f and g be as in (ii) in the proof of Theorem 1. By ( $\textbf{A3}$ ), for some large $M>0$ to be determined, for all $x\in\mathbb{N}_0$ ,
provided (C3), (C4), or (C5) holds. The rest of the proof is the same as that of Theorem 1(ii).
5.3. Proof of Theorem 3
Let $h_A=\mathbb{P}_{Y_0}(\tau_A<\infty)$ be the hitting probability [Reference Norris36]. In particular, $h_A$ is called the absorption probability if A is a closed communicating class. To verify conditions for certain absorption requires the following property for hitting probabilities. For any set $A\subseteq\mathcal{Y}$ , we write $h_A(i)$ for $h_A$ , to emphasize the dependence of the hitting probability on the initial state $i\in\mathcal{Y}$ . In particular, if $A=\{x\}$ is a singleton, we simply write $h_x$ for $h_A$ .
Assume without loss of generality that $\mathcal{Y}=\mathbb{N}_0$ , and $\partial\subseteq\{0\}$ .
(i) We first show recurrence and transience. The idea of applying the classical semimartingale approach originates from [Reference Lamperti29]. It suffices to show recurrence and transience for the embedded discrete-time Markov chai