# A Scaling Approach to Stochastic Chemical Reaction Networks

Lucie Laurence1<EMAIL_ADDRESS>INRIA Paris, 2 rue Simone Iff, 75589
Paris Cedex 12, France  and  Philippe Robert<EMAIL_ADDRESS>http://www-rocq.inria.fr/who/Philippe.Robert

###### Abstract.

We investigate the asymptotic properties of Markov processes associated to
stochastic chemical reaction networks (CRNs) driven by the kinetics of the law
of mass action. Their transition rates exhibit a polynomial dependence on the
state variable, with possible discontinuities of the dynamics along the
boundary of the state space. As a natural choice to study stability properties
of CRNs, the scaling parameter considered in this paper is the norm of the
initial state. Compared to existing scalings of the literature, this scaling
does not change neither the topology of a CRN, nor its reactions constants.
Functional limit theorems with this scaling parameter can be used to prove
positive recurrence of the Markov process. This scaling approach also gives
interesting insights on the transient behavior of these networks, to describe
how multiple time scales drive the time evolution of their sample paths for
example. General stability criteria are presented as well as a possible
framework for scaling analyses. Several simple examples of CRNs are
investigated with this approach. A detailed stability and scaling analyses of
a CRN with slow and fast timescales is worked out.

1Supported by PhD grant of ENS-PSL

###### Contents

  1. 1 Introduction
  2. 2 Mathematical Models of CRNs
  3. 3 Classical Stability Results
  4. 4 A Stability Criterion
  5. 5 Scaling Methods
  6. 6 Binary CRN Networks
  7. 7 Agazzi and Mattingly’s CRN
  8. 8 A CRN with Slow and Fast Timescales
  9. References
  10. A Stability Results
  11. B Technical Proofs
  12. C General Triangular Topologies

##  1\. Introduction

This paper investigates the asymptotic properties of Markov processes
associated to chemical reaction networks (CRNs) driven by the law of mass
action kinetics. The stability of a CRN is the positive recurrence property of
its Markov process. The state space of these processes is a subset of
$\mathbb{N}^{m}$, where $m{\geq}1$ is the number of chemical species. An
important feature of these processes is that their transition rates exhibit a
polynomial dependence on the state variable. Another important characteristic
is that there are possible discontinuities of the dynamical behavior along
some boundaries of the state space. This is due to the fact that some chemical
reactions require a minimal number of copies of some chemical species to take
place. This is referred to as boundary effects. Due to the polynomial property
of the reaction rates, multiple timescales may drive the time evolution of
these processes. To investigate the property of positive recurrence, one
starts from a “large” state and tries to figure out how it can go back to a
neighborhood of the origin. As it will be seen in several examples, depending
on which component is large, the timescale to return to $0$ can be very
different, this is where discontinuities of the kinetics play an important
role.

### Deficiency Zero Theorem

Early models of these systems in chemical physics are deterministic dynamical
systems in a continuous state space, $\mathbb{R}_{+}^{m}$, expressed as the
solution of a specific set of Ordinary Differential Equations (ODEs) expressed
with polynomials of the state variable. This is the so-called law of mass
action, see [21]. This set of ODEs is closely related to a weighted oriented
graph structure on a convenient space, the space of complexes. A complex is a
set of chemical species involved in a chemical reaction. The weights on the
edges of this graph are the reaction rate constants for the transformation of
a complex into another complex. A striking and impressive result in this
context is the Deficiency Zero Theorem which states that, under convenient
topological conditions for this oriented graph, the ODEs have a unique fixed
point which is locally stable. The main references are [27], [26] and [19].
See also [22] and [18] for a general detailed overview. Section 3 gives a
quick presentation of these results. Precise definitions are in Section 2.

These continuous models with state space $\mathbb{R}_{+}^{m}$ have been used
in chemical physics to describe the time evolution of concentrations of $m$
chemical species in a large volume. For a discrete state space such as
$\mathbb{N}^{m}$, the state variable gives the number of copies of $m$
chemical species. In such a stochastic context, an important result, proved in
[6] shows that, under the topological assumptions of the deficiency zero
theorem, the associated Markov process has an invariant distribution which is
the product of $m$ exponential distributions whose parameters are given by the
fixed point of the dynamical system of the deficiency zero theorem. It is a
consequence of the fact that the invariant distribution of these systems
satisfies a set of detailed balance equations. This algebraic approach, which
consists in plugging a product formula solution in the equilibrium equations,
has been used in subsequent works in [12, 11] and [30] and references therein,
to other classes of CRNs. An interesting feature of these Markov processes is
that the existence and uniqueness of an invariant distribution is essentially
expressed in terms of the topological structure of the graph ${\cal R}$ and
not on its associated weights as long as they are positive.

The assumptions under which the results of the deficiency zero theorem and its
consequences hold cover a large class of CRNs, but they are clearly not
necessary for the existence and uniqueness either of a fixed point in the
continuous case, or of an invariant distribution in the discrete case.

### Analogies with Queueing Networks

This context is reminiscent of an analogous situation for queueing networks at
the end of 1970s. A queueing network can be described as a finite set of
processing units among which several flows of requests are circulating. A
request requires the processing capacity of a subset of nodes for some amount
of time before leaving the network. Several large classes of queueing networks
having product form invariant distributions were identified, like Jackson’s
networks or loss networks, see [33]. Outside these classes, little was known,
even on the existence of equilibrium distributions. In the 1990s, several
simple queueing networks, seemingly stable, were shown to be in fact unstable.
See [10]. At this occasion, convenient stability results using Lyapunov
functions and scaling ideas were developed. Some of these results could be
traced back to earlier Has’minskiĭ’s works in the 1960s on extensions of
deterministic stability results of the beginning of the twentieth century in a
stochastic context. See [25].

### Criteria for Stability Properties

When the Markov process $(X(t))$ associated to a CRN with $m$ chemical
species, with a state space included in $\mathbb{N}^{m}$, does not have an
explicit invariant distribution, like a product form, the positive recurrence
of $(X(t))$ can be, in principle, proved with a Lyapunov function $f_{0}$ for
its infinitesimal generator ${\cal Q}$. This is expressed with a relation of
the type ${\cal Q}(f_{0})(x){\leq}{-}\gamma$, for some $\gamma{>}0$, holding
outside a finite subset $F$ of the state space. Under convenient conditions,
this gives an inequality of the type,

(1) $\mathbb{E}_{x}(f_{0}(X(t_{0})))-f_{0}(x)\leq{-}\frac{\gamma}{2}t_{0},$

for all $x{\not\in}F$ and $t_{0}$ sufficiently small. There are many examples
of specific CRNs where this method works, using entropy-like Lyapunov
functions, see [4], [7], Anderson et al. [3, 5], and [45]. The mentioned above
boundary effects of CRNs complicate in general significantly the search of
such a $f_{0}$. See [2] where polynomial functions are built as possible
Lyapunov functions, but with a significant technical work.

There is a related result of the literature, see [20], which may simplify some
aspects of this problem in some cases. Its interest lies in the fact that the
process may be observed at a convenient random instant depending, possibly, on
the initial state. This amounts to have $t_{0}$ replaced by a random
(stopping) time $\tau$ in Relation (1) and $\mathbb{E}(\tau)$ in its right-
hand side. The stopping time $\tau$ may be chosen so that the state $X(\tau)$
is in some convenient subspace. In this way, it may avoid the usual, annoying,
fine tuning of parameters of a given class of potential Lyapunov functions,
depending on the state where the process is at time $t_{0}$. The problem is on
the definition of $\tau$ depending on the initial state, which is easier to
handle a priori. Additionally, quite often, it is enough to consider the norm
$\|\cdot\|$ of the state as a Lyapunov function. There are still complicated
situations, as the one considered in Section 8, but the search of a Lyapunov
function is in general simplified with this approach. Its potential usefulness
in the study of the stability of CRNs has perhaps not been fully realized. See
Sections 6.1, 7 and 8 for its use in a simple context. See also [1].

### Scaling Analysis with the Norm of the Initial State

If having the proof of the existence and uniqueness of an invariant
distribution is an interesting and important result, it does not necessarily
give a real insight on the dynamical behavior of these networks. Obviously,
there are many, very different, Markov processes having an invariant measure
which is a product of exponential distributions. So even for CRNs of the
deficiency zero, little is known on the transient characteristics of the
associated Markov processes.

The study of the transient behavior of CRNs is presented in our paper as the
investigation of the convergence in distribution of a family of the sample
paths of the type

$\left(\frac{X_{i}(g(x)t)}{\|x\|^{\alpha_{i}}},i{=}1,\ldots,m\right)$

when the norm of the initial state $x$, $\|x\|$, goes to infinity, for a set
of coefficients $\alpha{=}(\alpha_{i})$ and a function $g$, which can be of
the order of $1/\|x\|^{\beta}$, for some $\beta{\geq}0$. Several choices are
in general possible for $\alpha{\in}\mathbb{R}_{+}^{m}$ and $\beta{\geq}0$,
and this is not a simple task as it will be seen.

For queueing networks, this line of research has been developed via the
convergence in distribution of scaled version, for time and space, of the
sample paths of the corresponding Markov processes. For Jackson networks, a
natural scaling is $g(x){=}\|x\|$ and $(\alpha_{1}){=}(1)$. The timescale is
sped-up by the norm of the initial state, in this way it provides a “fluid”
picture of these networks. See [39]. It should be noted that there is
essentially one  function $g$ and one value for the $\alpha_{i}$. This is
mainly due to the fact that the associated Markov processes are partially
homogeneous random walk. CRNs exhibit a much more complex behavior for the
choices of $g$ and $(\alpha_{i})$.

We do believe that this is an interesting and important ingredient to use to
investigate CRNs, in particular to have a close look at the various
interacting multiple timescales of these networks. It also turns out that such
a study can be used in the analysis of stability of CRNs. Finally, compared
with classical scalings already used in the CRN literature, see [8], [31] and
[32], or Proposition 2, it should be noted that the scalings introduced in
this paper do not modify, the topology, the reaction rates, or the weights on
the graph ${\cal R}$, of the CRN. See Section 5. Sections 6.1 and 7 give
simple examples of such an approach.

### An Interesting Example

Stability criteria and possible scaling analyses are presented and discussed
in a general framework in Sections 4 and 5. As often, in a complicated context
such as Markov process associated to CRNs, interesting, specific, examples can
provide a validation of an approach but also motivate some of its further
developments. In Section 8, an illustration of our approach is given by the
investigation of the CRN

$\emptyset\mathrel{\mathop{\xrightleftharpoons[\kappa_{1}]{\kappa_{0}}}}S_{1}{+}S_{2},\hskip
28.45274ptpS_{1}{+}S_{2}\mathrel{\mathop{\xrightleftharpoons[\kappa_{3}]{\kappa_{2}}}}pS_{1}{+}2S_{2},$

with $p{\geq}2$. It has been introduced and discussed in [1] for $p{=}2$. Its
transient behavior exhibits several unusual scaling properties with respect to
the norm of the initial state, which we will describe quickly. This is the
main motivation of our detailed analysis of section 8. Such a behavior is
probably not an exception in the context of CRNs, but more likely quite
common. Up to now it has not been really investigated. See [14] for a related
model. Additionally, it gives an interesting example of the use of time change
arguments to derive scaling results.

This CRN has the deficiency zero property. The associated Markov process has
in particular a product form measure. We do not use this property in our
analysis. If external arrivals of chemical species $S_{1}$ and/or $S_{2}$ are
added to the CRN, the deficiency zero property breaks down, nevertheless our
approach can be carried out much in the same way. The positive recurrence
property is proved with a convenient use of Filonov’s criterion mentioned
above.

The associated Markov process $(X_{1}(t),X_{2}(t))$ has several boundary
behaviors: at least $p$ copies of $S_{1}$ are required for the second reaction
to occur, and negative jumps of $(X_{1}(t))$ can occur only when $(X_{2}(t))$
is not null. It leads in fact to a kind of bi-modal behavior for the CRN.

A scaling analysis of this CRN is achieved for two classes of initial states.

  1. (a)

An initial state of the form $(N,b)$, with $b{\in}\mathbb{N}$ fixed.
Our main result, Theorem 26 shows that there exists some $t_{\infty}{>}0$,
such that the convergence in distribution of processes

$\lim_{N\to+\infty}\left(\frac{X_{1}(Nt)}{N},t{<}t_{\infty}\right)=(C(t_{\infty}{-}t),t{<}t_{\infty}),$

holds for some constant $C{>}0$.

  2. (b)

If the initial state is of the form $(a,N)$, $a{<}p$.
For the convergence in distribution of the associated occupation measures, see
Definition 31, then

$\lim_{N\to+\infty}\left(\frac{X_{2}(N^{p-1}t)}{N}\right)=(V(t))$

where $(V(t))$ is an explosive Markov process on $(0,1]$, with a
multiplicative structure, almost surely converging to $0$.

Starting from a large initial state of norm $N$, the scaled versions of the
processes for the CRNs of Sections 6.1 and 7 exhibit a rapid convergence to
$0$, exponential in several cases, and on slowed down timescales such as
$(t/N)$ or $(t/\sqrt{N})$.

In this example, to decrease the norm of the process, one has to speed-up the
timescale by a factor $N$ in (a) and $N^{p-1}$ in (b) and the decay in (a) is
only linear with respect to time. This is somewhat unusual in the current
mathematical literature of CRNs. It has in fact some analogies with the
scaling limits of queueing networks for which the timescale has to be sped by
$N$ and the limiting points of scaled sample paths are deterministic piecewise
linear functions. See [10] and Chapter 9 of [40].

In the limit, the first order in case (a) is a linear function on a finite
interval. In case (b), the CRN is described with an explosive Markov process
on $(0,1)$.

### Outline of the Paper

The general organization is as follows. In a first part we give a presentation
of CRNs, and general results on criteria of positive recurrence and on scaling
results. The second part, Sections 6 7 and 8, is devoted to the investigation
of examples. General mathematical results, i.e. for a large class of
topologies, on the stability and scaling properties of CRNs are certainly
important, but given the complexity of the kinetics of CRNs, the analysis of
specific examples is still a major source of progress to develop and validate
new tools to investigate CRNs.

Section 2 introduces general definitions for CRNs and the associated dynamical
systems. Classical results of the literature, including the deficiency zero
theorem, are recalled in Section 3.

In Section 4, a stability criterion involving a stopping time is introduced.
It is related to specific scaling schemes to study stability of CRNs, this is
discussed in Section 5. As it will be seen in subsequent sections, a scaling
approach can also be used to study the transient behavior of CRNs, via
convergence results of their sample paths.

The rest of the paper is devoted to the analysis of examples of binary CRNs. A
binary CRN is a CRN for which at most two chemical species are involved for
each transition of the Markov process. In Section 6, some general properties
of binary CRNs are presented.

The special case of a simple triangle topology for a binary CRN is
investigated in section 6.1. This case is not difficult but it illustrates
several aspects of the approach proposed in this paper: a simple Lyapunov
function taken at a convenient stopping time, and the boundary behavior of
this CRN leads to three types of scaling, depending on the initial state,
giving three different regimes for the time evolution of the state of the
network.

Section 7 is devoted to the analysis of a CRN proposed in [2] for which a
Lyapunov function is constructed to prove the stability of a version of the
CRN. We show that the criterion of Section 4 can be used with a simple
function to prove this result. A scaling picture for the time evolution of
this CRN is also presented. In Section 8, we study the stability and some
transient behaviors of the CRN mentioned above. For the sake of completeness,
in Section B of the Appendix we present some of the classical technical
arguments used, repeatedly, in the proofs of convergence in distribution of
sequence of processes. It includes the proof of functional laws of large
numbers and of an averaging principle. In Section C of the Appendix, a general
result for the stability of triangular networks with arbitrary complexes at
the vertices, a generalization of CRNs of Section 6.1, is shown. This is an
analogue of Theorem 7.6.1 of [18] for star networks.

##  2\. Mathematical Models of CRNs

We first recall the formal definitions for CRNs.

######  Definition 1.

A chemical reaction network (CRN) with $m$ chemical species, $m{\geq}1$, is
defined by a triplet $({\cal S},{\cal C},{\cal R})$,

  * —

${\cal S}{=}\\{1,\ldots,m\\}$ is the set of species;

  * —

${\cal C}$, the set of complexes, is a finite subset of $\mathbb{N}^{m}$;

  * —

${\cal R}$, the set of chemical reactions, is a subset of ${\cal C}^{2}$.

A chemical species of type $j{\in}{\cal S}$ is represented as $S_{j}$. A
complex $y{\in}{\cal C}$, $y{=}(y_{j})$ is composed of $y_{j}$ _molecules_ of
species $j{\in}{\cal S}$, its _size_ is $\|y\|{=}y_{1}{+}\cdots{+}y_{m}$. It
is also described as

$y=\sum_{j=1}^{m}y_{j}S_{j}.$

A chemical reaction $r{\in}{\cal R}$,
$r{=}(y_{r}^{-},y_{r}^{+}){=}\left(\left(y_{r,i}^{-}\right),\left(y_{r,i}^{+}\right)\right)){\in}{\cal
C}^{2}$ is represented as

$\sum_{i=1}^{m}y_{r,i}^{-}S_{i}\rightharpoonup\sum_{i=1}^{m}y_{r,i}^{+}S_{i}.$

The CRN can be seen as an oriented graph, called the _reaction graph_ , whose
vertices are complexes, i.e. in $\cal{C}$ and whose set of directed edges is
$\cal{R}$.

The state of the CRN is given by a vector
$x{=}(x_{i},1{\leq}i{\leq}m){\in}\mathbb{N}^{m}$, for $1{\leq}i{\leq}m$,
$x_{i}$ is the number of copies of chemical species $S_{i}$. A chemical
reaction $r{=}(y_{r}^{-},y_{r}^{+})$ corresponds to the change of state, for
$x{=}(x_{i})$,

(2) $x\longrightarrow x{+}y_{r}^{+}{-}y_{r}^{-},$

provided that $y_{r,i}^{-}{\leq}x_{i}$ holds for $1{\leq}i{\leq}m$, i.e. there
are at least $y_{r,i}^{-}$ copies of chemical species of type $i$, for all
$i{\in}{\cal S}$, otherwise the reaction cannot happen.

The notation $\emptyset$ refers to the complex associated to the null vector
of $\mathbb{N}^{m}$, $\emptyset{=}(0)$. For $y{=}(y_{i}){\in}{\cal C}$, a
chemical reaction of the type $(\emptyset,y)$ represents an external source
creating $y_{i}$ copies of species $i$, for $i{=}1$,…, $m$. A chemical
reaction of the type $(y,\emptyset)$ consists in removing $y_{i}$ copies of
species $i$, for $i{=}1$,…, $m$, provided that there are sufficiently many
copies of each species.

$S_{1}{+}S_{2}$$2S_{1}$$S_{2}$$\emptyset$$\kappa_{12}$$\kappa_{1}$$\kappa_{2}$$\kappa_{0}$$\kappa_{20}$

Figure 1. Example of a CRN

For the CRN of Figure 1, we have

  * —

${\cal S}{=}\\{1,2\\}$;

  * —

${\cal C}{=}\\{\emptyset,S_{1}{+}S_{2},2S_{1},S_{2}\\}$.

  * —

${\cal
R}{=}\\{(\emptyset,S_{2}),(S_{2},\emptyset),(S_{2},S_{1}{+}S_{2}),(S_{1}{+}S_{2},2S_{1}),(2S_{1},S_{2})\\}$.

###  2.1. Notations

Throughout the paper, we will use the following notations. For
$z{=}(z_{i}){\in}\mathbb{N}^{m}$ we denote

(3) $\|z\|{\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}}\sum_{i=1}^{m}|z_{i}|,\text{ and
}\|z\|_{\infty}{\stackrel{{\scriptstyle\text{\rm def.}}}{{=}}}\max_{1\leq
i\leq m}|z_{i}|.$

The quantity $|{\cal R}|$ is the cardinality of ${\cal R}$.

If $\kappa{=}(\kappa_{r},r{\in}{\cal R}){\in}\mathbb{R}_{+}^{|{\cal R}|}$ and
$y{=}(y_{r},r{\in}{\cal R}){\in}{\cal C}^{|{\cal R}|}$

(4) $\kappa_{\max}{=}\max_{r{\in}{\cal R}}\kappa_{r},\quad
y^{\pm}_{\max}{=}\max_{r{\in}{\cal R}}\|y^{\pm}_{r}\|,$

and we define the generalized factorials, for $z{=}(z_{i}){\in}\mathbb{N}^{m}$
and $y{=}(y_{i}){\in}{\cal C}$,

(5) $z!\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\prod_{i=1}^{m}z_{i}!,\quad
z^{(y)}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\frac{z!}{(z{-}y)!}=\prod_{i=1}^{m}\frac{z_{i}!}{(z_{i}{-}y_{i})!},$

with the convention that $z^{(y)}{=}0$, if there exists some $i_{0}{\in}{\cal
S}$ such that $y_{i_{0}}{>}z_{i_{0}}$.

A real-valued function $(x(t))$ on $\mathbb{R}_{+}$ is càdlàg if it is right
continuous and it has left-limits everywhere on $\mathbb{R}_{+}$, for $t{>}0$,
$x(t{-})$ denotes the left limit at $t$.

If ${\cal R}$ is a positive Borelian measure on $\mathbb{R}_{+}^{2}$ and $A$
is a Borelian subset of $\mathbb{R}_{+}$, we will use the following notation
for the differential term

${\cal
R}(A,\mathop{}\mathopen{}\mathrm{d}t)=\int\mathbbm{1}_{\left\\{s{\in}A\right\\}}{\cal
R}(\mathop{}\mathopen{}\mathrm{d}s,\mathop{}\mathopen{}\mathrm{d}t).$

Throughout the paper, the random variables ${\cal P}$ or $({\cal P}_{r}$,
$r{\in}{\cal R})$ on $\mathbb{R}_{+}^{2}$ are, implicitly, independent Poisson
processes with intensity measure the Lebesgue measure on $\mathbb{R}_{+}^{2}$.
They are used in the stochastic differential equations associated to CRNs.

###  2.2. Markov Processes

The dynamical behavior of a CRN, i.e. the time evolution of the number of
copies or the concentration of each of the $m$ chemical species is governed by
the law of mass action. See [44], [35] for surveys on the law of mass action
and the historical reference [21]. A vector $\kappa{=}(\kappa_{r},r{\in}{\cal
R})$ of positive numbers is added to the parameters of the model. For
$r{\in}{\cal R}$, $\kappa_{r}$ is the “rate” of the reaction $r$. We now
describe it in a deterministic context.

#### Deterministic Models

A deterministic model of a CRN $({\cal S},{\cal C},{\cal R})$ whose dynamics
are governed by the law of mass action is represented by a dynamical system
$(x(t)){=}(x_{i}(t),i{=}1,\ldots,m)$ with state space $\mathbb{R}_{+}^{m}$,
satisfying the set of ODEs,

(6) $\dot{x}(t)=\sum_{r{=}(y_{r}^{-},y_{r}^{+}){\in}{\cal
R}}\kappa_{r}\left(\prod_{i=1}^{m}x_{i}(t)^{y_{r,i}^{-}}\right)\left(y_{r}^{+}{-}y_{r}^{-}\right),$

It should be noted that, due to a possible lack of a global Lipschitz
property, a solution of ODE (6) may be defined on a finite interval only. The
solution $(x(t))$ may blow-up in finite time, i.e. one of its coordinates may
converge to infinity when $t$ converges to some finite $t_{\infty}$ from
below. In this case, we will use the convention that $x(t)$ is defined as
${\dagger}$, for $t{\geq}t_{\infty}$. The state ${\dagger}$ is the “point at
infinity”, the state space is in fact
$\mathbb{R}_{+}^{m}{\cup}\\{{\dagger}\\}$.

The stability property of a deterministic CRN refers to the local/global
stability properties of fixed points of the dynamical system 6. For the CRN of
Figure 1, this gives the ODEs,

$\displaystyle\dot{x}_{1}(t)$
$\displaystyle=\kappa_{2}\,x_{2}(t){+}\kappa_{12}\,x_{1}(t)x_{2}(t){-}2\kappa_{1}\,x_{1}(t)^{2},$
$\displaystyle\dot{x}_{2}(t)$
$\displaystyle=\kappa_{0}{-}\kappa_{20}\,x_{2}(t){+}\kappa_{1}\,x_{1}(t)^{2}{-}\kappa_{12}\,x_{1}(t)x_{2}(t).$

In a stochastic context, a mathematical model of a CRN $({\cal S},{\cal
C},{\cal R})$ is represented by a continuous time Markov jump process
$(X(t)){=}(X_{i}(t),i{=}1,\ldots,m)$ with state space $\mathbb{N}^{m}$. The
associated $Q$-matrix is given by, for $x{\in}\mathbb{N}^{m}$ and
$r{=}(y_{r}^{-},y_{r}^{+}){\in}{\cal R}$, the transition

(7) $x\longrightarrow x{+}y_{r}^{+}{-}y_{r}^{-},\text{occurs at rate
}\kappa_{r}x^{(y_{r}^{-})}{=}\kappa_{r}\frac{x!}{(x-y_{r}^{-})!},$

with Notation (5). Throughout the paper, it will be assumed that the process
starts and lives in $\cal{E}_{0}$, a subset of $\mathbb{N}^{m}$ where it is
irreducible.

The stability property of a stochastic CRN refers in general to the positive
recurrence property of $(X(t))$ on ${\cal E}_{0}$ for a given positive vector
$(\kappa_{r},r{\in}{\cal R})$.

The functional operator ${\cal Q}(f)$ associated to its $Q$-matrix is defined
by, for $x{\in}\mathbb{N}^{m}$,

(8) ${\cal Q}(f)(x)=\sum_{r{\in}{\cal
R}}\kappa_{r}x^{(y_{r}^{-})}\left(f\left(x{+}y_{r}^{+}{-}y_{r}^{-}\right){-}f(x)\right),$

for any function $f$ with finite support on $\mathbb{N}^{m}$.

Section B.1 of the Appendix gives a representation of the Markov process as
the solution of a Stochastic Differential Equation (SDE). This representation
is especially useful to establish functional limit theorems.

For the CRN of Figure 1, the non-trivial components of the $Q$-matrix of the
associated Markov process $(X(t)){=}(X_{1}(t),X_{2}(t))$ are given by

$(x_{1},x_{2})\longrightarrow(x_{1},x_{2}){+}\begin{cases}(0,1),&\kappa_{0},\\\
(0,{-}1),&\kappa_{20}x_{2},\end{cases}\qquad(x_{1},x_{2}){+}\begin{cases}(1,0),&\kappa_{2}\,x_{2},\\\
({-}2,1),&\displaystyle\kappa_{1}\,x_{1}(x_{1}{-}1),\\\
(1,{-}1),&\kappa_{12}\,x_{1}x_{2}.\end{cases}$

###  An Important CRN: The $\mathbf{M/M/\infty}$ queue

This is a simple CRN with an external input and one chemical species,

$\emptyset\mathrel{\mathop{\xrightleftharpoons[\mu]{\lambda}}}S_{1}.$

The ${M/M/\infty}$ queue with input parameter $\lambda{\geq}0$ and output
parameter $\mu{>}0$ is a Markov process $(L(t))$ on $\mathbb{N}$ with
transition rates

$x\longrightarrow\begin{cases}x{+}1&\lambda\\\ x{-}1&\mu x.\end{cases}$

The invariant distribution of $(L(t))$ is Poisson with parameter
$\rho{=}\lambda/\mu$.

This fundamental process can be seen as a kind of discrete Ornstein-Uhlenbeck
process. It has a long history, it has been used in some early mathematical
models of telephone networks at the beginning of the twentieth century, see
[15], also in stochastic models of natural radioactivity in the 1950’s, see
[24] and it is the basic process of mathematical models communication networks
analyzed in the 1970’s, see [38]. See Chapter 6 of [40].

Technical results on this stochastic process turn out to be useful to
investigate the scaling properties of some CRNs and, as we will see, in the
construction of couplings used in our proofs. See Sections 7 and B.2.2 for
example.

### Remarks

  1. (a)

Boundary Behaviors.
An important difference with deterministic CRNs is that there are boundary
effects, i.e. discontinuities in the dynamics. The transition
$x{\to}x{+}y^{+}{-}y^{-}$, i.e. the reaction

$\sum_{i=1}^{m}y_{i}^{-}S_{i}\rightharpoonup\sum_{i=1}^{m}y_{i}^{+}S_{i},$

occurs with a positive rate only if there are at least $y_{i}^{-}$ copies of
species $i$ for all $1{\leq}i{\leq}m$, otherwise its rate is $0$. Note that
there is no such constraint in the deterministic setting, see Relation (6). In
practice, boundary effects complicates significantly the analysis of
stochastic CRNs.

  2. (b)

Blow-Up Phenomena.
As in the deterministic case, an important issue concerning CRNs is the fact
that the system can be explosive. The proof of stability criteria have to take
into account this (possible) phenomenon. To illustrate that, we consider the
CRN with one chemical species and one chemical reaction

$2S_{1}\stackrel{{\scriptstyle 1}}{{\rightharpoonup}}3S_{1}.$

In a deterministic setting, Relation (6) is in this case
$\dot{x}(t)=x(t)^{2}$, $t{\geq}0$ whose solution starting from $1$ is
$1/(1{-}t)$ blows-up at $t{=}1$. In a stochastic context, the process $(X(t))$
starting from $X(0){=}2$ reaches infinity in a finite (random) amount of time
with the same distribution as

$\sum_{n{\geq}2}\frac{E_{n}}{n(n-1)},$

where $(E_{n})$ is an i.i.d. sequence of exponential random variables with
parameter $1$.

As in Section 2.2, a point ${\dagger}$ is added to the state space. If the
sequence of successive jumps of the process $(X(t))$ converges to some
(possibly) finite random variable $T_{\infty}$, we will use the convention
$X(t){=}{\dagger}$ for $t{\geq}T_{\infty}$.

###  2.3. Deterministic CRNs as a Scaled Version of Stochastic CRNs

It is possible to see the set of ODEs (6) with the parameters $(\kappa_{r})$
as a first order approximation of a scaled stochastic model of a CRN with the
same triplet $({\cal S},{\cal C},{\cal R})$ but with a set of scaled reaction
rates

(9) $\left(\kappa_{r}^{N},r{\in}{\cal
R}\right)=\left(\frac{\kappa_{r}}{N^{\|y_{r}^{-}\|{-}1}},r{\in}{\cal
R}\right).$

We denote by $(X^{N}(t)){=}(X_{i}^{N}(t))$ the associated Markov jump process.

This ad hoc scaling is not homogeneous in the sense that it does not
correspond to a uniform change of timescale. The timescale of reaction
$r{\in}{\cal R}$ is “slowed down” by a factor $1/N^{\|y_{r}^{-}\|{-}1}$. In
this manner, provided that all its coordinates are of the order of $N$, the
rate of any chemical reaction is of the same order, $O(1)$ with respect to
$N$. Consequently, this scaling removes an important feature of CRNs, multi-
timescales, i.e. that some reactions may be much faster than others. The main
scaling considered in this paper is with respect to the size of the initial
state, with a possible change of timescale, but independent of the chemical
reaction considered, the topology and the vector $\kappa$ of reaction rates
are unchanged.

The following proposition is a well-known result of the CRN literature for the
scaling assumption of Relation (9). See [37] for example.

######  Proposition 2.

If $(X^{N}(t))$ is the Markov process with $Q$-matrix given by Relation (7)
with the reaction rates given by Relations (9), if the sequence of initial
states is such that

$\lim_{N\to+\infty}\left(\frac{X_{i}^{N}(0)}{N}\right)=x_{0}{=}(x_{0,i}){\in}\mathbb{R}_{+}^{m},$

then the convergence in distribution

$\lim_{N\to+\infty}\left(\frac{X_{i}^{N}(t)}{N},t{\in}(0,t_{\infty})\right)=(x_{i}(t),t{\in}(0,t_{\infty})),$

holds, where $(x_{i}(t))$ is the solution of the set of ODEs (6) starting at
$(x_{0})$ and $t_{\infty}$ is its blow-up instant,

$t_{\infty}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\lim_{K\to+\infty}\inf\\{t{:}\|x(t)\|_{\infty}{\geq}K\\},$

with the convention that $\inf\emptyset{=}{+}\infty$.

###### Proof.

We give a quick proof of this classical result. For $K{\geq}1$, define

$H^{N}_{K}=\inf\left\\{t{>}0{:}\|X^{N}(t)\|_{\infty}{\geq}KN\right\\}\text{
and }t_{K}=\inf\left\\{t{>}0{:}\|x(t)\|_{\infty}{\geq}K\right\\},$

and

$\left(\overline{X}_{K}^{N}(t)\right)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\frac{X_{i}^{N}(t{\wedge}H^{N}_{K})}{N}\right).$

On the time interval $(0,H^{N}_{K})$, there are at least $C_{K}N$ jumps of the
process $(X_{N}(t))$, with
$C_{K}{=}\lfloor(K{-}\|x_{0}\|)/y_{\max}^{+}\rfloor$. The maximal jump rate of
the process on $(0,H^{N}_{K})$ is bounded by $\lambda_{\max}N$ with
$\lambda_{\max}{=}\kappa_{\max}K^{y_{\max}^{-}}|{\cal R}|$. Hence,

$\lim_{N\to+\infty}\mathbb{P}\left(H^{N}_{K}{\leq}t\right)=0,$

holds for any $t{<}C_{K}/\lambda_{\max}$.

We now use the SDE formulation of Section B.1 of the appendix to represent the
Markov process $(X^{N}(t))$. Relation (61) gives the identity

(10) $\overline{X}^{N}(t){=}\overline{X}^{N}(0){+}\sum_{r{\in}{\cal
R}}\frac{M_{r}(t{\wedge}H^{N}_{K})}{N}{+}\kappa_{r}\left(y_{r}^{+}{-}y_{r}^{-}\right)\int_{0}^{t{\wedge}H^{N}_{K}}\frac{X^{N}(s)!}{N^{\|y_{r}^{-}\|}(X^{N}(s)-y_{r}^{-})!}\mathop{}\mathopen{}\mathrm{d}s,$

where $(M_{r}(t)){=}((M_{r,i}(t)))$, $r{\in}{\cal R}$, is the set of
martingales defined by Relation (64). Relation (65) gives, for
$1{\leq}i{\leq}m$,

$\left\langle\sum_{r{\in}{\cal
R}}\frac{M_{r,i}}{N}\right\rangle(t{\wedge}H^{N}_{K})=\sum_{r{\in}{\cal
R}}\frac{\left\langle M_{r,i}\right\rangle}{N^{2}}(t{\wedge}H^{N}_{K})\\\
=\frac{1}{N}\sum_{r{\in}{\cal
R}}\left(y_{r,i}^{+}{-}y_{r,i}^{-}\right)^{2}\kappa_{r}\int_{0}^{t{\wedge}H^{N}_{K}}\frac{X(s)!}{N^{\|y_{r}^{-}\|}(X(s){-}y_{r}^{-})!}\mathop{}\mathopen{}\mathrm{d}s\\\
\leq\frac{1}{N}\sum_{r{\in}{\cal
R}}\left(y_{r,i}^{+}{-}y_{r,i}^{-}\right)^{2}\kappa_{r}K^{\|y_{r}^{-}\|}t,$

with Doob’s Inequality, we get that $(M_{r,i}(t{\wedge}H^{N}_{K})/N)$ is
converging in distribution to $0$, for all $1{\leq}i{\leq}m$. By using the
criterion of modulus of continuity, see [9], it is easy to show that the
sequence of processes $(\overline{X}_{K}^{N}(t))$ is tight for the convergence
in distribution.

For $t{<}t_{K}$, with high probability we have $H^{N}_{K}{\geq}t$, hence by
using Relation (10), we obtain that any limiting point of
$(\overline{X}_{K}^{N}(t),t{<}t_{K})$ satisfies the set of ODEs (6). In
particular the sequence $(\|X_{N}(t)\|,t{<}t_{K})$ is converging in
distribution to $(\|x(t)\|,t{<}t_{K})$. The proposition is proved. ∎

##  3\. Classical Stability Results

### Deterministic CRNs

We first state the classical results of [17] and [27] on the stability
properties of deterministic CRNs. See [18] for a broader picture of stability
results for CRNs. See also [22] for a quick, comprehensive, overview of these
results.

### Deficiency Zero Theorem

We briefly recall some definitions on CRNs.

  1. (a)

A linkage class is a connected component of the reaction graph. The quantity
$\ell$ denotes the number of linkage classes.

  2. (b)

The CRN is weakly reversible if every connected component of its reaction
graph is strongly connected.

  3. (c)

The Stoichiometric space $S$ is the vector subspace of $\mathbb{R}^{m}$
generated by $y_{r}^{+}{-}y_{r}^{+}$, $r{\in}{\cal R}$, its dimension is
denoted by $s$.

  4. (d)

The deficiency $\delta$ of the CRN is $|C{-}\ell{-}s|$, where $C$ is the total
number of complexes.

The main theorem can now be stated.

######  Theorem 3 (Feinberg (1979)).

Let $({\cal S},{\cal C},{\cal R})$ be a chemical reaction network with
deterministic mass action kinetics, if it is weakly reversible and with zero
deficiency, $\delta{=}0$, then there is exactly one equilibrium for the
dynamical system defined by Relation (6) in an irreducible subset of the state
space. Furthermore, in such a subset this equilibrium is locally
asymptotically stable.

This result shows that the CRN has a unique fixed point under a condition on
its topological structure, but independently of the set of coefficients
$\kappa{=}(\kappa_{r},r{\in}{\cal R})$.

###  3.1. Stochastic CRNs

Proposition (2) shows that a deterministic CRN can be seen as an asymptotic
stochastic CRN provided that the reaction rates are scaled conveniently as in
Relations (9). The interesting results of [6] show that the invariant
distribution of some stochastic CRNs can be expressed with the equilibrium of
a deterministic CRN. See also [5].

######  Theorem 4.

Let $(X(t))$ be the Markov process associated to a stochastic chemical
reaction network, irreducible on $\cal{E}_{0}$, whose $Q$-matrix is given by
Relation (8), if the CRN is weakly reversible and has a deficiency equal to
$0$, then, if $c{=}(c_{i})$ is the equilibrium point of Relation (6), the
positive measure on $\cal{E}_{0}$ defined by

(11) $\pi(x)=\prod_{i=1}^{m}\frac{c_{i}^{x_{i}}}{x_{i}!}$

is an invariant measure of $(X(t))$.

###  3.2. Analogies with Classical Queueing Networks

Theorem 4 states that the invariant distribution for the Markov processes of a
class of CRNs has a product form associated to the solution of a deterministic
set of equations.

A related result holds for a class of multi-dimensional Markov processes
associated to a class of queueing networks, referred to as Jackson Networks.
They can be described simply as follows. See also [6] for a discussion of
these aspects.

  1. (a)

There are $J$ sites for the location of particles. If
$x{=}(x_{j}){\in}\mathbb{N}^{J}$, for $1{\leq}j{\leq}m$, $x_{j}$ denotes the
number of particles at the site $j$.

  2. (b)

One of the particles at site $j$ leaves at rate $\mu_{j}{>}0$. It goes to site
$k{\in}\\{1,\ldots,J\\}$ with probability $p_{jk}$, or leave the network with
probability

$p_{j0}=1{-}\sum_{k=1}^{J}p_{jk}.$

  3. (c)

External particles arrive at the site $j$ at rate $\lambda_{j}{\geq}0$.

The matrix $P{=}(p_{jk})$ is the routing matrix, if its maximal eigenvalue is
${<}1$, then there exists a unique (non-negative) solution for the linear
system

(12) $y_{j}=\lambda_{j}{+}\sum_{k=1}^{J}y_{k}p_{kj}.$

The result due to [28],it shows that the sequence
$(\pi(x),x{\in}\mathbb{N}^{J})$ defined by

(13) $\pi(x)=\prod_{j=1}^{J}\left(\frac{y_{j}}{\mu_{j}}\right)^{x_{j}},\quad
x{=}(x_{j}){\in}\mathbb{N}^{J},$

is an invariant distribution. See also [33]. Relation (12) is the equivalent
of the equilibrium equation for the dynamical system defined by Relation (6)
and Relation (13) is the analog of Relation (11). Note that a condition has to
be added for the summability of $(\pi(x))$ defined by (13).

It is remarkable that in both cases, the invariant distribution $(\pi(x))$
satisfies a set of equations in the following way: for any $x{\in}{\cal S}$,
there exists a partition $({\cal A}_{p},p{\in}I)$ of the possible states of
the Markov process starting from $x$, such that, for any $p{\in}I$, the
relation

$\sum_{x^{\prime}{\in}{\cal
A}_{p}}\pi(x^{\prime})q(x^{\prime},x)=\pi(x)\sum_{x^{\prime}{\in}{\cal
A}_{p}}q(x,x^{\prime})$

holds. These relations are usually referred to as local balance equations. By
summing these identities for $p{\in}I$, this gives that $(\pi(x))$ satisfies
the equilibrium equation,

$\sum_{x^{\prime}{\in}{\cal
S}:x^{\prime}{\neq}x}\pi(x^{\prime})q(x^{\prime},x)=\pi(x)\sum_{x^{\prime}{\in}{\cal
S}:x^{\prime}{\neq}x}q(x,x^{\prime}).$

  1. (a)

For Jackson networks, there are $J{+}1$ subsets of transitions associated to
an internal arrival/departure at node $j$, $j{=}1,\ldots,J$ and a subset
related external arrival/definitive departure of particles. See Chapter 3 of
[33].

  2. (b)

For CRNs, the subsets are indexed by complexes, each one corresponding to
transitions involving a given complex. See Theorem 3.2 of [6].

##  4\. A Stability Criterion

The goal of this section is of formulating a criteria of positive recurrence
for Markov processes associated to CRNs. In this domain the literature is
quite rich for discrete time Markov chains. For continuous time processes,
these results can be used by studying the embedded Markov chain. This approach
can be cumbersome in a multi-dimensional state space, when there are many
possible transitions starting from a given state. Another specific
characteristic of continuous time case is the explosion phenomenon mentioned
above. If there is explosion, clearly the Markov process is transient, the
main problem is that it leads to complications in the design of a criterion of
positive recurrence in continuous time, since the process can die in finite
time.

In this section, $(X(t))$ is an irreducible Markov process on a subset ${\cal
E}_{0}$ of $\mathbb{N}^{m}$ whose $Q$ matrix is given by Relation (7). It is
assumed that on the probability space, there is a semi-group of shift
operators $(\theta_{t})$ so that the relation
$X(t{+}s){=}X(t){\circ}\theta_{s}$ holds almost surely for all $t$,
$s{\geq}0$. The natural filtration of $(X(t))$ is denoted by $({\cal F}_{t})$.
See Chapter I of [43] for a general presentation of this formalism for Markov
processes.

We define

(14) $t_{1}=\inf\\{s{>}0:X(s)\neq X(s{-})\\},$

the first instant of jump of $(X(t))$ and for $n{\geq}1$,

$t_{n+1}=\inf\\{s{>}t_{n}:X(s)\neq
X(s{-})\\}=t_{n}{+}t_{1}{\circ}\theta_{t_{n}},$

the sequence of the instants of successive jumps of the process. The Markov
process is non-explosive if and only if the sequence $(t_{n})$ is almost
surely converging to infinity.

######  Definition 5.

An energy function $f$ on ${\cal E}_{0}$ is a non-negative function such that,
for all $K{>}0$, the set $\\{x{\in}{\cal E}_{0}:f(x){\leq}K\\}$ is finite.

By convention the value of an energy function at the point at infinity
${\dagger}$ is ${+}\infty$. The classical energy functions used in a CRN
framework can be linear functions or an entropy function, for $x{=}(x_{i})$,

$f(x)=\sum_{i=1}^{m}a_{i}x_{i},\quad f(x)=\sum_{i=1}^{m}(x_{i}\ln
x_{i}{-}x_{i}{+}1).$

The following theorem is an adaptation to continuous time Markov processes of
a simple result for discrete time Markov chains, see [20], and another version
in the context of queueing networks, Theorem 8.13 of [40]. An additional
difficulty for CRNs in the statement of a stability criterion is the
possibility of explosion of the associated stochastic processes.

######  Theorem 6.

Let $(X(t))$ be an irreducible Markov process on ${\cal E}_{0}$ associated to
a CRN network with $Q$-matrix (7), if there exist

  1. (a)

an integrable stopping time $\tau$ such that $\tau{\geq}t_{1}{\wedge}\eta$,
for a constant $\eta{>}0$ and $t_{1}$ is the first jump of $(X(t))$, Relation
(14);

  2. (b)

an energy function $f$ on ${\cal E}_{0}$ and constants $K$ and $\gamma{>}0$
such that the relation

(15)
$\mathbb{E}_{x}\left(f(X(\tau))\right){-}f(x)\leq{-}\gamma\mathbb{E}_{x}(\tau),$

holds for all $x{{\in}}{\cal E}_{0}$ such that $f(x){\geq}K$,

then $(X(t))$ is a positive recurrent Markov process.

A function $f$ satisfying Condition (15) is usually referred to as a Lyapunov
function.

Note that if equation (15) holds, $\tau<T_{\infty}$ with $T_{\infty}$ the time
of explosion (possibly equal to $+\infty$), and the proof of the theorem shows
in particular that the process cannot explode.

As a consequence we state the classical positive recurrence criterion
involving the $Q$-matrix $Q$ of the Markov process.

######  Corollary 7.

Let $(X(t))$ be an irreducible Markov process on ${\cal E}_{0}$ associated to
a CRN network with $Q$-matrix $Q$ defined by Relation (7), if there exists an
energy function $f$ on ${\cal S}$ and $\gamma$, $K{>}0$ such that

(16) $Q(f)(x){\leq}{-}\gamma,\text{ if }f(x){\geq}K,$

then $(X(t))$ is positive recurrent.

###### Proof.

For $\tau=t_{1}$, and $x$ such that $f(x){\geq}K$, we have

$\mathbb{E}_{x}\left(f(X(t_{1}))-f(x)\right)=\mathbb{E}_{x}\left(\int_{0}^{t_{1}}Q(f)(X(s))\mathop{}\mathopen{}\mathrm{d}s\right)=Q(f)(x)\mathbb{E}_{x}(t_{1})\leq-\gamma\mathbb{E}_{x}(t_{1}),$

then we conclude with Theorem 6. ∎

The classical Condition (16) states that if the initial state has a large
level of energy, the next step decreases in average the energy by some fixed
positive quantity. There are various versions of this type of result. One of
them evaluates the expected value of $f(X(t_{0}))$ at a deterministic function
$t_{0}(x)$ of the initial state, see Proposition 4.5 of [10]. Finding a global
Lyapunov function may be cumbersome in general. It may imply to partition the
state space, for the location of $X(t_{0})$, to build and to glue piecewise,
local, Lyapunov functions. See [2] for a typical example. If the stopping time
$\tau$ of Theorem 6 is well chosen, the location of $X(\tau)$ is not a
concern. For this theorem, one has to partition the state space from the point
of view of the initial state, to define $\tau$, instead of the state at time
$t_{0}$. This is in general much simpler to handle.

Starting from a state $x$ of high energy does not necessarily lead quickly to
a lower energy level, i.e. $Q(f)(x)$ is not necessarily negative. Instead, it
may happen that one may have to wait for some amount of time, the quantity
$\tau$, before the energy can significantly decrease. The “price” to pay for
waiting is that the decay needed is not anymore some fixed negative quantity,
but, as it can be expected, a negative constant proportional to
$\mathbb{E}_{x}(\tau)$, the mean value of the waiting time.

Remark.
For the somewhat specific condition $\tau{\geq}t_{1}{\wedge}\eta$, it should
be noted that, at least some condition is required, since the trivial
$\tau{\equiv}0$ would work otherwise. In practice, it is quite natural that
the “convenient” stopping time $\tau$ should be greater than $t_{1}$, the
first instant when there is a change of state for the Markov process. There
are, however, situations when it is enough to consider a deterministic
stopping time $\tau$, when scaling limits of sample paths are used. Some of
our examples below exhibit this feature.

  1. (a)

If $\tau{\equiv}\eta{>}0$.
Condition (15) is just the classical Foster-Lyapunov condition for the
discrete Markov chain $(M_{n}){=}(X(n\eta))$, i.e. the process $(X(t))$ on the
discrete timescale $n{\mapsto}n\eta$,

(17) $\mathbb{E}_{x}(f(M_{1})){-}f(x)\leq{-}\gamma\eta,\quad\text{ for
}x{\in}{\cal S}_{0}\text{ such that }f(x){>}K.$

See [23], [10], and [36].

  2. (b)

If $\tau{=}t_{1}$.
Condition (15) is equivalent to condition (16), see the proof of corollary 7.
Most of stability analyses of CRN networks use this condition.

###### Proof of Theorem 6.

The proof uses essentially the same arguments as in the proof of Theorem 8.6
of [40]. The only specific difficulty lies in the fact that, a priori, the
Markov process can be explosive. Define the sequence of induced stopping times
$(s_{n})$ by induction, by $s_{1}{=}\tau$ and

(18) $s_{n+1}=s_{n}{+}\tau{\circ}\theta_{s_{n}}.$

By using the strong Markov property of $(X(t))$, see Theorem (9.4), Section
III.9 of [41], $(s_{n})$ is a non-decreasing sequence of stopping times. From
our assumption on $\tau$, we get that $(s_{n})$ is almost surely an increasing
sequence, i.e. $s_{n}{<}s_{n+1}$ for all $n{\geq}1$.

We define

$T_{K}=\inf\\{s{\geq}0:f(X(s))\leq K\\},\text{ and
}\nu\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\inf\\{n{\geq}0:f(X(s_{n}))\leq K\\},$

then, clearly $T_{K}{\leq}s_{\nu}$. Let, for $n{\geq}1$,

$Z_{n}\stackrel{{\scriptstyle\text{\rm def.}}}{{=}}f(X(s_{n}))+\gamma s_{n},$

then, with Relation (15) and the strong Markov property of $(X(t))$ for the
stopping time $s_{n}$, we obtain the relation

(19) $\mathbb{E}\left(Z_{n+1}\mid{\cal
F}_{s_{n}}\right)=Z_{n}{+}\mathbb{E}_{X(s_{n})}\left(f(X(\tau)){+}\gamma\tau\right){-}f(X(s_{n}))\leq
Z_{n},$

on the event $\\{\nu{>}n\\}$.

The process $(Z_{n\wedge\nu})$ is therefore a non-negative super-martingale,
in particular it is converging almost surely to a finite limit.

First, lets show that $\nu$ is almost surely finite. The almost sure
convergence of $(Z_{n\wedge\nu})$ gives that, almost surely on the event
$\\{\nu{=}{+}\infty\\}$, the sequence $(s_{n})$ is converging to a finite
limit. In particular the increments $(s_{n+1}{-}s_{n})$ are less than $\eta/2$
after some finite index $n_{0}$. Our assumption, $\tau{\geq}t_{1}{\wedge}\eta$
implies therefore that, for $n{>}n_{0}$,

(20) $s_{n+1}{\geq}s_{n}{+}t_{1}{\circ}\theta_{s_{n}},$

and, by induction, $s_{n}{\geq}t_{n-n_{0}}{\circ}\theta_{s_{n_{0}}}$ holds for
$n{>}n_{0}$.

Since $X(s_{n}){=}X(u_{n})$, where $u_{n}{\leq}s_{n}$ is the last jump instant
of $(X(t))$ before $s_{n}$, we have $u_{n}{<}u_{n+1}$ by Relation (20). Hence,
almost surely, on the event $\\{\nu{=}{+}\infty\\}$, the Markov process
$(X(t))$ explodes in finite time, so that

$\limsup_{n\to+\infty}f(X(u_{n}))=\limsup_{n\to+\infty}f(X(s_{n}))={+}\infty.$

Since $(Z_{n{\wedge}\nu})$ is converging almost surely to a finite limit, this
implies that the random variable $\nu$ is almost surely finite.

By integrating Relation (19) we get that

$\gamma\mathbb{E}_{x}(s_{n\wedge\nu})\leq\mathbb{E}_{x}\left(Z_{n\wedge\nu}\right)\leq\mathbb{E}_{x}(Z_{0})=f(x),$

the monotone convergence theorem gives the relation

$\mathbb{E}(T_{K})\leq\mathbb{E}(s_{\nu})\leq\frac{f(x)}{\gamma},$

we conclude with Proposition 33 of the appendix. The theorem is proved. ∎

##  5\. Scaling Methods

The scaling approaches presented in this section aim at providing a first
order description of the time evolution of a CRN, to investigate its stability
properties and its transient behavior.

### Classical Scalings

Proposition 2 gives an example of such a first order description proves such a
limiting result, but it is done with a significantly modified CRN, via a
specific scaling of reaction rates. If the coordinates of $(X_{j}(0))$ are of
the order of $N$ large, the constant $(\kappa_{r})$ are scaled by a power of
$N$ so that all reaction rates are of the same order.

Another, somewhat related, scaling has been considered by T. Kurtz and co-
authors to investigate CRNs with a hierarchy of timescales. Again the
constants $(\kappa_{r})$ of some reactions are assumed to have a convenient
polynomial growth with respect to a scaling parameter $N$. See [8], [31] and
[32].

The scaling considered in this paper does not modify the internal dynamics of
the CRN, the $\kappa$’s are not scaled in particular.

### Scaling with the Norm of the Initial State

The scaling considered is with respect to the norm of the initial state of the
Markov process. As we will see this scaling is inspired by Theorem 6 and it
can be used to prove stability results. At the same time it may give a first
order description of how the CRN behaves starting from a “large” state. This
scaling has been used repeatedly to investigate queueing networks. See [10]
and Chapter 9 of [40].

The following proposition is at the origin of the scaling approach to
investigate positive recurrence of CRNs.

######  Proposition 8.

Let $(X(t))$ be a Markov process associated to a CRN network with $Q$-matrix
defined by Relation (7), if there are an energy function $f$ on ${\cal
E}_{0}$, constants $\varepsilon$ and $\eta{>}0$, and an integrable stopping
time $\tau{\geq}t_{1}{\wedge}\eta$ such that the relations

$\limsup_{\begin{subarray}{c}x{\in}{\cal E}_{0}\\\
f(x)\to+\infty\end{subarray}}\frac{\mathbb{E}_{x}(f(X(\tau)))}{f(x)}\leq
1{-}\varepsilon\text{ and }C_{0}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\limsup_{\begin{subarray}{c}x{\in}{\cal E}_{0}\\\
f(x)\to+\infty\end{subarray}}\frac{\mathbb{E}_{x}(\tau)}{f(x)}<{+}\infty$

hold, then the Markov process $(X(t))$ is positive recurrent.

###### Proof.

One can find $K{>}0$ such that if $x{\in}{\cal E}_{0}$ satisfies $f(x){\geq}K$
then

$\mathbb{E}_{x}(f(X(\tau))){-}f(x)\leq-\frac{\varepsilon}{2}f(x)\leq-\frac{\varepsilon}{4C_{0}}\mathbb{E}_{x}(\tau),$

Theorem 6 concludes the proof. ∎

If we take $f$ as the norm $\|{\cdot}\|$, the above proposition suggests the
introduction of the scaled processes, for $x{\in}{\cal E}_{0}$,

$\left(\overline{X}_{x}(t)\right)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\frac{X\left(tg(x)\right)}{\|x\|},t{\geq}0\right),$

where $g$ is a positive function on ${\cal S}$.

Let $g$ such that $g(x){\leq}C_{0}\|x\|$, $x{\in}{\cal S}$, for some
$C_{0}{>}0$. Assume that, as $\|x\|$ goes to infinity, the family of processes
$(\overline{X}_{x}(t))$ converges in distribution to some deterministic
process $(\overline{x}(t))$, and that there exists some
$t_{0}{\in}\mathbb{R}_{+}$ such that $\|\overline{x}(t_{0})\|{\leq}\eta$, for
$\eta{\in}[0,1)$. The above proposition can be used to prove the positive
recurrence of the Markov process. Under appropriate uniform integrability
conditions, take $\tau{=}t_{0}g(x)$, then the relation
$\mathbb{E}(\|X(t_{0}g(x))\|)/\|x\|{<}\eta_{1}{\in}(0,1)$ holds if $\|x\|$ is
sufficiently large. See Corollary 9 for a related formal statement. The
stability of a CRN can be investigated by analyzing the stability of the
dynamical system $(\overline{x}(t))$, with an initial point on the unit sphere
of $\mathbb{R}_{+}^{m}$. If at some fixed instant $t_{0}$, it is strictly
inside the unit ball, then, modulo some technical conditions, the Markov
process $(X(t))$ is positive recurrent.

This is the general motivation for using a scaling approach in this context.
In practice, it may be more complicated, because of polynomial reaction rates
and boundary behaviors already mentioned, and of different timescales
involved. Still, it gives a kind of simple recipe of how to investigate these
questions in a first step. In our experience, it has the merit of
concentrating the efforts of the investigation on difficult aspects (like
boundary effects) of a given Markov process associated to a CRN.

The second motivation is that, additionally, the asymptotic evolution of
$\left(\overline{X}_{x}(t)\right)$ gives a first order description of the time
evolution of the CRN. In addition to a positive recurrence property, this may
give insight on the dynamical behavior of these networks.

######  Corollary 9.

Let $(X(t))$ be a Markov process associated to a CRN network $({\cal S},{\cal
C},{\cal R})$ with parameters $(\kappa_{r})$, if there exist
$\varepsilon{>}0$, $t_{0}{>}0$ and $\gamma{\geq}{-}1$ such that

$\limsup_{\|x\|\to+\infty}\frac{1}{\|x\|}\mathbb{E}_{x}\left(\left\|X\left({t_{0}}/{\|x\|^{\gamma}}\right)\right\|\right)\leq
1{-}\varepsilon,$

then $(X(t))$ is positive recurrent.

###### Proof.

The only thing lacking is condition (a) of Theorem 6. For $x{\in}{\cal
E}_{0}$, we define $\tau{=}t_{1}{\vee}(t_{0}/\|x\|^{\gamma})$, $\tau$ is
clearly a stopping time and

$\frac{1}{\|x\|}\mathbb{E}_{x}\left(\left\|X\left(\tau\right)\right\|\right)\leq\frac{1}{\|x\|}\mathbb{E}_{x}\left(\|X\left(t_{1}\right)\|\mathbbm{1}_{\left\\{t_{0}/\|x\|^{\gamma}<t_{1}\right\\}}\right)\\\
{+}\frac{1}{\|x\|}\mathbb{E}_{x}\left(\|X\left({t_{0}}\left/{\|x\|^{\gamma}}\right.\right)\|\mathbbm{1}_{\left\\{t_{0}/\|x\|^{\gamma}\geq
t_{1}\right\\}}\right).$

With Notation (4), since $\|X(t_{1})\|{\leq}\|x\|{+}y^{+}_{\max}$, we obtain
the relation

$\displaystyle\frac{1}{\|x\|}\mathbb{E}_{x}\left(\left\|X\left(\tau\right)\right\|\right)$
$\displaystyle\leq\frac{y^{+}_{\max}}{\|x\|}{+}\frac{1}{\|x\|}\mathbb{E}_{x}\left(\|x\|\mathbbm{1}_{\left\\{t_{0}/\|x\|^{\gamma}<t_{1}\right\\}}\right)$
$\displaystyle\hskip
85.35826pt+\frac{1}{\|x\|}\mathbb{E}_{x}\left(\|X\left({t_{0}}\left/{\|x\|^{\gamma}}\right.\right)\|\mathbbm{1}_{\left\\{t_{0}/\|x\|^{\gamma}\geq
t_{1}\right\\}}\right)$
$\displaystyle=\frac{y^{+}_{\max}}{\|x\|}+\frac{1}{\|x\|}\mathbb{E}_{x}\left(\left\|X\left(t_{0}\left/\|x\|^{\gamma}\right)\right.\right\|\right)\leq\frac{y^{+}_{\max}}{\|x\|}+(1{-}\varepsilon)$

since $(X(t))$ is $x$ on the time interval $(0,t_{1})$. We conclude with
Proposition 8. ∎

The scaling of the corollary gives only a simplified picture of a real CRN. In
general, the timescale used may depend on the initial state. Several examples
are investigated in Sections 6, 7 and 8.

There is however an interesting, global, scaling of this kind. Assuming that
$\|x\|$, the norm of the initial state is large, the rate of the fastest
reaction is at most of the order of $\|x\|^{y_{\max}^{-}}$, $y_{\max}^{-}$ is
the size of the largest complex initiating a reaction, see Definition (4). For
this reason, when looking at the scaled process, the timescale
$t{\mapsto}t/\|x\|^{y_{\max}^{-}{-}1}$ is quite natural.

######  Proposition 10.

If $(X^{N}(t))$ is the Markov process associated to a CRN network $({\cal
S},{\cal C},{\cal R})$ with parameters $(\kappa_{r})$ and whose initial state
is such that

$\lim_{N\to+\infty}\left(\frac{X_{i}^{N}(0)}{N}\right)=x_{0}{=}(x_{0,i}){\in}\mathbb{R}_{+}^{m},$

then there exists $t_{\infty}{>}0$, such that when $N$ goes to infinity the
family of processes

$\left(\overline{X}^{N}(t)\right)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\frac{1}{N}X\left({t}\left/{N^{y_{\max}^{-}{-}1}}\right),t{<}t_{\infty}\right)\right.,$

is converging in distribution to $(\ell(t),t{<}t_{\infty})$, the solution of
the ODE

(21) $\mathop{}\mathopen{}\mathrm{d}\ell(t)=\sum_{r\in{\cal
R},\,\|y^{-}_{r}\|=y^{-}_{\max}}\kappa_{r}\ell(t)^{y^{-}_{r}}(y^{+}_{r}{-}y^{-}_{r})$

on $(0,t_{\infty})$, with $y_{\max}^{-}$ defined by Relation (4).

This scaling has the effect of considering the same CRN but with a set of
reactions reduced to $\\{r{=}(y_{r}^{-},y_{r}^{+}){\in}{\cal
R}{:}\|y^{-}_{r}\|{=}y_{\max}^{-}\\}$.

###### Proof.

We proceed in an analogous way as in the proof of Proposition 2. The SDE
formulation of Section B.1 of the appendix to represent the Markov process
$(X^{N}(t))$ is used. Relation (61) gives for $t{\geq}0$,

(22) $\overline{X}^{N}(t)=\overline{X}^{N}(0)+\sum_{r{\in}{\cal
R}}\overline{M}^{N}_{r}(t)\\\ {+}\sum_{r{\in}{\cal
R}}\kappa_{r}\left(y_{r}^{+}{-}y_{r}^{-}\right)\int_{0}^{t}\frac{X^{N}({s}/{N^{y_{\max}^{-}{-}1}})!}{N^{y_{\max}^{-}}(X^{N}({s}/{N^{y_{\max}^{-}{-}1}}){-}y_{r}^{-})!}\mathop{}\mathopen{}\mathrm{d}s,$

where, for $r{\in}{\cal R}$, $(\overline{M}^{N}_{r}(t))$ is a martingale.
Stopping the process at $H^{N}_{K}$ as in the proof of Proposition 2, we
remark that for $r{=}(y_{r}^{-},y_{r}^{+}){\in}{\cal R}$ such that
$\|y^{-}_{r}\|<y_{\max}^{-}$,

$\frac{X^{N}(t)!}{N^{y_{\max}^{-}}(X^{N}(t){-}y_{r}^{-})!}\leq\frac{K^{y_{\max}^{-}-1}}{N},$

so that the corresponding terms in Relation (22) vanish, as processes, when
$N$ gets large. We can show easily that taking $t_{\infty}$ small enough,
$\lim_{N\to+\infty}\mathbb{P}\left(H^{N}_{K}{\leq}t\right)=0$, which allows us
to conclude. ∎

##  6\. Binary CRN Networks

In this section, we investigate simple examples of CRNs with complexes whose
size is at most $2$.

######  Definition 11.

A CRN network $({\cal S},{\cal C},{\cal R})$ is binary if any complex
$y{\in}{\cal C}$ is composed of at most two molecules, i.e. $\|y\|{\leq}2$.

The set of complexes can be represented as ${\cal
C}{=}J_{0}{\cup}J_{1}{\cup}J_{2}$, where, for $i{\in}\\{0,1,2\\}$, the subset
$J_{i}$ is the set of complexes of size $i$, note that $J_{i}$ can be empty.
If $y{\in}J_{1}$, with a slight abuse of notation, it will be represented as
$y{=}S_{y}$ for some $S_{y}{\in}{\cal S}$. Similarly,
$y{=}S_{y}^{1}+S_{y}^{2}$, for $S_{y}^{1},S_{y}^{2}{\in}{\cal S}$ when
$y{\in}J_{2}$.

######  Proposition 12.

If $(X(t))$ is the Markov process associated to a binary CRN network then,
almost surely, it is not explosive. Furthermore, for $T{>}0$, the family of
random variables in $\mathbb{R}_{+}^{2}$

$\left(X_{x}^{*}(T),x{\in}{\cal S}\right)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\frac{1}{\|x\|}\sup_{s\leq T}X(s),X(0){=}x{\in}{\cal
S}\right)$

is uniformly integrable.

###### Proof.

The SDEs associated to the CRN are given by

$\mathop{}\mathopen{}\mathrm{d}X(t)=\sum_{r\in{\cal
R}}(y_{r}^{+}{-}y_{r}^{-}){\cal
P}_{r}\left(\left(0,\kappa_{r}\frac{X_{r}!}{(X_{r}{-}y_{r}^{-})!}\right),\mathop{}\mathopen{}\mathrm{d}t\right),$

with our convention that ${\cal P}_{r}$, $r{\in}{\cal R}$, are independent
Poisson processes on $\mathbb{R}_{+}^{2}$ with intensity measure
$\mathop{}\mathopen{}\mathrm{d}s{\otimes}\mathop{}\mathopen{}\mathrm{d}t$ on
$\mathbb{R}_{+}^{2}$. See Section B.1 of the appendix. The binary condition
implies that for $r\in{\cal R}$ with $y_{r}^{-}{\in}J_{2}$, then

$\|y_{r}^{+}\|{-}\|y_{r}^{-}\|\leq 0.$

If $|{\cal R}|$ denotes the cardinality of ${\cal R}$, it is not difficult to
construct a coupling so that $\|X(t)\|{\leq}Z(t)$ holds, where $(Z(t))$ is the
solution of the SDE

$\mathop{}\mathopen{}\mathrm{d}Z(t)=2{\cal P}\left(\left(0,\kappa_{\max}|{\cal
R}|Z(t{-})\right),\mathop{}\mathopen{}\mathrm{d}t\right),$

with $Z(0){=}\|X(0)\|{=}\|x\|$. The process $(Z(t))$ is simply a pure birth
branching process with rate $\kappa_{\max}|{\cal R}|$. It is almost surely
non-explosive. For $0{\leq}t{\leq}T$, we have

$Z(t)=\|x\|+M_{Z}(t)+2\kappa_{\max}|{\cal
R}|\int_{0}^{t}Z(s)\mathop{}\mathopen{}\mathrm{d}s,$

where $(M_{Z}(t))$ is a martingale, with

$\left\langle M_{Z}\right\rangle(t)=4\kappa_{\max}|{\cal
R}|\int_{0}^{t}Z(s)\mathop{}\mathopen{}\mathrm{d}s.$

It is easily seen that $\mathbb{E}(Z(t)){=}\|x\|\exp(2\kappa_{\max}|{\cal
R}|t)$. If

$Z_{x}^{*}(t)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\sup_{s{\leq}t}\left(\frac{Z(s)}{\|x\|}\right),$

then,

$\frac{1}{9}Z_{x}^{*}(t)^{2}\leq 1{+}\frac{1}{\|x\|^{2}}\sup_{s\leq
T}|M_{Z}(s)|^{2}{+}4T(\kappa_{\max}|{\cal
R}|)^{2}\int_{0}^{t}Z_{x}^{*}(s)^{2}\mathop{}\mathopen{}\mathrm{d}s.$

Doob’s Inequality gives the inequality

$\mathbb{E}\left(\sup_{s\leq T}|M_{Z}(s)|^{2}\right)\leq
4\mathbb{E}(M_{Z}(T)^{2})\\\ =16\kappa_{\max}|{\cal
R}|\int_{0}^{T}\mathbb{E}(Z(s))\mathop{}\mathopen{}\mathrm{d}s\leq
32\kappa_{\max}|{\cal R}|T\|x\|e^{2\kappa_{\max}|{\cal R}|T}.$

and, with Gronwall’s Inequality, we obtain

$\sup_{x\in{\cal
S},x{\neq}0}\mathbb{E}\left(X_{x}^{*}(T)^{2}\right)\leq\sup_{x\in{\cal
S},x{\neq}0}\mathbb{E}\left(Z_{x}^{*}(T)^{2}\right)<{+}\infty.$

The family of random variable $(X_{x}^{*}(T))$ is uniformly integrable. The
proposition is proved. ∎

######  Proposition 13.

If $(X(t)$ is the Markov process associated to a binary CRN network then the
family of random variables

$T\left(\overline{X}(t)\right)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\frac{1}{\|x\|}X(t/\|x\|)\right),$

is tight when $x{=}X(0){\in}{\cal S}$ goes to infinity and any of its limiting
point $(\ell(t))$ is a continuous process satisfying the ODE

(23) $\dot{\ell}(t)=\sum_{\begin{subarray}{c}r\in{\cal R}\\\
y_{r}^{-}{\in}J_{2}\end{subarray}}\kappa_{r}\ell_{s_{y_{r}^{-},1}}(t)\ell_{s_{y_{r}^{-},2}}(t)(y_{r}^{+}{-}y_{r}^{-}).$

###### Proof.

This is a simple consequence of Proposition 10. ∎

The timescale $(t/\|x\|)$ and the space scale $1/\|x\|$ are valid for all
binary CRNs from the point of view of tightness properties. It does not mean
that they are the only ones, or the most meaningful. As it will be seen in
Section 6.1, depending on the type of initial state, it may happen that the
timescales $(t/\sqrt{\|x\|})$ or $(t)$ and the space scales $1/\sqrt{\|x\|}$
or $1$ are appropriate for the analysis of the asymptotic behavior of the time
evolution of the CRN. The timescale $(t/\|x\|)$ is well-suited when there are
complexes of size two and when the associated chemical species are all in
“large” number, of the order of $\|x\|$. Otherwise, it may be too slow to
change the state of the CRN, so that a faster timescale has to be used.

###  6.1. Triangular Binary Networks

We now consider a binary CRN with two chemical species, $m{=}2$, and three
distinct complexes $C_{i}$, $i{=}1$, $2$, $3$, of size ${\leq}2$ and the set
of routes is

${\cal R}{=}\\{(C_{1},C_{2}),(C_{2},C_{3}),(C_{3},C_{1})\\}.$

The purpose of this section is essentially pedagogical, to show, in a simple
setting, how the ideas of Sections 4 and 5 can be used in practice, on
stability and scaling properties.

As a side result, in Section C of the appendix, a proof of the positive
recurrence of a general class of triangle topologies is given: an arbitrary
set of chemical species, an arbitrary set of three complexes, i.e. whose sizes
are general, and a set of reactions ${\cal R}$ containing the three ones from
above. Proposition 35 of the appendix for triangle topologies can be seen as
the analogue of Theorem 7.6.1 of [18] for star networks. To the best of our
knowledge there are few such stability results with arbitrary complexes. We
have not been able to generalize this proof to a CRN with more than three
complexes.

$S_{2}$$S_{1}{+}S_{2}$$S_{1}$$\kappa_{2}$$\kappa_{12}$$\kappa_{1}$ (a)
Triangle T1

$S_{2}$$S_{1}{+}S_{2}$$S_{1}$$\emptyset$$\kappa_{2}$$\kappa_{02}$$\kappa_{01}$$\kappa_{12}$$\kappa_{1}$
(b) Triangle T2

Figure 2. Triangular CRNs

Note that the CRN T1 of Figure 2 is weakly reversible, with deficiency zero,
hence [6] shows that its associated Markov process is positive recurrent with
invariant distribution

(24) $\pi(x,y)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\frac{1}{Z}\frac{\rho_{1}^{x}}{x!}\frac{\rho_{2}^{y}}{y!},\quad(x,y){\in}{\cal
S},$

where ${\cal S}{\subset}\mathbb{N}^{2}$ is the state space, ${\cal
S}{=}\mathbb{N}^{2}{\setminus}\\{(0,0),(1,0)\\}$, $Z$ is the normalization
constant and $\rho_{1}{=}{\kappa_{2}}/{\kappa_{12}}$,
$\rho_{2}{=}{\kappa_{1}}/{\kappa_{12}}$.

The stability analysis of this CRN is revisited by using the criterion of
Theorem 6 and one of its consequences. This is, of course, not a new result
for T1, but this analysis can be done without extra-cost when each chemical
species has, in addition, an external source, as in the CRN T2 of Figure 2.
This CRN is not anymore weakly reversible, the results of [6] cannot be
applied.

The associated Markov process is denoted by
$(X^{N}(t)){=}(X^{N}_{1}(t),X^{N}_{2}(t))$. We denote by
$x^{N}{=}(x_{1}^{N},x_{2}^{N})$ the initial state, which is of norm $N$,
$N{=}x_{1}^{N}{+}x_{2}^{N}$, it is assumed that

(25)
$\lim_{N\to+\infty}\left(\frac{x_{1}^{N}}{N},\frac{x_{2}^{N}}{N}\right)=(\alpha_{1},1{-}\alpha_{1}),$

with $\alpha_{1}{\in}[0,1]$.

The scalings consider three types of regions of $\mathbb{N}^{2}$ for the
initial state: when the order of magnitude of the two coordinates are
respectively of the order of $(N,N)$, $(O(\sqrt{N}),N)$, or $(N,O(1))$. It is
shown that starting from a “large” state, three timescales play a role
depending on the asymptotic behavior of the initial state:

  1. (a)

$t\mapsto t/N$, when both components of the initial state are of the order of
$N$, i.e. when $0{<}\alpha_{1}{<}1$;

  2. (b)

$t\mapsto t/\sqrt{N}$, when $\alpha_{1}{=}0$ and $x^{N}_{1}$ is at most of the
order of $\sqrt{N}$.

  3. (c)

$t\mapsto t$, when $\alpha_{1}{=}1$ and $x_{2}^{N}$ is bounded by some
constant $K$.

The boundary effects mentioned in Section 2 play a role in case c), the second
coordinate remains in the neighborhood of the origin essentially.

For each of the three regimes, the scaled norm of the state is decreasing to
$0$, which is helpful to prove positive recurrence. The limit results show
additionally that the orders of magnitude in $N$ of both coordinates do not
change. In other words the space scale is natural and not the consequence of a
specific choice of the timescale. The following proposition gives a formal
statement of these assertions.

######  Proposition 14 (Scaling Analysis).

Under the assumptions (25) on the initial state of the CRN T1 of Figure 2,

  1. (a)

if $\alpha_{1}{>}0$, then for the convergence in distribution,

(26)
$\lim_{N\to+\infty}\left(\frac{X^{N}_{1}(t/N)}{N},\frac{X^{N}_{2}(t/N)}{N}\right)=(x_{a,1}(t),x_{a,2}(t)),$

where
$(x_{a,1}(t),x_{a,2}(t)){=}(\alpha_{1},(1-\alpha_{1})\exp(-\kappa_{12}\alpha_{1}t))$.

  2. (b)

If $\alpha_{1}{=}0$ and

$\lim_{N\to+\infty}\frac{x^{N}_{1}}{\sqrt{N}}=\beta{\in}\mathbb{R}_{+},$

then, for the convergence in distribution

(27)
$\lim_{N\to+\infty}\left(\frac{X^{N}_{1}(t/\sqrt{N})}{\sqrt{N}},\frac{X^{N}_{2}(t/\sqrt{N})}{N}\right)=(x_{b,1}(t),x_{b,2}(t)),$

where $(x_{b,1}(t),x_{b,2}(t))$ is the solution of the ODE

(28)
$\dot{x}_{b,1}(t)=\kappa_{2}x_{b,2}(t),\quad\dot{x}_{b,2}(t)={-}\kappa_{12}x_{b,1}(t)x_{b,2}(t),$

with $(x_{b,1}(0),x_{b,2}(0)){=}(\beta,1)$.

  3. (c)

If the initial state is $x^{N}{=}(N,k)$, for $k{\in}\mathbb{N}$, then, for the
convergence in distribution,

(29)
$\lim_{N\to+\infty}\left(\frac{X^{N}_{1}(t)}{N}\right)=(x_{c,1}(t))\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(e^{-\kappa_{1}t}\right).$

###### Proof.

See Section B.2 of the appendix. ∎

Stability Properties. We consider the CRN T2 of Figure 2, with external inputs
for both species. Here theorem 4 doesn’t apply, since the CRN is not weakly
reversible. We want to show a stability result on this CRN. It is not
difficult to see that the results of Proposition 14 holds in this case too
since on a finite time interval the number of external arrivals is almost
surely finite and independent of $N$.

######  Proposition 15.

The Markov process associated to the CRN T2 of Figure 2 is positive recurrent.

###### Proof.

Theorem 6 is used. We have to define the stopping time depending on the
initial state. Proposition 14 does not give a full partitioning of the
possible “large” initial states, some additional work has to be done. We
ignore the external arrivals, it is easily seen that similar arguments can be
used. With high probability, the stopping times chosen are smaller than the
instant of the first external arrival. There are three cases.

  1. (1)

The initial state is such that

$\lim_{N\rightarrow+\infty}\frac{x_{1}^{N}}{\sqrt{N}}\geq 1\quad\text{and
}\quad x_{2}^{N}\geq 1.$

We take $t^{N}_{1}$, the instant of the first jump of $(X_{N}(t))$. Elementary
calculations give the relation

$\limsup_{N\to+\infty}\mathbb{E}_{x^{N}}(\|X^{N}(t^{N}_{1})\|)-\|x^{N}\|\leq{-}\frac{\kappa_{12}}{\kappa_{12}+\kappa_{1}},\text{
and }\lim_{N\to+\infty}\sqrt{N}\mathbb{E}\left(t^{N}_{1}\right)=0.$

  2. (2)

The initial state is $x^{N}{=}(N,0)$.
Let, for $k_{0}{\in}\mathbb{N}$, $\tau_{k_{0}}$ be the instant when the
$(k_{0}{+}1)$th element $S_{1}$ is transformed into an $S_{2}$. The norm of
the process decreases when some of these molecules of $S_{2}$ disappear with
reaction $S_{1}{+}S_{2}{\rightharpoonup}S_{1}$. The probability for a molecule
of $S_{2}$ to be destroyed before a new transformation of $S_{1}$ into $S_{2}$
is lower bounded by $p{=}\kappa_{12}/(\kappa_{1}{+}\kappa_{12})$, therefore in
average there are more than $pk_{0}$ molecules of $S_{2}$ killed before
$\tau_{k_{0}}$. The reaction $S_{2}{\rightharpoonup}S_{1}{+}S_{2}$ could also
create some molecules during the time interval $[0,\tau_{k_{0}}]$. Its rate
being bounded by $\kappa_{2}k_{0}$, and noting that
$\mathbb{E}(\tau_{k_{0}}){\leq}k_{0}/(\kappa_{1}(N{-}k_{0}))$, it is not
difficult to show that there exists a constant $C_{0}$ such that the relation

$\mathbb{E}\left(\|X^{N}(\tau_{k_{0}})\|\right)\leq
N{-}k_{0}p{+}k_{0}\frac{C_{0}}{N{-}k_{0}},$

holds.

  3. (3)

The initial state is $x^{N}{=}(x^{N}_{1},N)$ with,

$\lim_{N\to+\infty}\frac{x^{N}_{1}}{\sqrt{N}}\leq 1.$

We will use the result (b) of Proposition 14 and its proof in the appendix
(for the convergence of the averages), to obtain the convergence

$\lim_{N\to+\infty}\frac{1}{N}\mathbb{E}_{x^{N}}(\|X(1/\sqrt{N})\|)=x_{b,2}(1){<}1.$

We only have to use Theorem 6 to conclude the proof of the proposition. ∎

##  7\. Agazzi and Mattingly’s CRN

In this section, we study the chemical reaction network introduced by [2],

(30)
$\emptyset\mathrel{\mathop{\xrightharpoonup{\kappa_{1}}}}S_{1}{+}S_{2},\quad
S_{2}\mathrel{\mathop{\xrightharpoonup{\kappa_{2}}}}\emptyset,\quad
pS_{1}{+}qS_{2}\mathrel{\mathop{\xrightharpoonup{\kappa_{3}}}}(q{+}1)S_{2}\mathrel{\mathop{\xrightharpoonup{\kappa_{4}}}}qS_{2},$

for $p$, $q{\in}\mathbb{N}$, $p>2$ and $q\geq 2$.

The continuous time Markov jump process $(X(t))=(X_{1}(t),X_{2}(t))$
associated to CRN (30) has a $Q$-matrix given by, for $x{\in}\mathbb{N}^{2}$,

$x\longrightarrow x{+}\begin{cases}e_{1}{+}e_{2}&\kappa_{1},\\\
{-}pe_{1}{+}e_{2}&\kappa_{3}x_{1}^{(p)}x^{(q)}_{2},\\\
{-}e_{2}&\kappa_{2}x_{2}{+}\kappa_{4}x^{(q+1)}_{2},\end{cases}$

where $e_{1}$, $e_{2}$ are the unit vectors of $\mathbb{N}^{2}$. This process
is clearly irreducible on $\mathbb{N}^{2}$, and non explosive since
$\emptyset{\rightharpoonup}A{+}B$ is the only reaction increasing the number
of molecules. See the proof of Proposition 12 for example. The fact that only
this reaction increases the norm of the state suggests that the proof of the
positive recurrence should not be an issue.

To prove this positive recurrence, see the proposition below, [2] use several
energy functions on $\mathbb{N}^{2}$, which are polynomial functions in
$x_{1}$ and $x_{2}$. The main technical difficulty is of gluing these
functions in order to have a global Lyapunov function for which the classical
Forster-Lyapunov theorem can be used. Note that there are also interesting
null recurrence and transience properties in this reference.

######  Proposition 16.

If $p{>}2$ and $q{\geq}2$, then the Markov process associated to the CRN (30)
is positive recurrent.

###### Proof.

Theorem 6 is used with a simple energy function, the norm
$\|x\|{=}x_{1}{+}x_{2}$ of the state $x{=}(x_{1},x_{2}){\in}\mathbb{N}^{2}$.
If the norm of the initial state is large enough, then the expected value of
the norm of the process taken at a convenient stopping time will be smaller,
so that Condition (15) of Theorem 6 holds.

### Step 1

As before, for $n{\geq}1$, $t_{n}$ denotes the instant of the $n$th jump.

$\mathbb{E}_{x}\left(\|X(t_{1})\|{-}\|x\|\right)=\left(2\kappa_{1}{-}\kappa_{2}x_{2}{-}(p{-}1)\kappa_{3}x_{1}^{(p)}x_{2}^{(q)}{-}\kappa_{4}x_{2}^{(q+1)}\right)\mathbb{E}_{x}\left[t_{1}\right],$

and, clearly, $\mathbb{E}_{x}(t_{1}){\leq}{1}/{\kappa_{1}}$.

If either $x_{2}{\geq}K_{1}{=}1{+}2\kappa_{1}/\kappa_{2}$ or
$q{\leq}x_{2}{<}K_{1}$ and
$x_{1}{\geq}K_{2}{=}1{+}2\kappa_{1}/((p{-}1)\kappa_{3}q!)$, then

$\mathbb{E}_{x}\left(\|X(t_{1})\|{-}\|x\|\right)\leq{-}\gamma\mathbb{E}_{x}(t_{1}),$

for some $\gamma{>}0$. Condition (15) holds for this set of initial states.

### Step 2

Now we consider initial states of the form $x_{0}^{N}{=}(N,b)$ with $b{<}q$
and $N$ large. The third and fourth reactions cannot occur until the instant

$\tau_{1}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\inf\\{t{>}0:X_{2}(t){\geq}q\\}.$

Until time $\tau_{1}$, the process $(X_{2}(t))$ has the sample path $(L(t))$
of an $M/M/\infty$ queue, see Section 2.2, with arrival rate $\kappa_{1}$ and
service rate $\kappa_{2}$. At time $\tau_{1}$ the state of the process has the
same distribution as the random variable

$(N{+}{\cal N}_{\kappa_{0}}(0,\tau_{1}),q),$

where ${\cal N}_{\kappa_{0}}$ is a Poisson process with rate $\kappa_{0}$.
Clearly $\tau_{1}$ is integrable as well as the random variable ${\cal
N}_{\kappa_{0}}(0,\tau_{1})$. We have also
$X_{1}(\tau_{1}{\wedge}t){=}N{+}\cal{N}_{\kappa_{1}}^{1}([0,t{\wedge}\tau_{1}])$,
so $\mathbb{E}_{(N,b)}[X_{1}(\tau_{1})]\leq N{+}\kappa_{1}C_{1}$, for some
constant $C_{1}$.

To summarize, starting from the initial state $x_{0}^{N}{=}(N,b)$ with
$b{<}q$, the quantities $\mathbb{E}_{x_{0}^{N}}(\tau_{1})$ and
$\mathbb{E}_{x_{0}^{N}}(X_{1}(\tau_{1})){-}N$ are bounded by a constant. We
are thus left to study the following case.

### Step 3

The initial state is $x_{0}^{N}{=}(N,q)$ with $N$ large.

As long as $X_{2}(t){\geq}q$, the third reaction is active, $p$ copies of
$S_{1}$ are removed and a copy of $S_{2}$ is created. Initially its rate is of
the order of $N^{p}$, the fastest reaction rate by far. We define $\nu$ as the
number of jumps before another reaction takes place.

$\nu\stackrel{{\scriptstyle\text{\rm def.}}}{{=}}\inf\\{n\geq
1:X(t_{n}){-}X(t_{n-1}){\neq}(-p,1)\\},$

$\mathbb{P}(\nu{>}k)=\prod_{i=0}^{k-1}\left(1{-}\frac{\kappa_{1}{+}\kappa_{2}(q{+}i){+}\kappa_{4}(q{+}i)^{(q+1)}}{\kappa_{3}(N{-}pi)^{(p)}(q{+}i)^{(q)}{+}\kappa_{1}{+}\kappa_{2}(q{+}i){+}\kappa_{4}(q{+}i)^{(q+1)}}\right),$

with the convention that $q^{(q+1)}{=}0$. For $i{\geq}1$,

$\frac{\kappa_{1}{+}\kappa_{2}(q{+}i){+}\kappa_{4}(q{+}i)^{(q+1)}}{\kappa_{3}(N{-}pi)^{(p)}(q{+}i)^{(q)}{+}\kappa_{1}{+}\kappa_{2}(q{+}i){+}\kappa_{4}(q{+}i)^{(q+1)}}\\\
\leq\frac{(\kappa_{1}{+}\kappa_{2}(q{+}i))(q{+}i)^{-(q+1)}{+}\kappa_{4}}{\kappa_{3}(N{-}pi)^{(p)}/i{+}(\kappa_{1}{+}\kappa_{2}(q{+}i))(q{+}i)^{-(q+1)}{+}\kappa_{4}}\leq\frac{iC_{0}}{(N{-}pi)^{(p)}{+}iC_{0}},$

for some apropriate constant $C_{0}{>}0$. Hence, if we fix
$0{<}\delta{<}1/2p$,

$\mathbb{E}_{x_{0}^{N}}(\nu)\geq\delta N\mathbb{P}(\nu{>}\delta N)\geq\delta
N\left(1{-}\frac{\delta NC_{0}}{(N{-}p\lfloor\delta N\rfloor)^{(p)}{+}\delta
NC_{0}}\right)^{\lfloor\delta N\rfloor},$

so that, since $p{>}2$,

(31) $\liminf_{N\to+\infty}\frac{1}{N}\mathbb{E}_{x_{0}^{N}}(\nu)\geq\delta.$

We define $\tau_{2}{=}t_{\nu}$, obviously

$\mathbb{E}_{x_{0}^{N}}(\tau_{2})\leq\frac{1}{\kappa_{1}},$

and we have

$\mathbb{E}_{x_{0}^{N}}(\|X(\tau_{2})\|{-}\|x_{0}^{N}\|)\leq(1{-}p)\mathbb{E}_{x_{0}^{N}}(\nu){+}2\leq-\gamma
N,$

for some $\gamma{>}0$ if $N$ is sufficiently large, using Relation (31).
Consequently is easy to see that there is a convenient constant $K$ such that
Condition (15) holds for this set of initial states and the stopping time
$\tau_{2}$, and also for the initial states of Step 2 and the stopping time
$\tau_{1}{+}\tau_{2}{\circ}\theta_{\tau_{1}}$. The proposition is proved. ∎

### A Scaling Picture

The key argument of the proof of the positive recurrence is somewhat hidden
behind an estimate of the expected value of the hitting time $\nu$ in Step 3.
It is not difficult to figure out that, starting from the state $(N,q)$, the
“right” timescale is $t{\mapsto}t/N^{p+q-1}$. In this section we sketch a
scaling argument to describe in more detail how the norm of the state goes to
$0$. It could also give an alternative way to handle Step 3.

Define the Markov jump process $(Z_{N}(t))=(Z^{N}_{1}(t),Z^{N}_{2}(t))$
corresponding to the last two reactions of the CRN network (30). Its
$Q$-matrix is given by, for $z{\in}\mathbb{N}^{2}$,

(32) $z\longrightarrow
z{+}\begin{cases}{-}pe_{1}{+}e_{2}&\kappa_{3}z_{1}^{(p)}z_{2}^{(q)},\\\
{-}e_{2}&\kappa_{4}z_{2}^{(q+1)},\end{cases}$

with initial state $(N,q)$. The scaling results of this section are obtained
for this process. It is not difficult to show that they also hold for the CRN
network (30) since the discarded reactions are on a much slower timescale.

Define the Markov jump process $(Y_{N}(t))=(Y^{N}_{1}(t),Y^{N}_{2}(t))$ whose
$Q$-matrix is given by, for $y{\in}\mathbb{N}^{2}$,

$y\longrightarrow y{+}\begin{cases}{-}pe_{1}{+}e_{2}&\kappa_{3}y_{1}^{(p)},\\\
{-}e_{2}&\kappa_{4}(y_{2}{-}q),\end{cases}$

with the same initial state. If $p{\geq}2$, with the same arguments as in the
proof of Proposition 14 (see Section B.2 of the appendix), it is not difficult
to show the convergence in distribution

(33)
$\lim_{N\to+\infty}\left(\frac{1}{N}\left(Y^{N}_{1},Y^{N}_{2}\right)\left(\frac{t}{N^{p-1}}\right)\right)=\left(y_{1}(t),y_{2}(t)\right)\\\
\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\frac{1}{\sqrt[p{-}1]{p(p{-}1)\kappa_{3}t{+}1}},\frac{1{-}y_{1}(t)}{p}\right)$

From this convergence we obtain that for any $\eta{\in}(0,1/p)$, the hitting
time $H^{N}_{Y}(\eta)$ of $\lfloor\eta N\rfloor$ by $(Y^{N}_{2}(t))$ is such
that $(N^{p-1}H^{N}_{Y}(\eta))$ converges in distribution to some constant.

For $t{\geq}0$, define the stopping time

$\tau^{N}_{t}=\inf\left\\{s{>}0:\int_{0}^{s}\frac{1}{Y^{N}_{2}(u)^{(q)}}\mathop{}\mathopen{}\mathrm{d}u\geq
t\right\\},$

and $(\widetilde{Z}^{N}(t)){=}(Y^{N}(\tau^{N}_{t}))$, then it is easy to check
that $(\widetilde{Z}^{N}(t))$ is a Markov process whose $Q$-matrix is given by
Relation (32). See Section III.21 of [41] for example. Consequently,
$(\widetilde{Z}^{N}(t))$ has the same distribution as $(Z^{N}(t))$.

######  Proposition 17.

If $p$, $q{\geq}2$, $(X^{N}(0)){=}(\lfloor\delta
N\rfloor,\lfloor(1{-}\delta)N/p\rfloor)$, for some $\delta{\in}(0,1)$, then
for the convergence in distribution

$\lim_{N\to+\infty}\left(\frac{1}{N}X^{N}\left(\frac{t}{N^{p+q-1}}\right)\right)=(x_{1}(t),x_{2}(t)){=}\left(\left(y_{1},\frac{1{-}y_{1}}{p}\right)\left(\phi^{-1}(t)\right)\right),$

with

$(y_{1}(t))=\left(\frac{\delta}{\sqrt[p{-}1]{p(p{-}1)\delta^{p-1}\kappa_{3}t{+}1}}\right)\text{
and }\phi(t)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\int_{0}^{t}\frac{p^{q}}{(1{-}y_{1}(s))^{q}}\mathop{}\mathopen{}\mathrm{d}s.$

###### Proof.

As mentioned above, from this initial state and this timescale, the processes
$(Z^{N}(t))$ and $(X^{N}(t))$ have the same asymptotic behavior for values of
the order of $N$. The proof uses the convergence (33) and the time-change
argument described above. ∎

The above proposition shows that on a convenient timescale, both coordinates
of $(X^{N}(t))$ are of the order of $N$. The scaled version of the first one
is converging to $0$, while the second component is increasing.

If $Y^{N}(0){=}(\lfloor\delta N\rfloor,\lfloor(1{-}\delta)N/p\rfloor)$, for
some $\delta{>}0$, let

$H^{N}=\inf\\{t{>}0:Y^{N}_{1}(t)\leq\sqrt[p]{N}\\}.$

By writing the evolution of $(Y^{N}(t))$ in terms of an SDE like Relation
(B.1), one easily obtains,

$\mathbb{E}(Y^{N}_{1}(H^{N}{\wedge}t))=\lfloor\delta N\rfloor{-}p\kappa_{3}\
E\left(\int_{0}^{H^{N}{\wedge}t}Y^{N}_{1}(s)^{(p)}\mathop{}\mathopen{}\mathrm{d}s\right)\\\
\leq\lfloor\delta
N\rfloor{-}p\kappa_{3}\left(\sqrt[p]{N}\right)^{(p)}E\left(H^{N}{\wedge}t\right),$

hence, by using the monotone convergence theorem, we obtain that

$E\left(H^{N}\right)\leq\frac{\lfloor\delta
N\rfloor}{p\kappa_{3}(\sqrt[p]{N})^{(p)}},\text{ hence,
}\sup_{N}E\left(H^{N}_{K}\right)<{+}\infty,$

since $p{\geq}2$. It is easily seen that the same property holds for
$(X^{N}_{1}(t))$.

To finish the description of the return path to $(0,0)$, we can assume
therefore that $X^{N}(0){=}(\lfloor\sqrt[p]{N}\rfloor,N)$. It is not difficult
to see that the reaction
$(q{+}1)S_{2}\mathrel{\mathop{\xrightharpoonup{\kappa_{4}}}}qS_{2}$ is driving
the evolution as long as $(X^{N}_{2}(t))$ is “large” since $(X^{N}_{1}(t))$
cannot grow significantly on the corresponding timescale. More formally, also
with the same arguments as in Section 14, the convergence in distribution

$\lim_{N\to+\infty}\left(\frac{1}{N}\left(X^{N}_{1},X^{N}_{2}\right)\left(\frac{t}{N^{q}}\right)\right)=\left(0,\frac{1}{\sqrt[q]{1{+}\kappa_{4}qt}}\right)$

holds.

##  8\. A CRN with Slow and Fast Timescales

In this section, the positive recurrence and scaling properties of the
following interesting CRN are investigated

(34)
$\emptyset\mathrel{\mathop{\xrightleftharpoons[\kappa_{1}]{\kappa_{0}}}}S_{1}{+}S_{2},\hskip
28.45274ptpS_{1}{+}S_{2}\mathrel{\mathop{\xrightleftharpoons[\kappa_{3}]{\kappa_{2}}}}pS_{1}{+}2S_{2},$

with $p{\geq}2$. It has been discussed in [1] for $p{=}2$.

This CRN exhibits several distinctive features of chemical reaction networks.
It provides an important example of a non-trivial CRN for which the results
and ideas of Section 5 and Section 4 can be used together with time change
arguments, and all of that within a (quite) limited technical framework. For
this reason, we have chosen to develop completely the technical arguments
used. The results obtained are interesting in their own right in fact.

Section 8.1 investigates the positive recurrence properties. It is an occasion
to have an other look at the choice of a Lyapunov function in view of
Condition 15 of Theorem 6. Section 8.2 considers the limiting behavior of the
sample paths of the CRN with a large initial state close to one of the axes.
As it can be expected, in both sections boundary effects play a very important
role: the second reaction cannot occur if there are less than $p$ copies of
$S_{1}$, and if the number of copies of $S_{2}$ is zero, only external
arrivals change the state of the CRN.

The Markov process $(X(t))=(X_{1}(t),X_{2}(t))$ associated to this CRN has a
$Q$-matrix $Q$ given by, for $x{\in}\mathbb{N}^{2}$,

$x\longrightarrow x{+}\begin{cases}e_{1}{+}e_{2}&\kappa_{0},\\\
-e_{1}{-}e_{2}&\kappa_{1}x_{1}x_{2},\end{cases}x\longrightarrow
x{+}\begin{cases}e_{2}&\kappa_{2}x_{1}^{(p)}x_{2},\\\
{-}e_{2}&\kappa_{3}x_{1}^{(p)}x_{2}^{(2)},\end{cases}$

where $e_{1}$, $e_{2}$ are the unit vectors of $\mathbb{N}^{2}$.

By using the SDE formulation of Section B.1 of the appendix, the associated
Markov process can be represented by the solution
$(X(t)){=}(X_{1}(t),X_{2}(t))$ of the SDE

(35) $\begin{cases}\mathop{}\mathopen{}\mathrm{d}X_{1}(t)&={\cal
P}_{X,0}((0,\kappa_{0}),\mathop{}\mathopen{}\mathrm{d}t){-}{\cal
P}_{X,1}((0,\kappa_{1}X_{1}X_{2}(t{-})),\mathop{}\mathopen{}\mathrm{d}t),\\\
\mathop{}\mathopen{}\mathrm{d}X_{2}(t)&={\cal
P}_{X,0}((0,\kappa_{0}),\mathop{}\mathopen{}\mathrm{d}t){-}{\cal
P}_{X,1}((0,\kappa_{1}X_{1}X_{2}(t{-})),\mathop{}\mathopen{}\mathrm{d}t)\\\
&\hskip 42.67912pt{+}{\cal
P}_{X,2}\left(\left(0,\kappa_{2}X_{1}^{(p)}X_{2}(t{-})\right),\mathop{}\mathopen{}\mathrm{d}t\right)\\\
&\hskip 42.67912pt{-}{\cal
P}_{X,3}\left(\left(0,\kappa_{3}X_{1}^{(p)}X_{2}^{(2)}(t{-})\right),\mathop{}\mathopen{}\mathrm{d}t\right),\end{cases}$

where ${\cal P}_{X,i}$, $i{\in}\\{0,1,2,3\\}$, are fixed independent Poisson
processes on $\mathbb{R}_{+}^{2}$ with intensity measure
$\mathop{}\mathopen{}\mathrm{d}s{\otimes}\mathop{}\mathopen{}\mathrm{d}t$. A
notation of this kind ${\cal P}_{A}$ or ${\cal P}_{A,i}$ will be used for
several $A$ in the following, with the same assumptions on the distribution
and the independence properties.

A slow return to $0$. The reactions of the second linkage class of this CRN
need $p$ copies of $S_{1}$ to be active. If the initial state is $(0,N)$,
copies of $S_{1}$ are created at rate $\kappa_{0}$, but they are removed
quickly at a rate greater than $\kappa_{1}N$. The first instant when $p$
copies of $S_{1}$ are present has an average of the order of $N^{p-1}$. See
Lemma 19. At that instant, the number of $S_{2}$ species is $N{+}p$, and the
second coordinate can then decrease, quickly in fact. The network exhibits a
kind of bi-modal behavior due to this boundary condition.

Starting from the initial state $x{=}(0,N)$, the time to decrease $(X_{2}(t))$
by an amount of the order of $N$ has thus an average of the order of
$\|N\|^{p-1}$. When $p{>}2$ and if we take the usual norm $\|\cdot\|$ as a
Lyapunov function, this results is at odds with one of the conditions of
positive recurrence criterion of Proposition 8. This problem could in fact be
fixed at the cost of some annoying technicalities. Our approach will be of
taking another simple, and somewhat natural, Lyapunov function. See Section
8.1. An initial state of the form $(N,0)$ leads also to another interesting
boundary behavior.

###  8.1. Positive Recurrence

######  Proposition 18.

The Markov process $(X(t))$ is positive recurrent.

Theorem 6 is used to prove this property. The proof is not difficult but it
has to be handled with some care. We will introduce two auxiliary processes
with which the process $(X_{N}(t))$ can be decomposed. One describes the
process when the first coordinate is below $p$ and the other when the second
coordinate is greater or equal to $1$. This representation gives a more formal
description of the bi-modal behavior mentioned above. Additionally, it will
turn out to be helpful to establish the scaling properties of this CRN in
Section 8.2. For $x{=}(x_{1},x_{2}){\in}\mathbb{N}^{2}$, we introduce

(36) $f_{p}(x){=}x_{1}{+}x_{2}^{p},$

$f_{p}$ will be our Lyapunov function. The strategy is of analyzing separately
the two boundary behaviors. The first one is essentially associated with the
initial state $(0,N)$ which we have already seen. The other case is for an
initial state of the form $(N,0)$, the problem here is of having the second
coordinate positive sufficiently often so that the first reaction can decrease
significantly the first coordinate.

####  8.1.1. Large Initial State along the Horizontal Axis

In this section it is assumed that the initial state is
$x(0){=}(x^{0}_{1},b)$, where $b{\in}\mathbb{N}$ is fixed and $x^{0}_{1}$ is
“large”. Without loss of generality one can assume $b{>}0$, otherwise nothing
happens until an external arrival.

As long as the second coordinate of $(X(t))$ is non-null the transitions
associated to ${\cal P}_{X,i}$, $i{=}2$, $3$ occur at a fast rate. When
$(X_{2}(t))$ is $0$, only one chemical reaction may happen, external arrivals
and at a “slow” rate $\kappa_{0}$.

We define by induction the non-increasing sequence $(T_{k})$ as follows,
$T_{0}{=}0$, and

$T_{k+1}=\inf\\{t{>}T_{k}:X_{1}(t){-}X_{1}(t{-}){=}{-}1\\}.$

The variables $(T_{k})$ are stopping times for the underlying filtration
$({\cal F}_{t})$ defined as in the appendix, see Relation (62).

For $t{>}0$, by using the fact that the Poisson process ${\cal P}_{X,i}$,
$i{=}1$, $2$, $3$ are independent and $(X_{2}(t))$ is greater than $1$ until
$T_{1}$ at least, we have

$\mathbb{P}(T_{1}{\geq}t)\leq\mathbb{E}\left(\exp\left({-}\kappa_{1}x_{1}^{0}\int_{0}^{t}X_{2}(s)\mathop{}\mathopen{}\mathrm{d}s\right)\right)\leq\exp\left({-}\kappa_{1}x^{0}_{1}t\right),$

hence $\mathbb{E}(T_{1}){\leq}1/(\kappa_{1}x^{0}_{1})$. Similarly, with the
strong Markov property, for $1{\leq}k{<}x^{0}_{1}$,

$\mathbb{E}(T_{k+1}{-}T_{k}){\leq}\frac{1}{\kappa_{0}}{+}\frac{1}{\kappa_{1}(x^{0}_{1}{-}k)},$

the additional term $1/\kappa_{0}$ comes from the fact that $X_{2}(T_{k})$ can
be zero, so that one has to wait for an exponential time with parameter
$\kappa_{0}$ to restart the CRN.

For $n_{0}{\geq}1$, we have seen that the random variable $T_{n_{0}}$ is
stochastically bounded by the sum of $2n_{0}$ i.i.d. exponentially distributed
random variables with some positive rate, hence

$C_{0}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\sup_{x_{1}^{0}>n_{0}}\mathbb{E}_{x^{0}_{1}}(T_{n_{0}})<{+}\infty$

Let ${\cal E}_{1}$ be the event when ${\cal P}_{X,1}$ has a jump before ${\cal
P}_{X,0}$ in SDE (35), then

$\mathbb{P}({\cal
E}_{1}^{c}){\leq}\frac{\kappa_{0}}{\kappa_{1}x^{0}_{1}{+}\kappa_{0}}.$

Similarly, for $k{\geq}2$, ${\cal E}_{k}$ is a subset of the event ${\cal
E}_{k-1}$ for which ${\cal P}_{X,1}$ has a jump before ${\cal P}_{X,0}$ after
the first time after $T_{k}$ when $(X_{2}(t))$ is greater than $1$, then

(37) $\mathbb{P}_{x_{1}^{0}}({\cal
E}_{k}^{c}){\leq}\sum_{i=0}^{k-1}\frac{\kappa_{0}}{\kappa_{1}(x^{0}_{1}{-}i){+}\kappa_{0}}{\leq}\frac{\kappa_{0}k}{\kappa_{1}(x^{0}_{1}{-}k){+}\kappa_{0}}.$

Let $s_{1}$ be the first instant of jump of ${\cal
P}_{X,0}((0,\kappa_{0}){\times}(0,t])$. From $t{=}0$, as long as the point
process ${\cal P}_{X,0}$, does not jump in SDE (35), that is, on the time
interval $[0,s_{1}]$, up to a change of time scale $t{\mapsto}X_{1}X_{2}(t)$,
the process $(X_{1}(t),X_{2}(t))$ has the same sequence of visited states as
the solution $(Y(t))$ of the SDE

(38) $\begin{cases}\mathop{}\mathopen{}\mathrm{d}Y_{1}(t)&={-}{\cal
P}_{X,1}((0,\kappa_{1}),\mathop{}\mathopen{}\mathrm{d}t),\\\
\mathop{}\mathopen{}\mathrm{d}Y_{2}(t)&={-}{\cal
P}_{X,1}((0,\kappa_{1})),\mathop{}\mathopen{}\mathrm{d}t)\\\ &\hskip
42.67912pt{+}{\cal
P}_{Y,2}\left(\left(0,\kappa_{2}Y_{1}(t{-})^{(p)-1}\right),\mathop{}\mathopen{}\mathrm{d}t\right)\\\
&\hskip 42.67912pt{-}{\cal
P}_{Y,3}\left(\left(0,\kappa_{3}Y_{1}^{(p)-1}(Y_{2}(t{-}){-}1)^{+}\right),\mathop{}\mathopen{}\mathrm{d}t\right),\end{cases}$

with the same initial state and the slight abuse of notation
$y^{(p){-}1}{=}y^{(p)}/y$. In particular if $u_{1}$ is the first instant when
$(Y_{1}(t))$ has a downward jump, an independent exponential random variable
with parameter $\kappa_{1}$, then the relation $Y_{2}(u_{1}){=}X_{2}(T_{1})$
holds on the event $\\{T_{X,1}{\leq}s_{1}\\}$.

From $t{=}0$, as long as the first coordinate of $(Y_{1}(t))$ does not change,
the second component $(Y_{2}(t))$ has the same distribution as
$(L_{b}((x_{1}^{0})^{(p)-1}t))$, where $(L_{b}(t))$ is a birth and death
process with birth rate $\kappa_{2}$ and death rate $\kappa_{3}(x{-}1)$ in
state $x{\geq}1$ and initial state $b$. It is easily seen that it is a
positive recurrent Markov process and that $(\mathbb{E}(L_{b}(t)^{p}))$ is a
bounded function. Consequently,

(39) $\sup_{x(0)}E\left(X_{2}(T_{1})^{p}\right)\leq C_{1}{<}{+}\infty,$

by induction, the same result holds for $T_{n_{0}}$ for a convenient constant
$C_{1}$.

Note that if $X_{2}(T_{1}{-}){=}1$, the next reaction happening after $T_{1}$
will be $\emptyset\rightharpoonup S_{1}+S_{2}$, and therefore the jump
downward of $X_{1}$ will be canceled. A decrease at time $T_{1}$ of
$(X_{1}^{N}(t))$ is sustainable if it happens when $X_{2}(T_{1}{-})\geq 2$
i.e. if $L_{b}((x_{1})^{(p)-1}u_{1}){\neq}0$. It is not difficult to construct
a coupling with $(L_{0}(t))$, a birth and death process starting at $0$, such
that $L_{b}(t){\geq}L_{0}(t)$ holds for all $t{\geq}0$. The convergence of
$(L_{0}(t))$ to equilibrium gives the existence of $K_{0}{\geq}0$ and
$\eta_{0}{>}0$ such that if $x^{0}_{1}{\geq}K_{0}$ then
$\mathbb{P}(L_{0}((x_{1})^{(p)-1}u_{1}){>}0){\geq}\eta_{0}$.

We can gather these results, and the stochastic bound on $T_{n_{0}}$, to get
the relations

$\displaystyle E_{x(0)}$
$\displaystyle(f_{p}(X(T_{n_{0}}))){-}f_{p}(x(0))\leq{-}{n_{0}}\eta_{0}\mathbb{P}_{x(0)}({\cal
E}_{n_{0}})$ $\displaystyle+\mathbb{E}_{x(0)}\left({\cal
P}_{X,0}\left((0,\kappa_{0}){\times}(0,T_{n_{0}}]\mathbbm{1}_{\left\\{{\cal
E}_{n_{0}}^{c}\right\\}}\right)\right)+E\left(X_{2}(T_{n_{0}})^{p}\right){-}b^{p}$
$\displaystyle\leq{-}\eta_{0}n_{0}+n_{0}\eta_{0}\mathbb{P}_{x(0)}({\cal
E}_{n_{0}}^{c}){+}\kappa_{0}C_{0}{+}C_{1}.$

One first choose $n_{0}$ so that $n_{0}{>}3(\kappa_{0}C_{0}{+}C_{1})/\eta_{0}$
and then with Relation (37), $K_{1}{\geq}K_{0}$ such that
$n_{0}\eta_{0}\mathbb{P}_{K_{1}}({\cal
E}_{k}^{c}){<}(\kappa_{0}C_{0}{+}C_{1})$. We obtain therefore that if
$x^{0}_{1}{>}K_{1}$, then

(40)
$\mathbb{E}_{x(0)}\left(f_{p}\left(X(T_{n_{0}})\right){-}f_{p}(x(0))\right)\leq-\delta,$

for some $\delta{>}0$ and
$\sup(\mathbb{E}_{x(0)}(T_{n_{0}}):x_{1}{\geq}K){<}{+}\infty$. Relation (40)
shows that Condition (15) of Theorem (6) is satisfied for our Lyapunov
function $f_{p}$ and stopping time $T_{n_{0}}$ for the initial state of the
form $(x^{0}_{1},b)$.

####  8.1.2. Initial State with a Large Second Component

In this section it is assumed that the initial state is $x(0){=}(a,x^{0}_{2})$
with $a{<}p$ and $x^{0}_{2}$ is large. We note that, as long as $(X_{1}(t))$
is strictly below $p$, the two coordinates experience the same jumps, the
quantity $(X_{2}(t){-}X_{1}(t))$ does not change. For this reason, for
$x{\geq}0$ and $k{\leq}p-1$, we introduce a process $(Z(k,x,t))$ which will be
used to express $(X(t))$ when its first coordinate is less than $p{-}1$. It is
the solution of the SDE

(41) $\mathop{}\mathopen{}\mathrm{d}Z(k,x^{0}_{2},t)={\cal
N}_{\kappa_{0}}(\mathop{}\mathopen{}\mathrm{d}t){-}{\cal
P}_{Z}((0,\kappa_{1}Z(k,x^{0}_{2},t{-})(x^{0}_{2}{-}k{+}Z(k,x^{0}_{2},t{-})))),\mathop{}\mathopen{}\mathrm{d}t),\\\
$

with $Z(k,x^{0}_{2},0){=}k$ and ${\cal P}_{Z}$ is a Poisson process on
$\mathbb{R}_{+}^{2}$.

Setting for $z<p$

$T_{Z}(z,x^{0}_{2})\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\inf\\{t{>}0:Z(z,x^{0}_{2},t){=}p\\},$

if $X(0){=}(0,x^{0}_{2})$, then it is easily seen that the relation

$(X(t{\wedge}T_{Z}(0,x^{0}_{2}))){\stackrel{{\scriptstyle\text{\rm
dist.}}}{{=}}}(Z(0,x^{0}_{2},t{\wedge}T_{Z}(0,x^{0}_{2})),x^{0}_{2}{+}Z(0,x^{0}_{2},t{\wedge}T_{Z}(0,x^{0}_{2})))$

holds by checking the jump rates.

We define, for $x{=}(x_{1},x_{2}){\in}\mathbb{N}^{2}$,

$\lambda(x)=\kappa_{0}{+}\kappa_{1}x_{1}x_{2}{+}\kappa_{2}x_{1}^{(p)}x_{2}{+}\kappa_{3}x_{1}^{(p)}x_{2}^{(2)},$

it is the total jump rate of $(X(t))$ in state $x$.

######  Lemma 19.

For $x_{1}^{0}{\geq}\kappa_{0}/(\kappa_{1}p)$,

$\limsup_{x^{0}_{2}\rightarrow+\infty}\frac{\mathbb{E}(T_{Z}(0,x^{0}_{2}))}{(x_{2}^{0})^{p-1}}\leq
C_{2},$

for some constant $C_{2}$.

###### Proof.

A simple coupling shows that the process $(Z(0,x,t)$ stopped at time
$T_{Z}(0,x)$ is lower bounded by a birth and death process $(U(t))$ starting
at $0$ with, in state $x$, a birth rate $\kappa_{0}$ and a death rate
$a_{1}{=}\kappa_{1}p(x+p)$. Denote by $H$ the hitting time of $p$ by $(U(t))$,
then it is easily seen, that, for $0{<}k{<}p$,

$(\mathbb{E}_{k}(H){-}\mathbb{E}_{k+1}(H))=\frac{a_{1}}{\kappa_{0}}(\mathbb{E}_{k-1}(H){-}\mathbb{E}_{k}(H))+\frac{1}{\kappa_{0}},$

with $\mathbb{E}_{0}(H){-}\mathbb{E}_{1}(H){=}1/\kappa_{0}$.In particular
$\mathbb{E}(T_{Z}(0,x)){\leq}\mathbb{E}_{0}(H)$. We derive the desired
inequality directly from this relation. ∎

  1. (a)

If $x_{1}{\geq}p$.
Define

$C_{1}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\sup_{x_{2}{\geq}1}\left(\frac{(x_{2}{+}p)^{(p)}{-}(x_{2})^{(p)}}{x_{2}^{p-1}}\right)<{+}\infty$

and

$\tau_{1}\stackrel{{\scriptstyle\text{\rm def.}}}{{=}}\inf\\{t{>}0:\Delta
X_{1}(t){+}\Delta X_{2}(t)\neq{-}1\\},$

where $\Delta X_{i}(t){=}X_{i}(t){-}X_{i}(t{-})$, for $i{\in}\\{1,2\\}$ and
$t{\geq}0$. The variable $\tau_{1}$ is the first instant when a reaction other
than $pS_{1}{+}2S_{2}\mathrel{\mathop{\rightharpoonup}}pS_{1}{+}S_{2}$ occurs.

For $1{\leq}k_{0}{<}x_{2}$, then

$\mathbb{P}_{x(0)}(X_{2}(\tau_{1})\leq
x_{2}{-}k_{0}-1)\geq\prod_{i=0}^{k_{0}}\frac{\kappa_{3}x_{1}^{(p)}(x_{2}{-}i)^{(2)}}{\lambda((x_{1},x_{2}{-}i))}\geq
p_{k_{0}}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\prod_{i=0}^{{k_{0}}}\frac{\kappa_{3}p^{(p)}(x_{2}{-}i)^{(2)}}{\lambda((p,x_{2}{-}i))}$

and there exists $K_{k_{0}}{\geq}k_{0}$ such that if $x_{2}{\geq}K_{k_{0}}$,
then

$(x_{2}{-}k_{0})^{(p)}{-}(x_{2})^{(p)}\leq{-}\frac{k_{0}}{2}x_{2}^{p-1}\text{
and }p_{k_{0}}\geq\frac{1}{2},$

from these relations, we obtain the inequality

(42)
$\mathbb{E}_{x(0)}\left(f_{p}\left(X(\tau_{1})\right){-}f_{p}(x)\right)\leq
1{+}\left((x_{2}{-}{k_{0}})^{(p)}{-}x_{2}^{(p)}\right)p_{k_{0}}{+}\left((x_{2}{+}1)^{(p)}{-}x_{2}^{(p)}\right)\\\
\leq\left(-\frac{k_{0}}{4}{+}1{+}C_{1}\right)x_{2}^{p-1}.$

We choose $k_{0}{=}\lceil 4(3{+}2C_{1})\rceil$, hence, for
$x_{2}{\geq}K_{k_{0}}$ the relation

$\mathbb{E}_{x(0)}\left(f_{p}\left(X(\tau_{1})\right){-}f_{p}(x)\right)\leq{-}2x_{2}^{p-1}.$

holds, and note that $\mathbb{E}(\tau_{1}){\leq}1/\kappa_{0}$.

  2. (b)

If $x_{1}{\leq}p{-}1$.
Define

$\tau_{0}=\inf\\{t{>}0:X_{1}(t){\geq}p\\},$

When $x_{1}=0$, the variable $\tau_{0}$ has the same distribution as
$T_{Z}(0,x_{2})$, otherwise it is easily seen that
$\mathbb{E}_{x(0)}[\tau_{0}]\leq\mathbb{E}[T_{Z}(0,x_{2})]$. Lemma 19 gives
therefore a constant $C_{2}>0$ so that

$\sup_{x_{2}\geq
K_{k_{0}}}\left(\frac{\mathbb{E}_{x}(\tau_{0})}{x_{2}^{p-1}}\right)<C_{2}.$

The state of the process at time $\tau_{0}$ is
$X(\tau_{0}){=}(p,x_{2}+(p-x_{1}))$, in particular

$\mathbb{E}_{x(0)}\left(f_{p}(X(\tau_{0})){-}f_{p}(x)\right)\leq
p{+}C_{1}x_{2}^{p-1},$

and at that instant, we are in case a).

The convenient stopping time is defined as
$\tau_{2}{\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}}\tau_{0}{+}\tau_{1}(\theta_{\tau_{0}})$. With $k_{0}$ and
$K_{k_{0}}$ as before, if $x_{2}{\geq}K_{k_{0}}$, by using Relation (42), we
obtain the inequality

$\mathbb{E}_{x}\left[f(X(\tau_{2}))-f(x)\right]\leq
p{+}C_{1}x_{2}^{p-1}{+}\mathbb{E}_{(p,x_{2}+(p-x_{1}))}\left[f(X(\tau_{1}))-f(x)\right]\\\
\leq
p{+}C_{1}x_{2}^{p-1}{+}\left(-\frac{k_{0}}{4}{+}1{+}C_{1}\right)\left(x_{2}+(p-x_{1})\right)^{p-1}\leq{-}x_{2}^{p-1}$

holds.

###### Proof of Proposition 18.

Theorem 6 can be used as a consequence of a), b), and Relation (40). ∎

###  8.2. A Scaling Picture

We investigate the scaling properties of $(X_{N}(t))$ when the initial state
is of the form $(N,0)$ or $(0,N)$ essentially. In the first case, an averaging
principle is proved on a convenient timescale. A time change argument is an
important ingredient to derive the main limiting result. In the second case,
the time evolution of the second coordinate of the process is non-trivial only
on “small” time intervals but with a “large” number of jumps, of the order of
$N$. This accumulation of jumps has the consequence that the convergence of
the scaled process cannot hold with the classical Skorohod topology on ${\cal
D}(\mathbb{R}_{+})$. There are better suited topologies to handle this kind of
situation. To keep the presentation simple, we have chosen to work with the
weaker topology in the space of random measures, for the occupation measures
of the sequence of scaled processes.

####  8.2.1. Horizontal Axis

For $N{\geq}1$, the initial state is $(x^{N}_{1},b)$, $b{\in}\mathbb{N}$ is
fixed, it is assumed that

(43) $\lim_{N\to+\infty}\frac{x_{1}^{N}}{N}=\alpha_{1}{>}0.$

When the process $(X_{2}(t))$ hits $0$, it happens only for a jump of ${\cal
P}_{X,1}$, all reactions but one are inactive. One has to wait for a jump of
${\cal N}_{\kappa_{0}}$ to restart the activity of the CRN.

We introduce the process $(Y_{N}(t)){=}(Y_{1}^{N}(t),Y_{2}^{N}(t))$, solution
of the SDE,

(44) $\begin{cases}\mathop{}\mathopen{}\mathrm{d}Y^{N}_{1}(t)&={\cal
N}_{\kappa_{0}}(\mathop{}\mathopen{}\mathrm{d}t){-}\mathbbm{1}_{\left\\{Y^{N}_{2}(t{-}){>}1\right\\}}{\cal
P}_{Y,1}((0,\kappa_{1}Y^{N}_{1}Y^{N}_{2}(t{-})),\mathop{}\mathopen{}\mathrm{d}t),\\\
\mathop{}\mathopen{}\mathrm{d}Y^{N}_{2}(t)&={\cal
N}_{\kappa_{0}}(\mathop{}\mathopen{}\mathrm{d}t){-}\mathbbm{1}_{\left\\{Y^{N}_{2}(t{-}){>}1\right\\}}{\cal
P}_{Y,1}((0,\kappa_{1}Y^{N}_{1}Y^{N}_{2}(t{-})),\mathop{}\mathopen{}\mathrm{d}t)\\\
&\hskip 42.67912pt{+}{\cal
P}_{Y,2}\left(\left(0,\kappa_{2}(Y^{N}_{1})^{(p)}Y^{N}_{2}(t{-})\right),\mathop{}\mathopen{}\mathrm{d}t\right)\\\
&\hskip 42.67912pt{-}{\cal
P}_{Y,3}\left(\left(0,\kappa_{3}(Y^{N}_{1})^{(p)}Y^{N}_{2}(t{-})^{(2)}\right),\mathop{}\mathopen{}\mathrm{d}t\right),\end{cases}$

with initial condition $(Y_{1}^{N}(0),Y_{2}^{N}(0)){=}(x_{1}^{N},b)$.

The process $(Y_{N}(t))$ behaves as $(X(t))$ except that its second coordinate
cannot be $0$, the associated transition is excluded. In state $(x,1)$ for
$(X(t))$, if the Poisson process ${\cal P}_{X,1}$ “rings” in Relation (35),
the state becomes $(x{-}1,0)$. It stays in this state for a duration which is
exponentially distributed with parameter $\kappa_{0}$ after which the state of
$(X(t))$ is back to $(x,1)$. These time intervals during which $(X_{2}(t))$ is
$0$ are, in some sense, “wiped out” to give $(Y_{N}(t))$. This can be
expressed rigorously via a time change argument. See Chapter 6 of [16] for
example.

Now the strategy to obtain a scaling result for $(X^{N}_{1}(t)$ is of
establishing a limit result for $(Y_{N}(t))$ and, with an appropriate change
of timescale, express the process $(X^{N}_{1}(t))$ as a “nice” functional of
$(Y_{N}(t))$.

Define

$\left(\overline{Y}_{1}^{N}(t)\right)=\left(\frac{Y_{1}^{N}(t)}{N}\right)\text{
and }\left\langle\mu_{N},f\right\rangle\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\int_{0}^{+\infty}f(s,Y_{2}^{N}(s))\mathop{}\mathopen{}\mathrm{d}s,$

if $f$ is a function on $\mathbb{R}_{+}{\times}\mathbb{N}$ with compact
support, $\mu_{N}$ is the occupation measure associated to $(Y^{N}_{2}(t))$.
See [34].

######  Proposition 20.

The sequence $(\mu_{N},(\overline{Y}_{1}^{N}(t)))$ is converging in
distribution to a limit $(\mu_{\infty},(y_{\infty}(t))$ defined by

$\left\langle\mu_{\infty},f\right\rangle=\int_{\mathbb{R}_{+}\times\mathbb{N}}f(s,x)\pi_{Y}(\mathop{}\mathopen{}\mathrm{d}x)\mathop{}\mathopen{}\mathrm{d}s,$

if $f{\in}{\cal C}_{c}(\mathbb{R}_{+}{\times}\mathbb{N})$, the function
$(y_{\infty}(t))$ is given by

(45)
$y_{\infty}(t)=\alpha_{1}\exp\left(-\frac{\kappa_{1}\kappa_{2}}{\kappa_{3}}t\right)\text{
for }t{\geq}0,$

and $\pi_{Y}$ is the distribution on $\mathbb{N}{\setminus}\\{0\\}$ defined
by, for $x{\geq}1$,

$\pi_{Y}(x)=\frac{1}{x!}\left(\frac{\kappa_{2}}{\kappa_{3}}\right)^{x}\frac{1}{e^{\kappa_{2}/\kappa_{3}}{-}1}.$

###### Proof.

The proof is quite standard. Because of the term $Y^{N}_{1}(t)^{(p)}$ in the
SDE of the process $(Y^{N}_{2}(t))$, the only (small) difficulty is to take
care of the fact that $(Y^{N}_{1}(t))$ has to be of the order of $N$,
otherwise $(Y^{N}_{2}(t))$ may not be a “fast” process. We give a sketch of
this part of the proof.

Let $a$, $b{\in}\mathbb{R}_{+}$ such that $0{<}a{<}\alpha_{1}{<}b$, and

$S_{N}=\inf\left\\{t{>}0,\overline{X}^{N}_{1}(t){\not\in}(a,b)\right\\}.$

Let $(L(t))$ a birth and death process on $\mathbb{N}$, when in state
$y{\geq}1$, its birth rate is $\beta y$ and the death rate is $\delta
y(y{-}1)$, with $\beta{=}(\kappa_{0}{+}\kappa_{2}b^{p})$ and
$\delta{=}\kappa_{3}a^{p}$. Its invariant distribution is a Poisson
distribution with parameter $\beta/\delta$ conditioned to be greater or equal
to $1$.

If $N$ is sufficiently large, we can construct a coupling of $(Y^{N}_{2}(t))$
and $(L(t))$, with $L(0){=}Y^{N}_{2}(0)$ and such that the relation

$Y^{N}_{2}(t){\leq}L(N^{p}t)$

holds for $t{\in}[0,S_{N})$.

For $t{>}0$,

$\frac{Y^{N}_{1}(t)}{N}\geq\frac{x_{1}^{N}}{N}{-}\kappa_{1}\int_{0}^{t}Y_{1}^{N}(s)Y_{2}^{N}(s)\,\mathop{}\mathopen{}\mathrm{d}s{-}M_{Y}^{N}(t),$

where $(M_{Y}^{N}(t))$ is the martingale given by

$\left(\frac{1}{N}\int_{0}^{t}\mathbbm{1}_{\left\\{Y^{N}_{2}(s{-}){>}1\right\\}}\left[{\cal
P}_{Y,1}((0,\kappa_{1}Y^{N}_{1}(s{-})Y^{N}_{2}(s{-})),\mathop{}\mathopen{}\mathrm{d}s){-}\kappa_{1}Y^{N}_{1}(s)Y^{N}_{2}(s))\mathop{}\mathopen{}\mathrm{d}s\right]\right),$

we have

(46)
$\frac{Y^{N}_{1}(t{\wedge}S_{N})}{N}\geq\frac{x_{1}^{N}}{N}{-}\kappa_{1}b\int_{0}^{t}L(N^{p}s),\mathop{}\mathopen{}\mathrm{d}s{+}M_{Y}^{N}(t{\wedge}S_{N}),$

and

$\left\langle
M_{Y}^{N}\right\rangle(t{\wedge}S_{N})\leq\frac{b}{N}\int_{0}^{t}L(N^{p}s)\mathop{}\mathopen{}\mathrm{d}s.$

By the ergodic theorem applied to $(L(t))$, almost surely

$\lim_{N\to+\infty}\int_{0}^{t}L(N^{p}s)\mathop{}\mathopen{}\mathrm{d}s=\lim_{N\to+\infty}\int_{0}^{t}\mathbb{E}(L(N^{p}s))\mathop{}\mathopen{}\mathrm{d}s\\\
=\lim_{N\to+\infty}\frac{1}{N^{p}}\int_{0}^{N^{p}t}L(s),\mathop{}\mathopen{}\mathrm{d}s=\frac{\beta}{\delta}\frac{t}{1{-}\exp({-}\beta/\delta)}.$

We deduce that $(M_{Y}^{N}(t),t\leq\eta)$ is converging in distribution to $0$
by Doob’s Inequality and, with Relation (46), that there exists $\eta{>}0$
such that

(47) $\lim_{N\to+\infty}\mathbb{P}(S_{N}{>}\eta)=1.$

For $\varepsilon{>}0$ and $K{>}0$,

$\mathbb{E}(\mu_{N}([0,\eta]{\times}[K,{+}\infty]))\leq\mathbb{E}\left(\int_{0}^{\eta{\wedge}S_{N}}\mathbbm{1}_{\left\\{Y^{N}_{2}(s)\geq
K\right\\}}\mathop{}\mathopen{}\mathrm{d}s\right){+}\eta\mathbb{P}(S_{N}\leq\eta)\\\
\leq\mathbb{E}\left(\int_{0}^{\eta}\mathbbm{1}_{\left\\{L(N^{p}s)\geq
K\right\\}}\right){+}\eta\mathbb{P}(S_{N}\leq\eta),$

again with the ergodic theorem and Relation (47), there exists some $N_{0}$
and $K{>}0$ such that
$\mathbb{E}(\mu_{N}([0,\eta]{\times}[K,{+}\infty])){\leq}\varepsilon$. Lemma
1.3 of [34] shows that the sequence of random measures $(\mu_{N})$ on
$\mathbb{R}_{+}{\times}\mathbb{N}$ restricted to $[0,\eta]{\times}\mathbb{N}$
is tight.

From there and in the same way as in Section B.2.2, it is not difficult to
conclude the proof of the proposition, on $[0,\eta]$ and extend by induction
this result on the time interval $[0,k\eta]$, for any $k{\geq}1$. ∎

Let ${\cal N}$ be a Poisson process on $\mathbb{R}_{+}^{3}$, independent of
the $\cal{P}_{Y}$ whose intensity measure is
$\mathop{}\mathopen{}\mathrm{d}s{\otimes}\mathop{}\mathopen{}\mathrm{d}t{\otimes}\kappa_{0}\exp(-\kappa_{0}a)\mathop{}\mathopen{}\mathrm{d}a$.
Recall that such a point process has the same distribution as

$\sum_{n{\geq}0}\delta_{(s_{n},t_{n},E_{n})},$

where $(s_{n})$ and $(t_{n})$ are independent Poisson processes on
$\mathbb{R}_{+}$ with rate $1$ independent of the i.i.d. sequence $(E_{n})$ of
exponential random variables with parameter $1$. See Chapter 1 of [40].

######  Definition 21 (Time Change).

Define the process $(A_{N}(t))$ by

$A_{N}(t)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(t+\int_{[0,t]\times\mathbb{R}_{+}}a\mathbbm{1}_{\left\\{Y^{N}_{2}(s{-})=1\right\\}}{\cal
N}((0,\kappa_{1}Y^{N}_{1}(s{-})Y^{N}_{2}(s{-})),\mathop{}\mathopen{}\mathrm{d}s,\mathop{}\mathopen{}\mathrm{d}a)\right),$

and its associated inverse function as

$B_{N}(t)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\inf\left\\{s>0:A_{N}(s)\geq t\right\\}.$

The instants of jump of $(A_{N}(t))$ correspond to the case when
$(Y^{N}_{2}(t))$ could switch from $1$ to $0$ for the dynamic of
$(X^{N}_{2}(t))$ and the size of the jump is the duration of time when
$(X^{N}_{2}(t))$ stays at $0$, its distribution is exponential with parameter
$\kappa_{0}$.

The process $(A_{N}(t))$ gives in fact the correct timescale to construct the
process $(X_{N}(t))$ with the process $(Y_{N}(t))$. We define the process
$(\widetilde{X}_{N}(t))$ on $\mathbb{N}^{2}$ by, for $t{\geq}0$,

(48) $\begin{cases}\widetilde{X}_{N}(A_{N}(t))=Y_{N}(t),\\\
\displaystyle\left(\widetilde{X}^{N}_{1}(u),\widetilde{X}^{N}_{2}(u)\right)=\left(Y^{N}_{1}(t{-}){-}1,0\right),\vspace{3mm}u{\in}[A_{N}(t-),A_{N}(t)).\end{cases}$

If $t$ is a jump instant of $(A_{N}(t))$, the process does not change on the
time interval $[A_{N}(t-),A_{N}(t))$. In this way, $(\widetilde{X}_{N}(t))$ is
defined on $\mathbb{R}_{+}$.

######  Lemma 22.

For $t{>}0$, then $A_{N}(B_{N}(t)){=}t$ if $t$ is not in an interval
$[A_{N}(u-),A_{N}(u))$ for some $u{>}0$, and the relation

$\sup_{t\geq 0}|\widetilde{X}_{N}(t){-}Y_{N}(B_{N}(t))|\leq 1$

holds.

###### Proof.

This is easily seen by an induction on the time intervals
$[A_{N}(s_{n}),A_{N}(s_{n+1}))$, $n{\geq}0$, where $(s_{n})$ is the sequence
on instants of jump of $(A_{N}(t))$, with the convention that $s_{0}{=}0$. ∎

######  Proposition 23.

The processes $(X_{N}(t))$ and $(\widetilde{X}_{N}(t))$ have the same
distribution.

###### Proof.

The proof is standard. The Markov property of $(\widetilde{X}_{N}(t))$ is a
consequence of the Markov property of $(Y_{N}(t))$ and the strong Markov
property of the Poisson process ${\cal N}$. It is easily checked that the
$Q$-matrices of $(X_{N}(t))$ and $(\widetilde{X}_{N}(t))$ are the same. ∎

######  Proposition 24.

For the convergence in distribution,

(49)
$\lim_{N\rightarrow+\infty}\left(\frac{A_{N}(t)}{N}\right)=(a(t))\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\alpha_{1}\frac{1}{\kappa_{0}(e^{\kappa_{2}/\kappa_{3}}-1)}\left(1{-}\exp\left({-}\frac{\kappa_{1}\kappa_{2}}{\kappa_{3}}t\right)\right)\right).$

###### Proof.

By using the fact that, for $0{\leq}u{\leq}T$,
$Y^{N}_{1}(u){\leq}x_{1}^{N}{+}{\cal P}_{Y,1}((0,\kappa_{0}){\times}(0,T])$
holds, the sequence of processes

$\left(\frac{1}{N}\int_{0}^{t}\frac{\kappa_{1}}{\kappa_{0}}\mathbbm{1}_{\left\\{Y^{N}_{2}(u){=}1\right\\}}Y^{N}_{1}(u)\mathop{}\mathopen{}\mathrm{d}u\right)$

is tight by the criterion of modulus of continuity. See Theorem 7.3 of [9] for
example. Proposition 20 shows that its limiting point is necessarily $(a(t))$.

We note that the process

$(M_{A,N}(t))=\left(\frac{1}{N}\left(A_{N}((0,t])-\frac{\kappa_{1}}{\kappa_{0}}\int_{0}^{t}\mathbbm{1}_{\left\\{Y^{N}_{2}(u){=}1\right\\}}Y^{N}_{1}(u)\mathop{}\mathopen{}\mathrm{d}u\right)\right),$

it is a square integrable martingale whose predictable increasing process is

$\left(\left\langle
M_{A,N}\right\rangle(t)\right)=\left(\frac{\kappa_{1}}{\kappa_{0}N}\int_{0}^{t}\mathbbm{1}_{\left\\{Y^{N}_{2}(u){=}1\right\\}}\frac{Y^{N}_{1}(u)}{N}\mathop{}\mathopen{}\mathrm{d}u\right).$

The martingale is vanishing as $N$ gets large by Doob’s Inequality. The
proposition is proved. ∎

Proposition 20 establishes a limit result for the sequence of processes
$(Y^{N}_{1}(t)/N)$. In our construction of $(X^{N}_{1}(t))$, time intervals,
whose lengths have an exponential distribution, are inserted. During these
time intervals the coordinates do not change. To have a non-trivial limit
result for $(X^{N}_{1}(t)/N)$, the timescale of the process has clearly to be
sped-up. It turns out that the convenient timescale for this is $(Nt)$, this
is a consequence of the convergence in distribution of $(A_{N}(t)/N)$
established in Proposition 24.

######  Proposition 25.

For the convergence in distribution, the relation

(50)
$\lim_{N\rightarrow+\infty}\left(B^{N}(Nt),t{<}t_{\infty}\right)=(a^{-1}(t))\\\
=\left({-}\frac{\kappa_{3}}{\kappa_{1}\kappa_{2}}\ln\left(\frac{\alpha_{1}{-}\kappa_{0}(e^{\kappa_{2}/\kappa_{3}}{-}1)t}{\alpha_{1}}\right),t{<}t_{\infty}\right),$

holds, where $(a(t))$ is defined in Proposition 24 and

$t_{\infty}=\frac{\alpha_{1}}{\kappa_{0}(\exp(\kappa_{2}/\kappa_{3}){-}1)}.$

###### Proof.

Note that both $(A_{N}(t))$ and $(B_{N}(t))$ are non-decreasing processes and
that the relation $A_{N}(B_{N}(t)){\geq}t$ holds for all $t{\geq}0$.

We are establishing the tightness property with the criterion of the modulus
of continuity. The constants $\varepsilon{>}0$, $\eta{>}0$ are fixed. For
$0{<}T{<}t_{\infty}$ we can choose $K{>}0$ sufficiently large so that
$a(K){>}T$ and we define

$h_{K}{=}\inf_{s{\leq}K}\left(a(s{+}\eta){-}a(s)\right),$

clearly $h_{K}{>}0$. By definition of $(B_{N}(t))$, we have

$\mathbb{P}(B_{N}(NT){\geq}K)=\mathbb{P}\left(\frac{A_{N}(K)}{N}\leq
T\right).$

The convergence of Proposition 24 shows that there exists $N_{0}$ such that if
$N{\geq}N_{0}$, the right-hand side of the last relation is less that
$\varepsilon$ and that

(51)
$\mathbb{P}\left(\sup_{0{\leq}u{\leq}K}\left|\frac{A_{N}(u{+}\eta){-}A_{N}(u)}{N}{-}(a(u{+}\eta){-}a(u))\right|\geq\frac{h_{K}}{2}\right)\leq\varepsilon$

holds.

For $\eta{>}0$, and $0{\leq}s{\leq}t{\leq}T$, if
$B_{N}(Nt){-}B_{N}(Ns){\geq}\eta$ holds, then

$A_{N}\left(B_{N}(Ns){+}\eta\right){-}A_{N}\left(B_{N}(Ns)\right)\leq
N(t{-}s)$

then, if $\delta{\leq}h_{K}/4$, for $N{\geq}N_{0}$,

$\mathbb{P}\left(\sup_{\begin{subarray}{c}0{\leq}s\leq t{\leq}T\\\
t-s\leq\delta_{0}\end{subarray}}\left|B_{N}(Nt){-}B_{N}(Ns)\right|{\geq}\eta\right)\\\
\leq\varepsilon{+}\mathbb{P}\left(\inf_{0{\leq}u{\leq}K}\left(\frac{A_{N}\left(u{+}\eta\right)}{N}{-}\frac{A_{N}\left(u\right)}{N}\right)\leq\frac{h_{K}}{4}\right)\leq
2\varepsilon,$

by Relation (51). The sequence of processes $(B_{N}(Nt))$ is therefore tight
and any of its limiting points is a continuous process. The convergence of
Proposition 24 shows that a limiting point has the same finite marginals as
the right-hand side of Relation (50). The proposition is proved. ∎

######  Theorem 26.

If $(X_{N}(t)){=}(X^{N}_{1}(t),X^{N}_{2}(t))$ is the Markov process associated
to the CRN (34) whose initial state is $(x^{N}_{1},b){\in}\mathbb{N}^{2}$,
$b{\in}\mathbb{N}$ and

$\lim_{N\to+\infty}{x_{1}^{N}}/{N}=\alpha_{1}{>}0,$

then, the convergence in distribution

$\lim_{N\to+\infty}\left(\frac{X^{N}_{1}(Nt)}{N},t{<}t_{\infty}\right)=\left(\kappa_{0}(e^{\kappa_{2}/\kappa_{3}}{-}1)(t_{\infty}{-}t),t{<}t_{\infty}\right).$

holds, with
$t_{\infty}{=}\alpha_{1}/(\kappa_{0}(\exp(\kappa_{2}/\kappa_{3}){-}1))$.

###### Proof.

Proposition 24 and 25 show that the sequence of processes

$\left(\left(\frac{Y_{N}(t)}{N},t{>}0\right),\left(B_{N}(Nt),t{<}t_{\infty}\right)\right)$

is converging in distribution to
$((y_{\infty}(t)),(a^{-1}(t),t{<}t_{\infty}))$. Consequently, the relation

$\lim_{N\to+\infty}\left(\frac{Y_{N}(B_{N}(Nt))}{N},t{<}t_{\infty}\right)=\left(y_{\infty}\left(a^{-1}(t)\right),t{<}t_{\infty}\right)$

holds for the convergence in distribution. We conclude the proof of the
proposition by using Lemma 22. ∎

####  8.2.2. Vertical Axis

For $N{\geq}1$, the initial state is $x_{N}(0){=}(a,x^{N}_{2})$, it is assumed
that $a{<}p$ and

(52) $\lim_{N\to+\infty}\frac{x_{2}^{N}}{N}=1.$

As seen in Section 8.1.2 when the first coordinate is strictly less than $p$,
with a second coordinate of the order of $N$, it takes an amount of time of
the order of $N^{p-1}$ for the process $(X_{1}^{N}(t))$ to hit $p$. See Lemma
19. In a second, short phase, a decrease of the second coordinate takes place
before returning below $p$. We now establish two limiting results

######  Lemma 27.

If $(Z(z,N,t))$ is the solution of the SDE (41) with initial state $z{<}p$,
and $T_{Z}(z,N)$ is its hitting time of $p$ then, the sequence
$(T_{Z}(z,N)/N^{p-1})$ converges in distribution to an exponential random
variable with parameter

(53) $r_{1}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\frac{\kappa_{0}}{(p{-}1)!}\left(\frac{\kappa_{0}}{\kappa_{1}}\right)^{p-1}.$

###### Proof.

The proof is standard. It can be done by induction on $p{\geq}2$ with the help
of the strong Markov property of $(Z(z,N,t))$ for example. ∎

We now study the phase during which $(X^{N}_{1}(t))$ is greater or equal to
$p$. Define $(T^{N}_{k})$ the non-decreasing sequence of stopping time as
follows, $T^{N}_{0}{=}0$ and, for $k{\geq}0$,

(54)
$T^{N}_{k+1}=\inf\\{t{\geq}T^{N}_{k}:X^{N}_{1}(t){=}p{-}1,X^{N}_{1}(t{-})=p\\}.$

######  Proposition 28 (Decay of $(X^{N}_{2}(t))$).

Under Assumption 52 for the initial condition, for the convergence in
distribution

$\lim_{N\to+\infty}\left(\frac{X^{N}_{2}(T^{N}_{1})}{X^{N}_{2}(0)},\frac{T^{N}_{1}}{X^{N}_{2}(0)^{p-1}}\right)\stackrel{{\scriptstyle\text{\rm
dist.}}}{{=}}\left(U^{\delta_{1}},E_{1}\right),$

where $U$ is a uniform random variable on $[0,1]$, independent of $E_{1}$ an
exponential random variable with parameter $r_{1}$ defined by Relation (53),
and

(55) $\delta_{1}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\frac{\kappa_{3}(p{-}1)!}{\kappa_{1}}.$

###### Proof.

Let $H_{N}$ be the hitting time of $p$ for $(X^{N}_{1}(t))$, $H_{N}$ has the
same distribution as $T_{Z}(k,x_{2}^{N})$. Its asymptotic behavior is given by
Lemma 27. Since the reactions
$pS_{1}{+}S_{2}{\rightleftharpoons}pS_{1}{+}2S_{2}$ are inactive on the time
interval $[0,H_{N}]$, we have
$X^{N}_{2}(H_{N}){=}x^{N}_{2}{+}p{-}a\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}x_{2}^{\tau}$. Let $\tau_{N}$ such that $H_{N}{+}\tau_{N}$ is the
first instant when $(X^{N}_{1}(t))$ returns to $p{-}1$. With the strong Markov
property the time origin is translated to $H_{N}$, it is enough to study the
asymptotic behavior of $X^{N}_{2}(\tau_{N})$ starting from $p$.

It is not difficult to see that, with high probability external arrivals do
not play a role during the time interval $[0,\tau_{N})$ simply because the
other reaction rates are of the order of $N$ or $N^{2}$. We will ignore them.
We can therefore assume that $X^{N}_{1}(s)=p$ until $\tau_{N}$. In the same
way it is easily seen that the sequence of random variables $(N\tau_{N})$ is
tight.

After time $0$, the reaction $x{\to}x{-}e_{2}$ occurs until time $\tau_{1,N}$
when one of the other reactions happens. Since we are interested at the final
value $X_{2}^{N}(\tau_{1,N})$, modulo a time change, it is equivalent to look
at the Markov process with $Q$-matrix

$x\longrightarrow x{+}\begin{cases}-e_{1}{-}e_{2}&\kappa_{1},\\\
e_{2}&\kappa_{2}(p{-}1)!,\\\
{-}e_{2}&\kappa_{3}(p{-}1)!(x_{2}{-}1)^{+},\end{cases}$

When $x{\to}x{-}e_{1}{-}e_{2}$ or $x{\to}x{+}e_{2}$ occurs, i.e. after an
exponentially random variable $F_{1}$ with parameter
$\kappa_{1}{+}\kappa_{2}(p{-}1)!$, the state of $(X_{2}(t))$ at this instant
is

$X^{N}_{2}(\tau_{1,N}{-})\stackrel{{\scriptstyle\text{\rm
dist.}}}{{=}}\sum_{i{=}1}^{x_{2}^{\tau}}\mathbbm{1}_{\left\\{E_{i}{\geq}F_{1}\right\\}},$

where $(E_{i})$ is an i.i.d. sequence of exponential random variables with
parameter $\kappa_{3}(p{-}1)!$, and
$\left|X^{N}_{2}(\tau_{1,N})-X^{N}_{2}(\tau_{1,N}{-})\right|\leq 1$. For the
convergence in distribution,

$\lim_{N\to+\infty}\frac{X^{N}_{2}(\tau_{1,N})}{X_{2}^{N}(0)}=\exp\left({-}\kappa_{3}(p{-}1)!F_{1}\right).$

The reaction $x{\to}x{-}e_{1}{-}e_{2}$ occurs at time $\tau_{1,N}$ with
probability $1{-}q_{1}$, with

$q_{1}{=}\frac{\kappa_{2}(p{-}1)!}{\kappa_{1}{+}\kappa_{2}(p{-}1)!},$

and in this case $\tau_{N}{=}\tau_{1,N}$. Otherwise, there is a new cycle of
length $\tau_{2,N}$ and that

$\lim_{N\to+\infty}\frac{X^{N}_{2}(\tau_{1,N}){+}\tau_{2,N})}{X_{2}^{N}(0)}=\exp\left({-}\kappa_{3}(p{-}1)!(F_{1}{+}F_{2})\right),$

where $(F_{i})$ is an i.i.d. sequence with the same distribution as $F_{1}$.
By induction we obtain the convergence in distribution

$\lim_{N\to+\infty}\frac{X^{N}_{2}(\tau_{N})}{X_{2}^{N}(0)}=\exp\left({-}\kappa_{3}(p{-}1)!\sum_{1}^{G}F_{i}\right),$

where $G$ is a random variable independent of $(F_{i})$ with a geometric
distribution with parameter $q_{1}$, $\mathbb{P}(G{\geq}n){=}q_{1}^{n-1}$ for
$n{\geq}1$. Trivial calculations gives the desired representation. ∎

In view of the last result it is natural to expect that the convergence of the
scaled process $(X^{N}_{2}(t/N^{p-1})/N$ to a Markov process with jumps. The
only problem is that, as we have seen in the last proof, there is a large
number of jumps downwards, of the order of $N$, on the time interval of length
$\tau_{N}$ of the previous proof. Event if $\tau_{N}$ is arbitrarily small
when $N$ gets large, there cannot be convergence in the sense of the classical
Skorohod topology. There are topologies on the space of càdlàg functions
${\cal D}(\mathbb{R}_{+})$ for which convergence in distribution may hold in
such a context. See [29] for example. For the sake of simplicity, we present a
convergence result formulated for a weaker topology expressed in terms of
occupation measures.

We now introduce a Markov process on $(0,1]$ as the plausible candidate for a
limiting point of $(X^{N}_{2}(t/X^{N}_{2}(0)^{p-1})/N$.

######  Definition 29.

The infinitesimal generator ${\cal A}$ of a Markov process $(U(t))$ on $(0,1]$
is defined by, for $f{\in}{\cal C}_{c}((0,1])$,

(56) ${\cal
A}(f)(x)=\frac{r_{1}}{x^{p-1}}\mathbb{E}\left(f\left(xU^{\delta_{1}}\right){-}f(x)\right),\quad
x{\in}(0,1],$

where $r_{1}$, $\delta_{1}$ are constants defined by Relations (53) and (55).

Analytically the operator ${\cal A}$ can be expressed as

${\cal
A}(f)(x)=\frac{r_{1}}{x^{p-1}}\int_{0}^{1}\left(f\left(xu^{\delta_{1}}\right){-}f(x)\right)\mathop{}\mathopen{}\mathrm{d}u,\quad
x{\in}(0,1].$

######  Proposition 30.

If $(U(t))$ is a Markov process on $(0,1]$ with infinitesimal generator ${\cal
A}$, then, with probability $1$, it is an explosive process converging to $0$.

###### Proof.

Assume that $U(0){=}\alpha{\in}(0,1]$. By induction, the sequence of states
visited by the process has the same distribution as $(V_{n})$ with, for
$n{\geq}0$,

$V_{n}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\alpha\exp\left({-}\delta_{1}\sum_{i=1}^{n}E_{i}\right),$

where $(E_{i})$ is an i.i.d. sequence of exponentially distributed random
variables with parameter $1$. The sequence of instants jumps has the same
distribution as

$\left(t^{V}_{n}\right)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\sum_{i=1}^{n}(V_{i-1})^{p-1}\frac{\phi_{i}}{r_{1}}\right),$

where $(\phi_{i})$ is an i.i.d. sequence of exponentially distributed random
variables with parameter $1$, independent of $(E_{i})$. It is easily seen that
$(V_{n})$ converges to $0$ almost surely and that the sequence
$(\mathbb{E}(t^{V}_{n}))$ has a finite limit. The proposition is proved. ∎

######  Definition 31 (Scaled occupation measure of $(X^{N}_{2}(t))$).

For $N{\geq}1$, $\Lambda_{N}$ is the random measure on
$\mathbb{R}_{+}{\times}(0,1]$ defined by, for $f{\in}{\cal
C}_{c}(\mathbb{R}_{+}{\times}(0,1])$,

(57)
$\left\langle\Lambda_{N},f\right\rangle=\frac{1}{N^{p-1}}\int_{0}^{+\infty}f\left(\frac{s}{N^{p-1}},\frac{X^{N}_{2}(s)}{N}\right)\mathop{}\mathopen{}\mathrm{d}s.$

We can now state our main scaling result for large initial states near the
vertical axis.

######  Theorem 32.

If $(X_{N}(t)$ is the Markov process associated to the CRN (34) whose initial
state is $(a,x^{N}_{2}){\in}\mathbb{N}^{2}$, $a{\leq}p{-}1$, and such that

$\lim_{N\to+\infty}{x_{2}^{N}}/{N}=\alpha{>}0,$

then the sequence $(\Lambda_{N})$ defined by Relation (57) converges in
distribution to $\Lambda$, the occupation measure of $(V(t))$ a Markov process
with infinitesimal generator ${\cal A}$ starting at $\alpha$, i.e. for
$f{\in}{\cal C}_{c}(\mathbb{R}_{+}{\times}(0,1])$,

$\left\langle\Lambda,f\right\rangle=\int_{0}^{+\infty}f(s,V(s))\mathop{}\mathopen{}\mathrm{d}s.$

###### Proof.

Without loss of generality, due to the multiplicative properties of the
convergence, see Proposition 28, we can take $\alpha{=}1$ and assume that
$X^{N}_{2}(0){=}N$. Recall that the Laplace transform of a random measure $G$
on $\mathbb{R}_{+}{\times}(0,1]$ is given by

${\cal L}_{G}(f)\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\mathbb{E}\left(\exp\left(-\left\langle
G,f\right\rangle\right)\right),$

for a non-negative function $f{\in}{\cal C}_{c}(\mathbb{R}_{+}{\times}(0,1])$.
See Section 3 of [13].

To prove the convergence in distribution of $(\Lambda_{N})$ to $\Lambda$, it
is enough to show that the convergence

$\lim_{N\to+\infty}{\cal L}_{\Lambda_{N}}(f)={\cal L}_{\Lambda}(f),$

holds for all non-negative functions $f{\in}{\cal
C}_{c}(\mathbb{R}_{+}{\times}(0,1])$. See Theorem 3.2.6 of [13] for example.

If $f{\in}{\cal C}_{c}(\mathbb{R}_{+}{\times}(0,1])$, its support is included
in some $[0,T]{\times}(\eta,1]$, for $\eta{>}0$ and $T{>}0$.

Let $(T^{N}_{k})$ the sequence of stopping times defined by Relation (54). The
Laplace transform of $\Lambda_{N}$ at $f$ is given by

(58) ${\cal
L}_{\Lambda_{N}}(f)=\mathbb{E}\left(\exp\left({-}\sum_{k{\geq}0}\int_{T^{N}_{k}/N^{p-1}}^{T^{N}_{k+1}/N^{p-1}}f\left(s,\frac{X^{N}_{2}\left(T^{N}_{k}\right)}{X^{N}_{2}(0)}\right)\mathop{}\mathopen{}\mathrm{d}s\right)\right).$

Let $(t^{V}_{k},V_{k})$ be the sequence of couples of instants of jumps and
its value of the Markov process $(V(t))$, as defined in the proof of
Proposition 30. For $\varepsilon{>}0$, there exists some $n_{0}$ such that

$\mathbb{P}\left(\alpha\prod_{i=1}^{n_{0}}V_{i}\geq\frac{\eta}{2}\right)\leq\varepsilon/2,$

holds, and, consequently,

(59) $\left|{\cal
L}_{\Lambda}(f)-\mathbb{E}\left(\exp\left({-}\sum_{k=0}^{n_{0}{-}1}\int_{t^{V}_{k}}^{t^{V}_{k+1}}f\left(s,V_{k}\right)\mathop{}\mathopen{}\mathrm{d}s\right)\right)\right|\leq\varepsilon.$

Proposition 28 shows that, for the convergence in distribution,

$\lim_{N\to+\infty}\left(\frac{X^{N}_{2}(T^{N}_{k+1})}{X^{N}_{2}(T^{N}_{k})},\frac{T^{N}_{k+1}{-}T^{N}_{k}}{X^{N}_{2}(T^{N}_{k})^{p-1}},k{\geq}0\right)=\left(U_{k}^{\delta_{1}},E_{k},k{\geq}0\right),$

where $(U_{k})$ and $(E_{k})$ are i.i.d. independent sequence of random
variables whose respective distributions are uniform on $[0,1]$, and
exponential with parameter $r_{1}$. Hence, there exists $N_{0}$ such that if
$N{\geq}N_{0}$, then

(60) $\begin{cases}\displaystyle\left|{\cal
L}_{\Lambda_{N}}(f)-\mathbb{E}\left(\exp\left({-}\sum_{k=0}^{n_{0}-1}\int_{T^{N}_{k}/N^{p-1}}^{T^{N}_{k+1}/N^{p-1}}f\left(s,\frac{X^{N}_{2}\left(T^{N}_{k}\right)}{X^{N}_{2}(0)}\right)\mathop{}\mathopen{}\mathrm{d}s\right)\right)\right|\leq
2\varepsilon,\\\
\displaystyle\mathbb{P}\left(\frac{X^{N}_{2}(T^{N}_{k+1})}{X^{N}_{2}(T^{N}_{k})}{\leq}1,\forall
k{\in}\\{0,\ldots,n_{0}\\}\right)\geq 1{-}\varepsilon\end{cases}$

Define, for $n{>}0$,

$(I^{N}_{n})\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\left(\sum_{k=0}^{n-1}\int_{T^{N}_{k}/N^{p-1}}^{T^{N}_{k+1}/N^{p-1}}f\left(s,\frac{X^{N}_{2}\left(T^{N}_{k}\right)}{X^{N}_{2}(0)}\right)\mathop{}\mathopen{}\mathrm{d}s\right)$

In views of Relations (59) and (60), all we have to do is to prove, for every
$n{>}0$, the convergence in law of $(I^{N}_{n})$ to

$I_{n}\stackrel{{\scriptstyle\text{\rm
def.}}}{{=}}\int_{0}^{t^{V}_{n}}f\left(s,V(s)\right)\mathop{}\mathopen{}\mathrm{d}s=\sum_{k=0}^{n-1}\int_{t^{V}_{k}}^{t^{V}_{k+1}}f\left(s,V_{k}\right)\mathop{}\mathopen{}\mathrm{d}s,$

as $N$ gets large.

We will prove by induction on $n{>}0$, the convergence in distribution

$\lim_{N\to+\infty}\left(I^{N}_{n},\left|\ln\left(\frac{X^{N}_{2}\left(T^{N}_{n}\right)}{X^{N}_{2}(0)}\right)\right|,\frac{T^{N}_{n}}{X^{N}_{2}(0)^{p-1}}\right)\\\
=\left(\int_{0}^{t^{V}_{n}}f\left(s,V(s)\right)\mathop{}\mathopen{}\mathrm{d}s,\left|\ln\left(V_{n}\right)\right|,t^{V}_{n}\right).$

We will show the convergence of the Laplace transform of the three random
variables taken at $(a,b,c)$, for $a$, $b$, $c{>}0$.

For $n=1$, this is direct consequence of Proposition 28. If it holds for
$n{\geq}1$, the strong Markov property of $(X^{N}(t))$ for the stopping time
$T^{N}_{n}$ gives the relation

$H_{N}(a,b,c)\stackrel{{\scriptstyle\text{\rm