Detailed balance, local detailed balance, and global potential for stochastic chemical reaction networks

TL;DR

This paper explores various forms of detailed balance in stochastic chemical reaction networks, establishing their relationships, implications for global potential functions, and providing conditions for stochastic detailed balance with applications to high-dimensional systems.

Contribution

It clarifies the equivalence and differences among four types of detailed balance conditions and introduces a new sufficient condition for stochastic detailed balance in complex networks.

Findings

Four types of detailed balance are equivalent under certain conditions.

Local detailed balance implies the existence of a global potential function.

A new sufficient condition for stochastic detailed balance is proposed, enabling the construction of multistable high-dimensional networks.

Abstract

Detailed balance of a chemical reaction network can be defined in several different ways. Here we investigate the relationship among four types of detailed balance conditions: deterministic, stochastic, local, and zero-order local detailed balance. We show that the four types of detailed balance are equivalent when different reactions lead to different species changes and are not equivalent when some different reactions lead to the same species change. Under the condition of local detailed balance, we further show that the system has a global potential defined over the whole space, which plays the central role in the large deviation theory and the Freidlin-Wentzell-type metastability theory of chemical reaction networks. Finally, we provide a new sufficient condition for stochastic detailed balance, which is applied to construct a class of high-dimensional chemical reaction networks that both satisfies stochastic detailed balance and displays multistability.

Figures

Equations

$$R_i:\sum _{j=1}^d{\nu }_i^j{S}_j\stackrel{k_i}{\to }\sum _{j=1}^d{\nu }_i^{\prime j}{S}_j,1\le i\le N,$$

$$\mathcal{S}=\{S_1,\dots ,S_d\},\mathcal{C}=\{{\nu }_i,{\nu }_i^{\prime }|i=1,\dots ,N\},\mathcal{R}=\{R_i|i=1,\dots ,N\}$$

$${\mathcal{R}}_p^{+}$$

$${\mathcal{R}}_p^{-}$$

$$R_{pl}^{+}:$$

$$R_{pl}^{-}:$$

$$E_V=\left\{\frac{n}{V}:n=(n_1,\dots ,n_d)\in {\mathbb{Z}}_{\ge 0}^d\right\}$$

$$q_{x,x+\frac{{\omega }_p}{V}}^V=\sum _{l=1}^{r_p}\frac{k_{pl}^{+}}{V^{|{\nu }_{pl}|-1}}\frac{(Vx)!}{(Vx-{\nu }_{pl})!},$$

$$q_{x,x-\frac{{\omega }_p}{V}}^V=\sum _{l=1}^{r_p}\frac{k_{pl}^{-}}{V^{|{\nu }_{pl}^{\prime }|-1}}\frac{(Vx)!}{(Vx-{\nu }_{pl}^{\prime })!},$$

$$q_{x,x}^V=-\sum _{p=1}^r\left(q_{x,x+\frac{{\omega }_p}{V}}^V+q_{x,x-\frac{{\omega }_p}{V}}^V\right),$$

$$\begin{array}{rl}\frac{\mathrm{d}{\mathbb{P}}_x^V(t)}{\mathrm{d}t}=& \sum _{p=1}^r{q}_{x-\frac{{\omega }_p}{V},x}^V{\mathbb{P}}_{x-\frac{{\omega }_p}{V}}^V(t)+\sum _{p=1}^r{q}_{x+\frac{{\omega }_p}{V},x}^V{\mathbb{P}}_{x+\frac{{\omega }_p}{V}}^V(t)\\ & -\sum _{p=1}^r\left(q_{x,x+\frac{{\omega }_p}{V}}^V+q_{x,x-\frac{{\omega }_p}{V}}^V\right){\mathbb{P}}_x^V(t),x\in E_V.\end{array}$$

$$\{\begin{array}{rl}\frac{\mathrm{d}x(t)}{\mathrm{d}t}& =\sum _{p=1}^r[f_p^{+}(x(t))-f_p^{-}(x(t))]{\omega }_p,\\ x(0)& =x_0,\end{array}$$

$$f_p^{+}(x)=\sum _{l=1}^{r_p}{k}_{pl}^{+}{x}^{{\nu }_{pl}},f_p^{-}(x)=\sum _{l=1}^{r_p}{k}_{pl}^{-}{x}^{{\nu }_{pl}^{\prime }},x\in {\mathbb{R}}_{\ge 0}^d,$$

$$\underset{V\to \infty }{\lim }{\mathbb{P}}_{x_0^V}(\underset{t\in [0,T]}{\sup }\Vert X^V(t)-x(t)\Vert \le \delta )=1,$$

$$H(x,\theta )=\sum _{p=1}^r\left[f_p^{+}(x)\left(e^{{\omega }_p\cdot \theta }-1\right)+f_p^{-}(x)\left(e^{-{\omega }_p\cdot \theta }-1\right)\right],x\in {\mathbb{R}}_{\ge 0}^d,\theta \in {\mathbb{R}}^d,$$

$$L(x,y)=\underset{\theta \in {\mathbb{R}}^d}{\sup }(\theta \cdot y-H(x,\theta )),x\in {\mathbb{R}}_{\ge 0}^d,y\in {\mathbb{R}}^d.$$

$$I_{x_0,T}(\varphi )=\{\begin{array}{ll}{\int }_0^{T}L(\varphi (t),\dot{\varphi }(t))\mathrm{d}t,& \text{if ϕ is absolutely continuous and ϕ(0)=x0},\\ \infty ,& \text{otherwise}.\end{array}$$

$$\begin{array}{rl}\underset{V\to \infty }{\mathrm{liminf}}\frac{1}{V}\log {\mathbb{P}}_{x_0^V}(X^V(\cdot )\in {\Gamma }^o)& \ge -\underset{\varphi \in {\Gamma }^o}{\inf }{I}_{x_0,T}(\varphi ),\\ \underset{V\to \infty }{\mathrm{limsup}}\frac{1}{V}\log {\mathbb{P}}_{x_0^V}(X^V(\cdot )\in \bar{\Gamma })& \le -\underset{\varphi \in \bar{\Gamma }}{\inf }{I}_{x_0,T}(\varphi ),\end{array}$$

$$W(x_0,y)=\inf \{I_{x_0,T}(\varphi ):\varphi (0)=x_0,\varphi (T)=y,T\ge 0\},x_0,y\in {\mathbb{R}}_{\ge 0}^d.$$

$$\underset{V\to \infty }{\lim }\frac{1}{V}\log {\mathbb{E}}_{x_0^V}{\tau }_V=\underset{y\in \partial D}{\inf }W(c,y),$$

$$W(c,y_0)=\underset{y\in \partial D}{\inf }W(c,y),$$

$$\underset{V\to \infty }{\lim }{\mathbb{P}}_{x_0^V}(\Vert X^V({\tau }_V)-y_0\Vert <\delta )=1,$$

$$\underset{\delta \to 0}{\lim }\underset{V\to \infty }{\lim }\frac{1}{V}\log {\mathbb{P}}_{x_0^V}(\Vert X^V({\tau }_V)-z_0\Vert <\delta )=W(c,y_0)-W(c,z_0).$$

$$k_{pl}^{+}{c}^{{\nu }_{pl}}=k_{pl}^{-}{c}^{{\nu }_{pl}^{\prime }},\text{for any}1\le p\le r,1\le l\le r_p.$$

$${\pi }_x^V{q}_{x,y}^V={\pi }_y^V{q}_{y,x}^V,\text{for any}x,y\in E_V.$$

$$q_{x_1,x_2}^V{q}_{x_2,x_3}^V\cdots q_{x_n,x_1}^V=q_{x_2,x_1}^V{q}_{x_3,x_2}^V\cdots q_{x_1,x_n}^V$$

$$\sum _{p=1}^r{\xi }_p\log \frac{f_p^{+}(x)}{f_p^{-}(x)}=0,\text{for any}x\in {\mathbb{R}}_{>0}^d,$$

$${\omega }_q\cdot \nabla \left(\log \frac{f_p^{+}(x)}{f_p^{-}(x)}\right)={\omega }_p\cdot \nabla \left(\log \frac{f_q^{+}(x)}{f_q^{-}(x)}\right),\text{for any}x\in {\mathbb{R}}_{>0}^d.$$

$$\begin{array}{rl}& C:0\to \mathrm{sgn}({\xi }_1){\omega }_1\to \cdots \to {\xi }_1{\omega }_1\\ & \to {\xi }_1{\omega }_1+\mathrm{sgn}({\xi }_2){\omega }_2\to \cdots \to {\xi }_1{\omega }_1+{\xi }_2{\omega }_2\to \cdots \\ & \to {\xi }_1{\omega }_1+{\xi }_2{\omega }_2+\cdots +\mathrm{sgn}({\xi }_r){\omega }_r\to \cdots \to {\xi }_1{\omega }_1+{\xi }_2{\omega }_2+\cdots +{\xi }_r{\omega }_r=0,\end{array}$$

$$f\left(\frac{1}{V}\right):=\log \frac{q_{{\eta }_1,{\eta }_2}^V{q}_{{\eta }_2,{\eta }_3}^V\cdots q_{{\eta }_L,{\eta }_1}^V}{q_{{\eta }_2,{\eta }_1}^V{q}_{{\eta }_3,{\eta }_2}^V\cdots q_{{\eta }_1,{\eta }_L}^V}=0.$$

$$f(0):=\underset{V\to \infty }{\lim }f\left(\frac{1}{V}\right)=0\text{(zero-order information)},$$

$$f^{\prime }(0):=\underset{V\to \infty }{\lim }\frac{f\left(\frac{1}{V}\right)-f(0)}{\frac{1}{V}}=0\text{(first-order information)}.$$

Full text

Keywords: chemical reaction systems; microscopic reversibility; Kolmogorov’s cycle condition; quasi-potential; large deviation; metastability

1 Introduction

The mathematical theory of chemical reaction networks has attracted massive attention over the past two decades due to its wide applications in biology, chemistry, ecology, and epidemics [1]. If a reaction system is well mixed and the numbers of molecules are very large, random fluctuations can be ignored and the evolution of the concentrations of all chemical species can be modeled deterministically as a set of ordinary differential equations (ODEs) based on the law of mass action. If the chemical species are presented in low numbers, however, random fluctuations can no longer be ignored and the evolution of the system is usually modeled stochastically as a continuous-time Markov chain on a high-dimensional lattice, which is widely known as a density-dependent Markov chain [2]. The Kolmogorov forward equation for a density-dependent Markov chain is the well-known chemical master equation, which is first introduced by Delbrück [3]. At the center of the mathematical theory of chemical reaction networks is a limit theorem proved by Kurtz [4, 5, 6], which states that when the system size tends to infinity, the trajectories of the stochastic model of a reaction system will converge to those of the deterministic model over any compact time interval, whenever the initial condition converges. Thus far, stochastic reaction networks have served as a fundamental model for the single-cell stochastic gene expression dynamics of gene regulatory networks [7, 8, 9, 10, 11, 12, 13, 14]. Recently, the limit theorem of Kurtz has been generalized to stochastic gene regulatory networks with bursting dynamics [15, 16].

The limit theorem of Kurtz [5] can be viewed as the law of large numbers for stochastic reaction networks. The corresponding large deviation theory and the Freidlin-Wentzell-type metastability theory for stochastic reaction networks have also been studied by many authors [17, 18, 19, 20, 21, 22] and were rigorously established by Agazzi et al. under the mass action kinetics [23]. At the center of the metastability theory is a quantity called quasi-potential, which plays a crucial role in the analysis of the exit time and exit distribution from a basin of attraction, as well as the most probable transition paths between multiple attractors when the system size is large [24]. However, the quasi-potential is usually defined via an abstract variational expression and hence it is very difficult to obtain a general analytical expression for the quasi-potential.

There are two types of reaction networks that should be distinguished, those satisfy detailed balance and those violate detailed balance. In terms of physical chemistry, detailed balance is a fundamental thermodynamic constraint for closed systems. If there is no sustained energy supply, then a chemical system, when it reaches the steady state, must satisfy detailed balance [25]. In the modelling of many biochemical systems such as enzymes [26] and ion channels [27], detailed balance has become a basic requirement [28, 29]. However, in the literature, there are two different definitions of detailed balance for a chemical reaction network. From the deterministic perspective, detailed balance means that there is no net concentration flux between any pair of reversible reactions, in which case there is no chemical potential difference and thus the system is in chemical equilibrium. From the stochastic perspective, detailed balance means that there is no net probability flux between any pair of microstates on the high-dimensional nonnegative integer lattice, where each microstate is defined as the ordered tuple of concentrations of all chemical species. To distinguish between them, we refer to the former as deterministic detailed balance and the latter as stochastic detailed balance. Some authors believed that the two types of detailed balance are equivalent [30]. However, Joshi [31] pointed out recently that they are not equivalent; they are equivalent when different reactions lead to different species changes, while they are in general not equivalent for systems having two reactions that lead to the same species change — deterministic detailed balance implies stochastic detailed balance and the opposite is not true.

In this paper, in addition to deterministic and stochastic detailed balance, we propose a third type of detailed balance, which is called local detailed balance. This new type of detailed balance characterizes the local asymptotic behavior of a reaction network as the system size tends to infinity. We prove that the three types of detailed balance (deterministic, stochastic, and local) are equivalent when different reactions lead to different species changes, while local detailed balance is even weaker than the other two when some different reactions lead to the same species change — stochastic detailed balance implies local detailed balance and the opposite is not true. This is the first main result of the present paper. More importantly, under the condition of local detailed balance, we prove that a stochastic reaction network has a global potential that can be computed explicitly and concisely. The global potential reduces to the quasi-potential within each basin of attraction. In general, the quasi-potential is only defined within each basin of attraction. However, local detailed balance guarantees that the system has a global potential that can be defined over the whole space. This is the second main result of the present paper.

In [31], Joshi gave the sufficient and necessary condition for deterministic detailed balance. While the author also provided a weaker sufficient condition for stochastic detailed balance, it is difficult to apply it in practice since an infinite number of restrictions need to be verified. In this paper, we provide a simpler sufficient condition for stochastic detailed balance that is more applicable in practice. This new sufficient condition is imposed directly on rate constants and only a finite number of restrictions need to be verified. This is the third main result of the present paper.

The structure of this paper is organized as follows. In Section 2, we recall the basic concepts of chemical reaction networks and state some preliminary results. In Section 3, we reveal the relationship among four types of detailed balance and give some counterexamples. In Section 4, we show that a global potential exists for stochastic reaction networks satisfying local detailed balance and obtain the explicit expression of the global potential. The remaining sections are devoted to the detailed proofs of the main theorems.

2 Model and preliminary results

Let $\mathbb{Z}_{\geq 0},\mathbb{R}_{\geq 0}$, and $\mathbb{R}_{>0}$ denote the sets of nonnegative integers, nonnegative real numbers, and positive real numbers, respectively. Recall that a chemical reaction system is composed of a collection of chemical species $\{S_{1},\dots,S_{d}\}$ and a family of reactions

[TABLE]

where $\nu^{j}_{i},\nu^{\prime j}_{i}\in\mathbb{Z}_{\geq 0}$ are the molecule numbers of $S_{i}$ consumed and created in one instance of that reaction, respectively. For simplicity, we write $\nu_{i}=(\nu_{i}^{1},\dots,\nu_{i}^{d})$ and $\nu^{\prime}_{i}=(\nu^{\prime 1}_{i},\dots,\nu^{\prime d}_{i})$, which are called complexes. Then the reaction $R_{i}$ can be written more concisely as $\nu_{i}\rightarrow\nu^{\prime}_{i}$. Moreover, $\nu^{\prime}_{i}-\nu_{i}$ is called the reaction vector of $R_{i}$. Let

[TABLE]

denote the collections of chemical species, complexes, and reactions respectively. Then the ordered triple $\{\mathcal{S},\mathcal{C},\mathcal{R}\}$ is called a chemical reaction network.

A chemical reaction network is called reversible if for any reaction $R_{i}:\nu_{i}\rightarrow\nu^{\prime}_{i}\in\mathcal{R}$, there exists a reverse reaction $R^{-}_{i}:\nu^{\prime}_{i}\rightarrow\nu_{i}\in\mathcal{R}$ [31]. For any pair of reversible reactions $R_{i}$ and $R_{i}^{-}$, we say that $R_{i}$ is a forward reaction if $\nu_{i}<\nu^{\prime}_{i}$, where the symbol “$<$” is understood in the lexicographic order, namely $\nu_{i}<\nu^{\prime}_{i}$ if and only if $\nu_{i}^{j}<\nu^{\prime j}_{i}$ for the first $j$ at which $\nu_{i}^{j}$ and $\nu^{\prime j}_{i}$ differ; otherwise, $R_{i}$ is called a backward reaction. Throughout the paper, we assume that all reaction networks under consideration are reversible.

Most previous papers assumed that different reactions have different reaction vectors. However, in many reaction networks, multiple different reactions may have the same reaction vector. For example, the reaction $S_{1}\rightarrow S_{2}$ and the enzyme catalyzed reaction $S_{1}+E\rightarrow S_{2}+E$ may coexist in a biochemical reaction system with the latter having a larger reaction rate, where $E$ denotes an enzyme. To cover such systems, here we consider the more general case where each reaction vector may correspond to multiple different reactions. For convenience, we introduce the following definition, which plays an important role in the present paper.

Definition 2.1**.**

Two reactions are called equivalent if they have the same reaction vector.

From this definition, two equivalent reactions are either both forward or both backward. Following the notation in [31], let $V(\mathcal{R})=\{\nu^{\prime}_{i}-\nu_{i}|R_{i}\ \text{is a forward reaction}\}$ denote the collection of reaction vectors for forward reactions. Throughout the paper, the elements in $V(\mathcal{R})$ will be listed as $\omega_{1},\cdots,\omega_{r}$, where $1\leq r\leq N$. For any $\omega_{p}=(\omega_{p}^{1},\cdots,\omega_{p}^{d})\in V(\mathcal{R})$, we set

[TABLE]

Then we can relabel the elements in $\mathcal{R}^{+}_{p}$ and $\mathcal{R}^{-}_{p}$ as

[TABLE]

where $r_{p}$ represents the number of forward reactions with the same reaction vector $\omega_{p}$ and we call it the multiplicity of $\omega_{p}$. Here we mainly focus on the case of $r_{p}>1$ for some $1\leq p\leq r$. For any $1\leq p\leq r$ and $1\leq l\leq r_{p}$, we set $\nu_{pl}=(\nu^{1}_{pl},\dots,\nu^{d}_{pl})$ and $\nu^{\prime}_{pl}=(\nu^{\prime 1}_{pl},\dots,\nu^{\prime d}_{pl})$.

We first recall the stochastic model of reaction networks. For each $1\leq j\leq d$, let $N_{j}(t)$ denote the number of molecules of the chemical species $S_{j}$ at time $t$. Then the concentration of $S_{j}$ at time $t$ is given by $X^{V}_{j}(t)=N_{j}(t)/V$, where $V$ is the system size. Let $X^{V}(t)=(X^{V}_{1}(t),\dots,X^{V}_{d}(t))$ denote the concentration process of all chemical species. At the mesoscopic level, the process $\{X^{V}(t):t\geq 0\}$ can be modeled by a continuous-time Markov chain on the $d$-dimensional nonnegative integer lattice

[TABLE]

with transition rate matrix $Q^{V}=(q^{V}_{x,y})$ whose elements are defined as follows: for any $1\leq p\leq r$ and any $x\in E_{V}$,

[TABLE]

where $|\nu|=\sum_{j=1}^{d}\nu_{j}$ is called the order of the complex $\nu=(\nu_{1},\cdots,\nu_{d})$ and we write $x!=\prod_{j=1}^{d}x_{j}!$ for each vector $x=(x_{1},\cdots,x_{d})\in\mathbb{Z}^{d}_{\geq 0}$. Let $\mathbb{P}^{V}_{x}(t)$ denote the probability of observing state $x\in E_{V}$ at time $t$. Then the evolution of the stochastic model is governed by the following chemical master equation:

[TABLE]

We next recall the deterministic model of reaction networks. For each $1\leq j\leq d$, let $x_{j}(t)$ denote the concentration of the chemical species $S_{j}$ at time $t$. At the macroscopic level, the concentration process $x(t)=(x_{1}(t),\dots,x_{d}(t))$ of all chemical species can be modeled by the following ordinary differential equation with mass action kinetics:

[TABLE]

where $x_{0}$ is the initial concentration vector and

[TABLE]

where we write $x^{y}=\prod_{j=1}^{d}x_{j}^{y_{j}}$ for any vectors $x,y\in\mathbb{R}^{d}_{\geq 0}$. The relationship between the mesoscopic stochastic model and the macroscopic deterministic model is revealed by the following celebrated Kurtz theorem [4, 5]: For any $\delta,T>0$, whenever $x_{0}^{V}\in E_{V}$ and $x_{0}^{V}\rightarrow x_{0}\in\mathbb{R}^{d}_{\geq 0}$, then

[TABLE]

where $\mathbb{P}_{x_{0}^{V}}(\cdot)=\mathbb{P}(\cdot|X^{V}(0)=x_{0}^{V})$ and $\|x\|$ denotes Euclidean norm of $x\in\mathbb{R}^{d}$. This implies that as the system size tends to infinity, the trajectories of the stochastic model will converge to those of the deterministic model on any compact time interval, whenever the initial value converges.

The limit theorem in (3) can be viewed as the law of large numbers for the stochastic model. The corresponding large deviation principle was proved recently by Agazzi et al. [23, Theorem 1.6] and is stated as follows. The Hamiltonian of a stochastic reaction network is defined as

[TABLE]

where $x\cdot y=\sum_{j=1}^{d}x_{j}y_{j}$ denotes the usual scalar product on $\mathbb{R}^{d}$. The Lagrangian of a stochastic reaction network is then defined as the Legendre-Fenchel transform of the Hamiltonian with respect to the variable $\theta$, namely

[TABLE]

The Lagrangian is nonnegative because $L(x,y)\geq 0\cdot y-H(x,0)=0$. Moreover, it is not hard to prove that $L(x,y)=\infty$ for any $y\notin\mathrm{span}(V(\mathcal{R}))$. This is because any $y\notin\mathrm{span}(V(\mathcal{R}))$ can be decomposed uniquely as $y=y_{1}+y_{2}$, where $y_{1}\in\mathrm{span}(V(\mathcal{R}))$ and $0\neq y_{2}\in\mathrm{span}(V(\mathcal{R}))^{\perp}$. Thus for any $K>0$, we have $L(x,y)\geq Ky_{2}\cdot y-H(x,Ky_{2})=K\|y_{2}\|^{2}$, where we have used the fact that $H(x,\theta)=0$ for any $\theta\in\mathrm{span}(V(\mathcal{R}))^{\perp}$. Since $K$ is arbitrarily chosen, we conclude that $L(x,y)=\infty$.

To proceed, let $D_{[0,T]}(\mathbb{R}^{d}_{\geq 0})$ denote the space of càdlàg functions $\phi:[0,T]\rightarrow\mathbb{R}^{d}_{\geq 0}$ equipped with the topology of uniform convergence. For any $x_{0}\in\mathbb{R}^{d}_{\geq 0}$, let $I_{x_{0},T}:D_{[0,T]}(\mathbb{R}^{d}_{\geq 0})\rightarrow[0,\infty]$ be the function defined as

[TABLE]

Using the properties of the Lagrangian, it is easy to see that $I_{x_{0},T}(\phi)=\infty$ if there exists $0\leq t\leq T$ such that $\phi(t)\notin x_{0}+\mathrm{span}(V(\mathcal{R}))$. With these notation, Agazzi et al. [23, Theorem 1.6] proved the following result: provided that the network is asiphonic and strongly endotactic (see [23, Definitions 1.8 and 1.9] for detailed definitions), for any $x^{V}_{0}\in E_{V}$ and $x^{V}_{0}\rightarrow x_{0}$, the law of the process $\{X^{V}(t):t\in[0,T]\}$ with $X^{V}(0)=x^{V}_{0}$ satisfies a large deviation principle with rate $V$ and good rate function $I_{x_{0},T}$. The large deviation principle means that for any measurable set $\Gamma\subset D_{[0,T]}(\mathbb{R}^{d}_{\geq 0})$, we have

[TABLE]

where $\Gamma^{o}$ and $\bar{\Gamma}$ denote the interior and closure of $\Gamma$, respectively. Combining (3) and (6), it is easy to see that $I_{x_{0},T}(x)=0$, where $x=x(t)$ is the solution of the deterministic model (1). The rate function can be used to define the following quasi-potential:

[TABLE]

Intuitively, $W(x_{0},y)$ represents the “cost” for the stochastic reaction network to move from $x_{0}$ to $y$. It is easy to see that the quasi-potential is nonnegative and jointly continuous in $x_{0}$ and $y$ [24]. Using the properties of the Lagrangian, it is easy to see that $W(x_{0},y)=\infty$ if $y\notin x_{0}+\mathrm{span}(V(\mathcal{R}))$.

Agazzi et al. [23, Theorem 1.15] also deals with the Freidlin-Wentzell-type metastability theory for chemical reaction networks, where the quasi-potential plays a central role. For simplicity, we consider the case where the domain under consideration contains only one stable equilibrium point. Specifically, we assume that the following four conditions are satisfied:

(a) Let $D$ be a bounded open domain in $\mathbb{R}^{d}_{\geq 0}$ with a piecewise $C^{2}$ boundary $\partial D$.

(b) Let $c\in D$ be an asymptotically stable equilibrium point of the deterministic model (1).

(c) The set $\bar{D}=D\cup\partial D$ is attracted to $c$, which means that whenever $x_{0}\in\bar{D}$, the solution of the deterministic model (1) starting from $x_{0}$ satisfies $x(t)\in D$ for each $t>0$ and $x(t)\rightarrow c$ as $t\rightarrow\infty$.

(d) There exists a ball $B\subset\bar{D}$ such that for any $x\in B$ and $y\in\bar{D}$, the set $\bar{D}$ contains the line segment between $x$ and $y$.

It is easy to check that Assumptions A.3 and A.4 in [23] are satisfied under these conditions. Then the Kurtz theorem implies that when $V$ is sufficiently large, the trajectory of the stochastic model will stay in the domain $D$ over any compact time interval with overwhelming probability. However, it is still possible for the system to escape from $D$. The mean exit time from $D$ has the following asymptotic behavior:

[TABLE]

where $\tau_{V}=\inf\{t\geq 0:X^{V}(t)\notin D\}$ denotes the exit time of $X^{V}$ from $D$. Moreover, if there is a unique $y_{0}\in\partial D$ such that

[TABLE]

then for any $\delta>0$, the exit position from $D$ has the following asymptotic behavior:

[TABLE]

and for any $\delta>0$ and $z_{0}\in\partial D$, the exit distribution from $D$ has the following asymptotic behavior:

[TABLE]

Intuitively, when $V$ is sufficiently large, the stochastic model will escape from $D$ around a particular point $y_{0}\in\partial D$ at which the quasi-potential restricted to $\partial D$ attains its minimum.

3 Detailed balance for chemical reaction networks

In this section, we investigate the relationship among different types of detailed balance conditions for chemical reaction networks. Before stating our results, we first recall the definitions of deterministic and stochastic detailed balance for chemical reaction networks [31].

Definition 3.1**.**

We say that a reaction network satisfies deterministic detailed balance (or reaction network detailed balance [31]) if there exists $c\in\mathbb{R}^{d}_{>0}$ such that

[TABLE]

Here $c$ is called a chemical equilibrium state of the reaction network.

Clearly, any chemical equilibrium state $c$ is also an equilibrium point of the deterministic model (1) and thus it is also called a detailed balanced equilibrium point. It has been shown that for mass action kinetics, if one positive equilibrium point of the deterministic model is detailed balanced, then every positive equilibrium point is detailed balanced [31, 32].

Definition 3.2**.**

We say that a reaction network satisfies stochastic detailed balance (or Markov chain detailed balance [31]) if for any $V>0$, there exists a probability measure $\pi^{V}=(\pi^{V}_{x})$ on $E_{V}$ such that

[TABLE]

Note that here we do not require that $\pi^{V}$ is a probability distribution.

A simple method of verifying stochastic detailed balance is to use the Kolmogorov criterion [33], which states that a reaction network satisfies stochastic detailed balance if and only if for any $V>0$, the transition rates satisfy the following Kolmogorov cycle condition:

[TABLE]

for any finite number of states $x_{1},\dots,x_{n}\in E_{V}$. In other words, the Kolmogorov criterion states that a reaction network satisfies stochastic detailed balance if and only if for any $V>0$, the product of the transition rates of the stochastic model along any cycle is equal to that along the reversed cycle.

Besides deterministic and stochastic detailed balance, we introduce another type of detailed balance which is defined as follows. This new type of detailed balance will play an important role in constructing the global potential of a chemical reaction network.

Definition 3.3**.**

(i) We say that a reaction network satisfies zero-order local detailed balance if for any integers $\xi_{1},\dots,\xi_{r}$ satisfying $\sum_{p=1}^{r}\xi_{p}\omega_{p}=0$, we have

[TABLE]

where $f^{+}_{p}(x)$ and $f^{-}_{p}(x)$ are the functions defined in (2).

(ii) We say that a reaction network satisfies first-order local detailed balance if for any $1\leq p,q\leq r$ with $p\neq q$, we have

[TABLE]

(iii) We say that a reaction network satisfies local detailed balance if it satisfies both zero-order and first-order local detailed balance.

Remark 3.4**.**

The ideas behind the above definition are explained as follows. For any integers $\xi_{1},\dots,\xi_{r}$ satisfying $\sum_{p=1}^{r}\xi_{p}\omega_{p}=0$, we can construct a cycle $C$ in the integer lattice $\mathbb{Z}^{d}$, which is given by

[TABLE]

where $\mathrm{sgn}(x)$ is the sign function which takes the value of $1$ if $x>0$, takes the value of [math] if $x=0$, and takes the value of $-1$ if $x<0$. Obviously, for any $x^{V}\in E_{V}$ and $x^{V}\rightarrow x\in\mathbb{R}^{d}_{>0}$, the cycle $C$ can induce a cycle $x^{V}+C/V$ in $E_{V}$ around $x$, which becomes increasingly smaller as $V$ increases. For convenience, let $\eta_{1}\rightarrow\eta_{2}\rightarrow\dots\rightarrow\eta_{L}\rightarrow\eta_{1}$ denote the induced cycle in $E_{V}$, where $L=\sum_{p=1}^{r}|\xi_{p}|$ is the number of transitions in the cycle. If a reaction network satisfies stochastic detailed balance, then it follows from Kolmogorov’s cycle condition that

[TABLE]

Note that the left-hand side of this equality is a function of $1/V$. Since $f(1/V)=0$ for all $V>0$, we have

[TABLE]

and

[TABLE]

Roughly speaking, the condition (9) extracts the zero-order information of the equality (11) as $V\rightarrow\infty$ and the condition (10) extracts the first-order information of the equality (11) as $V\rightarrow\infty$. Since the induced cycle becomes increasingly smaller as $V$ increases, (9) and (10) actually contain the zero-order and first-order local information of detailed balance around $x\in\mathbb{R}^{d}_{>0}$, respectively.

It is a well-known result that deterministic detailed balance implies stochastic detailed balance for a chemical reaction network [31, Theorem 5.9]. The following theorem reveals the relationship between stochastic and local detailed balance.

Theorem 3.5**.**

If a reaction network satisfies stochastic detailed balance, then it also satisfies local detailed balance. In other words, stochastic detailed balance implies local detailed balance.

Proof.

The proof of the theorem will be given in Section 5. ∎

The next corollary follows from Theorem 3.5 and [31, Theorem 5.9] immediately.

Corollary 3.6**.**

For a chemical reaction network, the following statements hold:

(a) Deterministic detailed balance implies stochastic detailed balance.

(b) Stochastic detailed balance implies local detailed balance.

(c) Local detailed balance implies zero-order local detailed balance.

The above corollary reveals the inclusion relationship among four types of detailed balance: deterministic, stochastic, local, and zero-order local detailed balance, as illustrated in Fig. 1. Deterministic detailed balance is the strongest and zero-order local detailed balance is the weakest. The following proposition reveals when the four types of detailed balance are equivalent.

Proposition 3.7**.**

If a chemical network has no equivalent reactions, then the following statements are equivalent:

(a) The network satisfies deterministic detailed balance.

(b) The network satisfies stochastic detailed balance.

(c) The network satisfies local detailed balance.

(d) The network satisfies zero-order local detailed balance.

Proof.

By Corollary 3.6, we only need to prove that (d) implies (a). If the network satisfies zero-order local detailed balance, for any integers $\xi_{1},\cdots,\xi_{r}$ satisfying $\sum_{p=1}^{r}\xi_{p}\omega_{p}=0$, we have

[TABLE]

Since the network has no equivalent reactions, we have $r_{p}=1$ for any $1\leq p\leq r$. This shows that

[TABLE]

where $\log x=(\log x_{1},\cdots,\log x_{d})$. Combining the above two equations shows that

[TABLE]

Thus we conclude that for any integers $\xi_{1},\cdots,\xi_{r}$ satisfying $\sum_{p=1}^{r}\xi_{p}(\nu_{p1}^{\prime}-\nu_{p1})=0$, we have

[TABLE]

This is exactly the so-called Wegscheider cycle condition, which is widely known as the sufficient and necessary condition for deterministic detailed balance [34, Proposition 1]. Therefore, we have proved that (d) implies (a). ∎

We have seen that if a chemical network has no equivalent reactions, then the four types of detailed balance are equivalent. For chemical networks having equivalent reactions, however, the four types of detailed balance are no longer equivalent, as can be seen from the following three counterexamples.

The following example [35, 31] gives a reaction network that satisfies stochastic detailed balance but violates deterministic detailed balance.

Example 3.8**.**

Consider the following well-known Schlögl model [35]:

[TABLE]

The stochastic model of this reaction network is a one-dimensional birth-death process and thus must satisfy stochastic detailed balance. Moreover, it is easy to check that deterministic detailed balance is satisfied if and only if $k^{+}_{11}/k^{-}_{11}=k^{+}_{12}/k^{-}_{12}$ [35]. In other words, if $k^{+}_{11}/k^{-}_{11}\neq k^{+}_{12}/k^{-}_{12}$, then deterministic detailed balance is violated.

The following example gives a reaction network that satisfies local detailed balance but violates stochastic detailed balance.

Example 3.9**.**

Consider the following chemical reaction system:

[TABLE]

By definition, the forward reactions are given by

[TABLE]

and the backward reactions are given by

[TABLE]

It is easy to see that the first two forward reactions have the same reaction vector $\omega_{1}=1$ and the last three forward reactions also have the same reaction vector $\omega_{2}=2$. The multiplicities of the two reaction vectors are given by $r_{1}=2$ and $r_{2}=3$, respectively.

We first prove that the system satisfies local detailed balance. Clearly, the two reaction vectors are linearly related by $\xi_{1}\omega_{1}+\xi_{2}\omega_{2}=0$ with $\xi_{1}=2$ and $\xi_{2}=-1$. It is easy to check that

[TABLE]

Therefore, we have

[TABLE]

which shows that zero-order local detailed balance is satisfied. Moreover, we have

[TABLE]

which shows that first-order local detailed balance is also satisfied.

We next prove that the system violates stochastic detailed balance. For each $V>0$, consider the following cycle in $E_{V}$:

[TABLE]

The transition rates along this cycle and its reversed cycle are given by

[TABLE]

Direct computation shows that

[TABLE]

It is easy to check that the left-hand side of this equation is equal to 1 if and only if $k_{1}^{+}/k^{-}_{1}=k^{+}_{2}/k^{-}_{2}$, which means that stochastic detailed balance is violated if $k_{1}^{+}/k^{-}_{1}\neq k^{+}_{2}/k^{-}_{2}$.

The following example gives a reaction network that satisfies zero-order local detailed balance but violates local detailed balance.

Example 3.10**.**

Consider the following chemical reaction system:

[TABLE]

By definition, the forward reactions are given by

[TABLE]

and the backward reactions are given by

[TABLE]

The first two forward reactions have the same reaction vector $\omega_{1}=(1,0)$, the next two forward reactions have the same reaction vector $\omega_{2}=(0,1)$, and the reaction vector of the last forward reaction is given by $\omega_{3}=(1,1)$. The multiplicities of the three reaction vectors are given by $r_{1}=r_{2}=2$ and $r_{3}=1$, respectively.

We first prove that the system satisfies zero-order local detailed balance. Clearly, the three reaction vectors are linearly related by $\xi_{1}\omega_{1}+\xi_{2}\omega_{2}+\xi_{3}\omega_{3}=0$ with $\xi_{1}=\xi_{2}=1$ and $\xi_{3}=-1$. It is easy to check that

[TABLE]

Therefore, we have

[TABLE]

which shows that zero-order local detailed balance is satisfied. On the other hand, it is easy to check that

[TABLE]

Clearly, the left-hand sides of the above two equations are equal if and only if $k_{1}^{+}/k^{-}_{1}=k^{+}_{2}/k^{-}_{2}$. In other words, if $k_{1}^{+}/k^{-}_{1}\neq k^{+}_{2}/k^{-}_{2}$, then first-order local detailed balance is violated and thus the system does not satisfy local detailed balance.

In [31, Theorem 4.2], Joshi gave the following sufficient and necessary condition for deterministic detailed balance.

Theorem 3.11**.**

[31, Theorem 4.2] A reaction network satisfies deterministic detailed balance if and only if the following two conditions are satisfied:

(a) For any $1\leq p\leq r$, the rate constants of the reactions in $\mathcal{R}^{+}_{p}$ and $\mathcal{R}^{-}_{p}$ satisfy

[TABLE]

(b) For any integers $\xi_{1},\cdots,\xi_{r}$ satisfying $\sum_{p=1}^{r}\xi_{p}\omega_{p}=0$, we have

[TABLE]

While Joshi [31, Theorem 5.14] also gave a sufficient condition for stochastic detailed balance, it is difficult to apply it in practice since an infinite number of restrictions need to be verified. Here we give a simpler sufficient condition for stochastic detailed balance that is more applicable in practice. To state our sufficient condition, we need the following definition.

Definition 3.12**.**

For each $1\leq p\leq r$, we say that the reaction vector $\omega_{p}$ satisfies the orthogonality condition if

[TABLE]

It is easy to see that if $r_{p}=1$, then the orthogonality condition is automatically satisfied for $\omega_{p}$. The following theorem provides a new sufficient condition for stochastic detailed balance. This sufficient condition is imposed directly on rate constants and only a finite number of restrictions need to be verified.

Theorem 3.13**.**

Suppose that the reaction vectors $\omega_{1},\cdots,\omega_{r}$ are linearly independent. Suppose that for each $1\leq p\leq r$, either one of the following two conditions is satisfied:

(a) The reaction vector $\omega_{p}$ satisfies the orthogonality condition.

(b) The rate constants of the reactions in $\mathcal{R}^{+}_{p}$ and $\mathcal{R}^{-}_{p}$ satisfy

[TABLE]

Then the reaction network satisfies stochastic detailed balance.

Proof.

The proof of this theorem will be given in Section 6. ∎

From Theorem 3.11, if the reaction vectors $\omega_{1},\cdots,\omega_{r}$ are linearly independent, then (13) holds trivially and hence the reaction network satisfies deterministic detailed balance if and only if (14) is satisfied for all $1\leq p\leq r$.

We next use Theorem 3.13 to construct more examples of chemical reaction networks that satisfy stochastic detailed balance but violate deterministic detailed balance.

Example 3.14**.**

Consider the following chemical reaction system:

[TABLE]

By definition, the forward reactions are given by

[TABLE]

and the backward reactions are given by

[TABLE]

The first two forward reactions have the same reaction vector $\omega_{1}=(1,0)$ and the last two forward reactions have the same reaction vector $\omega_{2}=(1,-1)$. The multiplicities of the two reaction vectors are given by $r_{1}=r_{2}=2$, respectively, and the two reaction vectors are linearly independent. Moreover, we have

[TABLE]

It is easy to check that $(\nu_{22}^{j}-\nu^{j}_{21})\omega_{1}^{j}=0$ for $j=1,2$. This shows that the orthogonality condition is satisfied for the reaction vector $\omega_{2}$. By Theorem 3.13, the reaction network satisfies stochastic detailed balance if the condition (b) holds for $p=1$, namely $k^{+}_{11}/k^{-}_{11}=k^{+}_{12}/k^{-}_{12}$. Thus, if $k^{+}_{11}/k^{-}_{11}=k^{+}_{12}/k^{-}_{12}$ but $k^{+}_{21}/k^{-}_{21}\neq k^{+}_{22}/k^{-}_{22}$, then the system satisfies stochastic detailed balance but violates deterministic detailed balance.

4 Global potential for chemical reaction networks

In this section, we investigate the quasi-potential of chemical reaction networks, which plays a central role in the Freidlin-Wentzell-type metastability theory. The quasi-potential defined in (7) has two important features: (i) it is locally defined within each basin of attraction and in general cannot be globally defined over the whole space and (ii) it is defined in the variational form which is usually too complicated to be computed explicitly. The following theorem shows that, under the condition of local detailed balance, the quasi-potential not only can be defined globally over the whole space but also has an explicit and concise expression.

In the following, we always assume that the stochastic model $X^{V}$ starts from some $x^{V}_{0}\in E_{V}$ which satisfies $x^{V}_{0}\rightarrow x_{0}\in\mathbb{R}^{d}_{\geq 0}$ as $V\rightarrow\infty$. Recall that the critical points of a function $U\in C^{1}(\mathbb{R}^{d}_{>0})$ are those points in $\mathbb{R}^{d}_{>0}$ at which $\nabla U$ vanish.

Theorem 4.1**.**

Suppose that a reaction network satisfies local detailed balance. Let $\{\omega_{i_{1}},\cdots,\omega_{i_{m}}\}$ be an arbitrary basis of ${\rm span}(V(\mathcal{R}))$ and let $M=(\omega_{i_{1}}^{T},\cdots,\omega^{T}_{i_{m}})$ be a $d\times m$ matrix, where $m$ is the dimension of ${\rm span}(V(\mathcal{R}))$. Let $F:\mathbb{R}^{d}_{>0}\rightarrow\mathbb{R}^{d}$ be a vector field defined as

[TABLE]

Then the following five statements hold:

(a) The definition of the vector field $F$ is independent of the choice of the basis of ${\rm span}(V(\mathcal{R}))$. In addition, for any $1\leq p\leq r$, we have

[TABLE]

(b) The vector field $F$ has a potential function $U\in C^{\infty}(\mathbb{R}^{d}_{>0})$, namely

[TABLE]

(c) The potential function $U$ satisfies

[TABLE]

where $x=x(t)$ is the solution of the deterministic model (1), and the equality holds if and only if the deterministic model starts from any one of its equilibrium points.

(d) The critical points of $U$ within $\mathbb{R}^{d}_{>0}\cap(x_{0}+\mathrm{span}(V(\mathcal{R})))$ are also the equilibrium points of the deterministic model (1).

(e) Let $c\in\mathbb{R}^{d}_{>0}\cap(x_{0}+\mathrm{span}(V(\mathcal{R})))$ be an equilibrium point of the deterministic model (1). If $y\in\mathbb{R}^{d}_{>0}$ is attracted to $c$ for the deterministic model (1), then

[TABLE]

where $W$ is the quasi-potential defined in (7).

Proof.

The proof of this theorem will be given in Section 7. ∎

It is easy to see that if two potential functions $U_{1},U_{2}\in C^{\infty}(\mathbb{R}^{d}_{>0})$ both satisfy (17), namely

[TABLE]

then $U_{1}$ and $U_{2}$ must coincide with each other (up to a constant) on $\mathbb{R}^{d}_{>0}\cap(x_{0}+\mathrm{span}(V(\mathcal{R})))$.

Definition 4.2**.**

The potential function $U$ introduced in Theorem 4.1 is called the global potential of a reaction network.

Remark 4.3**.**

In the classic Freidlin-Wentzell theory for randomly perturbed dynamical systems [36], it was shown that when detailed balance is satisfied, the system has a global potential that can be computed explicitly [36, Chapter 4, Theorem 3.1]. The above theorem is actually the counterpart of this result for stochastic reaction networks. It shows that under the condition of local detailed balance, we are able to construct a global potential of the reaction network, which is exactly the same as the quasi-potential (up to a constant) within each basin of attraction.

Remark 4.4**.**

Combining Theorems 3.5 and 4.1, we can see that stochastic detailed balance ensures the existence of a global potential. From the deterministic perspective, the vector field $F$ defined in (15) can be understood as the chemical force of the deterministic model. Recall that $F$ has a potential function if and only if the line integral of $F$ (the work exerted by the chemical force) along any smooth closed curve is zero. Similarly, from the stochastic perspective, the chemical force of the stochastic model between any two states $x,y\in E_{V}$ is defined as $\log q^{V}_{x,y}/q^{V}_{y,x}$ [37, Theorem 2.5]. Then the Kolmogorov cycle condition guarantees that the work exerted by the chemical force along any cycle is zero, namely

[TABLE]

for any cycle $x_{1}\rightarrow x_{2}\rightarrow\cdots\rightarrow x_{n}\rightarrow x_{1}$. This intuitively explains why stochastic detailed balance implies the existence of a global potential.

The following corollary shows that the global potential $U$ satisfies the time-independent Hamilton-Jacobi equation.

Corollary 4.5**.**

Suppose that a reaction network satisfies local detailed balance. Then we have

[TABLE]

where $H(x,\theta)$ is the Hamiltonian defined in (4) and $F(x)$ is the vector field defined in $\eqref{vectorfield}$. In particular, we have

[TABLE]

where $U(x)$ is the global potential introduced in Theorem 4.1.

Proof.

Since the network satisfies local detailed balance, it follows from Theorem 4.1(a) that

[TABLE]

The rest of the proof follows immediately from Theorem 4.1(b). ∎

If a reaction network satisfies deterministic detailed balance, then any detailed balanced equilibrium point $c\in\mathbb{R}^{d}_{>0}$ of the deterministic model (1) must also be complex balanced (see [30] for the detailed definition of this concept). Every complex balanced equilibrium point $c=(c_{1},\dots,c_{d})\in\mathbb{R}^{d}_{>0}$ can be used to construct a similar potential

[TABLE]

which turns out to be a Lyapunov function of the deterministic model [38, 39]. The following theorem reveals the relationship between the global potential $U$ introduced in Theorem 4.1 and the potential $\tilde{U}$ defined in (19).

Theorem 4.6**.**

Suppose that a reaction network satisfies deterministic detailed balance with detailed balanced equilibrium point $c\in\mathbb{R}^{d}_{>0}$. Let $U$ be the global potential of the reaction network and let $\tilde{U}$ be the Lyapunov function defined in (19). Then $U$ coincides with $\tilde{U}$ (up to a constant) on $\mathbb{R}^{d}_{>0}\cap(x_{0}+{\rm span}(V(\mathcal{R})))$.

Proof.

Since $c\in\mathbb{R}^{d}_{>0}$ is a detailed balanced equilibrium point, it follows from (8) that

[TABLE]

where $\log c=(\log c_{1},\dots,\log c_{d})$. This shows that for any $x\in\mathbb{R}^{d}_{>0}$,

[TABLE]

Moreover, it follows from Theorem 4.1(a),(b) and (12) that for any $x\in\mathbb{R}^{d}_{>0}\cap(x_{0}+{\rm span}(V(\mathcal{R})))$,

[TABLE]

Combining the above two equations yields

[TABLE]

Then the function $H=\tilde{U}-U$ must satisfy

[TABLE]

For any $x,y\in\mathbb{R}^{d}_{>0}\cap(x_{0}+{\rm span}(V(\mathcal{R})))$, let $\phi:[0,1]\rightarrow\mathbb{R}^{d}_{>0}$ be an arbitrary smooth curve satisfying

[TABLE]

Then we have

[TABLE]

where we have used the fact that $\dot{\phi}(\cdot)\in{\rm span}(V(\mathcal{R}))$. This implies the desired result. ∎

From the above theorem, $U$ must coincide with $\tilde{U}$ on $\mathbb{R}^{d}_{>0}\cap(x_{0}+{\rm span}(V(\mathcal{R})))$ if the reaction network satisfies deterministic detailed balance. Moreover, Theorem 4.1(b) shows that $U$ is the potential function of the vector field $F$. However, the following counterexample shows that $\tilde{U}$ is not necessarily the potential function of the vector field $F$, even the reaction network satisfies deterministic detailed balance.

Example 4.7**.**

Consider the following chemical reaction system:

[TABLE]

Clearly, the system satisfies deterministic detailed balance and $V(\mathcal{R}))=\{(1,-1)\}$. It is easy to check that the vector field $F$ is given by

[TABLE]

Suppose that $x_{0}=(1,0)$. Then there is a unique detailed balanced equilibrium point

[TABLE]

Then the Lyapunov function $\tilde{U}$ associated with the equilibrium point $c$ is given by

[TABLE]

This shows that

[TABLE]

It is easy to check that $-\nabla\tilde{U}\neq F$ on $\mathbb{R}^{2}_{>0}\cap(x_{0}+{\rm span}(V(\mathcal{R})))$ unless $(x_{1},x_{2})=c$.

We next give an example showing the application of the results in this section.

Example 4.8**.**

We revisit the chemical reaction system given in Example 3.14. Recall that the system satisfies stochastic detailed balance when $k^{+}_{11}/k^{-}_{11}=k^{+}_{12}/k^{-}_{12}$. Here we assume that this condition is satisfied. Under this condition, it follows from Theorem 3.5 that the system also satisfies local detailed balance. Note that the reaction vectors in $V(\mathcal{R})$ are given by $\omega_{1}=(1,0)$ and $\omega_{2}=(1,-1)$, which are linearly independent. Therefore, the matrix $M$ is given by

[TABLE]

and thus the vector field $F$ is given by

[TABLE]

where we have used the condition $k^{+}_{11}/k^{-}_{11}=k^{+}_{12}/k^{-}_{12}$. It then follows from Theorem 4.1(b) that the system has a global potential which satisfies

[TABLE]

Integrating the above equation gives the following explicit expression of the global potential:

[TABLE]

where $C$ is a constant which can be chosen so that the minimum of $U$ is zero. By Theorem 4.1(d), any critical point $c=(c_{1},c_{2})\in\mathbb{R}^{2}_{>0}$ of the global potential $U$ must satisfy $\nabla U(c)=0$, namely

[TABLE]

Note that $c_{2}$ satisfies a cubic equation, which is capable of having three distinct positive real roots. In this case, the global potential $U$ has two local minimum points $c_{A}$ and $c_{C}$ and one saddle point $c_{B}$, as illustrated in Fig. 2. It is easy to see that $c_{A}$ and $c_{C}$ are stable equilibrium points of the deterministic model (1) and $c_{B}$ is an unstable equilibrium point. Therefore, by applying Theorem 3.5, we have constructed a high-dimensional chemical reaction network that both satisfies stochastic detailed balance and displays multistability, i.e. multiple attractors for the deterministic model.

When $V$ is large, a stochastic reaction network can transition between multiple attractors with low probability events. In analogy to the classic Freidlin-Wentzell theory [36], if $x_{0}$ is in the basin of attraction of $c_{A}$, then for any $\delta>0$, it follows from [23, Theorem 1.15] that the transition time between attractors has the following asymptotic behavior:

[TABLE]

where $\sigma(B(c_{C},\delta))$ denotes the hitting time of the ball centered at $c_{C}$ with radius $\delta$. By Theorem 4.1(e), for any $y\in\mathbb{R}^{2}_{>0}$ staying in the basin of attraction of $c_{A}$, we have

[TABLE]

Taking $y\rightarrow c_{B}$ in the above equation and applying the continuity of the quasi-potential finally yield

[TABLE]

Note that the right-hand side of this equation is independent of $c_{C}$. This is because the trajectory of the stochastic model will be attracted by $c_{C}$ with very fast speed once it has escaped from the basin of attraction of $c_{A}$. Therefore, the hitting time of $B(c_{C},\delta)$ is mainly determined by the exit time from the basin of attraction of $c_{A}$.

Remark 4.9**.**

The reaction networks in Examples 3.8 and 3.14 are called multistable systems since they are capable of producing multiple positive equilibrium points of the deterministic model. So far, many results have been obtained to identify whether a deterministic reaction network admits multiple equilibrium points [40, 41, 42]. Here we mainly focus on the stochastic model and our results can be applied to investigate the stochastic transitions between multiple attractors.

5 Proof of Theorem 3.5

Proof of Theorem 3.5.

By the definition of the transition rates of the stochastic model, it is easy to see that for any $1\leq p\leq r$, $x^{V}\in E_{V}$, and $x^{V}\rightarrow x\in\mathbb{R}^{d}_{>0}$, we have

[TABLE]

This clearly shows that

[TABLE]

To prove that the system satisfies local detailed balance, we first prove that zero-order local detailed balance is satisfied. For any sufficiently large $V$ and any integers $\xi_{1},\dots,\xi_{r}$ satisfying $\sum_{p=1}^{r}\xi_{p}\omega_{p}=0$, we construct the following cycle in $E_{V}$:

[TABLE]

where $\mathrm{sgn}(x)$ is the sign function which takes the value of $1$ if $x>0$, takes the value of [math] if $x=0$, and takes the value of $-1$ if $x<0$. Applying the Kolmogorov cycle condition to this cycle yields

[TABLE]

Taking the limit of $V\rightarrow\infty$ in this equation and applying (21) yield

[TABLE]

Taking logarithms on both sides yields

[TABLE]

which shows that the system satisfies zero-order local detailed balance.

We next prove that first-order local detailed balance is satisfied. To this end, for any sufficiently large $V$ and any $1\leq p,q\leq r$ with $p\neq q$, we consider the following parallelogram cycle in $E_{V}$:

[TABLE]

Applying the Kolmogorov cycle condition to this cycle yields

[TABLE]

Taking logarithms on both sides of this equation yields

[TABLE]

By the mean value theorem, it is not hard to prove that

[TABLE]

Combining the above two equations yields

[TABLE]

which shows that first-order local detailed balance is also satisfied. ∎

6 Proof of Theorem 3.13

To prove that a reaction network satisfies stochastic detailed balance, it suffices to prove that for any $V>0$ and any cycle $\eta_{1}\rightarrow\eta_{2}\rightarrow\cdots\rightarrow\eta_{L}\rightarrow\eta_{1}$ in $E_{V}$, the Kolmogorov cycle condition

[TABLE]

is satisfied. To this end, we need the following lemma.

Lemma 6.1**.**

Under the conditions in Theorem 3.13, for any $x\in E_{V}$ and any integers $\xi_{1},\cdots,\xi_{r}$, whenever $q^{V}_{x,x+\frac{\omega_{p}}{V}}>0$ and $q^{V}_{\tilde{x},\tilde{x}+\frac{\omega_{p}}{V}}>0$ for some $1\leq p\leq r$, we have

[TABLE]

where

[TABLE]

Proof.

We first discuss the relationship between the reactions that can occur at $x$ and the reactions that can occur at $\tilde{x}$. For each $1\leq p\leq r$, let

[TABLE]

be the family of reactions that belong to $\mathcal{R}^{+}_{p}$ and can occur at $x$. We claim that if $\omega_{p}$ satisfies the orthogonality condition, $\mathcal{R}^{+}_{p}(x)\neq\varnothing$, and $\mathcal{R}^{+}_{p}(\tilde{x})\neq\varnothing$, then $\mathcal{R}^{+}_{p}(x)=\mathcal{R}^{+}_{p}(\tilde{x})$.

To prove this, set

[TABLE]

Since $\mathcal{R}^{+}_{p}(x)\neq\varnothing$, there exists $1\leq l_{1}\leq r_{p}$ such that $Vx_{j}\geq\nu^{j}_{pl_{1}}$ for all $1\leq j\leq d$. It then follows from the definition of $J_{p}$ that $Vx_{j}\geq\nu^{j}_{pl}$ for all $j\in J_{p}$ and $1\leq l\leq r_{p}$. Similarly, since $\mathcal{R}^{+}_{p}(\tilde{x})\neq\varnothing$, we conclude that $V\tilde{x}_{j}\geq\nu^{j}_{pl}$ for all $j\in J_{p}$ and $1\leq l\leq r_{p}$. On the other hand, if $\omega_{p}$ satisfies the orthogonality condition, then $\omega^{j}_{q}=0$ for all $j\notin J_{p}$ and $q\neq p$. This shows that

[TABLE]

for all $j\notin J_{p}$. Therefore, for all $1\leq j\leq d$ and $1\leq l\leq r_{p}$, we have proved that $Vx_{j}\geq\nu^{j}_{pl}$ holds if and only if $V\tilde{x}_{j}\geq\nu^{j}_{pl}$ holds. This clearly shows that $\mathcal{R}^{+}_{p}(x)=\mathcal{R}^{+}_{p}(\tilde{x})$.

We are now in a position to prove (22). For each $1\leq p\leq r$, note that

[TABLE]

where we have used the fact that $\omega_{p}-\nu^{\prime}_{pl}=-\nu_{pl}$ for any $1\leq p\leq r$ and $1\leq l\leq r_{p}$. Since $q^{V}_{x,x+\frac{\omega_{p}}{V}}>0$ and $q^{V}_{\tilde{x},\tilde{x}+\frac{\omega_{p}}{V}}>0$, we have $\mathcal{R}^{+}_{p}(x)\neq\varnothing$ and $\mathcal{R}^{+}_{p}(\tilde{x})\neq\varnothing$. This shows that $\mathcal{R}^{+}_{p}(x)=\mathcal{R}^{+}_{p}(\tilde{x})$. Similarly, we have

[TABLE]

To prove (22), we only need to prove

[TABLE]

Since $\nu^{\prime}_{pl}=\nu_{pl}+\omega_{p}$ for any $1\leq p\leq r$ and $1\leq l\leq r_{p}$, we only need to prove

[TABLE]

If $\omega_{p}$ does not satisfy the orthogonality condition, then (14) must hold. In this case, it is easy to see that

[TABLE]

If $\omega_{p}$ satisfies the orthogonality condition, then we have

[TABLE]

where the second and third equalities in (24) follow from the fact that $\nu^{j}_{pl_{1}}=\nu^{j}_{pl_{2}}$ for any $1\leq l_{1},l_{2}\leq r_{p}$ and $j\in J_{p}$ and the fact that $\tilde{x}_{j}=x_{j}$ for any $j\notin J_{p}$. Therefore, we have proved (23). This completes the proof. ∎

We are now in a position to prove Theorem 3.13.

Proof of Theorem 3.13.

Let $\eta_{1}\rightarrow\eta_{2}\rightarrow\cdots\rightarrow\eta_{L}\rightarrow\eta_{1}$ be an arbitrary cycle in $E_{V}$. We first prove that there is a one-to-one correspondence between the transitions in the cycle. Since $\omega_{1},\cdots,\omega_{r}$ are linearly independent, for each $1\leq p\leq r$, the number of transitions resulting from the reaction vector $\omega_{p}$ in the cycle must be equal to the number of transitions resulting from the reaction vector $-\omega_{p}$. Otherwise, the reaction vector $\omega_{p}$ can be linearly expressed by other reaction vectors in $V(\mathcal{R})$. This contradicts the fact that the elements in $V(\mathcal{R})$ are linearly independent.

We then pair the transitions resulting from the reaction vector $\omega_{p}$ with the transitions resulting from the reaction vector $-\omega_{p}$ in the following manner. First, we project the cycle onto the one-dimensional line spanned by $\omega_{p}$, as illustrated in Fig. 3. Note that the projection mentioned here means oblique projection rather than orthogonal projection (suppose that $\omega_{1},\cdots,\omega_{r}$ are linearly independent vectors; if a vector $\omega$ can be linearly expressed by $\omega_{1},\cdots,\omega_{r}$ as $\omega=a_{1}\omega_{1}+\cdots+a_{r}\omega_{r}$, then the (oblique) projection of $\omega$ onto the direction of $\omega_{p}$ is simply defined as $a_{p}\omega_{p}$). The image of the projection onto the one-dimensional line is still a cycle and only the transitions resulting from $\pm\omega_{p}$ exist in the projected cycle. Second, we decompose the transitions in the projected cycle into multiple floors and we pair each transition with one of its reversed transitions located on the same floor, as depicted in Fig. 3. This is always possible since the number of transitions resulting from $\omega_{p}$ on each floor must be equal to that resulting from $-\omega_{p}$. Based on this method, each transition

[TABLE]

can be paired with a transition

[TABLE]

for some integers $\xi_{1},\cdots,\xi_{r}$. Obviously, this method establishes a one-to-one correspondence between the transitions in the cycle. Applying Lemma 6.1 to the paired transitions yields

[TABLE]

This shows that

[TABLE]

Due to the one-to-one correspondence between the transitions in the cycle, we finally obtain

[TABLE]

Thus we have proved that the Kolmogorov cycle condition is satisfied for each cycle, which implies that the system satisfies stochastic detailed balance.

∎

7 Proof of Theorem 4.1

To prove Theorem 4.1, we need some lemmas. The following lemma is exactly Theorem 4.1(a).

Lemma 7.1**.**

Suppose that a reaction network satisfies zero-order local detailed balance. Then the vector field $F$ defined in (15) is independent of the choice of the basis of ${\rm span}(V(\mathcal{R}))$. Moreover, for any $1\leq p\leq r$, we have

[TABLE]

Proof of Lemma 7.1.

Let $\{\omega_{i_{1}},\dots,\omega_{i_{m}}\}$ and $\{\omega_{j_{1}},\dots,\omega_{j_{m}}\}$ be two arbitrary bases of ${\rm span}(V(\mathcal{R}))$ and let $M=(\omega^{T}_{i_{1}},\dots,\omega^{T}_{i_{m}})$ and $\tilde{M}=(\omega^{T}_{j_{1}},\dots,\omega^{T}_{j_{m}})$ be two matrices. Since the columns of $M$, as well as the columns of $\tilde{M}$, are linearly independent, there exists an invertible matrix $A=(a_{pk})\in M_{m\times m}(\mathbb{\mathbb{R}})$ such that $\tilde{M}=MA$, where $M_{m\times m}(\mathbb{\mathbb{R}})$ denotes the set of $m\times m$ matrices whose components are all real numbers. Since $\mathrm{rank}(M)=m$, the matrix $M$ must have an invertible $m\times m$ submatrix. Since the entries of $M$ and $\tilde{M}$ are all rational numbers, using Cramer’s rule of computing the inverse matrix, it is easy to see that the entries of $A$ are all rational numbers. Note that each column of $\tilde{M}$ can be linearly expressed by the columns of $M$ as

[TABLE]

Since the system satisfies zero-order local detailed balance and since $a_{pk}\in\mathbb{Q}$ for all $1\leq p,k\leq m$, we obtain

[TABLE]

Since $\mathrm{rank}(A)=\mathrm{rank}(A^{T}A)$ for an arbitrary real matrix $A$ [43, Section 3.5], it is easy to see that the matrices $M^{T}M$ and $\tilde{M}^{T}\tilde{M}$ are both invertible. Thus we obtain

[TABLE]

This implies that the definition of the vector field $F$ is independent of the choice of the basis of ${\rm span}(V(\mathcal{R}))$.

On the other hand, since $\{\omega_{i_{1}},\cdots,\omega_{i_{m}}\}$ is a basis of ${\rm span}(V(\mathcal{R}))$, for each $1\leq p\leq r$, the reaction vector $\omega_{p}$ can be linearly expressed by $\omega_{i_{1}},\cdots,\omega_{i_{m}}$ as

[TABLE]

where $a=(a_{1},\cdots,a_{m})\in\mathbb{R}^{m}$. This shows that

[TABLE]

Since the system satisfies zero-order local detailed balance and since the entries of $(M^{T}M)^{-1}M^{T}\omega_{p}^{T}$ are all rational numbers, we immediately obtain

[TABLE]

This completes the proof. ∎

The following lemma is exactly Theorem 4.1(b).

Lemma 7.2**.**

Suppose that a reaction network satisfies local detailed balance. Then the vector field $F$ has a potential function $U\in C^{\infty}(\mathbb{R}^{d}_{>0})$, namely

[TABLE]

Proof.

For any $x\in x_{0}+{\rm span}(V(\mathcal{R}))$, there exists $y=(y_{1},\cdots,y_{m})\in\mathbb{R}^{m}$ such that

[TABLE]

Since $\{\omega_{i_{1}},\cdots,\omega_{i_{m}}\}$ is a basis of ${\rm span}(V(\mathcal{R}))$, this equation establishes a one-to-one correspondence between $x\in\mathbb{R}^{d}_{>0}\cap(x_{0}+{\rm span}(V(\mathcal{R})))$ and $y\in E$, where $E$ is some convex open subset of $\mathbb{R}^{m}$. We denote this correspondence by

[TABLE]

The inverse of this affine transformation is then given by

[TABLE]

For convenience, we introduce a function $h:E\rightarrow\mathbb{R}^{m}$ as

[TABLE]

It is then easy to check that

[TABLE]

Since the system satisfies first-order local detailed balance, we immediately obtain

[TABLE]

To proceed, we construct a smooth differential 1-form on $E$ as

[TABLE]

It then follows from (26) that

[TABLE]

This shows that the differential form $\omega$ is closed. Since $E$ is a convex open subset of $\mathbb{R}^{m}$, it then follows from the Poincaré lemma [44, Theorem 17.14] that the $k$th de Rham cohomology of $E$ is vanishing for each $k\geq 1$, which means that $\omega$ is exact. In other words, there exists a function $\tilde{U}\in C^{\infty}(E)$ such that

[TABLE]

This clearly shows that $h=-\nabla\tilde{U}$. Since $\tilde{U}$ is a smooth function on the open set $E\subset\mathbb{R}^{m}$, we can extend it to a function $\tilde{U}\in C^{\infty}(\mathbb{R}^{m})$. To proceed, we define the potential function $U\in C^{\infty}(\mathbb{R}^{d}_{>0})$ as

[TABLE]

For any $x\in\mathbb{R}^{d}_{>0}\cap(x_{0}+{\rm span}(V(\mathcal{R})))$, straightforward computations show that

[TABLE]

This completes the proof. ∎

The following lemma is exactly Theorem 4.1(c),(d).

Lemma 7.3**.**

Suppose that a reaction network satisfies local detailed balance. Then

[TABLE]

where $x=x(t)$ is the solution of the deterministic model (1), and the equality holds if and only if the deterministic model starts from any one of its equilibrium points. Moreover, the critical points of $U$ within $\mathbb{R}^{d}_{>0}\cap(x_{0}+\mathrm{span}(V(\mathcal{R})))$ are also the equilibrium points of the deterministic model.

Proof.

Let $x=x(t)$ be the solution of the deterministic model (1). It then follows from Lemma 7.1 that

[TABLE]

where the equality holds for some $t\geq 0$ if and only if $f_{p}^{+}(x(t))=f^{-}_{p}(x(t))$ for all $1\leq p\leq r$, which implies that $x(t)$ is an equilibrium point of the deterministic model (1). If the deterministic model does not start from any one of its equilibrium points, then $x(t)$ is not an equilibrium point and thus the equality in (27) cannot be attained. Finally, if $c\in\mathbb{R}^{d}_{>0}\cap(x_{0}+\mathrm{span}(V(\mathcal{R})))$ is a critical point of $U$, then we have $F(c)=-\nabla U(c)=0$. It then follows from (16) that $f^{+}_{p}(c)=f^{-}_{p}(c)$ for each $1\leq p\leq r$, which shows that $c$ is an equilibrium point of the deterministic model. ∎

The above lemma shows that the asymptotic stability of equilibrium points of the deterministic model (1) can be analyzed with the aid of the potential function $U$, according to the classic Lyapunov stability criterion [45, Chapter 30].

Lemma 7.4**.**

Let $\{\omega_{i_{1}},\cdots,\omega_{i_{m}}\}$ be an arbitrary basis of ${\rm span}(V(\mathcal{R}))$ and let $M=(\omega_{i_{1}}^{T},\cdots,\omega^{T}_{i_{m}})$ be a $d\times m$ matrix. Then for any $x\in\mathbb{R}^{d}_{>0}$ and $y\in{\rm span}(V(\mathcal{R}))$, we have

[TABLE]

where $A=M(M^{T}M)^{-1}$ and $L(x,y)$ is the Lagrangian defined in (5). Furthermore, there exists a unique $\theta=\theta_{0}\in\mathbb{R}^{m}$ such that the maximum is attained.

Proof.

Let $\{v_{1},\cdots,v_{d-m}\}$ be an arbitrary basis of ${\rm span}(V(\mathcal{R}))^{\perp}$ and let

[TABLE]

be two matrices. We next prove that the square matrix $B$ is invertible. To this end, we consider the system of linear equations $Bz^{T}=0$. Since $Bz^{T}=0$, we have $M^{T}z^{T}=0$ and $N^{T}z^{T}=0$. Since the columns of $M$ and $N$ constitute a basis of $\mathbb{R}^{d}$, we conclude that $z=0$. Thus the system of linear equations $Bz^{T}=0$ has only the zero solution, which implies that $B$ is invertible.

For any $x\in\mathbb{R}^{d}_{>0}$ and $y\in{\rm span}(V(\mathcal{R}))$, we have

[TABLE]

where we have used the fact that $B$ is invertible. Since $y\in{\rm span}(V(\mathcal{R}))$, there exist $k_{1},\cdots,k_{m}\in\mathbb{R}$ such that $y=k_{1}\omega_{i_{1}}+\dots+k_{m}\omega_{i_{m}}$. Since the columns of $M$ and the columns of $N$ are orthogonal, we have

[TABLE]

Similarly, we can prove that the last $d-m$ components of $\omega_{p}B^{T}$ are [math] for all $1\leq p\leq r$. Thus the supremum in the last equality of (28) can be taken over the first $m$ components of $\theta$, namely

[TABLE]

where

[TABLE]

and $A=M(M^{T}M)^{-1}$. Direct computations show that the Hessian matrix of the function $G(x,y,\theta)$ with respect to $\theta$ is given by

[TABLE]

where

[TABLE]

Since $A^{T}(\omega_{i_{1}}^{T},\dots,\omega^{T}_{i_{m}})=(M^{T}M)^{-1}M^{T}M=I$ and $g_{p}(x,\theta)>0$ for any $x\in\mathbb{R}^{d}_{>0}$, it is easy to check that the Hessian matrix is negative definite. Therefore, $G(x,y,\theta)$ is a strictly concave function with respect to $\theta$ [46, Corollary 3.8.6]. We next prove that

[TABLE]

Since $yA=(k_{1},\cdots,k_{m})$, for any $\theta=(\theta_{1},\cdots,\theta_{m})$, we have

[TABLE]

Therefore, (30) follows directly from the fact that exponential functions grow much faster than linear functions. Since $G(x,y,\theta)$ is a strictly concave function with respect to $\theta$ and since (30) holds, it follows that $G(x,y,\theta)$ must attain its maximum at a unique $\theta=\theta_{0}\in\mathbb{R}^{m}$ [47, Chapter B, Theorem 4.1.1]. ∎

For any absolutely continuous trajectory $\phi:[0,T]\rightarrow\mathbb{R}^{d}_{>0}$ satisfying

[TABLE]

we define

[TABLE]

to be the line integral of the vector field $F$ along the trajectory $\phi$. If the system satisfies local detailed balance, then the vector field $F$ has a potential function $U$. In this case, $S(\phi)$ can be represented by the potential function $U$ as

[TABLE]

which only depends on the endpoints of the trajectory $\phi$ and thus is “path-independent”. Moreover, we define the reversed trajectory of $\phi$ as

[TABLE]

Lemma 7.5**.**

Suppose that a reaction network satisfies local detailed balance. Then for any absolutely continuous trajectory $\phi:[0,T]\rightarrow\mathbb{R}^{d}_{>0}$ satisfying

[TABLE]

we have

[TABLE]

where $\phi^{-}$ is the reversed trajectory of $\phi$.

Proof of Lemma 7.5.

Since $\phi(\cdot)\in x_{0}+\mathrm{span}(V(\mathcal{R}))$, we have $\dot{\phi}(\cdot)\in\mathrm{span}(V(\mathcal{R}))$. It then follows from Lemma 7.4 that for each $0\leq t\leq T$, there exists a unique $\theta=\theta_{1}(t)\in\mathbb{R}^{m}$ such that the following maximum is attained:

[TABLE]

where $G(x,y,\theta)$ is the function defined in (29). Using Lemma 7.4 again shows that for each $0\leq t\leq T$, there exists a unique $\theta=\theta_{2}(t)\in\mathbb{R}^{m}$ such that the following maximum is attained:

[TABLE]

Since the maximum in (32) is attained at $\theta=\theta_{1}(t)$, taking the derivatives of $G(\phi(t),\dot{\phi}(t),\theta)$ with respect to $\theta$ and evaluating at $\theta=\theta_{1}(t)$ yields

[TABLE]

where $A=M(M^{T}M)^{-1}$. By the proof of Lemma 7.4, $G(\phi(t),\dot{\phi}(t),\theta)$ is a strictly concave function with respect to $\theta$. This shows that $\theta=\theta_{1}(t)$ is the unique solution of the following equation

[TABLE]

Similarly, since the maximum in (33) is attained at $\theta=\theta_{2}(t)$, we obtain

[TABLE]

To proceed, let

[TABLE]

It then follows from (16) that

[TABLE]

This clearly shows that

[TABLE]

which implies that

[TABLE]

Substituting this equation into (35), it is easy to check that $\theta=-\theta_{2}(t)-f(\phi(t))$ is also a solution of the equation (34). By the uniqueness of the solution of the equation (34), we immediately obtain

[TABLE]

It thus follows from (32), (33), and (36) that

[TABLE]

For the reversed trajectory $\phi^{-}$, note that

[TABLE]

Finally, we obtain

[TABLE]

This completes the proof. ∎

The following lemma is exactly Theorem 4.1(e).

Lemma 7.6**.**

Suppose that a reaction network satisfies local detailed balance. Let $c\in\mathbb{R}^{d}_{>0}\cap(x_{0}+\mathrm{span}(V(\mathcal{R})))$ be an equilibrium point of the deterministic model (1). If $y\in\mathbb{R}^{d}_{>0}$ is attracted to $c$ for the deterministic model (1), then

[TABLE]

Proof.

Let $\phi:[0,T]\rightarrow\mathbb{R}^{d}_{>0}$ be an arbitrary absolutely continuous trajectory satisfying

[TABLE]

It thus follows from Lemma 7.5 that

[TABLE]

where $\phi^{-}$ is the reversed trajectory of $\phi$. Taking the infimum over $\phi$ on both sides of this equation yields

[TABLE]

On the other hand, let $\phi_{y}(t)$ denote the trajectory of the deterministic model (1) starting from $y$. In addition, let

[TABLE]

denote the reversed trajectory of $\phi_{y}$ over the interval $[0,T]$. Applying Lemma 7.5 again shows that

[TABLE]

Moreover, let

[TABLE]

be an absolutely continuous trajectory from $c$ to $\phi_{y}(T)$. Recall that for any $y\in{\rm span}(V(\mathcal{R}))$, we have

[TABLE]

where $A=M(M^{T}M)^{-1}$. Since $y$ is attracted to $c$, for any $\epsilon>0$, we have $\|\phi_{y}(T)-c\|\leq\epsilon$ when $T$ is sufficiently large. For convenience, set

[TABLE]

For any $\theta=(\theta_{1},\cdots,\theta_{m})\in\mathbb{R}^{m}$, whenever $\|x-c\|\leq\epsilon$, we have

[TABLE]

Thus for any $y=k_{1}\omega_{i_{1}}+\dots+k_{m}\omega_{i_{m}}\in\mathrm{span}(V(\mathcal{R}))$, whenever $\|x-c\|\leq\epsilon$, we have

[TABLE]

where

[TABLE]

It is easy to see that the function $f^{\star}$ is continuous on $\mathrm{span}(V(\mathcal{R}))$. Since $\|\dot{\zeta_{T}}(t)\|\equiv 1$ and $f^{\star}$ is bounded on compact sets, there exists a constant $C_{2}>0$ such that for any $T\geq 0$ and $0\leq t\leq\|\phi_{y}(T)-c\|$,

[TABLE]

Thus when $T$ is sufficiently large, it follows from (38) that for any $0\leq t\leq\|\phi_{y}(T)-c\|$,

[TABLE]

where we have used the fact that $\|\zeta_{T}(t)-c\|\leq\epsilon$ when $T$ is sufficiently large. Finally, we obtain

[TABLE]

Combining $\zeta_{T}$ and $\psi_{T}$, we obtain an absolutely continuous trajectory from $c$ to $y$. Therefore, we have

[TABLE]

Since $y$ is attracted to $c$, taking $T\rightarrow\infty$ in the above equation yields

[TABLE]

where we have used the arbitrariness of $\epsilon$. Finally, the desired result follows from (37) and (39). ∎

Acknowledgement

We thank Prof. Vadim A. Malyshev and Sergey A. Pirogov for stimulating discussions via email, and thank Prof. Hao Ge and Dr. Xiao Jin for providing some useful references. We are also grateful to the anonymous referees for their valuable comments and suggestions which helped us greatly in improving the quality of this paper. C. J. acknowledges support from the NSAF grant in National Natural Science Foundation of China with grant No. U1930402. D.-Q. Jiang is supported by National Natural Science Foundation of China with grant No. 11871079.