Bifurcation analysis of a reaction-diffusion-advection predator-prey system with delay

1.   Introduction

In this paper we prove the existence of a stationary distribution and the non-explosivity of stochastically modeled reaction networks with a specific form of non-mass action kinetics (see (3.2)), and generalize two results related to reaction networks with mass action kinetics to the non-mass action setting.

In ^[1], it was shown that if a deterministically modeled reaction network (with deterministic mass action kinetics) is complex balanced with complex balanced equilibrium $c\in \mathbb R_{\ge 0}^m$, then the associated stochastic model (with stochastic mass action kinetics and the same choice of rate constants) admits the stationary distribution

It was also shown in ^[1] that if the intensity functions of the stochastic model are not of mass action, but of the form

where $\theta_i: \mathbb{Z}\to\mathbb{R}_{\geq 0}$ and $\theta_i(0) = 0$, then the stochastic model admits the stationary measure

where $c$ is the complex balanced equilibrium of the associated deterministic mass action model with rate constants $\{\kappa_k\}$. Related work characterizing stationary distributions for stochastic models with non-mass action kinetics can be found in ^[4]. Later, we will show that $\pi$ in (1.2) is summable under some mild growth conditions on $\theta_i$, thus it can be normalized to provide a stationary distribution.

The main result in ^[1] related to mass action models (i.e. the characterization of the stationary distribution for complex balanced models given in (1.1)) has been the starting point for a number of research endeavors. In particular,

1. The main result of ^[5] gives a condition for a stochastically modeled reaction network to be non-explosive. Applying this result for complex balanced networks, whose stationary distributions are given by (1.1), it can be concluded that for any initial distribution, complex balanced stochastic mass-action systems are non-explosive.

2. The converse of the main result of ^[1] related to mass action models was proven in ^[3]. That is, it was shown that if (1.1) is the stationary distribution of a stochastically modeled reaction network with mass action kinetics, then the associated deterministic model with mass action kinetics is complex balanced.

3. In ^[2], it was shown that a particular scaling limit of the stationary distribution (1.1) is a well known Lyapunov function. This result has been useful in the study of large deviations related to stochastic models of chemical systems ^[6,7].

In this paper we generalize the three results above to the case of stochastically modeled reaction networks with non-mass action kinetics.

The outline of the remainder of the paper is as follows. In Section 2, we review the relevant definitions and mathematical models for chemical reaction networks. In Section 3, we state the previous results on reaction networks that were briefly described above. Finally, in Section 4 we provide our main results, which provide the existence of a stationary distribution and the non-explosivity of stochastically modeled reaction networks with non-mass action kinetics and generalize the previous findings on mass action systems to the non-mass action settings.

2.   Chemical reaction networks

2.1.   Reaction networks and key definitions

We consider a set of $m$ species $\{S_1, S_2, .., S_m\}$ undergoing a finite number of reaction types enumerated by $k \in \{1, \dots, K\}$. For the $k^{th}$ reaction type we denote by $y_k$ and $y_k' \in \mathbb{Z}^m_{\geq 0}$ the vectors representing the number of molecules of each species consumed and created in one instance of the reaction. For example, the reaction $S_1+S_2\to 2S_2$ has $y_k = (1, 1)$ and $y_k' = (0, 2)$, if the system has 2 species. The vectors $y_k$ and $y_k'$ are called complexes of the system. $y_k$ is called the source complex and $y_k'$ is called the product complex. A complex can be both a source complex and a product complex.

Definition 2.1. Let $\mathcal{S} = \{S_1, ..., S_m\}$, $\mathcal{C} = \cup_k\{y_k, y_k'\}$, and $\mathcal{R} = \cup_k\{y_k\to y_k'\}$ be the sets of species, complexes and reactions respectively. The triple $\{\mathcal{S, C, R}\}$ is called a reaction network.

To each reaction network $\{\mathcal{S, C, R}\}$, there is a unique directed graph constructed as follows. The nodes of the graph are the complexes. A directed edge is placed from $y_k$ to $y_k'$ if and only if there is a reaction $y_k \to y_k'$. Each connected component is called a linkage class of the graph. We denote by $\ell$ the number of linkage class.

Definition 2.2. A reaction network $\{\mathcal{S, C, R}\}$ is called weakly reversible if the associated directed graph is strongly connected.

Definition 2.3. The linear subspace $S = \text{span}_k\{y_k'-y_k\}$ generated by all reaction vectors is called the stoichiometric subspace of the network. For $c\in \mathbb{R}^m_{\geq 0}$ we say $c+S = \{x\in \mathbb{R}^m |x = c+s$ for some $s\in S\}$ is a stoichiometric compatibility class, $(c+S)\cap \mathbb{R}^m_{\geq 0}$ is a non-negative stoichiometric compatibility class, and $(c+S)\cap \mathbb{R}^m_{ > 0}$ is a positive stoichiometric compatibility class. Denote $dim(S) = s$.

Finally, we provide the definition of the deficiency of a reaction network ^[15].

Definition 2.4. The deficiency of a chemical reaction network $\{\mathcal{S, C, R}\}$ is $\delta = |\mathcal{C}|-\ell-s$, where $|\mathcal{C}|$ is the number of complexes, $\ell$ is the number of linkage classes, and $s$ is the dimension of the stoichiometric subspace of the network.

2.2.   Dynamical system models

2.2.1.   Stochastic model

The most common stochastic model for a reaction network $\{\mathcal{S, C, R}\}$ treats the system as a continuous time Markov chain whose state at time $t$, $X(t)\in \mathbb Z^m_{\ge 0}$, is a vector giving the number of molecules of each species present with each reaction modeled as a possible transition for the chain. The model for the $k^{th}$ reaction is determined by the source and product complexes of the reaction, and a function $\lambda_k$ of the state that gives the transition intensity, or rate, at which the reaction occurs. In the biological and chemical literature, transition intensities are referred to as propensities.

Given that the $k^{th}$ reaction happens at time $t$, the state is updated by the addition of the reaction vector $y_k'-y_k$,

A common choice for the intensity functions $\lambda_k$ is to assume the system satisfies the stochastic version of mass action kinetics. In this case, the functions have the form

where $\kappa_k > 0$ is called the rate constant. Under the assumption of mass action kinetics and a non-negative initial condition, it follows that the dynamics of the system is confined to the particular non-negative stoichiometric compatibility class determined by the initial value $X(0)$, namely $X(t)\in (X(0)+S)\cap \mathbb{R}^m_{\geq 0}$.

Simple book-keeping implies that $X(t)$ satisfies

where $R_k(t)$ gives the number of times reaction $k$ has occurred by time $t$. Kurtz showed that $X$ can be represented as the solution to the stochastic equation

where the $Y_k$ are independent unit-rate Poisson process ^[8].

Another way to characterize the models of interest is via Kolmogorov's forward equation, termed the chemical master equation in the biology and chemistry literature, which describes how the distribution of the process changes in time. Letting $p_\mu(x, t)$ give the probability that $X(t) = x$ assuming an initial distribution of $\mu$, the forward equation is

Constant solutions to the forward equation, i.e. those satisfying

are stationary measures for the process, and if they are summable they can be normalized to give a stationary distribution. Assuming the associated stochastic model is non-explosive, stationary distributions characterize the long-time behavior of the stochastically modeled system.

Finding stationary distributions is in general a difficult task. However, as discussed in the Introduction and in later sections, an explicit form for $\pi(x)$ can be found when the associated deterministic mass action model has a complex balance equilibrium.

2.2.2.   Deterministic model

Under an appropriate scaling limit (the classical scaling detailed in Section 3.2.1) the continuous time Markov chain model of the previous section becomes

where

where $\kappa_k > 0$ is the rate constant, and where for two vectors $u, v\in \mathbb{R}^m_{\geq 0}$ we denote $u^v = \prod_i u_i^{v_i}$ with the convention $0^0 = 1$. Later, we will also utilize the notation $uv$ for the vector whose $i$th component is $u_iv_i$.

We say that the deterministic system (2.3) has deterministic mass action kinetics if the rate functions $f_k$ have the form (2.4). The system 2.3 is equivalent to the system of ODEs

The trajectory with initial condition $x_0$ is confined to the non-negative stoichiometric compatibility class $(x_0+S)\cap \mathbb{R}^m_{\geq 0}$.

Some mass action systems have complex balanced equilibria ^[25,13], which has been shown to play an important role in many biological mechanisms ^[9,10,11,12]. An equilibrium point $c$ is said to be complex balanced if for all $z\in \mathcal{C}$, we have

where the sum on the left is over reactions for which $z$ is the product complex and the sum on the right is over reactions for which $z$ is the source complex.

In ^[13] it was shown that if there exists a complex balanced equilibrium $c\in \mathbb{R}^m_{\geq 0}$ then

1. There is one, and only one, positive equilibrium point in each positive stoichiometric compatibility class.

2. Each such equilibrium point is complex balanced.

3. Each such complex balanced equilibrium point is locally asymptotically stable relative to its stoichiometric compatibility class.

In ^[24], a proof is presented showing global stability relative to the stoichiometric compatibility class. Because of the above, we say that a system is complex balanced if it admits a complex balanced equilibrium.

The condition of a reaction network being complex balanced can be checked more conveniently by using a classical result ^[14]. See also ^[15,16].

Theorem 2.1. If the reaction system is weakly reversible and has a deficiency of zero, then for any choice of rate constants $\{\kappa_k\}$ the deterministically modeled system with mass action kinetics is complex balanced.

3.   Previous results

3.1.   Product-form stationary distribution

As mentioned in the previous section, complex balanced systems deserve special attention, and thus characterizations of complex balanced systems are of interest. The following Theorem in ^[1] provides an explicit form for the stationary distribution of complex balanced systems.

Theorem 3.1. Let $\{\mathcal{S, R, C}\}$ be a reaction network. Suppose that when modeled deterministically with mass action kinetics and rate constants $\{\kappa_k\}$ the system is complex balanced with a complex balanced equilibrium $c\in \mathbb{R}^m_{\geq 0}$. Then the stochastically modeled system with intensities (2.1), with the same rate constants $\{\kappa_k\}$, admits the stationary distribution

See also ^[5], which shows that these systems are non-explosive, implying $\pi$ yields the limiting distributions of the process.

In the case when the stochastic model does not have mass action kinetics, ^[1] also provides an extended result. In particular, ^[1] considers generalized intensity functions as mentioned in several past papers ^[17,18],

where $\kappa_k$ are positive rate constants, $\theta_i: \mathbb{Z} \to \mathbb{R}_{\geq 0}$, and $\theta_i(x) = 0$ if $x\leq 0$. The functions $\theta_i$ should be thought of as the "rate of association" of the $i^{th}$ species ^[17]. For a system with intensity functions (3.2), the product form stationary distribution is quite similar to the one in Theorem 3.1.

Theorem 3.2. Let $\{\mathcal{S, R, C}\}$ be a reaction network. Suppose that when modeled deterministically with mass action kinetics and rate constants $\{\kappa_k\}$ the system is complex balanced with a complex balanced equilibrium $c\in \mathbb{R}^m_{\geq 0}$. Then the stochastically modeled system with general intensity functions (3.2), with the same rate constants $\{\kappa_k\}$, admits the stationary measure

In Section 4, we will show that $\pi$ in (3.3) is summable under some mild growth condition on $\theta_i$, and thus it can be normalized to a stationary distribution.

The result in Theorem 3.1 is a cornerstone for later research in both application and theory. It has been used in various studies in applied and computational biology such as tumor growth and invasion ^[19], phosphorylation systems ^[9], studies of the chemical master equation ^[26], etc. It is also the background for many theoretical developments in the field of chemical reaction networks ^[2,3]. Interestingly, it has been proven that the converse is also true.

Theorem 3.3. (^[3]). Let $\{\mathcal{S, R, C}\}$ be a reaction network and consider the stochastically modeled system with rate constants $\{\kappa_k\}$ and mass action kinetics (2.1). Suppose that for some $c \in \mathbb{R}^m_{\ge 0}$ the stationary distribution for the stochastic model is (3.1). Then $c$ is a complex balanced equilibrium for the associated deterministic model with mass action kinetics and rate constants $\{\kappa_k\}$.

3.2.   Non-equilibrium potential and Lyapunov functions

Another interesting result comes from the scaling behavior of the stationary distribution for complex balanced system. We first provide a key definition.

Definition 3.1. Let $\pi$ be a probability distribution on a countable set $\Gamma$ such that $\pi(x) > 0$ for all $x\in\Gamma$. The non-equilibrium potential of the distribution $\pi$ is the function $\phi_\pi: \Gamma \to \mathbb{R}$, defined by

In ^[2] it was shown that under an appropriate scaling, the limit of the non-equilibrium potential of the stationary distribution of a complex balanced system converges to a certain well-known Lyapunov function.

Definition 3.2. Let $E \subset \mathbb{R}^m_{\geq 0}$ be an open subset of $\mathbb{R}^m_{\geq 0}$ and let $f: \mathbb{R}^m_{\geq 0} \to \mathbb{R}$. A function $\mathcal{V}: E\to \mathbb{R}$ is called a Lyapunov function for the system $\dot x = f(x)$ at $x_0 \in E$ if $x_0$ is an equilibrium point for $f$, that is $f(x_0) = 0$, and

1. $\mathcal{V}(x) > 0$ for all $x\neq x_0, x\in E$ and $\mathcal{V}(x_0) = 0$.

2. $\nabla\mathcal{V}(x)\cdot f(x)\leq 0$, for all $x\in E$, with equality if and only if $x = x_0$, where $\nabla\mathcal{V}$ denotes the gradient of $\mathcal{V}$.

In particular, the non-equilibrium potential of the stationary distribution in (3.1) converges to the usual Lyapunov function of Chemical Reaction Network Theory

The next section discusses the scaling in which the convergence happens. It is called the classical scaling in the literature.

3.2.1.   The classical scaling

Here we present a brief introduction to the classical scaling. For more detailed discussions, see ^[20,21,22].

Let $|y_k| = \sum_i y_{ki}$ and let $V$ be the volume of the system times Avogadro's number. Suppose $\{\kappa_k\}$ are the rate constants for the stochastic model. We defined the scaled rate constants as follows

and denote the scaled intensity function for the stochastic model by

Note that if $x \in \mathbb Z^m_{\ge 0}$ gives the counts of the different species, then $\tilde x : = V^{-1} x$ gives the concentrations in moles per unit volume. Then, by standard arguments

where the final equality defines $\lambda_k$ and justifies the definition of deterministic mass action kinetics.

Denote the stochastic process determining the counts by $X^V(t)$, then normalizing the original process $X^V$ by V and defining $\bar X^V: = \frac{X^V}{V}$ gives us

where we are utilizing the representation (2.2). Since the law of large numbers for the Poisson process implies $V^{-1}Y(Vu)\approx u$, we may conclude that a good approximation to the process $\bar X^V$ is the function $x = x(t)$ defined as the solution to the ODE

which is exactly (2.5).

A corollary of Theorem 3.1 gives us the stationary distribution for the classically scaled system.

Theorem 3.4. Let $\{\mathcal{S, R, C}\}$ be a reaction network. Suppose that when modeled deterministically with mass action kinetics and rate constants $\{\kappa_k\}$ the system is complex balanced with a complex balanced equilibrium $c\in \mathbb{R}^m_{\geq 0}$. For $V > 0$, let $\{\kappa_k^V\}$ satisfy (3.5). Then the stochastically modeled system on $\mathbb Z^d_{\ge 0}$ with rate constants $\{\kappa_k^V\}$ and intensity functions (3.6) admits the stationary distribution

An immediate implication of Theorem 3.4 is that a stationary distribution for the scaled model $\bar X^V$ is

3.2.2.   Convergence of the non-equilibrium potential

The main finding in ^[2] is concerned with the scaling limit of the stationary distribution $\tilde \pi^V$ of (3.8).

Theorem 3.5. Let $\{S, C, R\}$ be a reaction network and let $\{\kappa_k\}$ be a choice of rate constants. Suppose that, modeled deterministically, the system is complex balanced. For $V > 0$, let $\{\kappa^V_k\}$ be related to $\{ \kappa_k\}$ via (3.5). Fix a sequence of points $\tilde x^V\in \frac1V \mathbb Z^d_{\ge 0}$ for which $\lim_{V\to \infty}\tilde {x}^V = \tilde x \in \mathbb R^d_{ > 0}$. Further let $c$ be the unique complex balanced equilibrium within the positive stoichiometric compatibility class of $\tilde x$.

Let $\pi^V$ be given by (3.7) and let $\tilde \pi^V$ be as in (3.8), then

where $\mathcal V$ is the Lyapunov function for the ODE model satisfying (3.4).

4.   Main results for networks with general kinetics

In this section, we first show that for stochastically modeled reaction networks with non-mass action kinetics defined via (3.2) whose associated mass action system is complex balanced, the stationary measure (3.3) can be normalized to yield a stationary distribution. We further show that these stochastic models are non-explosive. We then extend Theorems 3.3 and 3.5 from Section 3 to the non-mass action case.

4.1.   Existence of a stationary distribution and non-explosivity of non-mass action systems

We begin with a theorem proving that the stochastic models considered in Theorem 3.2 are positive recurrent when only mild growth conditions are placed on the functions $\theta_i$.

Theorem 4.1. Let $\{S, C, R\}$ be a reaction network with rate constants $\{\kappa_k\}$. Suppose that when modeled deterministically, the associated mass action system is complex balanced with equilibrium $c\in \mathbb{R}^m_{ > 0}$. Suppose that $\theta_i$ and $\lambda_k$ satisfy the conditions in and around (3.2). Moreover, suppose that for each $i$ we have $\lim_{x\to \infty} \theta_i(x) = \infty$. Then,

1. the measure $\pi$ given in (3.3) is summable over $\mathbb Z^m_{\ge 0}$, and a stationary distribution exists for the stochastically modeled process, and moreover

2. the stochastically modeled process is non-explosive.

Proof. We first show that $\pi$ is summable over $\mathbb Z^m_{\ge 0}$. We have

so long as each sum in the final expression is finite. Thus it is sufficient to prove that $\sum\limits_{x\in \mathbb Z_{\ge 0}} \frac{c_i^{x}}{\theta_i(1)\cdots\theta_i(x)}$ is finite for each $i$. By the ratio test

where the last equality is due to the assumption that $\lim_{x\to \infty} \theta_i(x) = \infty$. Hence the sum is convergent.

We turn to showing that the process is non-explosive. From ^[5], to show that the process is non-explosive, it is sufficient to show

From (4.10) and (4.12), we need to show

Let $s_i = \max_k \{y_{ki}\}$, where the max is over all source complexes, and let $R$ be the number of reactions. Let $n_i > s_i$ be such that $\theta_i(x_i) > 1, \cdots, \theta_i(x_i-s_i+1) > 1$ for all $x_i > n_i$. Then

where $C = R\max_{i = 1}^m\{c_i^{s_i}\}$, and the last inequality follows from part 1. Thus the process is non-explosive.

4.2.   Generalization of Theorem 3.3

We are set to provide the next theorem of the current paper, which is the converse statement of Theorem 3.2 and generalizes Theorem 3.3 from ^[3]. In the theorem below, we assume $\lim_{x\to \infty} \theta_i(x) = \infty$ for each $i$. In Corollary 1, we generalize the result to allow $\lim_{x\to \infty} \theta_i(x) \in \{0, \infty\}$ for each $i$.

Theorem 4.2. Let $\{\mathcal{S, R, C}\}$ be a reaction network and consider the stochastically modeled system with rate constants $\{\kappa_k\}$ and intensity functions (3.2). Suppose that $\lim_{x\to \infty} \theta_i(x) = \infty$ for each $i = 1, \dots, m$ and that for some $c \in \mathbb{R}^m_{\ge 0}$ a stationary measure for the stochastic model satisfies (3.3). Then $c$ is a complex balanced equilibrium for the associated deterministic model with mass action kinetics and rate constants $\{\kappa_k\}$.

Proof. By assumption, we have that $\pi$ satisfies

Plugging (3.2) and (3.3) into this equation yields

Canceling and moving terms when necessary, we have

Enumerating the reaction on the right by their product complexes, and the reactions on the left by their source complexes, the equation above becomes

Since the above holds for all $x\in \mathbb{Z}^m_{\ge 0}$, the two sides are equal as functions. Hence, if the functions in the set

are linearly independent, then we must have

which is the condition for the associated mass action system to be complex balanced.

Thus it remains to show that the functions in the set (4.1) are linearly independent. We will prove that the functions are linearly independent by induction on the number of species, and begin with the base case $m = 1$, which we provide as a lemma.

Lemma 4.1. When $m = 1$ (i.e. the system has only one species) the functions in the set (4.1) are linearly independent.

Proof of Lemma 4.1. Let $\mathcal{C} = \{z_1, \dots, z_n\}$ ordered so that $z_i < z_{i+1}$ for each $i = 1, \dots, n-1$. Suppose, in order to find a contradiction, the functions in the set (4.1) are linearly dependent. Then there exist $\alpha_1, \cdots, \alpha_r \in \mathbb R$ with $r\le n$ and $\alpha_r \ne 0$, such that

Let $M = \frac{|\alpha_1|}{|\alpha_r|}+\cdots+\frac{|\alpha_{r-1}|}{|\alpha_r|}$. Since $\theta(x)\to\infty$, as $x\to \infty$, we can find an $N > 0$ such that $\forall x > N$, we have $\theta(x-z_r+1) > M$ and $\theta(x), \dots, \theta(x-z_r+1)\geq 1$. In this case,

This contradicts (4.2). Therefore, (4.1) must be linearly independent.

We turn to the inductive step. Thus, we now assume that functions of the form (4.1) for distinct complexes $z$ are linearly independent when there are $m-1$ species. We must show that this implies linear independence when there are $m$ species.

Enumerate the complexes as $\mathcal{C} = \{z^1, z^2, \ldots, z^n\}$. Suppose that there are $\alpha_1, \dots, \alpha_n$ for which

We will show that each $\alpha_i = 0$.

First note that we can not have $z^1_i = z^2_i = \ldots = z^n_i$ for each $i = 1, \dots, m$, for otherwise all the complexes are the same. Thus, and without loss of generality, we assume that not all of the $z_1^k$ are equal. In particular, we will assume that $z_1^1, \dots, z_1^n$ consists of $p$ distinct values with $2\le p \le n$. We will also assume that the complexes are ordered so that the first $r_1$ terms of $z_1^k$ are the same, the second $r_2$ terms are the same, etc. That is,

We now consider the left hand side of (4.3) as a function of $x_1$ alone. For $j = 1, \dots, p$, we define

By Lemma 4.1, the functions $f_j,$ for $j = 1, \dots, p$, are linearly independent. Combining similar terms in (4.3) we have

From the independence of the $f_j$, it must be the case that each bracketed term above is zero.

Without loss of generality, we just consider the first bracketed term in (4.5):

which we know is equal to zero. The goal now is to apply our inductive hypothesis to conclude that each of $\alpha_1, \dots, \alpha_{r_1}$ is equal to zero.

For each of $k = 1, \dots, r_1$, we let $\tilde z^k = (z^k_2, \dots, z^k_m)$. Then each term in the sum (4.6) is a function on $\mathbb R^{m-1}$ of the general form (4.1) with new complexes $\tilde z^k \in \mathbb R^{m-1}$. To use the inductive hypothesis, we must argue that the $\tilde z^k$ are distinct. Consider, for example, the first two terms: $\tilde z^1$ and $\tilde z^2$. By (4.4), we know that $z^1_1 = z^2_1$; that is, the coefficient of species 1 for the two complexes are the same. If we also had $\tilde z^1 = \tilde z^2$, then all the coefficients of the species would be the same for the two complexes, contradicting the fact that they are distinct complexes (i.e. $z^1 \ne z^2$). Hence, it must be that $\tilde z^1 \ne \tilde z^2$. Thus, by the inductive hypothesis, all the terms of the sum (4.6) are linearly independent, and $\alpha_1 = \cdots = \alpha_{r_1} = 0$. Repeating this argument for the other bracketed terms completes the proof.

We have proven the independence of (4.1) in all cases, which completes the proof of the theorem.

The following relaxes the condition in Theorem 4.2 that the limit of the functions $\theta_i$ must be infinity.

Corollary 1. Let $\{\mathcal{S, R, C}\}$ be a reaction network and consider the stochastically modeled system with rate constants $\{\kappa_k\}$ and intensity functions (3.2). Suppose that $\lim_{x\to \infty} \theta_i(x) \in\{0, \infty\}$ for each $i = 1, \dots, m$ and that for some $c \in \mathbb{R}^m_{\ge 0}$ a stationary measure for the stochastic model satisfies (3.3). Suppose further that $\theta_i(x) > 0$ for $x$ large enough. Then $c$ is a complex balanced equilibrium for the associated deterministic model with mass action kinetics and rate constants $\{\kappa_k\}$.

Proof. Without loss of generality, assume $\lim_{x\to \infty} \theta_i(x) = 0$ for $i\le \ell$ and $\lim_{x\to \infty} \theta_i(x) = \infty$ for $i \ge \ell + 1.$ The proof is the same as that of Theorem 4.2 in that we must prove the linear independence of the functions in (4.1). Let

For $x$ large enough that $\theta_i(x) > 0$, let $\phi_i(x) = \frac{1}{\theta_i(x)}$ for each $i \le \ell$. Then we have $\lim_{x\to \infty} \phi_i(x) = \infty$. Now (4.7) becomes

Let $w_k = \max_{1\le j \le n}\{z_k^j\}$. Then from (4.8) we have

This is similar to the set-up of Theorem 4.2 (since each $\phi_i(x) \to \infty$, as $x \to \infty$) and we can conclude $\alpha_1 = \ldots = \alpha_n = 0$ and complete the proof.

4.3.   Generalization of Theorem 3.5

This section is concerned with the convergence of the non-equilibrium potential of the stationary distribution of systems with general kinetics, under some appropriate scaling. In particular, we would like to have a similar result as Theorem 3.5 for the case of general kinetics. One difficulty that arises is that the classical scaling is not, in general, appropriate for our purposes. This is illustrated by the example below.

Example 1. Consider the reaction network with one species A and reactions given by

with the intensity function given by (3.2), where the rate constants are $\kappa_{\emptyset\to A} = \kappa_{A\to\emptyset} = 1$ and $\theta(x) = x^2$. Consider the process under the classical scaling. We obtain a stationary distribution for the scaled model $\tilde X^V$ from Theorem 3.4 and (3.8):

We consider the limiting behavior of the non-equilibrium potential $-\frac{1}{V}\ln(\tilde\pi^V(\tilde x^V))$, as $V\to \infty$. Using Stirling's approximation,

We need to estimate $M$ when $V\to \infty$. From Lemma in the Appendix, we have

for some constants $a, b \in \mathbb R$. Thus

Clearly, $\lim_{V\to \infty}-\frac{1}{V}\ln(\tilde\pi^V(\tilde x^V)) = \infty$, and we do not have convergence of the non-equilibrium potential under the classical scaling.

With the above example in mind, we provide an alternative scaling.

4.3.1.   The modified scaling

Define $|y_k| = \sum_i y_{ki}$ and let $V$ be a scaling parameter. For each reaction $y_k\to y_k'$ let $\kappa_k$ be a positive parameter. We now define the rate constant for $y_k\to y_k'$ as

where the parameter $d$ is a vector to be chosen (they will depend upon the limiting values $\lim_{x\to \infty} \theta_i(x)$). Note that the classical scaling is the case when $d = (1, 1, \ldots, 1)$. Then we define the scaled intensity function

where, as usual, $\theta_i: \mathbb{Z} \to \mathbb{R}_{\geq 0}$, and $\theta_i(x) = 0$ if $x\leq 0$.

Theorem 4.3. Let $\{S, C, R\}$ be a reaction network with rate constants $\{\kappa_k\}$. Suppose that when modeled deterministically, the associated mass action system is complex balance with equilibrium $c\in \mathbb{R}^m_{ > 0}$. For some V, let $\{\kappa_k^V\}$ be related to the $\{\kappa_k\}$ via (4.9). Then the stochastically modeled system with scaled intensity function (4.10) has stationary measure

If (4.11) is summable, then a normalizing constant $M$ can be found so that

is a stationary distribution.

Proof. The proof is similar to that of Theorem 3.1, as found in ^[1], except care must be taken to ensure that the terms associated with the scaling parameter $V$ cancel appropriately.

Let $X^V$ be the process associated with the intensities (4.10) and let $\tilde X^V = V^{-1} X^V$ be the scaled process. By Theorem 4.3, we have that for $x^V \in \frac1V \mathbb Z^m_{\ge 0}$

is a stationary measure for the scaled process. If (4.13) is summable, which is ensured by Theorem 4.1 so long as $\theta_i(x) \to \infty$ as $x\to \infty$, then

is a stationary distribution. In the next section, we consider the the limiting behavior of $-\frac1V \ln(\tilde \pi^V(\tilde x^V))$ as $V \to \infty$ for a class of $\theta_i$.

4.3.2.   Limiting behavior of $-\frac1V \ln(\tilde \pi^V(\tilde x^V))$

We make the following assumption on the functions $\theta_i$.

Assumption 1. We assume that (ⅰ) $\theta_i: \mathbb{Z} \to \mathbb{R}_{\geq 0}$, (ⅱ) $\theta_i(x) = 0$ if $x\leq 0$, and (ⅲ) there exists $d, A\in \mathbb{R}^m_{ > 0}$ such that $\lim_{x_i\to \infty} \frac{\theta_i(x_i)}{x_i^{d_i}} = A_i$ for each $i$.

Roughly speaking, this class of functions act like power functions when $x$ is large. We will utilize functions satisfying Assumption 1 to build intensity functions as in (4.10). We will show that if the deterministic mass action system is complex balanced, then the limiting behavior of the scaled non-equilibrium potential of the stochastically modeled system with intensities (4.10) is a Lyapunov function for the ODE system

where we recall that $Ax^d$ is the vector with $i$th component $A_ix_i^{d_i}$. This result therefore generalizes Theorem 3.5 (which is Theorem 8 in ^[2]).

Lemma 4.2. Let $\{S, C, R\}$ be a reaction network with rate constants $\{\kappa_k\}$. Suppose that when modeled deterministically, the associated mass action system is complex balanced with equilibrium $c\in \mathbb{R}^m_{ > 0}$. Let $d, A\in \mathbb R^m_{ > 0}$. Then the system (4.15) is complex balanced with equilibrium vector $\tilde c$ satisfying

Proof. The proof consists of verifying that for each $z\in \mathcal C$,

where the sum on the left consists of those reactions with source complex $z$ and the sum on the right consists of those with product complex $z$. This is immediate from the definition of $\tilde c$.

We now turn to the scaled models, and prove that the properly scaled non-equilibrium potential converges to a Lyapunov function for the ODE system (4.15).

Theorem 4.4. Let $\{S, C, R\}$ be a reaction network with rate constants $\{\kappa_k\}$. Suppose that when modeled deterministically, the associated mass action system is complex balanced with equilibrium $c\in \mathbb{R}^m_{ > 0}$.

Fix $d, A\in \mathbb R^m_{ > 0}$ and let $\theta_i$ be a choice of functions satisfying Assumption 1. For $V > 0$ and the $d > 0$ already selected, let $\{\kappa_k^V\}$ be related to $\{\kappa_k\}$ as in (4.9) and let the intensity functions for the stochastically modeled system be (4.10). Let $\tilde \pi^V$ be the stationary distribution for the scaled process guaranteed to exist by Theorems 4.3 and 4.1 and given by (4.12).

Fix a sequence of points $\tilde x^V \in \tfrac1V \mathbb Z^m_{\ge 0}$ for which $\lim_{V\to \infty} \tilde x^V = \tilde x \in \mathbb Z^m_{ > 0}$. Then

where $\mathcal V$ is defined by the final equality, and moreover $\mathcal V$ is a Lyapunov function for the ODE system (4.15).

Note that by taking $d = (1, .., 1)$ and $A = (1, .., 1)$, the limit of the $\theta_i$ in Assumption 1 is simply mass action kinetics. Hence, the main result in ^[2] is contained within the above theorem.

Proof. Using (4.12) and (4.14) we have

We analyze the limiting behavior of the different pieces of the last expression.

1. We begin with the first term

where $M$ is defined using (4.12). In Lemma in the appendix we show that as $V\to\infty$

where by $a_V\sim b_V$, as $V\to \infty$, we mean $\lim_{V\to \infty} (a_V - b_V) = 0$. We may then apply Lemma to (4.17) to conclude there are constants $a, b$ such that

Taking the limit $V\to \infty$, we see that only the first term remains, which yields

2. We use Stirling's approximation with the middle terms

Taking the limit $V \to \infty$, and noting that $\tilde x^V \to \tilde x$, we have

3. We turn to the final term. By using an argument similar to (4.17), there is a constant $C > 0$ for which

Taking the limit $V \to \infty$, and noting that $\tilde x^V \to \tilde x$, we have

Combining the three parts, we conclude (4.16) holds. The fact that the limit is a Lyapunov function is proven in Lemma 4.3 below.

Lemma 4.3. The function given by (4.16),

is a Lyapunov function for the system (4.15).

Proof. We have

Let $f$ be the right-hand side of (4.15) and recall that $c$ is a complex balanced equilibrium of the mass action model. We have

where we used the inequality $a(\ln{b}-\ln{a})\leq b-a$ and the last equation holds because $c$ is the complex balanced equilibrium for the mass-action system.

Acknowledgments

We gratefully acknowledge financial support from Army Research Office grant W911NF-18-1-0324.

Conflict of interest

The authros declare no conflict of interest in this paper.