1. Introduction
Stochastic biological models based on continuous-time Markov chains (CTMCs) are commonly used to model complex cellular behaviours, gene expression, and the evolution of DNA [Reference Ewens19]. Noise is an inherent extrinsic and intrinsic property of such biological systems and cannot be ignored without compromising the conclusions and the accuracy of the models.
In many cases it is reasonable to expect well-behaved biological systems, and thus also well-behaved stochastic models. Therefore, in these cases, it is natural to assume the modelling CTMC is ergodic, that is, that there exists a unique stationary distribution which describes the system in the long run. In other cases, for example for population processes without immigration, the population eventually goes extinct almost surely, and thus the ergodic stationary distribution is trivially the Dirac delta measure at zero. In these cases, it makes sense to study the quasi-stationary distribution (QSD), that is, the long-time behaviour of the process before extinction (usually called the Q-process) [Reference Méléard and Villemonais34]. Jointly, stationary distributions and QSDs are referred to as limit distributions in the present paper.
The stationary distribution (provided it exists) is generally difficult to state in explicit form, except in a few cases. If the underlying stochastic process has a detailed balanced structure, then the stationary distribution takes a product form; see [Reference Anderson, Craciun and Kurtz2, Reference Cappelletti and Wiuf9, Reference Gelenbe and Mitrani23, Reference Kelly30, Reference Mairesse and Nguyen32, Reference Whittle46] for various network models, [Reference Méléard and Villemonais34] for birth–death processes (BDPs), and [Reference Economou18, Reference Gelenbe22, Reference Hoessly and Mazza28, Reference Neuts38, Reference Ramaswami42] for generalizations of such processes. QSDs with explicit expressions appear in even rarer cases [Reference Van Doorn45].
While an explicit expression might not be known in general, less suffices in many cases. For example, if an expression for the tail distribution is known, then the existence and relative sizes of moments could be assessed from the decay rate of the tail distribution. Additionally, relative recurrence times could be assessed for stationary distributions of CTMCs.
With this in mind, we establish results for the tail behaviour of stationary distributions and QSDs, provided such exist. In particular, we concentrate on CTMCs on the non-negative integers $\mathbb{N}_0$ with asymptotic power law transition rate functions (to be made precise). Our approach is based on generic identities derived from the flux-balance equation [Reference Kelly30] for limit distributions and stationary measures (Theorems 1 and 2). The identity for stationary distributions/measures might be seen as a difference equation, which has order one less than the difference equation obtained directly from the master equation. More interestingly, for any given state x, the left-hand side of the identity consists of terms evaluated in states $\ge x$ only, while the right-hand side of the identity consists of terms evaluated in states $<x$ only, and all terms have non-negative coefficients. For BDPs, the identities coincide with the classical recursive expressions for stationary distributions and QSDs.
Furthermore, the identities allow us to study the tail behaviour of limit distributions, provided they exist, and to characterize their forms (Theorems 4 and 3). Specifically, in Section 5, for CTMCs with transition rate functions that are asymptotically power laws, we show that there are only three regimes: the decay follows either (i) a Conley–Maxwell–Poisson distribution (light-tailed), (ii) a geometric distribution, or (iii) a sub-exponential distribution. Similar trichotomy results appear in the literature in other contexts [Reference Asmussen5, Reference Maulik and Zwart33], but, to our knowledge, only for stationary distributions. Our approach is based on repeated use of the identities we establish, combined with certain combinatorial identities (e.g., Lemma 1).
Importantly, we successfully obtain QSD tail asymptotics. To the best of our knowledge, no similar classification results on QSD tail asymptotics have been established, despite the fact that the Lyapunov function approach has been used frequently to establish ergodicity of QSDs [Reference Champagnat and Villemonais12]. A superficial reason may be the inherent difference between stationary distributions and QSDs: stationary distributions are equilibria of the master equation while QSDs are not. A deeper explanation may be that the drift of CTMCs (or discrete-time Markov chains (DTMCs)) may not directly relate to the decay rate of QSDs. For an absorbing CTMC, known conditions for exponential ergodicity of stationary distributions are not even sufficient to establish uniqueness of QSDs [Reference Van Doorn45].
This difference is also reflected in our results, where an extra condition is required to establish QSD tail asymptotics (Theorems 3 and 4). Our novel approach successfully addresses tail asymptotics of both stationary distributions and QSDs at the same time, based on the similarity of the algebraic equations they satisfy (Theorems 1 and 2).
We apply our main results to biochemical reaction networks, a general single-cell stochastic gene expression model, an extended class of branching processes, and stochastic population processes with bursty reproduction, none of which are BDPs.
The Lyapunov function approach is widely taken as a standard method to prove ergodicity of Markov processes [Reference Meyn and Tweedie36]. Additionally, the approach has been used to obtain tail asymptotics of stationary distributions, assuming ergodicity [Reference Aspandiiarov and Iasnogorodski6, Reference Bertsimas, Gamarnik and Tsitsiklis8, Reference Denisov, Korshunov and Wachtel16, Reference Denisov, Korshunov and Wachtel17, Reference Menshikov and Popov35, Reference Xu, Hansen and Wiuf47]. In contrast, our results do not require ergodicity or any other finite moment condition. For DTMCs, ergodicity follows from the existence of a stationary distribution on an irreducible set [Reference Miller37]; in contrast, this is not true for CTMCs [Reference Miller37]. In particular, for explosive Markov chains, potentially with more than one stationary distribution, the cited results fail, whereas our results also hold in this case.
The trichotomy pattern which we observe for the tail asymptotics is not surprising, as it has already been observed for BDPs [Reference Takine44], as well as for processes on continuous state spaces, such as the Lindley process [Reference Asmussen5] and the exponential functional of a Lévy process drifting to infinity [Reference Maulik and Zwart33]. The techniques applied in these papers do not seem applicable in our setting.
It is noteworthy that the identities we establish could be used to calculate the limit distribution recursively, up to an error term that depends on only a few generating terms ( $\pi(0)$ in the case of BDPs) and the truncation point of the limit distribution. The error term is given by the tail distribution, and thus the decay rate of the error term can be inferred from the present work. Approximation of limit distributions will be pursued in a subsequent paper. The main results of this paper, together with the approach, can be extended to DTMCs in a straightforward way.
2. Preliminaries
2.1. Sets and functions
Denote the sets of real numbers, positive real numbers, integers, positive integers, and non-negative integers by $\mathbb{R}$ , $\mathbb{R}_+$ , $\mathbb{Z}$ , $\mathbb{N}$ , and $\mathbb{N}_0$ , respectively. For $m, n\in\mathbb{N}$ , let $\mathbb{R}^{m\times n}$ denote the set of m-by-n matrices over $\mathbb{R}$ . Furthermore, for any set B, let $\#B$ denote its cardinality and $\unicode{x1d7d9}_B$ the corresponding indicator function. For $b\in\mathbb{R}$ , $A\subseteq\mathbb{R}$ , let $bA=\{ba\,:\, a\in A\}$ and $A+b=\{a+b\,:\, a\in A\}$ . Given $A\subseteq\mathbb{R}$ , let $\inf A$ and $\sup A$ denote the infimum and supremum of the set, respectively. By convention, $\sup A=-\infty$ and $\inf A=+\infty$ if $A=\varnothing$ ; $\sup A=+\infty$ if A is unbounded from above; and $\inf A=-\infty$ if A is unbounded from below. For $x, y\in\mathbb{N}_0$ with $x\ge y$ , let $x^{\underline{y}}=\frac{x!}{(x-y)!}$ .
Let f and g be non-negative functions on an unbounded set $A\subseteq\mathbb{N}_0$ . We write $f(x)\lesssim g(x)$ if there exist $C,N>0$ such that
that is, $f(x)=\mathrm{O}(g(x))$ since f is non-negative. (Here O refers to the standard big-O notation.) The function f is said to be asymptotic power law (APL) if there exists $r_1\in\mathbb{R}$ such that $\lim_{x\to\infty}\frac{f(x)}{x^{r_1}}=a$ exists and is finite. Hence $r_1=\lim_{x\to\infty}\frac{\log f(x)}{\log x}$ . An APL function f is called hierarchical (APLH) on A with $(r_1,r_2,r_3)$ if there further exist $r_2, r_3$ with $r_2+1\ge r_1>r_2>r_3\ge r_1-2$ , and $a>0$ , $b\in\mathbb{R}$ , such that for all large $x\in A$ ,
The requirement $r_2+1\ge r_1$ and $r_3\ge r_1-2$ comes from the analysis in Sections 6–7, where asymptotic Taylor expansions of functions involve the powers of the first few leading terms. Here $r_1$ , $r_2$ , and $r_3$ are called the first, second, and third power of f, respectively. All rational functions, polynomials, and real analytic APL functions are APLH. Not all APL functions are APLH; e.g., $f(x)=(1+(\!\log(x+1))^{-1})x$ on $\mathbb{N}$ is not APLH.
Given a function f which is APLH, the first power in the expansion is uniquely determined, while the other two powers are not. In other words, given the asymptotic expansion (1), there may exist a family of APLH functions admitting the same asymptotic expansion. Let $r_2^*$ and $r_3^*$ be the infima over all $(r_2,r_3)\in\mathbb{R}^2$ such that f is APLH on A with $(r_1,r_2,r_3)$ . For convention, we always choose the minimal powers $\big(r_1,r_2^*,r_3^*\big)$ , whenever f is APLH with $\big(r_1,r_2^*,r_3^*\big)$ . As an example, $f(x)=x^2+3x+4$ is APLH on $\mathbb{N}_0$ with $(2,r_2,r_3)$ for any $2>r_2>r_3\ge 1$ (in which case $b=0$ ) or $1=r_2>r_3\ge 0$ ( $b=3$ ). In this case, f is APLH on $\mathbb{N}_0$ with minimal powers $\big(r_1,r_2^*,r_3^*\big)=(2,1,0)$ . In contrast, take $f(x)=x+x^{1/3}\log x$ . Then f is APLH on $\mathbb{N}_0$ with $(1,r_2,r_3)$ for any $r_2>r_3>1/3$ ( $b=0$ ). In this case, $r_1=1$ and $r_2^*=r_3^*=1/3$ , but f is not APLH on $\mathbb{N}$ with $(1,1/3,1/3)$ . For any real analytic APLH function f on $\mathbb{N}_0$ , f is APLH on $\mathbb{N}_0$ with $\big(r_1,r_2^*,r_3^*\big)$ , where $r_1=\lim_{x\to\infty}\frac{\log f(x)}{\log x}$ , $r_2^*=r_1-1$ , and $r_3^*=r_1-2$ .
2.2. Markov chains
Let $(Y_t\,:\, t\ge 0)$ (or $Y_t$ for short) be a CTMC with state space $\mathcal{Y}\subseteq\mathbb{N}_0$ and transition rate matrix $Q=(q(x,y))_{x,y\in\mathcal{Y}}$ ; in particular, each entry is finite. Recall that 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 Norris39]. Assume $\partial\subsetneq\mathcal{Y}$ is a (possibly empty) finite closed absorbing set, which is to the left of $\partial^{\textsf{c}}=\mathcal{Y}\setminus\partial$ . Here, the relative position of $\partial$ to $\partial^{\textsf{c}}$ ensures that the only way for the Markov chain to end in an absorbing state is by jumping from a transient state backward to an absorbing state (this property is used in Proposition 1 below). Furthermore, define the set of jump vectors
and let $\Omega_{\pm}=\{\omega\in\Omega\,:\, {\text{sgn}}(\omega)=\pm1\}$ be the sets of forward and backward jump vectors, respectively.
For any probability distribution $\mu$ on $\partial^{\textsf{c}}$ , define
where $\mathbb{P}_x$ denotes the probability measure of $Y_t$ with initial condition $Y_0=x\in\mathcal{Y}$ . Any positive measure $\mu$ on a set $A\subseteq\mathbb{N}_0$ can be extended naturally to a positive measure on $\mathbb{N}_0$ with no mass outside A, $\mu(\mathbb{N}_0\setminus\! A)=0$ .
A (probability) measure $\pi$ on $\mathcal{Y}$ is a stationary measure (distribution) of $Y_t$ if it is a non-negative equilibrium of the so-called master equation [Reference Gillespie24]:
(Here and elsewhere, functions defined on $\mathcal{Y}$ are put to zero when evaluated at $x\not\in\mathcal{Y}\subseteq\mathbb{Z}$ .) Any stationary distribution $\pi$ of $Y_t$ satisfies
for all $t\ge0$ , where $2^{\mathcal{Y}}$ is the power set of $\mathcal{Y}$ [Reference Champagnat and Villemonais12].
Let $\tau_{\partial}=\inf\{t>0\,:\, Y_t\in\partial\}$ be the entrance time of $Y_t$ into the absorbing set $\partial$ . We say $Y_t$ admits certain absorption if $\tau_{\partial}<\infty$ almost surely for all $Y_0\in\partial^{\textsf{c}}$ . Moreover, the process associated with $Y_t$ conditioned never to be absorbed is called a Q-process [Reference Champagnat and Villemonais11]. A probability measure $\nu$ on $\partial^{\textsf{c}}$ is a quasi-stationary distribution (QSD) of $Y_t$ if for all $t\ge0$ ,
3. Identities for limit distributions
Let $\omega_*=\text{gcd}(\Omega)$ be the (unique) positive greatest common divisor of $\Omega$ . Define the scaled largest positive and negative jump, respectively, as
Furthermore, for $j\in\{\omega_-,\ldots,\omega_++1\}$ , define
Hence, $\varnothing=A_{\omega_-}\subseteq A_j\subseteq A_{j+1}\subseteq A_0=\Omega_-$ for $\omega_-< j<0$ , and $\varnothing=A_{\omega_++1}\subseteq A_{j+1}\subseteq A_j\subseteq A_1=\Omega_+$ for $1< j< \omega_+$ .
The following classical result provides a necessary condition for QSDs.
Proposition 1. ([Reference Collet, Martínez and San Martín15].) Assume $\partial\neq\varnothing$ . Let $\nu$ be a QSD of $Y_t$ on $\partial^{\textsf{c}}$ . Then for $x\in\mathbb{N}_0\!\setminus\!\partial$ ,
where
is finite.
Proof. For $\partial=\{0\}$ and $x\in\partial^{\textsf{c}}$ , the identity follows from [Reference Collet, Martínez and San Martín15]. The same argument applies if $\partial^{\textsf{c}}\not=\{0\}$ . For $x\in\mathbb{N}_0\setminus\mathcal{Y}$ , the identity holds trivially (with both sides being zero), which can be argued similarly to the proof of Theorem 1.
Before stating a generic identity for stationary distributions which originates from the master equation of the CTMC, we provide an example.
Example 1. Consider a CTMC on $\mathbb{N}_0$ with $\Omega=\{1,-2\}$ and transition rate functions
where $\kappa_i$ ( $i=1,2,3$ ) are positive constants. Such rate functions can be associated with a weakly reversible stochastic reaction network (see Example 10 in Section 5). This CTMC is ergodic on $\mathbb{N}$ and there exists a stationary distribution $\pi$ , which solves the master equation
where $a(x)=\kappa_3x^{\underline{3}}$ and $b(x)=(\kappa_1x+\kappa_2x^{\underline{2}})$ . One can show by induction that the terms in (4) can be separated so that those with $\pi(y)$ for $y\ge x$ and those with $\pi(y)$ for $y<x$ are on different sides of the equality, in such a way that all coefficients are positive. Naive rearrangement of the terms does not suffice. We find
In addition, to our surprise, we observe that (4) is a third-order difference equation, while (5) is a second-order difference equation. Hence, the identity (5), though equivalent to (4), helps simplify the relationship among all terms of the stationary distribution.
The following generic identities generalize the observations in Example 1, and provide an equivalent definition of a stationary distribution, perhaps in a more handy form. The identities are the so-called flux-balance equation [Reference Kelly30, Lemma 1.4] specialized to one dimension. We emphasize that stationary measures are not necessary finite, while stationary distributions and QSDs are probability distributions.
Theorem 1. The following statements are equivalent:
-
(1) $\pi$ is a stationary measure (stationary distribution) of $Y_t$ on $\mathcal{Y}$ ;
-
(2) $\pi$ is a positive measure (probability distribution) on $\mathcal{Y}$ satisfying, for $x\in\mathbb{N}_0$ ,
\begin{multline*}\sum_{\omega\in\Omega_-}\sum_{j=\omega\omega^{-1}_*+1}^{0}q(x-j\omega_*, x-j\omega_*+\omega)\pi(x-j\omega_*)\\=\sum_{\omega\in\Omega_+}\!\sum_{j=1}^{\omega\omega^{-1}_*}q(x-j\omega_*,x-j\omega_*+\omega)\pi(x-j\omega_*)<\infty;\end{multline*} -
(3) $\pi$ is a positive measure (probability distribution) on $\mathcal{Y}$ satisfying, for $x\in\mathbb{N}_0$ ,
\begin{multline*}\sum_{j=\omega_-+1}^{0}\pi\left(x-j\omega_*\right)\sum_{\omega\in A_j}q(x-j\omega_*,x-j\omega_*+\omega)\\= \sum_{j=1}^{\omega_+}\pi\left(x-j\omega_*\right)\sum_{\omega\in A_j}q(x-j\omega_*,x-j\omega_*+\omega)<\infty.\end{multline*}
Proof. We only prove the equivalent representations for stationary measures; the equivalent representations for stationary distributions then follow. We assume without loss of generality that $\omega_*=1$ .
Recall from [Reference Kelly30, Lemma 1.4] that $\pi$ is a stationary measure satisfying (2) if and only if for any $A\subseteq\mathcal{Y}$ ,
where $A^c=\mathcal{Y}\setminus A$ . (In [Reference Kelly30], the identity (6) is proved for finite $\mathcal{Y}$ . It can be shown by induction that it also holds for countable $\mathcal{Y}$ .) Let $\pi(x)=0$ for $x\notin \mathcal{Y}$ and $q(x,y)=0$ for $(x,y)\notin\mathcal{Y}^2$ ; then (6) holds for $A\subseteq\mathcal{Y}$ if and only if for $A\subseteq\mathbb{Z}$ ,
In particular, for any $x\in \mathbb{N}_0$ , let $A=\{y\in\mathbb{Z}\,:\, y\ge x\}$ . Then (7) implies that
Conversely, (8) implies (2) (i.e., $\pi$ is a stationary measure), which is shown by subtracting from both sides of (8) the same equations but with x replaced by $x+1$ . This is possible since both sides are finite. We now compare (8) with Theorem 1(3). Recall that $A_{\omega_-} = A_{\omega_++1} = \varnothing$ . In the latter equation, let $y=x-j$ and $z=x-j+\omega$ ; then Theorem 1(3) follows from (8). Furthermore, let $j=x-y$ and $\omega=z-y$ , then (8) is obtained from Theorem 1(3). Observe that Theorem 1(2) and Theorem 1(3) are equivalent by Fubini’s theorem.
A special form of Theorem 1 under more assumptions has been stated in the context of stochastic reaction networks [Reference Hansen26, Proposition 5.4.9].
For any positive measure $\mu$ on $\mathbb{N}_0$ , let
be the tail distribution (or simply the tail) of $\mu$ .
The following identities for QSDs also present equivalent definitions of the latter.
Theorem 2. Assume $\partial\neq\varnothing$ . Then the following statements are equivalent:
-
(1) $\nu$ is a QSD of $Y_t$ on $\partial^{\textsf{c}}$ ;
-
(2) $\nu$ is a probability measure on $\partial^{\textsf{c}}$ , and for $x\in\mathbb{N}_0\setminus\partial$ ,
\begin{align*}&\sum_{\omega\in\Omega_-}\sum_{j=\omega\omega_*^{-1}+1}^{0}q(x-j\omega_*,x-j\omega_*+\omega)\nu(x-j\omega_*)=\theta_{\nu}T_{\nu}(x)\\&+\sum_{\omega\in\Omega_+}\!\sum_{j=1}^{\omega\omega_*^{-1}}q(x-j\omega_*,(x-j\omega_*+\omega)\nu(x-j\omega_*)<\infty;\end{align*} -
(3) $\nu$ is a probability measure on $\partial^{\textsf{c}}$ , and for $x\in\mathbb{N}_0\setminus\partial$ ,
\begin{align*}&\sum_{j=\omega_-+1}^{0}\nu\left(x-j\omega_*\right)\sum_{\omega\in A_j}q(x-j\omega_*,x-j\omega_*+\omega)=\theta_{\nu}T_{\nu}(x)\\&+\sum_{j=1}^{\omega_+}\nu\left(x-j\omega_*\right)\sum_{\omega\in A_j}q(x-j\omega_*,x-j\omega_*+\omega)<\infty.\end{align*}
Proof. The proof is similar to that of Theorem 1 and is thus omitted.
If a CTMC jumps unidirectionally (e.g., it is a pure birth or a pure death process), then all stationary measures, if such exist, are concentrated on absorbing states [Reference Xu, Hansen and Wiuf48]. In contrast, an absorbed pure birth process has no QSDs. However, an absorbed pure death process may admit one or more QSDs, as illustrated below.
Example 2. Consider a pure death process $Y_t$ on $\mathbb{N}_0$ with linear death rates $d_j=dj$ for $j\in\mathbb{N}_0$ . Let $0<a\le 1$ , and define $\nu$ as follows:
where $\Gamma(\cdot)$ is the gamma function. Then $\nu$ is a QSD of $Y_t$ on $\mathbb{N}$ with $\partial=\{0\}$ . Hence, there is a family of QSDs fulfilling Theorem 2.
Formulae for stationary distributions and QSDs of BDPs follow directly from Theorems 1 and 2.
Corollary 1. ([Reference Anderson4, Reference Cavender10].) (i) Let $Y_t$ be a BDP on $\mathbb{N}_0$ with birth and death rates $b_j$ and $d_j$ , respectively, such that $b_{j-1}>0$ and $d_j>0$ for all $j\in\mathbb{N}$ . If $\pi$ is a stationary distribution for $Y_t$ , then
(ii) Let $Y_t$ be a BDP on $\mathbb{N}_0$ with birth and death rates $b_j$ and $d_j$ , respectively, such that $b_0=0$ , and $b_j>0$ and $d_j>0$ for all $j\in\mathbb{N}$ . Then a probability distribution $\nu$ on $\mathbb{N}$ is a QSD trapped into 0 for $Y_t$ if and only if
Proof. Here $\Omega=\{-1,1\}$ , $\omega_*=\omega_+=1$ , $\omega_-=-1$ , and $\mathcal{Y}=\mathbb{N}_0$ . Moreover, $q(j,j-1)=d_j$ and $q(j,j+1)=b_j$ for $j\in\mathbb{N}$ .
(i) $\partial=\varnothing$ . Since $\pi$ is a stationary distribution on $\mathbb{N}_0$ , it follows from Theorem 1 that
Hence the conclusion is obtained by induction.
(ii) $\partial=\{0\}$ and $\partial^{\textsf{c}}=\mathbb{N}$ . It follows from Theorem 2 that a probability measure $\nu$ is a QSD on $\mathbb{N}$ if and only if
that is, (9) holds.
Regarding the tail distributions, we have the following identities.
Corollary 2. Assume $\Omega$ is finite and $\partial=\varnothing$ . Let $\pi$ be a stationary distribution of $Y_t$ on $\mathcal{Y}$ . Then, for $x\in\mathbb{N}_0$ ,
where $A_j$ is defined in (3).
Proof. Assume without loss of generality that $\omega_*=1$ and $0\in\mathcal{Y}$ . The left-hand side of the equation in Theorem 3(3) is
while the right-hand side equals
which together yield the desired identity.
Corollary 3. Assume $\Omega$ is finite, $\partial\neq\varnothing$ , and let $\nu$ be a QSD of $Y_t$ on $\partial^{\textsf{c}}$ . Then for all $x\in\mathbb{N}_0\setminus\partial$ ,
where $A_j$ is defined in (3).
Proof. Similar to that of Corollary 2.
4. Asymptotic tails of limit distributions
To establish the asymptotic tails of limit distributions, we assume the following:
( $\textbf{A1}$ ) $\#\Omega<\infty$ .
( $\textbf{A2}$ ) $\mathcal{Y}$ is unbounded, and for $\omega\in\Omega$ , $q(x,x+\omega)=a_{\omega}x^{R_{\omega}^1}+b_{\omega}x^{R_{\omega}^2}+\mathrm{O}\big(x^{R_{\omega}^3}\big)$ is an APLH function in x on $\mathcal{Y}$ with $\big(R_{\omega}^1,R_{\omega}^2,R_{\omega}^3\big)$ for some constants $a_{\omega}, b_{\omega}$ . The APLH function is assumed strictly positive for all large $x\in\mathcal{Y}$ .
( $\textbf{A3}$ ) $\partial^{\textsf{c}}$ is irreducible.
Assumption ( $\textbf{A1}$ ) guarantees that the chain has bounded jumps only, which enables us to use the identities established in the previous section to estimate the tails. Assumption ( $\textbf{A2}$ ) is common in applications, and ensures that $\partial^{\textsf{c}}$ is unbounded too, as $\partial$ is finite by assumption. In particular, ( $\textbf{A2}$ ) is satisfied provided the following assumption holds:
$(\textbf{A2})^{\prime}$ For $\omega\in\Omega$ , $q(x,x+\omega)$ is a strictly positive polynomial for all large $x\in\mathcal{Y}$ .
Assumption ( $\textbf{A3}$ ) is assumed to avoid non-essential technicalities. Moreover, ( $\textbf{A3}$ ) means that either $Y_t$ is irreducible or the conditional process of $Y_t$ before entering $\partial$ is irreducible. This assumption is satisfied for many known one-dimensional infinite CTMCs modelling biological processes (e.g., for population processes). In addition, ( $\textbf{A3}$ ) implies that $\Omega_+\neq\varnothing$ and $\Omega_-\neq\varnothing$ (otherwise there are no non-singleton communicating classes).
The following parameters are well-defined and finite. Let
and define
where
(These values are only used to define $\sigma_1$ , $\sigma_2$ .)
Hence $R^1_-\ge R_{\omega}^2\ge R_-\,-1$ , for some $\omega\in\Omega_-$ with $R_{\omega}^1=R_-$ . Similarly, $R^1_+\ge R_+-1$ , $R^2_-\ge R_-\,-2$ , and $R^2_+\ge R_+-2$ . This implies that $0<\sigma_1\le1$ and $\sigma_1<\sigma_2\le2$ . If all transition rate functions are real analytic, then by convention, $R_{\omega}^2=R_{\omega}^1-1,\ R_{\omega}^3=R_{\omega}^1-2$ , for all $\omega\in\Omega$ , and hence $\sigma_1=1$ and $\sigma_2=2$ . Furthermore, let
The parameters also admit limit-free representations:
The form with the limit emphasizes that the parameters are coefficients of monomials of certain leading degrees. Furthermore, define
Note that $\alpha\le0$ implies $R_-\ge R_+$ . Moreover, we have $\alpha_+, \alpha_->0$ and $0<\delta\le\Delta$ , by ( $\textbf{A3}$ ). Furthermore, $\alpha$ , $\alpha_-$ , $\alpha_+$ , $\beta$ , and $\vartheta$ do not depend on the choice of second and third powers of the transition rate functions, whereas $\sigma_1$ , $\sigma_2$ , $\gamma$ , $\Delta$ , and $\delta$ do depend on the powers.
To state the results on the asymptotic tails of limit distributions, we classify probability distributions into the following classes.
Let $\mathcal{P}$ be the set of probability distributions on A. For $a,b>0$ , define
where $T_{\mu}(x)$ is the tail distribution of a probability measure $\mu$ , and $\mathrm{o}(\cdot)$ refers to the standard little-o notation. Furthermore, define
The sets $\mathcal{P}_{a}^{i+}$ , $i=1,2,3$ , are decreasing in a, while $\mathcal{P}_{a}^{i-}$ , $i=1,2,3$ , are increasing in a. Similarly, $\mathcal{P}_{a,b}^{2+}$ is decreasing in both a and b, while $\mathcal{P}_{a,b}^{2-}$ is increasing in both a and b. Moreover, it is readily verified that
The probability distributions in $\mathcal{P}^{2+}_{1}\cap\mathcal{P}^{2-}_{1}$ decay as fast as exponential distributions and are therefore exponential-like. Similarly, those in $\mathcal{P}^{1+}$ are super-exponential; those in $\mathcal{P}^{2-}_{<1}$ are sub-exponential, and in particular those in $\mathcal{P}^{3+}\cap\mathcal{P}^{3-}$ are power-like [Reference Johnson, Kemp and Kotz29].
The Conley–Maxwell–Poisson (CMP) distribution on $\mathbb{N}_0$ with parameter $(a,b)\in\mathbb{R}_+^2$ has probability mass function given by the following [Reference Johnson, Kemp and Kotz29]:
In particular, ${\textsf{CMP}}_{a,1}$ is a Poisson distribution. For every probability distribution $\mu\in\mathcal{P}^{1+}\cap\mathcal{P}^{1-}$ , there exist $(a_1,b_1), (a_2,b_2)\in\mathbb{R}_+^2$ such that
Conversely, every CMP distribution is an element in $\mathcal{P}^{1+}\cap\mathcal{P}^{1-}$ , and hence super-exponential. The zeta distribution on $\mathbb{N}$ with parameter $a>1$ has probability mass function given by the following [Reference Johnson, Kemp and Kotz29]:
where $\zeta(a)=\sum_{i=1}^{\infty}i^{-a}$ is the Riemann zeta function of a. For every probability distribution $\mu\in\mathcal{P}^{3+}\cap\mathcal{P}^{3-}$ , there exist $a_1, a_2>1$ such that
Conversely, every zeta distribution is an element in $\mathcal{P}^{3+}\cap\mathcal{P}^{3-}$ , and hence sub-exponential.
We first provide the tail asymptotics for QSDs. We point out that parameter conditions for the existence and ergodicity of QSDs are given in [Reference Xu, Hansen and Wiuf48].
Theorem 3. Assume ${(\textbf{A1})}$ – ${(\textbf{A3})}$ and $\partial\neq\varnothing$ . Assume there exists a QSD $\nu$ of $Y_t$ on $\partial^{\textsf{c}}$ . Then $\alpha\le0\le R$ . Furthermore, we have the following:
-
If $R=0$ , then $\alpha_-\ge\theta_{\nu}$ . If in addition $R_-=R_+$ , then $\alpha\le-\theta_{\nu}$ , and if $R_->R_+$ , then $\beta\ge\theta_{\nu}$ .
-
If $R_-=R_+>0$ and $\alpha=0$ , then $R\ge\sigma_1$ .
Moreover, if $R_->R_+$ , we have the following:
-
(i) If $R=0$ and $\beta>\theta_{\nu}$ , then $\nu\in\mathcal{P}^{1-}_{(R_-\,-R_+)\omega_*^{-1}}$ .
-
(ii) If $R=0$ , $\beta=\theta_{\nu}$ , and $R_-\,-R_+\le1$ , then $\nu\in\mathcal{P}^{2-}_{1}$ .
-
(iii) If $R=0$ , $\beta=\theta_{\nu}$ , and $R_-\,-R_+>1$ , then $\nu\in\mathcal{P}^{2-}_{R_-\,-R_+-1}$ .
-
(iv) If $R>0$ , then $\nu\in\mathcal{P}^{1-}_{(R_-\,-R_+)\omega_*^{-1}}$ . If in addition $R>1$ , then $\nu\in\mathcal{P}^{1+}_{(R_-\,-R_+)(\omega_+\omega_*)^{-1}}$ .
If $R_+=R_-$ , we have the following:
-
(v) If $R>0$ and $\alpha<0$ , then $\nu\in\mathcal{P}^{2-}_{1}$ . If in addition $R>1$ , then $\nu\in\mathcal{P}^{2+}_{1}$ .
-
(vi) If $R>0$ and $\alpha=0$ , we have the following:
-
– If $R=\sigma_1<1$ , then $\nu\in\mathcal{P}^{2-}_{1-R}$ .
-
– If $R\ge\sigma_1=1$ , then $\nu\in\mathcal{P}^{3-}$ .
-
– If $\min\{1,R\}>\sigma_1$ , then $\nu\in\mathcal{P}^{2-}_{1<}$ ,
-
-
(vii) If $R=0$ , we have the following:
-
– If $\alpha+\theta_{\nu}=0$ and $\sigma_1<1$ , then $\nu\in\mathcal{P}^{2-}_{1-\sigma_1}$ .
-
– If $\alpha+\theta_{\nu}=0$ and $\sigma_1=1$ , then $\nu\in\mathcal{P}^{3-}$ ,
-
– If $\alpha+\theta_{\nu}<0$ , then $\nu\in\mathcal{P}^{2-}_{1}$ .
-
Furthermore, we have the following:
-
(viii) If $R=1$ , then $\nu\in\mathcal{P}^{3+}_{\theta_{\nu}\alpha_+^{-1}}$ .
-
(ix) If $0<R<1$ , then $\nu\in\mathcal{P}^{2+}_{1-R}$ .
-
(x) If $R=0$ and $\alpha_->\theta_{\nu}$ , then $\nu\in\mathcal{P}^{2+}_{1}$ .
-
(xi) If $R=0$ and $\alpha_-=\theta_{\nu}$ , then $\nu\in\mathcal{P}^{1+}_{-\sigma_1\omega_-^{-1}}$ .
Corollary 4. Assume ${(\textbf{A1})}$ – ${(\textbf{A3})}$ and $\partial\neq\varnothing$ . No QSDs have a tail decaying faster than a CMP distribution. Moreover, any QSD, if such exists, is super-exponential if $R_->\max\{1,R_+\}$ or (xi) holds, exponential-like if (a) $R_-=R_+=0$ and $\alpha+\theta_{\nu}<0$ or (b) $R_-=R_+>1$ and $\alpha<0$ holds, and sub-exponential if (vi) or $R_-=R_+=\alpha+\theta_{\nu}=0$ holds; in particular, it decays no faster than a power-like distribution if $R\ge \sigma_1=1$ and $\alpha=0$ .
Analogously, we further characterize the tails of the stationary distributions. Parameter conditions for the existence and ergodicity of stationary distributions are given in [Reference Xu, Hansen and Wiuf48].
Theorem 4. Assume ${(\textbf{A1})}$ – ${(\textbf{A2})}$ and $\partial=\varnothing$ . Assume there exists a stationary distribution $\pi$ of $Y_t$ on $\mathcal{Y}$ with unbounded support. Then $\alpha\le0$ , and in particular, when $\alpha=0$ ,
-
if $\Delta=0$ then $\sigma_1<1$ ;
-
if $\sigma_1=1$ then $\Delta>1$ .
Moreover, we have the following:
-
(i) If $R_->R_+$ , then $\pi\in\mathcal{P}^{1+}_{(R_-\,-R_+)(\omega_+\omega_*)^{-1}}\cap\mathcal{P}^{1-}_{(R_-\,-R_+)\omega_*^{-1}}$ .
-
(ii) If $R_-=R_+$ and $\alpha<0$ , then $\pi\in\mathcal{P}^{2+}_{1}\cap\mathcal{P}^{2-}_{1}$ .
-
(iii) If $\alpha=0$ , $\Delta>0$ , and $\sigma_1<1$ , then $\pi\in\mathcal{P}^{2+}_{1-\sigma_1}\cap\mathcal{P}^{2-}_{1-\sigma_1}$ .
-
(iv) If $\alpha=0$ and $\sigma_1=1$ , then $\pi\in\mathcal{P}^{3+}_{\Delta-1}$ . In particular, if in addition (iv) $^\prime$ $\delta>1$ , then $\pi\in\mathcal{P}^{3+}_{\Delta-1}\cap\mathcal{P}^{3-}_{\delta-1}$ .
-
(v) If $\alpha=0$ , $\Delta=0$ , and $\sigma_2<1$ , then $\pi\in\mathcal{P}^{2-}_{1-\sigma_2}$ .
-
(vi) If $\alpha=0$ , $\Delta=0$ , and $\sigma_2\ge1$ , then $\pi\in\mathcal{P}^{3-}$ .
As a consequence, a trichotomy regarding the tails of the stationary distributions can be derived.
Corollary 5. Assume ${(\textbf{A1})}$ – ${(\textbf{A2})}$ and $\partial=\varnothing$ . Any stationary distribution of $Y_t$ on $\mathcal{Y}$ with unbounded support, if such exists, is
Corollary 6. Assume ${(\textbf{A1})}$ , ${(\textbf{A2})^{\prime}}$ , ${(\textbf{A3})}$ , $\partial=\varnothing$ , $R\ge 3$ , and $(R-1)\vartheta-\alpha_+\le0$ . Any stationary distribution of $Y_t$ on $\mathcal{Y}$ with unbounded support, if such exists, is ergodic. In particular, $Y_t$ is non-explosive.
Proof. By ${(\textbf{A2})^{\prime}}$ , $R\in\mathbb{N}_0$ and $\sigma_1=1$ . By [Reference Xu, Hansen and Wiuf48, Theorem 1], under ${(\textbf{A3})}$ , $Y_t$ is explosive if either (1) $R\ge2$ and $\alpha>0$ , or (2) $R\ge3$ , $\alpha=0$ , and $\gamma-\vartheta>0$ . By Theorem 4, $\alpha\le0$ . If a stationary distribution exists and if $Y_t$ is non-explosive, then under ${(\textbf{A3})}$ , $Y_t$ is is positive recurrent and the stationary distribution is unique and ergodic [Reference Norris39]. When $\alpha=0$ and $R\ge3$ , it follows from Theorem 4 that $\Delta>1$ , i.e., $\gamma-\vartheta<(R-1)\vartheta-\alpha_+\le0$ as assumed. Hence, $Y_t$ is always non-explosive and thus the stationary distribution is ergodic.
We make the following remarks:
-
The estimate of the tail does not depend on the choice ( $R_{\omega}^2,\ R_{\omega}^3$ ) of the transition rate functions when $\alpha<0$ , whereas it may depend on this choice when $\alpha=0$ . In this case, the larger $\sigma_1$ and $\sigma_2$ are, the sharper the results are.
-
Generically, no limit distributions (in the cases covered) can decay faster than a CMP distribution or slower than a zeta distribution.
-
The unique gap case in Theorem 4 is $\alpha=0$ , $\sigma_1<1$ , and $\Delta=0$ .
-
Assume $(\textbf{A2})^{\prime}$ . By Corollary 6, if the chain $Y_t$ is explosive and a stationary distribution exists, then $\alpha=0$ and $R\ge3$ .
-
Although not stated explicitly in Theorem 4, the tail asymptotics of a stationary distribution of a BDP (Cases (i)–(iii) and (iv) $^\prime$ ) is sharp up to the leading order, in comparison with Proposition 3. Similarly, when $R_->R_+$ , the tail asymptotics is sharp up to the leading order for upwardly skip-free processes (Case (i)) [Reference Anderson4]. In comparison with the sharp results provided in Proposition 4, the results obtained in Theorem 4 capture the leading tail asymptotics of a stationary distribution, e.g., in Case (iii).
-
The assumption that $R>1$ in Corollary 4 is crucial. Indeed, as Examples 5 and 6 illustrate, when $R=1$ and $\alpha<0$ , the QSD may still exist and has either geometric or zeta-like tail. This means $\alpha<0$ is not sufficient for any QSD to have an exponential tail. It remains to see whether a QSD with CMP tail may exist when $R=1$ . Moreover, we emphasize that the conditions $R>1$ and $\alpha<0$ ensure the existence of a unique ergodic QSD assuming $(\textbf{A2})^{\prime}$ [Reference Xu, Hansen and Wiuf47]. Hence, such ergodic QSDs are not heavy-tailed.
-
Let $Y_t$ be a CTMC and $\widetilde{Y}_t$ a BDP on the same state space. If the critical parameters are the same ( $\alpha$ , $\beta$ , R, etc.), then the two processes share the same qualitative characterization in terms of existence of limit distributions, ergodicity, explosivity, etc. [Reference Xu, Hansen and Wiuf47]. One might conjecture that the decay of limit distributions (provided such exist) is also the same and takes the sharp form induced by the BDP. In other words, one might conjecture that the tail asymptotics of a general CTMC coincide with that of a BDP representative.
In the following, we illustrate and elaborate on the results by example. The examples have real analytic APLH transition rate functions, and thus $\sigma_1=1$ and $\sigma_2=2$ .
Assumption ${(\textbf{A1})}$ is crucial for Theorem 4.
Example 3. Consider a model of stochastic gene regulatory expression [Reference Shahrezaei and Swain43], given by $\Omega=\{-1\}\cup\mathbb{N}$ and
where $a>0$ , and $b_j\ge0$ for all $j\in\mathbb{N}_0$ . Here the backward jump $-1$ represents the degradation of mRNA with unity degradation rate, and the forward jumps $j\in\mathbb{N}$ account for bursty production of mRNA with transcription rate a and burst size distribution $(b_j)_{j\in\mathbb{N}_0}$ . When $b_j=(1-\delta)\delta^j$ for $0<\delta<1$ , the stationary distribution is the negative binomial distribution [Reference Shahrezaei and Swain43]:
When $a=1$ , $\pi$ is also geometric. Theorem 4, if it did apply, would seem to suggest $T_{\pi}$ decays like a CMP distribution, since $1=R_->R_+=0$ . The technical reason behind this is that the proof applies Corollary 2, which requires ${(\textbf{A1})}$ . It does not seem possible to directly extend the result in Theorem 4 to CTMCs with unbounded jumps.
Example 4. Consider a BDP with birth and death rates
where $a>0$ , $b\in\mathbb{N}$ , and the S(i,j) denote the Stirling numbers of the second kind [Reference Abramowitz and Stegun1]. Here $R_+=0<R_-=b$ , and $\alpha=-S(b,b)=-1<0$ . This BDP has an ergodic stationary distribution on $\mathbb{N}_0$ [Reference Xu, Hansen and Wiuf47], and the unique stationary distribution is $\pi={\textsf{CMP}}_{a,b}$ .
By Theorem 3(v), the tail of a QSD decays no faster than an exponential distribution when $\alpha<0$ and $0\le R_-=R_+\le1$ , which is also confirmed by the examples below.
Example 5. Consider the linear BDP on $\mathbb{N}_0$ with birth and death rates
where b and d are positive constants [Reference O’Neill40]. For this process, a QSD $\nu$ is
Hence $T_{\nu}$ decays as fast as the zeta distribution with parameter 2. Here $\alpha=-d/2<0$ , and $R=R_-=R_+=1$ .
Example 6. Consider the linear BDP on $\mathbb{N}_0$ with $b_j=b j$ and $d_j=d_1 j$ with $0<b<d_1$ and $j\in\mathbb{N}$ [Reference Van Doorn45]. For this process, a QSD $\nu$ is
a geometric distribution. Here $\alpha=b-d_1<0$ , and $R=R_-=R_+=1$ . By Parts (iv) and (viii) of Theorem 3, the tail of the QSD decays no faster than a CMP distribution and no more slowly than a zeta distribution, if $R_+<R_-=1$ , which is also confirmed by the example below.
Example 7. Consider a BDP on $\mathbb{N}_0$ given by
Here $R=R_-=1>R_+=0$ , $\alpha=-1$ , and $\sigma_1=1$ , $\sigma_2=2$ . Using the same Lyapunov function constructed in the proof of [Reference Xu, Hansen and Wiuf47, Theorem 4.4], it can be shown that there exists a uniquely ergodic QSD. By Theorem 3(iv), the tail of the QSD decays no more slowly than a CMP distribution. Indeed, the QSD is given by
The tail of the QSD decays like a Poisson distribution.
Example 8. Consider a BDP on $\mathbb{N}_0$ given by
Here $R=R_-=R_+=2>1$ , and $\alpha=0$ . Corollary 4 states that any QSD (if it exists) is heavy-tailed and its tail decays no faster than a zeta-like distribution. Indeed, a QSD of the process is given by
Example 9. Consider a quadratic BDP on $\mathbb{N}_0$ , given by
Then $R_-=R_+=R=2>1$ , $\alpha=-1/2$ . Hence there exists a uniquely ergodic QSD [Reference Xu, Hansen and Wiuf47]. By Corollary 4, this QSD decays exponentially. Indeed, the QSD is given by
5. Applications
In this section, we apply the results on asymptotic tails to diverse models in biology. We emphasize that for all models/applications, the transition rate functions are real analytic APLH on a subset of $\mathbb{N}_0$ , and thus $\sigma_1=1$ and $\sigma_2=2$ , which we will not further mention explicitly.
5.1. Biochemical reaction networks
In this section, we apply the results of Section 4 to some examples of stochastic reaction networks (SRNs) with mass-action kinetics. These are used to describe interactions of constituent molecular species with many applications in systems biology, biochemistry, genetics, and beyond [Reference Gardiner21, Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani41]. An SRN with mass-action kinetics is a CTMC on $\mathbb{N}_0^d$ ( $d\ge 1$ ) encoded by a labelled directed graph [Reference Anderson, Kurtz, Koeppl and Setti3]. We concentrate on SRNs on $\mathbb{N}_0$ with one species (S). In this case the graph is composed of reactions (edges) of the form $n\text{S}\longrightarrow^{\kappa} m\text{S}$ , $n,m\in\mathbb{N}_0$ (n molecules of species S is converted into m molecules of the same species), encoding a jump from x to $x+m-n$ with propensity $q(x,x+m-n)=\kappa x(x-1)\ldots(x-n+1)$ , $\kappa>0$ . Note that multiple reactions might result in the same jump vector.
In general little is known about the stationary distributions of a reaction network, let alone the QSDs, provided either such exist [Reference Anderson, Craciun and Kurtz2, Reference Gupta, Briat and Khammash25, Reference Hansen and Wiuf27, Reference Hoessly and Mazza28]. Special cases include complex balanced networks (in arbitrary dimension) which have Poisson product-form distributions [Reference Anderson, Craciun and Kurtz2, Reference Cappelletti and Wiuf9], reaction networks that are also BDPs, and reaction networks with irreducible components, each with a finite number of states.
Example 10. To show how general the results are, we consider two SRNs, neither of which is a BDP.
(i) Consider a reaction network with a strongly connected graph [Reference Xu, Hansen and Wiuf47]:
For this reaction network, $\Omega=\{1,-2\}$ , and
Hence, $\alpha=-\kappa_3<0$ . It is known that there exists a unique exponentially ergodic stationary distribution $\pi$ on $\mathbb{N}$ [Reference Xu, Hansen and Wiuf47]. (The state 0 is neutral [Reference Xu, Hansen and Wiuf48].) By Theorem 4, $\pi\in\mathcal{P}^{1+}_1\cap\mathcal{P}^{1-}_1$ . Hence $\pi$ is light-tailed and $T_{\pi}$ decays as fast as a Poisson distribution, since $\omega_+=\omega_*=1$ , $R_+=2$ , and $R_-=3$ . However, the stationary distribution is generally not Poisson. If $\kappa_2^2=\kappa_1\kappa_3$ , then the reaction network is complex-balanced; hence the stationary distribution is Poisson [Reference Anderson, Craciun and Kurtz2]. If the parameter identity is not fulfilled, then the distribution cannot be Poisson in this case [Reference Cappelletti and Wiuf9].
(ii) Consider a similar reaction network including direct degradation of S [Reference Xu, Hansen and Wiuf47]:
The threshold parameters are the same as in (i), and it follows from [Reference Xu, Hansen and Wiuf47] that the reaction network has a uniformly exponentially ergodic QSD $\nu$ . By Theorem 1, $T_{\nu}$ decays like a CMP distribution.
Example 11. The following bursty Schlögl model was proposed in [Reference Falk, Mendler and Drossel20]:
where $j\in\mathbb{N}$ . When $j=1$ , it reduces to the classical Schlögl model.
The associated process has a unique ergodic stationary distribution $\pi$ on $\mathbb{N}_0$ [Reference Xu, Hansen and Wiuf47]. Bifurcation with respect to patterns of the ergodic stationary distribution is discussed in [Reference Falk, Mendler and Drossel20], based on a diffusion approximation in terms of the Fokker–Planck equation. Using the results established in this paper, the tail distribution can be characterized rigorously. In fact $\pi\in\mathcal{P}^{1+}_{j^{-1}}\cap\mathcal{P}^{1-}_1$ . Hence $\pi$ is light-tailed and $T_{\pi}$ decays like a CMP distribution.
Proof. We have $\Omega=\{-1,1\}\cup\{j\}$ . The ergodicity follows from [Reference Xu, Hansen and Wiuf47, Subsection 4.3] as a special case. It is straightforward to verify that ${(\textbf{A1})}$ – ${(\textbf{A3})}$ are all satisfied. Moreover, $R_+=2$ and $R_-=3$ . The conclusion then follows from Theorem 4.
Example 12. Consider the following one-species SRN:
In this case, $\omega_*=3$ , and the ambient space $\mathbb{N}_0$ is divided into disjoint irreducible sets $3\mathbb{N}$ , $3\mathbb{N}_0+1$ , and $3\mathbb{N}_0+2$ , as well as an absorbing state 0. By simple calculation, we have $R_-=3>R_+=1$ , and $\alpha = -3\kappa_2<0$ . Hence, this SRN is positive recurrent on the sets $3\mathbb{N}_0+1$ and $3\mathbb{N}_0+2$ , while it admits a unique QSD on $3\mathbb{N}$ [Reference Xu, Hansen and Wiuf48, Theorem 7]. According to Theorem 4, the tails of the stationary distributions decay like CMP on $3\mathbb{N}_0+1$ and $3\mathbb{N}_0+2$ , and according to Theorem 3, the tail of the QSD also decays as CMP on $3\mathbb{N}$ . Since the transition rate functions take a common form (polynomial) on all three irreducible sets, the parameters are the same on all three irreducible sets (open and closed), and consequently, the tail asymptotics are the same on the two closed irreducible sets. The same would hold true for the open irreducible sets, if there were more than one open irreducible set.
Example 13. Consider the following one-species S-system modelling a gene regulatory network [Reference Chowdhury, Chetty and Evans14]:
with the generalized mass action kinetics (GMAK)
where $\kappa_1, \kappa_2>0$ are the reaction rate constants, and $\xi_2>\xi_1>0$ are the indices of GMAK.
By Stirling’s formula,
hence $q(x,x+1)$ is APLH with $(\xi_1,\xi_1-1,\xi_1-2)$ and $q(x,x-2)$ is APLH with $(\xi_2,\xi_2-1,\xi_2-2)$ . Then $R_-=\xi_2>R_+=\xi_1$ , $\omega_*=1$ , $\omega_-=-2$ , and $\omega_+=1$ . Using the same Lyapunov function constructed in the proof of [Reference Xu, Hansen and Wiuf47, Theorem 4.4], it can be shown that there exists a uniquely ergodic stationary distribution $\pi$ on $\mathbb{N}_0$ with support $\mathbb{N}$ . By Theorem 4, $\pi\in\mathcal{P}^{1+}_{\xi_3-\xi_1}\cap\mathcal{P}^{1-}_{\xi_3-\xi_1}$ .
5.2. An extended class of branching processes
Consider an extended class of branching processes on $\mathbb{N}_0$ [Reference Chen13] with transition rate matrix $Q=(q(x,y))_{x,y\in\mathbb{N}_0}$ given by
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:
( $\textbf{H1}$ ) $\mu(0)>0$ , $\mu(0)+\mu(1)<1$ .
( $\textbf{H2}$ ) $\sum_{y\in\mathbb{N}}q(0,y)y<\infty$ , $M=\sum_{k\in\mathbb{N}_0}k\mu(k)<\infty$ .
( $\textbf{H3}$ ) r(x) is a polynomial of degree $R\ge1$ for large x.
The tail asymptotics of infinite stationary measures in the null-recurrent case is investigated in [Reference Liu and Zhao31] under ( $\textbf{H1}$ )–( $\textbf{H2}$ ) for general r. Here we assume r is polynomial ( $\textbf{H3}$ ). The following is a consequence of the results of Section 5.
Theorem 5. Assume ${(\textbf{H1})}$ – ${(\textbf{H3})}$ , $Y_0\neq0$ , and that $\mu$ has finite support.
-
(i) Assume $q_0>0$ . Then there exists an ergodic stationary distribution $\pi$ on $\mathbb{N}_0$ if (i-1) $M<1$ or (i-2) $M=1$ and $R>1$ . Moreover, $T_{\pi}$ decays like a geometric distribution if (i-1) holds, and like a zeta distribution if (i-2) holds.
-
(ii) Assume $q_0=0$ . Then there exists an ergodic QSD $\nu$ on $\mathbb{N}$ if (ii-1) $M<1$ and $R>1$ or (ii-2) $M=1$ and $R>2$ . Moreover, $T_{\nu}$ decays like a geometric distribution if (ii-1) holds, while it decays no faster than a zeta distribution if (ii-2) holds.
Proof. For all $k\in\Omega$ , let
By ${(\textbf{H1})}$ , $\mu(k)>0$ for some $k\in\mathbb{N}$ . Hence, regardless of q(0), by positivity of r, ${\textrm(\textbf{A1})}$ – ${\textrm(\textbf{A3})}$ are satisfied with $\Omega_{-}=\{-1\}$ and $\Omega_{+}=\{j\in\mathbb{N}: j+1\in\text{supp}\,\mu \text{ or } q(0,j)>0\}$ . Let $r(x)=ax^R+bx^{R-1}+\mathrm{O}\big(x^{R-2}\big)$ with $a>0$ . It is straightforward to verify that $R_{+}=R_{-}=R$ , $\alpha=a(M-1)$ . The ergodicity follows from [Reference Xu, Hansen and Wiuf47], and the tail asymptotics follow from Theorems 4 and 3.
5.3. Stochastic population processes under bursty reproduction
Two stochastic population models with bursty reproduction are investigated in [Reference Be’er and Assaf7].
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 Assaf7], and $\mu$ is the burst size distribution.
Approximations of the mean time to extinction and QSD are discussed in [Reference Be’er and Assaf7] against various burst size distributions of finite mean (e.g., a Dirac measure, a Poisson distribution, a geometric distribution, a negative-binomial distribution). The existence of an ergodic QSD for this population model is established in [Reference Xu, Hansen and Wiuf47]. However, the tail of the QSD is not addressed therein.
Theorem 6. Assume $\mu$ has finite support. Let $\nu$ be the unique ergodic QSD on $\mathbb{N}$ trapped to zero for the Verhulst logistic model $Y_t$ . Then $T_{\nu}$ decays like a CMP distribution.
Proof. We have $\Omega=\{-1\}\cup\text{supp}\,\mu$ , $q(x,x-1)=\frac{c}{K}x^2+x$ , $q(x,x+k)=c\mu(k)x$ , for $k\in\text{supp}\,\mu$ and $x\in\mathbb{N}$ . Since $\mu$ has finite support, ${(\textbf{A1})}$ – ${(\textbf{A3})}$ are satisfied. Moreover, since $\text{supp}\,\mu\neq\varnothing$ , we have $R_-=2$ and $R_+=1$ . Again, the ergodicity result follows from [Reference Xu, Hansen and Wiuf47]. The tail asymptotics follow directly from Theorem 4.
In subsequent sections, we provide proofs of the main results in Section 4. Since the proof of Theorem 3 is based on that of Theorem 4, we begin by proving Theorem 4 in the next section.
6. Proof of Theorem 4
Let
Note that $\beta=\alpha_0$ . From Lemma 1, $\alpha_-=\omega_*\sum_{j=\omega_-+1}^0\alpha_{j}$ and $\alpha_+=\omega_*\sum_{j=1}^{\omega_+}\alpha_j$ . By ${(\textbf{A3})}$ ,
Since
we have
Since ( $\textbf{A1}$ )–( $\textbf{A2}$ ) imply that ${\cap}_{\omega\in\Omega}\ \{x\in\mathcal{Y}\,:\, q(x,x+\omega)=0\}$ is finite, it easily follows that both $\Omega_-\neq\varnothing$ and $\Omega_+\neq\varnothing$ since $\text{supp}\,\pi$ is unbounded. Hence $\alpha_-\ge\alpha_0>0$ , $\alpha_+\ge\alpha_1>0$ , and $-\infty<\omega_-<\omega_+<\infty$ .
For ease of exposition and without loss of generality, we assume throughout the proof that $\omega_*=1$ (recall the argument in the proof of Theorem 1). Hence $\mathbb{N}_0+b\subseteq\mathcal{Y}\subseteq\mathbb{N}_0$ for some $b\in\mathbb{N}_0$ by Proposition 2.
Most of the inequalities below are based on the identities in Theorem 1 and Corollary 2. Therefore, we use ‘LHS’ (‘RHS’) with a label in the subscript as shorthand for the left-hand side (right-hand side) of the equation with the given label.
The claims that $\Delta=0$ implies $\sigma_1<1$ and that $\sigma_1=1$ implies $\Delta>1$ are proved in Step I below.
We first show $\alpha\le0$ . Suppose for the sake of contradiction that $\alpha>0$ . Then either (1) $R_+>R_-$ or (2) $R_+=R_-$ and $\alpha_+>\alpha_-$ holds. Define the auxiliary function
From (11) it follows that
Let
Then there exist $N_3,\ N_4\in\mathbb{N}$ with $N_3>N_1, N_4$ such that for all $x\ge N_3$ ,