Constrained Langevin approximation for the Togashi-Kaneko model of autocatalytic reactions

1.   Introduction

In 2001, Togashi and Kaneko ^[1] introduced a simple model of autocatalytic reactions that displays a peculiar "switching behavior" in some regions of the parameter space. The system state switches between patterns where a few species are abundant and the remaining species are almost absent, demonstrating multi-stability at those patterns; see the right column of Figures 1 and 2 for some sample trajectories. Paraphrasing ^[2], it is believed that "the switching is triggered by a single molecule of a previously extinct species that drives the system to a different pattern through a sequence of quick reactions". The emergence of such multi-stability induced by the small number effect, called discreteness-induced transitions (DITs) in ^[1], has been observed in many complicated models in physics, biology and other scientific fields. For instance, it has been reported in catalytic chemical reactions ^[3,4,5,6], reaction-diffusion systems ^[7,8], gene regulatory networks ^[9,10], cancer tumor evolution ^[11], virus replication ^[12], ecology ^[13].

The widespread nature of DITs has attracted many theoretical studies on the model in ^[1] and its variants. These studies include analysis for the switching time ^[13,14,15], stationary distributions ^[2,16], separation of time scale ^[6,17] and multimodality ^[18,19]. In ^[14,15], the authors analyzed a mass-conserved reaction network, namely the autocatalytic reactions (1.3) below with $d\in\{2, 3\}$ species, together with mutations between species,

The total number $N$ of molecules among all species remains constant in time for this model, making it more amenable to analysis. For $d = 2$, by [14,Eq (7)], the mean time to move from one boundary state to the other is approximately

In ^[15], a similar model with $d = 3$ species (and more reactions) is studied, where a noise-induced reversal of chemical currents was observed as the total number of molecules decreased.

For the model in ^[1], henceforth called the TK model, mass is no longer conserved. The reaction network of the TK model consists of a cycle of autocatalytic reactions, together with inflow and outflow reactions. When there are $d$ species $\{A_i\}_{i = 1}^d$, the reaction network is

Suppose $X^i_t$ represents the number of molecules for $A_i$ at time $t$. It is standard to assume that the vector $X_t = (X^i_t)_{i = 1}^d\in \mathbb{Z}_+^d$ evolves according to a continuous-time Markov chain (CTMC) over time, with transition rates specified by mass-action kinetics. More precisely, we let $X^{d+1} = X^1$ by convention, and construct this stochastic process $X = (X_t)_{t\in\mathbb{R}_+}$ as the solution to the stochastic equation

where $\{e_i\}_{i = 1}^d$ is the standard basis of $\mathbb{R}^d$ and $\{\mathcal{N}^{\kappa}_i, \mathcal{N}^{\lambda}_i, \mathcal{N}^{\delta}_i\}_{i = 1}^d$ are independent Poisson processes with unit rate; see the monograph ^[20] for basic properties of this equation. Note that the sum of all coordinates $\sum_{i = 1}^dX^i$ is equal in distribution to the process $S$ that solves the following equation that no longer depends on $\kappa$:

where $\mathcal{N}_1$ and $\mathcal{N}_2$ are independent Poisson processes with unit rate. Therefore, the total mass $\sum_{i = 1}^dX^i_t$ is an immigration-death process that converges, as $t\to\infty$, to the Poisson distribution with mean $\frac{d\lambda}{\delta}$ exponentially fast with rate $\delta$ (see [21,Chapter 9]).

In [2,Theorem 4.1] it was shown that the CTMC described by (1.5)-(1.6) is positive recurrent and converges exponentially fast to its stationary distribution when $\delta > 0$ and $\lambda > 0$. For the special case when $\delta = \frac{d\, \kappa}{d-1}$, an explicit form of the stationary distribution is known in [2,Theorem 4.3]. Beyond these, not much is proven about (1.5)-(1.6) in general dimensions. Also, no asymptotic formula like (1.2) is known for (1.5)-(1.6), and it is not clear how to directly compare the model described by (1.3)-(1.4) with the model described by (1.3) and (1.1), because the reaction vector at the boundary of $\mathbb{Z}_+^d$ are different for the two models.

A standard model reduction technique is to look at the mean field approximation or the diffusion (Langevin) approximation of (1.5)-(1.6). In these approximations, it is assumed that the initial molecule count for each species is proportional to a parameter $V$ which is Avogadro's number times the volume of the vessel in which the reactions take place. However, neither of these two approximations is a good predictor: ordinary differential equations ODEs do not capture the DIT in (1.5)-(1.6) because the DIT for this case is due to the successive extinction and revival of species, and the mean-field approximation breaks down when the abundance of some species are not high. Central limit theorems and diffusion approximations ^[20,22] can capture the fluctuation around the deterministic ODE for stochastic chemical reaction networks, but they can have negative coordinates, which is unrealistic and poses a technical issue: the diffusion process may not remain well-defined when a coordinate becomes negative.

To address these issues, an obliquely reflected diffusion called constrained Langevin approximation (CLA) was proposed in ^[23,24] as a better diffusion approximation to the CTMC that arises from chemical reaction networks. This obliquely reflected diffusion process has state space in the positive orthant $\mathbb{R}_+^d$ (where $d$ is the number of species) and thus respects the constraint that chemical concentrations are never negative. Intuitively, a CLA behaves like a diffusion process inside the strictly positive orthant $\mathbb{R}^d_{ > 0}$ and reflects instantaneously at a boundary face in the direction specified by a vector field. Special care needs to be taken for reflection at the intersection of two or more faces. It was demonstrated in ^[23] through numerical studies that, in addition to having the correct support, the stationary distribution for CLA can capture the behavior of the CTMC more accurately than the usual diffusion approximation.

Existing analyses of the TK model and its variants are mostly restricted to models with a small number of species, and they do not cover the analysis of the CLA of the corresponding CTMC. In this paper, we analyze the CLA for the TK model (1.5)-(1.6) in general dimensions, and we perform a simulation study for both the CTMC (1.5)-(1.6) and the CLA.

Organization of this paper. In Section 2, we show that the CLA possesses the Feller property, is positive Harris recurrent and converges exponentially fast to the unique stationary distribution $\pi$. We also show that $\pi$ has finite moments and we characterize it in terms of an elliptic partial differential equation. The proofs of these results, heavily based on the Foster Lyapunov function approach, are presented in Section 4. Finally, Section 3 contains our simulation study for the TK model and its CLA in dimensions $d = 2, 3$ and higher. We demonstrate that, at least in dimensions $d = 2$ and 3, the CLA can capture both the stationary distribution and the expected transition time of the TK model when $V$ is large enough and is proportional to the initial molecule count. For $d = 2$, the switching time distribution of the CLA is less accurate when $V$ is small. Stochastic simulations for some trajectories and the stationary distribution of the TK model in dimension 6 has been carried out and are discussed in Section 3.3.

2.   Analytical results for the constrained Langevin approximation

In this section, we describe and analyze a CLA to the TK model (1.5)-(1.6). Precisely, the CLA is the strong solution to (2.2). In the three subsections below, we first establish the wellposedness of Eq (2.2) and the Feller property of the CLA. We then show that the CLA is positive Harris recurrent and exponentially ergodic. Finally, we characterize the stationary distribution $\pi$ and show that it has finite moments. The proofs of these results are in Subsection 4.3.

Under the classical scaling, the initial molecule counts are proportional to a scaling parameter $V$ and the rate constants are of order $\kappa = O(V^{-1})$, $\lambda = O(V)$ and $\delta = O(1)$ as $V\to\infty$. So we let

where $\kappa'$, $\lambda'$, and $\delta'$ are constants that will emerge in the mean-field approximation as $V\to\infty$. In chemical reactions, $V$ denotes Avogadro's number times the volume of the vessel in which all reactions take place.

Following the general method in ^[23,24], a CLA for the TK model (1.3) and (1.4) is described by the stochastic differential equation with reflection,

where $W$ is a $d$-dimensional Brownian motion, $b:\mathbb{R}_+^d\to \mathbb{R}_+^d$ and $\gamma:\partial\mathbb{R}_+^d\to \mathbb{R}_+^d$ are functions given by

where $\{e_k\}_{k = 1}^n$ is the standard basis in $\mathbb{R}^d$, $x = (x_k)_{k = 1}^d$ and $|\cdot|$ is the usual Euclidean norm. The function $\sigma$ is the $d\times d$-matrix-valued function on $\mathbb{R}_+^d$ given by $\sigma(x) = \sqrt{\Gamma(x)}$ where, by [23,Eq (29)],

where $e_{i, j}\in \mathbb{R}^{d\times d}$ is the matrix whose $(i, j)$-th entry is one and all other entries are zero. Note that the matrix $\Gamma(x)$ is symmetric and strictly positive definite (or uniformly elliptic) for all $x \in \mathbb{R}_+^d$; see Subsection 4.1 for details.

Remark 1. The square root matrix $\sigma$ of $\Gamma$ ($\sigma = \sqrt{\Gamma}$) in (2.2) is implicit and needs to be calculated in practice. Another, more explicit equation that gives the same process $Z$ in distribution is to use higher dimensional Brownian motion. Namely, we first let $\widetilde{W}$ be a 3-dimensional Brownian motion when $d = 2$, and $\widetilde{W}$ be a $2d$-dimensional Brownian motion for $d\ge 3$. We then replace ${\sigma}(Z^{(V)}_t)\, d{W}_t$ in (2.2) by ${\widetilde{\sigma}}(Z_t)\, d\widetilde{W}_t$, where $\widetilde{\sigma}$ is an explicit $d\times 2d$-matrix-valued function for $d\ge 3$. For example, for $d = 3$, $\widetilde{W}$ is a 6-dimensional Brownian motion and

For $d = 2$, we let

The matrix $\widetilde{\sigma}$ can be obtained from (1.5)-(1.6) as in [23,Section 3], under the classical scaling given by (2.1).

2.1.   Path-wise existence and uniqueness of the CLA

The path-wise existence of obliquely reflected diffusions can fail ^[25]. However, for (2.2), the path-wise existence of a solution ensured by [24,Theorem 6.1] or the argument in [26,Theorem 5.1], since the reflection angles are nice enough.

Following ^[27], we say a Markov process $X$ on $\mathbb{R}_+^d$ is Feller if the map $x \mapsto \mathbb{E}_x[f(X_t)]$ is bounded and continuous for any bounded continuous function $f: \mathbb{R}_+^d \to \mathbb{R}$. We cannot find any existing result that can tell us whether the solution to (2.2) is a Feller process or not. So we give a proof of the Feller property (see Proposition 3 in Section 4).

Theorem 1 (Path-wise solution and Feller property). For each $V\in (0, \infty)$ and initial condition $z_0\in \mathbb{R}_+^d$, there exists a unique path-wise solution to Eq (2.2). The solutions start from different points in $\mathbb{R}_+^d$ forming a family of Feller continuous strong Markov processes in $\mathbb{R}_+^d$.

Let $Z = Z^{(V)}$ be the solution to (2.2) from now on, and omit the superscript $V$ when there is no ambiguity. A solution $Z$ to (2.2) is a good approximation of $\frac{X}{V}$ on any compact time interval as $V\to\infty$ with $Z_0 = \frac{X_0}{V}$. More precisely, the CLA was obtained in ^[24] as the scaling limit of a sequence of jump-diffusion processes that are believed to be good approximations of the stochastic reaction network. These jump-diffusion processes behave like the standard Langevin approximation in the interior of the positive orthant and a rescaled version of the CTMC on the boundary of the orthant. Although a rigorous connection between the CLA and the CTMC is still missing, the simulation results in ^[23,24] demonstrated that the CLA is a remarkably good approximation of the CTMC.

Our proofs depend heavily on the Foster-Lyapunov function approach ^[27]. We consider the function $U:\mathbb{R}^d_+ \to \mathbb{R}_+$ defined by

where $|x|_1 = \sum_{i = 1}^d |x_i|$. Importantly, this function is compatible with the reflection field on the boundary (Lemma 1 in Section 4) and leads to Lyapunov inequalities (Lemma 2 in Section 4), which enables our stability analysis for $Z$ that is needed to establish the Feller property positive Harris recurrence, and exponential ergodicity of $Z$.

2.2.   Positive Harris recurrence and exponential ergodicity of the CLA

In [2,Theorem 4.1] it was shown that the discrete TK model (1.5)-(1.6) is positive Harris recurrent and converges exponentially fast to its stationary distribution. Here we obtain the analogous results for the CLA. Recall from [27,Sections 3 and 4] that $Z$ is called Harris recurrent if there exists a $\sigma$-finite measure $\mu$ on $\mathbb{R}^d_+$ such that whenever $\mu(A) > 0$, we have $\mathbb{P}_x(\tau_A < \infty) = 1$ for all $x\in \mathbb{R}_+^d$, where $\tau_A: = \inf\{t\in\mathbb{R}_+:\, Z_t\in A\}$ is the hitting time of a Borel set $A$. If a process $Z$ is Harris recurrent, then there exists a unique $\sigma$-finite invariant measure [28,Theorem 1]. If, furthermore, the invariant measure is finite, then $Z$ is called positive Harris recurrent.

Theorem 2 (Positive recurrence). The solution $Z$ to (2.2) is positive Harris recurrent and hence has a unique stationary distribution $\pi$. Furthermore, all moments of $\pi$ are finite.

Next, we consider the rate of convergence to stationarity. We say that $Z$ is $f$-exponentially ergodic for a function $f$ if the law of $Z_t$ converges to $\pi$ exponentially fast in the following sense: there exists a constant $\beta < 1$ and a function $B:\mathbb{R}_+^d\to \mathbb{R}_+$ such that

where the $\|\cdot\|_f$-norm is defined as $\|\mu\|_f : = \sup_{g:\, |g| \le f} \left|\int_{\mathbb{R}_+^d} g d\mu\right|$ and the supremum is taken over the space of Borel measurable functions $g$ on $\mathbb{R}_+^d$ with $|g(x)|\leq f(x)$ for all $x\in \mathbb{R}_+^d$.

Theorem 3 (Exponential ergodicity). The solution $Z$ to (2.2) is $f$-exponentially ergodic with $f = U+1$, where $U$ is defined in (2.5).

Theorem 3 implies that (2.6) remains true if we replace $\|\cdot\|_{f}$ by the total variation distance $\|\cdot \|_{\rm TV}$. This is because $\|\mu_1-\mu_2\|_{\rm TV}\leq \|\mu_1-\mu_2\|_{f}$ when $f\geq 1$.

Positive recurrence can fail for reflected diffusions on unbounded domains. For example, it fails for the reflected Brownian motion with a positive drift on $[0, \infty)$. Hence suitable conditions on the state-dependent coefficients and the reflection vector field are needed. In ^[29], the authors consider a reflected diffusion on a convex polyhedral cone $G\subset \mathbb{R}^d$ with the vertex at the origin and the reflection vector $v_i$ is assumed to be constant on each face $G_i$ of the cone. Let $\mathfrak{C}$ be the cone spanned by $\{-v_i\}_{i}$. The main result [29,Theorem 2.2] asserts that the reflected diffusion is positive recurrent and has a unique invariant distribution if there is a bounded set $A\subset G$ such that the vector $b(x)$ is in the cone $\mathfrak{C}$ and uniformly away from the boundary of $\mathfrak{C}$, for all $x\in G\setminus A$.

Unfortunately, the result in ^[29,30] cannot be applied to our CLA directly since the reflection vector field for $Z$ is state-dependent. On other hand, the papers ^[24,26,31] consider state-dependent reflection vector field on non-smooth domains, but these results are not concerned with positive recurrence. We shall prove Theorems 2 and 3 by the Foster-Lyapunov function approach in this paper.

2.3.   Characterization of the stationary distribution of the CLA

Characterization of stationary distributions of a general class of reflected diffusions is given in ^[31]. It was shown that a stationary distribution, should it exist, must satisfy an adjoint linear elliptic partial differential equation with oblique derivative boundary conditions. This equation is called the basic adjoint relationship for reflected Brownian motion in ^[32,33]. Here, in Subsection 4.3, we verify the conditions in [31,Theorems 2 and 3] and apply those results to our CLA.

Following ^[31], we let $C_c^2(\mathbb{R}^d_+)$ be the space of twice continuously differentiable functions on $\mathbb{R}_+^d$ with compact support. We consider the space of functions

where $C_c^2(\mathbb{R}^d_+)\oplus \mathbb{R}$ denotes the space of functions in $C_c^2(\mathbb{R}^d_+)$ plus a constant in $\mathbb{R}$. We also let ${\cal L}$ be the differential operator:

Proposition 1. A probability measure $\pi$ on $\mathbb{R}_+^d$ is a stationary distribution of $Z$ if and only if $\pi(\partial \mathbb{R}_+^d) = 0$ and $\int_{\mathbb{R}_+^d} {\cal L}f(x)\, \pi(dx) \le 0$ for all $f \in {\cal H}$.

A more explicit way to characterize the stationary distribution is through a partial differential equation as in [31,Theorem 3]. For this we recall that the adjoint operator ${\cal L}^*$ of ${\cal L}$ is

Proposition 2. Suppose that there exists a nonnegative integrable function $p \in {\cal C}^2(\mathbb{R}_+^d)$ that satisfies the following three relations:

1) ${\cal L}^*p(x) = 0$ for all $x \in \mathbb{R}_+^d$;

2) For each $i\in\{1, 2, \cdots, d\}$ and $x \in \left\{x_i = 0\right\}\cap \left\{x_j > 0, \, \forall j \in \{1, 2, \cdots, d\}\backslash{\{i\}}\right\}$,

3) For each $1 \le i\neq j \le d$ and $x \in \{x_i = 0\}\cap \{x_j = 0\} \cap \partial \mathbb{R}_+^d$, $p(x) = 0$.

Then the probability measure on $\mathbb{R}_+^d$ defined by

is a stationary distribution for the process $Z$.

The probability measure $\pi$ in Proposition 1 exists and is unique, by Theorem 2. We do not know if the function $p \in {\cal C}^2(\mathbb{R}_+^d)$ in Proposition 2 exists.

3.   Simulation study for the TK model and the CLA

In this section, we present simulation results for the CTMC (1.5)-(1.6) and its associated CLA (2.2). In particular, the dynamical properties of the CLA, including its stationary distribution $\pi$ (guaranteed in Theorem 2) and its finite time trajectory, are compared with those of the CTMC in various dimensions. To begin, in Figure 1 we show some sample trajectories of the CTMC $X$ that solves (1.5)-(1.6) in dimensions $d = 2, 3, 4$, under different choices of parameters. In the left column, trajectories were simulated with the parameters $\lambda = 4, \delta = 1/8$ and $\kappa = 1/32$ in (1.3) and (1.4), where the trajectory oscillates around the fixed line $y = 32$; In comparison, all trajectories in the right column were simulated with the parameters $\lambda = 1/2, \delta = 1/64$ and $\kappa = 1/32$, and all processes tend to spend most of the time on the boundary, i.e., some species are almost extinct while others are abundant. This switching behavior can also be observed in higher dimensions, as one can see in Figure 2 for $d = 5, 6$ under the conditions of $\lambda = 1/2, \delta = 1/64$ and $\kappa = 1/32$.

In Section 3.1, we simulate the CTMC and the CLA for $d = 2$. Our simulations suggest that the stationary distributions of the CTMC are well approximated by those of the CLA. The hitting time distributions of the CTMC are well approximated by those of the CLA when $V$ is large enough; this approximation is less accurate when $V$ is small. In Section 3.2, similar simulations and suggestions are obtained for $d = 3$; there, we also emphasize a discrepancy of the finite trajectory property of TK models between $d = 2$ and $d = 3$. In Section 3.3, we focus on higher-dimensional TK models. For example, we give a description of the switching behaviors between meta-stable patterns of the 6-dimensional CTMC.

Simulation schemes.All CTMCs are simulated by using the Gillespie algorithm ^[34]. Special care needs to be taken when simulating the CLA, when the trajectory is near the boundary. Here all CLAs are simulated via the modified Euler-Maruyama method proposed in ^[35]. An alternative simulation scheme for the CLA may also be developed by a suitable discrete version of the local time ^[36,37].

3.1.   Simulation results for the 2-dimensional TK model

For the dimension $d = 2$, our simulations indicate that the CLA (2.2) nicely captures the stationary distribution of the CTMC (1.5)-(1.6). Furthermore, we consider the time between the extinction events of the two species (called the switching time). We demonstrate that the switching time distribution of the CLA captures that of the CTMC when the volume $V$ is large enough, and this approximation is less accurate when $V$ is small.

For the 2-dimensional CTMC, explicit expression of the stationary distribution was derived in ^[2]. For the classical scaling (2.1) with

the system exhibits different stationary behavior for different choices of $D$: when $D > 2/V$, the stationary distribution is unimodal; whereas when $D = 2/V$, stationary distribution conditioning on the level sets $x+y = n$ is uniform; and when $D < 2/V$, the stationary distribution is heavily concentrated on both boundaries where one species is almost extinct. Such behavior can be visualized in the first row of Figure 3, where the stationary distributions are plotted for different choices of $D$.

In the second row of Figure 3, stationary densities, obtained via time averaging over long-time trajectories of the CLA (2.2), are plotted for the same parameters as the associated CTMC. In all three cases, the CLA in (2.2) accurately captures the stationary behavior of the CTMC in the first row of Figure 3.

For $D < 2/V$, the finite trajectory of the 2-dimensional TK model can be observed in Figure 1b, where the process spends most of its time on the boundary, i.e., one species is almost extinct. Hence, the switching time between boundaries is an important metric describing the finite-time dynamics of TK models. For simplicity, we will define the switching time as the first time that $X^2$ reaches 0, assuming the process starts initially with no $A_1$, i.e., $(X^1_0, X^2_0) = (0, 2V)$. Note that finite time properties, including switching time distributions, cannot be extracted from the stationary distributions.

Note that for both the CTMC and the CLA, the trajectory will inevitably visit between boundaries for all $V > 0$ (even for large $V$) because both stochastic models are positive recurrent. For a larger $V$, the trajectory typically fluctuates around an equilibrium point for a long time before hitting the boundary and visits between boundaries less often than that for a smaller $V$. In Figure 4a, we illustrate the switching time for a range of $V$ values that covers all three parameter regimes: $D < 2/V$, $D = 2/V$ and $D > 2/V$.

In an attempt to obtain an explicit formula for the switching time, we consider a 1-dimensional approximation of the 2-dimensional CLA. Roughly, we assume that the inflow and outflow reactions occur at a much slower rate than that of the autocatalytic reactions (i.e., $\lambda', \delta' \ll \kappa'$). Then the total mass evolves at a much slower time scale than that of the autocatalytic reactions (which are mass-conserving) so that $Z^1_t + Z^2_t\approx Z^1_0 + Z^2_0 = n$ for a long time period. We then approximate $Z^2_t$ by $n-Z^1_t$ and the equation for $Z^1$ is given by

where $L_t$ is the local time of $S_t$ on the boundary of the interval $[0, n]$ and $m$ is the inward normal vector at the boundary of the interval (so $m(0) = 1$ and $m(n) = -1$). Then we let $\tau^{[0, n)} = \inf\{t\geq 0:\, S_t \notin [0, n)\}$ be the first time the process $S_t$ hits $n$; then, the expectation $\mathbb{E}_x\left[\tau^{[0, n)}\right]$ solves the boundary value problem as in ^[38]:

where $\mathcal{L}$ is the generator of (3.1); namely, suppose $f \in C^2([0, n])$; then,

Hence the expected hitting time of $n$ can be expressed explicitly:

where

In Figure 4, the mean switching time is plotted against different parameters in (2.1), while fixing other parameters whose values can be found within the caption. The swithcing time was computed for the CTMC (1.5)-(1.6) as well as the CLA 2D given by (2.2), and the 1-dimensional reflected approximation in (3.1), which is termed as CLA 1D in Figure 4.

As $V\rightarrow \infty$ in Figure 4a, the switching time of the CTMC can be better approximated by CLA 2D in (2.2). Such a trend can also be observed when $\kappa < 1$ in Figure 4b, whereas the CLA 1D given by (3.1) becomes better approximations as $\kappa\rightarrow \infty$. Heuristically as $\kappa \rightarrow \infty$, the dynamics of the CTMC are dominated by the autocatalytic reactions of (1.3), which coincides with the motivation of the CLA 1D given by (3.1) where the total mass evolves at a much slower time scale. This approximation is less accurate for $V < 40$ in Figure 4a.

We compare the switching time distributions between the CTMC and the CLA via the histogram in Figure 5. Parameters of the simulation are given by $V = 64,$ and $\kappa' = 1$ in (2.1) for all four figures. In addition, in Figure 5a and 5b $\lambda' = \delta' = 1/32$ so that the CTMC in (1.5)-(1.6) possesses a uniform distribution when conditioned on the level sets $\{x+y = n\}$; in Figure 5a and 5b, $\lambda' = \delta' = 1/64$, so that the CTMC in (1.5)-(1.6) possesses a bimodal stationary distribution as mass is concentrated on the meta-stable states near the boundary. In both cases, the mean and variance of the switching time are quite close and the shape of the distributions is nicely recovered.

However, as $\lambda'$ or $\delta'\rightarrow 0$, the switching time estimates using CLA 2D (2.2) are no longer close to the CTMC, as shown in Figure 4c and d. In Figure 6, such cases were investigated for $V = 64, \lambda' = \delta' = 1/256$ and $\kappa' = 1$, where the switching time distribution of the CTMC (1.5)-(1.6) is approximated by the CLA 1D (3.1). This approximation yields a much better result than the CLA 2D in (2.2), as inflow and outflow reactions occur much less frequently than autocatalytic reactions. However, the approximation via the CLA 1D (3.1) still underestimates the mean and variance of the switching time for the exact CTMC model.

3.2.   Simulation result for the 3-dimensional TK model

In this section, we first explain the discrepancy between the finite time dynamics of 2-dimensional and 3-dimensional TK models. Then, we show via simulation that the CLA in (2.2) recovers the stationary distribution of the associated CTMC given by (1.5)-(1.6). Furthermore, we propose the cycling time, which is analogous to switching time when $d = 2$, and the distribution of cycling time of the CTMC can be well approximated via its associated CLA when $V$ is large, and this approximation is less accurate when $V$ is small.

The main contrast between $d = 2$ and $d = 3$ for the CTMC, as observed in Figure 1b and 1d, is that the transition time between meta-stable states for $d = 2$ is much longer than that for $d = 3$. A heuristic explanation for this is that for the 2-dimensional TK model, the two autocatalytic reactions given by (1.3) move in opposite directions while sharing the same intensity, hence making it difficult for the process to move between meta-stable states; whereas, for higher-dimensional TK models ($d\geq 3$), autocatalytic reactions move the process in a cyclic direction, namely, $A_1\rightarrow A_2 \rightarrow \cdots \rightarrow A_d \rightarrow A_1$. More specifically, when $A_1$ is abundant, the dynamics of $A_1$ are mainly driven by autocatalytic reactions $A_d+A_1\rightarrow 2A_1$ and $A_1+A_2\rightarrow 2A_2$. However, the firing of both reactions tends to decrease the rate of $A_d+A_1\rightarrow 2A_1$ and increase the rate of $A_1+A_2\rightarrow 2A_2$, leading to an imbalance toward the gaining of $A_2$; hence most of $A_1$ will be changed to $A_2$.

Similar to 2-dimensional TK models, positive recurrence has been established for $d = 3$ in ^[2]; however, the explicit form of the stationary distribution is only derived when $\delta = \frac{3}{2}\kappa$. In Figure 7, the stationary distribution of the CTMC given by (1.5)-(1.6) is plotted on the hyperplane $\{ x+y+z = 3V\}$ for $\delta > \frac{3}{2}\kappa$, $\delta = \frac{3}{2}\kappa$ and $\delta < \frac{3}{2}\kappa$ in the first row. Similar multi-modalities of stationary distributions are observed when $\delta < \frac{3}{2}\kappa$. In the second row of Figure 7, the stationary distribution of the associated CLA of (2.2) is plotted by using the densities near the hyperplane, namely $\left\{ |x+y+z-3|\leq \frac{1}{128}\right\}$. In all three cases, the CLA captures the stationary behavior of the exact CTMC model.

As discussed regarding the different transition times between meta-stable states for $d = 2$ and $d = 3$, the dominant species when $d = 3$ form a cycle, $A_1\rightarrow A_2\rightarrow A_3\rightarrow A_1$. Note that such cyclic behavior is not exclusive to parameter regimes when stationary distributions are concentrated near the boundaries. In Figure 8, trajectories of the 3-dimensional TK model are plotted for $V = 256, \; \kappa' = 1$ and $D = \lambda' = \delta'$ given by 1/32 and 3/512 respectively. The stationary distribution is unimodal or uniform when conditioning on the hyperplane $\{x+y+z = 3V\}$ respectively. In both cases, cyclic behavior still persists within the trajectory; hence, we propose cycling time as the analogous quantity for describing finite-time dynamics for $d\geq 3$, which can be defined as the second time that $X^3$ reaches peak abundance when the initial condition is given by $(X^1_0, X^2_0, X^3_0) = (0, 0, 3V)$. Note that the switching time is also well defined in 3-dimensional TK models; however, the switching could end up in different regions, i.e., if initially there are only species $A_3$, it could move to the boundary with only species $A_1$ as well as the boundary with only species $A_2$. In comparison, the cycling time is more straightforward and consistently captures the average behavior of the 3-dimensional TK model.

We will simulate and compare cycling time distribution for the CTMC given by (1.5)-(1.6) and its associated CLA given by (2.2). In Figure 9, the mean switching time is plotted against different parameters in (2.1)and we fixed the other parameters whose values can be found within the caption. Cycling time was computed for the CTMC as well as the CLA 3D in (2.2). Similar to the 2D results, the mean cycling time is better captured when $V\rightarrow \infty$, or when $\kappa\rightarrow 0$.

In Figure 10, the cycling time distribution for the CTMC is compared with its associated CLA (2.2). Simulation of the cycling time in Figure 10a and 10b was obtained by using $V = 256, \; \kappa' = 1$ and $\lambda' = \delta' = 1/32$, which corresponds to a unimodal stationary distribution. In this case, the cycling time distribution of the CTMC, as well as its mean and variance, is nicely recovered by its associated CLA.

In Figure 10c and 10d, we used the same $V$ and $\kappa'$ while choosing $\lambda' = \delta' = 3/512$, which corresponds to the stationary distribution being uniform when conditioning on the hyperplane. However, the cycling time distribution associated with CLA underestimates the mean and variance of the cycling time associated with the original CTMC.

3.3.   Simulation results for the general d-dimensional TK model

In this section, we present stochastic simulations for higher dimensional TK models. We first give a precise description of the 6-dimensional TK model in terms of two switchings: a slow switching between even and odd species, as well as fast switching between dominant regions near the boundary. We then investigate the effect of the dimension $d$ on the mean cycling time.

The dynamics of the 4-dimensional CTMC given by (1.5)-(1.6) are elaborated in the original publication of ^[1]. When the system volume $V$ is large, its trajectory fluctuates around an equilibrium point for a long time before hitting the boundary; when $V$ is small, random fluctuation dominates the system. Interesting behaviors emerge when $V$ resides in an intermediate range, where the extinction of species slowly switches between odd and even species in auto-catalytic reaction loops. More precisely, the model switches between states that are abundant in odd species and states that are abundant in even species, during which the other species are almost extinct. Moreover, within a temporal domain of abundant odd species, there are fast switches between species $A_1$ and $A_3$ with a large imbalance between the two species, either $X^1\gg X^3$, or $X^1\ll X^3$.

To describe the finite time dynamics for the $6-$dimensional TK model, we simulate trajectories of the CTMC given by (1.5)-(1.6), with the parameters $V = 64, \lambda' = \delta' = 1/256$ and $\; \kappa' = 1$. To investigate its switching behavior, we define the process $B(t)$ below that captures the disparity between odd species and even species:

with the initial condition $X^i_0 = V$ for all $i = 1, 2, \cdots 6$. Exponential convergence to the stationary distribution of the CTMC, as established in ^[2], guarantees that the denominator $6V$ is approximately the average total population when $T$ is large. A similar expression is considered in ^[1] for the 4-dimensional TK model to identify regions of $V$ where DITs persist.

In Figure 11a, the finite time trajectory of $B(t)$ is plotted, and $B(t)$ switches between regions near $B(t) = 1$ and $B(t) = -1$, which are regions with abundant odd or even species respectively. The stationary distribution of $B(t)$ is plotted in Figure 11b via time averaging, which also yields a bimodal distribution having peaks near $-1$ and $1$.

Next, we investigate the dynamics in between slow switches by analyzing trajectories with abundant odd species. In particular, the joint distribution of

conditioned on odd species being abundant, is plotted in Figure 12a under stationarity via time averaging. Throughout the simulation, the abundance condition was approximated by $B(t)\geq 0.95$. The joint distribution is concentrated within three boundary regions, which can be specified by $\Omega_i = \{ \rho_i \approx 0\}$ for $i = 1, 3, 5$, as each boundary implies that one odd species is almost extinct while two other odd species are abundant. Moreover, these regions are not symmetric with respect to the two dominant species, in the sense that $(\rho_1, \rho_3)$ is on average $(0.4, 0.6)$ when $\rho_5 \approx 0$. Such asymmetry exists since the states with more $A_1$ are sensitive with respect to the birth of $A_2$. More specifically, gaining $A_2$ (hence losing $A_1$) and losing $A_2$ (hence gaining $A_3$) is characterized by the autocatalytic reactions $A_1+A_2\rightarrow 2A_2$ and $A_2+A_3\rightarrow 2A_3$, the rates of which are completely determined by the relative counts of $A_1$ and $A_3$. For states with $\rho_1\geq \rho_3$, higher counts of $A_1$ speed up the gain of $A_2$; hence, there is the transition $A_1\rightarrow A_2 \rightarrow A_3$ into states with $\rho_1\leq \rho_3$; on the other hand, states with $\rho_1\leq \rho_3$ are likely to remain unchanged, since $A_2$ is more likely to be exhausted.

In addition to asymmetry, the dynamic of odd species, conditioned on odd species being abundant, moves between three dominant regions $\{\Omega_i\}_{i = 1, 3, 5}$ in a clockwise manner as plotted in Figure 12a, $\Omega_5\rightarrow \Omega_3 \rightarrow \Omega_1 \rightarrow \Omega_5$. More specifically in the region $\Omega_5$ where $A_5$ is almost extinct, a birth of $A_4$ would transition $A_3$ into $A_5$, which moves the process into the region $\Omega_3$. Births of all other species can not change the dominant molecules of $A_1$ and $A_3$. Similar movement applies recursively leading to a clockwise cycle in Figure 12a.

To support our claims, a sample trajectory of odd species is shown in Figure 12b, during which we only have abundant odd species. The process stays within the region $\Omega_5$ at time $t = 9000$, which switches into the region $\Omega_3$ at the time $t = 9010$, and then proceeds to $\Omega_1$ around $t = 9020$. At time $t = 9025$ the process returns to $\Omega_5$ and completes a clockwise cycle.

We summarize the dynamics for $d = 6$. Two types of switches can be utilized to describe the finite time dynamics of the $6-$dimensional TK model when DITs persist. In particular, abundant molecule species switch between odd and even on a slower time scale. In between these switches, the dynamics dominated by odd species will cycle between boundary regions as $\Omega_5\rightarrow \Omega_3 \rightarrow \Omega_1 \rightarrow \Omega_5$, or $\Omega_6\rightarrow \Omega_4 \rightarrow \Omega_2 \rightarrow \Omega_6$ when the dynamics are dominated by even species.

Last but not least, the mean cycling time is plotted against some dimension $d$ in Figure 13. As in previous sections, the cycling time was computed for each trajectory as the second time that $X^d_t$ reaches peak abundance (initially there are only species $d$), and the mean cycling times were averaged over independent samples. Despite longer cycles as $d$ increases for $d\geq 5$, the cycling time decreases as $d$ increases from 3 to 5. This is because if the birth of species $A_2$ occurs before the birth of species $A_1$, the switches from $A_d\rightarrow A_1$ and $A_1\rightarrow A_2$ would happen simultaneously, as can be observed in Figure 2. Such simultaneous switches speed up the cycling and hence reduce the cycling time, whereas these switches are more likely to occur in higher dimensions. However as $d$ increases further, such an effect is properly averaged out and the mean cycling time approximately increases linearly in $V$.

4.   Proofs

This section contains the proofs of the results stated in Section 2. The proofs will be collected in Subsection 4.3, after we develop some preliminary estimates for the CLA in Subsections 4.1 and 4.2.

4.1.   Estimates for the coefficients of CLA

Let $\Gamma, \sigma, b$ and $\gamma$ be defined as in (2.3) and (2.4), and note that $b$ is globally Lipschitz for the dimension $d = 2$, but only locally Lipschitz for $d\geq 3$ due to the terms $\sum_{k = 1}^d e_k(\kappa'(x_{k-1} - x_{k+1})x_k$. Precisely,

where $|x| = \sqrt{\sum_{i = 1}^dx_i^2}$ is the Euclidean norm of $x$, and $|x|_{\infty}: = \max_{1\leq i\leq d}|x_i|$.

Note also that $\Gamma(x)$ is symmetric and positive for all $x \in \mathbb{R}^d$. Precisely, for all $x \in \mathbb{R}^d_+$ and $\theta \in \mathbb{R}^d$,

We can show that $\sigma$ is locally Lipschitz and grows linearly. We note that

So the operator norm of the $d\times d$ matrix grows quadratically:

for some $C > 0$ by the equivalence of norms in $\mathbb{R}^n$. Then we have

Local Lipschitzness of $\sigma(\cdot)$ is given by (4.1) and the Powers-Stormer inequality [39,Lemma 4.2]:

where $\|\cdot\|_{HS}$ denotes the Hilbert-Schmidt norm. Since all norms are equivalent in a finite dimensional vector space and each entry in $\Gamma(\cdot)$ is a second order polynomial, $\Gamma(x)$ is locally Lipschitz; hence, $\sqrt{\Gamma(\cdot)} = \sigma(\cdot)$ is also locally Lipschitz.

4.2.   Lyapunov inequalities and moment estimates

Recall the Lyapunov function $U$ defined in (2.5) and the differential operator ${\cal L}$ in (2.8). We also write $U^p(x) : = \left(U(x)\right)^p$ for simplicity.

Lemma 1. Let $\gamma$ be as in (2.3); then, we have $\nabla U^p(x) \cdot \gamma(x) \le 0$ for all $x \in \mathbb{R}_+^d$ and $p\in\mathbb{N}$.

Proof. Observe that for all $x \in \mathbb{R}^d_+$,

where $\{e_i\}_{i = 1}^d$ is the standard basis in $\mathbb{R}^d$; since $\gamma(x) = \frac{b(x)}{|b(x)|}$, we only need to check $\nabla U^p(x) \cdot b(x) \le 0$:

Lemma 2 (Lyapunov inequalities). Let $U$ and ${\cal L}$ be defined in (2.5) and (2.8) respectively. Then for all $p\in\mathbb{N}$, there exist constants $c_p, c_p' > 0$ such that

for all $x \in \mathbb{R}_+^d$. Furthermore, there is a compact set $K \subset \mathbb{R}_d^+$ and $f: \mathbb{R}^d_+\to [1, \infty)$ such that

for some positive constants $c_1'$ and $c_2'$.

Proof. For $x = (x_1, \cdots, x_d) \in \mathbb{R}^d_+$ and $1 \le i, j \le d, x\in \mathbb{R}_+^d$, we have that

Applying the differential operator ${\cal L}$ to $U^p(x)$,

Let $c_p = \frac{p\delta'}{2}$ then

which is bounded above in $\mathbb{R}^d_+$ by some positive constant $c_p'$ due to the fact that $\left(\frac{2p-1}{V}\left(|x|_1 + \frac{d\lambda'}{\delta'}\right) - \frac{1}{2}\left(|x|_1 - \frac{d\lambda'}{\delta'}\right)^2\right)$ is quadratic with the leading coefficient being a negative number. This proves (4.2).

To prove (4.3), we let $p = 1$ and $f(x) : = U(x) + 1$ and let $c_1' = \frac{\delta'}{2}$; then,

By the basic facts of quadratic functions, ${\cal L}U(x) + c_1' f(x)$ is only positive on a compact set in $\mathbb{R}_+^d$; we call it $K$, and since it is also continuous, it is uniformly bounded by some $c_2' > 0$. So (4.3) holds.

Lemma 3. Let $x \in \mathbb{R}^d_+$ and $Z$ be the solution to (2.2) with $Z_0 = x$; then,

Proof. Suppose $Z$ is the process that solves the CLA (2.2) that starts at $x \in \mathbb{R}_+^d$ and let $U$ be defined as (2.5); then, by Ito's formula, we have

where $M_t$ is the local martingale term that has the following explicit expression:

where $(\sigma_{i, j})_{i, j = 1}^d = \sigma = \sqrt{\Gamma}$ is the dispersion matrix and $\{W^{(j)}\}_{j = 1}^d$ are independent 1-dimensional Brownian motions. By Lemma 1, the process $\int_0^t\nabla U(Z_s) \cdot \gamma(Z_s) dL_s$ is non-positive and decreasing since $\gamma(x) = \frac{b(x)}{\|b(x)\|}$. Therefore we have the following almost sure inequality for all $t \ge 0$:

Note that $\sum_{i, j = 1}^d \Gamma_{i, j}(x) = d\lambda' +\delta' |x|_1 = \delta'\sqrt{U(x)} +2 d\lambda'$ and $\delta' U(x) \ge 0$ for all $x$, therefore we can rewrite the above inequality as

For each $m\in\mathbb{N}$, we let $\{\tau_m'\}_{m = 1}^\infty$ be the family of stopping times defined by $\tau_m' = \inf\{t\geq 0: |Z_t|_1 \ge m\}$ and let $\tau_m = \tau'_m \wedge m$. Define the stopped process $\widetilde{Z}^{(m)}$ by $\widetilde{Z}^{(m)} = (\widetilde{Z}^{(m), i}_\cdot)_{1 \le i \le d} = (Z^{i}_{\cdot \wedge \tau_m})_{1 \le i \le d} = Z_{\cdot \wedge \tau_m}$. Let $M_t^{(m)} = M_{t\wedge \tau_m}$ and replace $Z_t$ and $M^t$ in (4.8) by $\widetilde{Z}^{(m)}_t$ and $M_t^{(m)}$ respectively; then raise both sides by a power of $2p$; by Jessen's inequality we get that, for all $T > 0$ and $t \in [0, T]$, there is a $C$ depending only on $T$ and $p$ such that

Now take the sup over time and expectation $\mathbb{E}_x$ to get

We wish to obtain an inequality of the Gronwall type; note that the square root function is concave on $\mathbb{R}_+$, so by Jensen's inequality we get

Denote $\langle M\rangle_\cdot$ as the quadratic variation of a process $M$. By the BDG inequality, there is an absolute constant $C > 0$ that depends only on $p$ such that

Furthermore, for all $T > 0$, there exist constants $C_1,$ and $C_2$ (depending on $T$ and $p$) such that for $t\in[0, T]$,

where the above positive constant $C_2$ depends only on $p$ and $T$ and is again independent of $m$. Finally, by Jensen's inequality again, we see that there is $C_0 > 0$ that depends on $T, p$ and $V$ such that the following inequality holds for all $t \in [0, T]$

We let $y(t): = \mathbb{E}_x\left[\sup_{s \in [0, t]} U^{2p}(Z_s^{(m)})\right]$ and combine the inequalities (4.9), (4.10) and (4.12) to get the following:

We let $\omega(s) = (\delta'^{2p} +(2d\lambda')^p)\sqrt{s} + (2d\lambda')^{2p} + \delta'^p s^\frac{3}{4}$, which is positive and strictly increasing on $[0, \infty)$. Therefore, $\Phi:[0, \infty)\to (0, \infty)$ defined by $\Phi(t) = \int_0^t \frac{1}{\omega(s)}ds$ is strictly increasing and continuous. Hence $\Phi^{-1}$ exists and is continuous. Also, note that $\mathbb{E}_x[\sup_{s \in [0, t]}U^{2p}(\widetilde{Z}^{(m)}_s)]$ is continuous in $t$, so by the Gronwall-type inequality [40,Theorem 4,p3] with $C_0U^{2p}(x)$ and $C_0$ in the place of $M$ and $\Psi$, we have

Note that the right hand side does not depend on $m$, so by taking $m\to\infty$ on the left hand side and invoking Fatou's lemma, we see that $\mathbb{E}_x[\sup_{s \in [0, t]}U^{2p}(Z_s)]$ is finite for all $t \ge 0$ which implies (4.5).

Proposition 3. The solutions to the (2.2) starting from different $x\in \mathbb{R}_+^d$ form a Feller process.

Proof. Let $\tau^x_M = \inf\{t \ge 0:|Z^x_t|_1 \ge M\}$ for $M > M^* = \frac{d \lambda'}{\delta'}+1$ as in [24,Eq (4.25)]; then, $\tau^x_M$ is a stopping time, where we note that $u$ in [24,Eq (4.25)] is equal to $(1, 1, \cdots, 1)$ in our case. We will show that for all $f \in C_b(\mathbb{R}^d_+)$ and $t \ge 0$ with $x \in \mathbb{R}_+^d$,

where $Z^x$ and $Z^y$ are strong solutions to the CLA (2.2) that start at $x, y \in \mathbb{R}_+^d$ respectively. This implies the Feller property since the function $x \mapsto \mathbb{E}[f(X_t^x)]$ is bounded.

Fixing $x \in \mathbb{R}^d_+$, we first note that for any $\epsilon > 0, \exists M_\epsilon > M^* > 0$ such that for all $M \ge M_\epsilon$ we have

Indeed, by the Markov inequality, we have

and by (4.14) and continuity of $\Phi$ in the proof of Lemma 3, the map $\mathbb{R}_+^d\ni y \mapsto\mathbb{E}\left[\sup_{s\in [0, t]}|Z^y_t|^4\right]$ is uniformly bounded on a compact set; hence, the numerator on the right-hand side is bounded by some $C_x > 0$ depending only on $x$. So the right-hand side goes to zero as $M \to \infty$.

Assume $y \in B(x, 1)$ and $M > M_\epsilon$; then, we have the following decomposition of (4.15):

where $Z^{x, (M)}_\cdot$ denotes the stopped process $Z^x_{t \wedge\tau_M^x}.$ Then following the proof of Theorem 6.1 in ^[24], with its modification to the proof of Theorem 5.1 of ^[26] to the stopped process, we obtain the following inequality similar to [24,Eq (6.4)]: for each $T > 0$, there is a constant $C > 0$ such that for $t\in[0, T]$,

By the Gronwall's inequality, as $B(x, 1) \ni y \to x$, the left hand side of (4.18) goes to zero. Therefore, from (4.17) we see that

and since $\epsilon > 0$ is arbitrary, we see that (4.15) holds.

4.3.   Proofs for the results in Section 2

Proof of Theorem 1. The reaction network of the TK model contains inflows and outflows of all species (1.4). Furthermore, the autocatalytic reactions given by (1.3) satisfy the mass-conserving/mass-dissipating assumption in [24,Definition 3.1(a)] with $u = (1, 1, \cdots, 1)\in \mathbb{R}^d$. Hence [24,Assumption 3.1] is satisfied. Strong uniqueness of the CLA follows from [24,Theorem 6.1], and weak existence follows from [24,Section 7]. Now strong existence and weak uniqueness follow from the Yamada-Watanabe-Engelbert theorem ^[41].

The Feller property is given by Proposition 3; hence, it also has strong Markov property by ^[41].

The following result from [27,Theorem 4.2] (see also [42,Theorem 2.2.12]) provides a condition on Lyapunov functions that guarantees the existence of a unique invariant distribution for Feller diffusions. Note that a skeleton chain of a Feller diffusion also possesses the Feller property. By [43,Proposition 2.2], every compact subset is petite for a Feller diffusion and all of its skeleton chains.

Theorem 4. ([27,Theorem 4.2]) Let $X$ be a Feller diffusion. If $U$ is a positive function such that for some positive constants $c_1, c_2 > 0$, a function $f:\mathbb{R}^d\to[1, \infty)$ and a compact set $K\subset \mathbb{R}^d_+$, $U$ is bounded on $K$ and the following inequality holds for $X$:

then the diffusion is positive Harris recurrent and there is an invariant probability measure $\pi$ for $X$; also, any invariant probability $\pi$ satisfies $\int_{\mathbb{R}_+^d} f(x)\pi(dx) < \infty$.

The following result from ^[27] (see also [42,Theorem 2.2.15]) provides a condition on Lyapunov functions that guarantees exponential ergodicity of Feller diffusions.

Theorem 5. ([27,Theorem 6.1]) Let $X$ be a Feller diffusion. Assume that there exists a norm-like function $U$ and constants $c_1 > 0$ and $c_2\in \mathbb{R}$ such that $X$ satisfies

for all $x\in\mathbb{R}^d_+$. Then $X$ is $f$-exponentially ergodic with $f = U+1$.

Conditions (4.19) and (4.20) are called (CD2) and (CD3) respectively in ^[27], and they are satisfied for our CLA (process $Z$) thanks to Lemma 2.

Proof of Theorem 2. Note that Proposition 3 implies that $Z$ is a Feller process. Let $U$ be the Lyapunov function defined as (2.5), and the inequality (4.3) implies that $U$ satisfies the inequality (4.19). Therefore, by Theorem 4, $Z$ is positive Harris recurrent and hence has a unique invariant probability measure [27,Section 4.1].

It remains to show that all moments of the stationary distribution $\pi$ are finite. By Ito's formula, Lemmas 1 and 2, we have that for each $p \in\mathbb{N}$, there exist some positive constants $c_p$ and $c_p'$, and the following inequality holds for $t\in\mathbb{R}_+$:

By rearranging terms and dividing by $t$ and $c_p$, it follows that

Now, let us define $U^p_M$ as the truncated function of $U^p$ at $M \ge 0$, that is,

Then $U^p_M$ is a bounded continuous function. Since $Z_t$ converges to $\pi$ in law, we have that

Therefore, by (4.22) and (4.21)

Now, take $M\to \infty$ on the left hand side and by the monotone convergence theorem, we have

Since the inequality holds for all $p \in\mathbb{N}$, we may conclude that all moments of $\pi$ are finite.

Proof of Theorem 3. Let $U$ be the Lyapunov function defined in (2.5). Inequality (Eq 4.2) says that $U$ satisfies the inequality (Eq 4.20); hence, we can apply Theorem 5 to conclude that $Z$ is $f$-exponentially ergodic with $f = U+1$.

Proof of Proposition 1. Following ^[31], we let $n(x) = \{\sum_{i \in {\cal I}(x)} \alpha_i e_i, \alpha_i > 0\}$ be the set of interior normal vectors to the domain $\mathbb{R}_+^d$ at $x \in \partial\mathbb{R}_+^d$, where ${\cal I}(x) = \{1 \le i \le d: x_i = 0\}$. Let

This definition is a bit different from that of [31,Eq (6)] where they define $d(x)$ as a set valued function, but since our reflection $\gamma(x)$ is well defined for all $x \in \partial \mathbb{R}_+^d$ including the non-smooth part, we can set it to be single valued. If $x = (x_1, \cdots, x_d) \in \partial \mathbb{R}_+^d$, then there is some $i$ between $1$ and $d$ such that $x_i = 0$; hence, $e_i \in n(x)$, so $b_i(x) = \lambda'$ and

Hence ${\cal V} : = \partial \mathbb{R}_+^d\backslash {\cal U} = \emptyset$.

We prove Proposition 1 by checking all of the conditions in [31,Theorem 2]: note that $\Gamma(x)$ is uniformly elliptic for all $x \in \mathbb{R}_+^d$ by (4.1), and that the reflection $\gamma$ is piece-wise $C^2(\partial \mathbb{R}_+^d)$. [31,Assumption 2] is satisfied since ${\cal V} = \emptyset$ by [31,Remark 3.4]. The wellposedness of the submartingale problem in the statement of [31,Theorem 2] is given by [44,Theorem 1] and Theorem 1. Now, all assumptions of [31,Theorem 2] are satisfied, which proves our statement.

Proof of Proposition 2. We prove the statement by checking all conditions in [31,Theorem 3]: by the proof of Proposition 1 we see that [31,Assumption 2] is satisfied and the corresponding submartingale problem is well posed. Furthermore, all entries of $\Gamma(\cdot)$ and $b(\cdot)$ are smooth since they are polynomials, so [31,Theorem 3] holds with $\mathbb{R}_+^d, \gamma(x), b(x), \Gamma(x)$ in the place of $G, d(x), b(x), a(x)$ in Eqs (12)–(16) of [31,p. 1341].

Acknowledgments

This research was initiated during the American Institute of Mathematics (AIM) workshop "Limits and control of stochastic reaction networks" in July 2021. The authors gratefully acknowledge the support of AIM and the organizers of the workshop. The authors are indebted for the stimulating discussion during the monthly TK group meetings with Lea Popovic, Ruth Williams, Grzegorz Rempala, Hye Won Kang, Enrico Bibbona, Siri Paola, Wasiur Khuda Bukhsh and Felipe Campos Vergara. This research was partially supported by NSF awards DMS 1855417 and DMS 2152103 and ONR grant TCRI N00014-19-S-B001 to W. T. Fan.

Conflict of interest

The authors declare no conflicts of interest.