# A partition method for bounding continuous-time Markov chain models of general reaction network

Guillaume Ballif, Laurent Pfeiffer, Jakob Ruess

<figure><img src="image_1.png" /></figure>

## To cite this version:

Guillaume Ballif, Laurent Pfeiffer, Jakob Ruess. A partition method for bounding continuous-time Markov chain models of general reaction network. 2025. hal-05319664

HAL Id: hal-05319664

https://hal.science/hal-05319664v1

Preprint submitted on 17 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A partition method for bounding continuous-time Markov chain models of general reaction network

Guillaume Ballif*†, Laurent Pfeiffer‡ and Jakob Ruess§

**Abstract:** In this work, we present a general method to establish properties of multi-dimensional continuous-time Markov chains representing stochastic reaction networks. This method consists of grouping states together (via a partition of the state space), then constructing two one-dimensional birth and death processes that lower and upper bound the initial process under simple assumptions on the infinitesimal generators of the processes. The construction of these bounding processes is based on coupling arguments and transport theory. The bounding processes are easy to analyse analytically and numerically and allow us to derive properties on the initial continuous-time Markov chain. We focus on two important properties: the behavior of the process at infinity through the existence of a stationary distribution and the error in truncating the state space to numerically solve the master equation describing the time evolution of the probability distribution of the process. We derive explicit formulas for constructing the optimal bounding processes for a given partition, making the method easy to use in practice. We finally discuss the importance of the choice of the partition to obtain relevant results and illustrate the method on two examples of chemical reaction network.

**Mathematics Subject Classification:** 60J27, 60G10, 92C40.

## 1 Introduction

Stochastic reaction networks have been widely used in many fields such as chemistry [5], biology [2, 19] and ecology, medicine [12], epidemiology [7], economics [10] to model complex and often multi-dimensional random phenomena. Mathematically, it is straightforward to represent such models as continuous-time Markov chains (CTMCs) but the multi-dimensional nature of the real process greatly complicates the study of properties of the corresponding CTMC (such as explosion in finite time and the behaviour when time tends to infinity) and the numerical simulation of trajectories. For example, the existence of a stationary distribution is a crucial property in the study of the long-term behaviour but whether or not a given CTMC admits a stationary distribution is often hard to say *a priori*. This question has been studied extensively and several methods have been developed to address it: coupling techniques [29], conditions on the semi-group or the study of the generator spectrum [8, 17] or existence of a Lyapunov function [14, 31]. However, the size of the state space and the assumptions to be verified make these methods difficult to apply systematically to multidimensional CTMCs. In some cases, structural properties of CTMCs allow one to obtain results on the existence or even the form of the stationary distribution. These properties may relate to the transition rates (complex balanced system [3, 4], mass-action kinetics and extension [6]), the structure of the reaction network (state space decomposition [27]) or the size of the state space (one-dimensional CTMC [26]). However, methods that exploit such properties cannot be applied systematically to multidimensional CTMCs due to the restrictive assumptions required.

From a numerical point of view, the time evolution of the probability distribution of a CTMC is described by a master equation (also called Kolmogorov equation) composed of an infinite system of ordinary differential equations. The size of the state space often makes solving the master equation numerically impossible due to the exponential increase in computation time with the state space dimension. Two classes of methods have been developed in the literature to circumvent this problem: Monte Carlo methods and approximation methods.

<footnote>*Inria Saclay, 91120 Palaiseau, France. Mail: guillaume.ballif@inria.fr</footnote>

<footnote>†Corresponding author</footnote>

<footnote>‡Université Paris-Saclay, CNRS, CentraleSupélec, Inria, Laboratoire des signaux et systèmes, and Fédération de Mathématiques de CentraleSupélec, 91190, Gif-sur-Yvette, France. Mail: laurent.pfeiffer@inria.fr</footnote>

<footnote>§Inria Saclay, 91120 Palaiseau, France. Mail: jakob.ruess@inria.fr</footnote>
Monte Carlo methods [22] are based on the simulation of a large number of trajectories of the process with the stochastic simulation algorithm (SSA) [22] also known as the Gillespie algorithm in order to reconstruct the probability distribution of the CTMC. Various improvements have been developed over the past few years to improve the precision and computation speed of the method [1, 23]. However, a large number of trajectories of the process is required to obtain good approximations of the probability distribution, which is very costly in terms of computing time.

The second type of methods aims to solve an approximation of the master equation [9, 18, 28]. One of the most commonly used methods is the truncation of the state space to replace numerical solving of the master equation by that of a linear system of ordinary differential equations [25, 32, 35] but in the case of a multi-dimensional CTMC, even after truncation the state space is often still much too large for numerical solution of this linear system to be efficient or even just feasible. Moreover, most of these methods do not provide an estimate of the error caused by this approximation that can be calculated in practice.

In this work, we present a general method for constructing one-dimensional bounding processes (in a certain sense) for any CTMC under some simple assumptions on the infinitesimal generator. We can then analyse the one-dimensional bounding processes to derive properties of the multi-dimensional original CTMC that would be hard to obtain directly. Concretely, we focus on two essential aspects when studying CTMC models: the behavior of the process at infinity through the existence of a stationary distribution and the quantification of the approximation error when the infinite state space is truncated. The method for constructing the bounding processes is based on an ordered partition of the state space (that we will refer to as classes) that allows us to reduce the dimension of the CTMC by studying its behavior only on classes. Coupling arguments and transport theory are used to derive these bounding results (upper and lower bounds).

Different methods have been developed in order to reduce the state space of a CTMC to be able to derive properties on this CTMC: aggregation of states (also called lumping) for Markov chain on a finite state space [11] or for some particular cases of Markov chains on an infinite state space [21, 30, 33], decomposition of the reaction network into subnetworks [24] or comparisons between processes through coupling [20]. Some of these methods assume the existence of a stationary distribution in order to deduce its properties (sensitivity to parameters and explicit form) [13].

To make our framework useful in practice, we derive explicit formulas for constructing the infinitesimal generator of the optimal bounding processes from assumptions on the generator of the original process. These explicit formulas are used to illustrate in an example the results that we can derive on the multi-dimensional original CTMC. The choice of optimal classes to partition the state space (i.e. leading to the best possible bounds) is also discussed by comparing the upper-bounding process and properties derived from a naive choice and those derived from a set of classes taking into account the structure of the reaction network.

The manuscript is organized as follows. Section 2 focuses on the construction of an upper-bounding process. The construction is done in two main steps. In Theorem 2.1, we construct an upper-bounding process based on an infinitesimal generator satisfying some monotonicity properties. We provide in Proposition 2.2 an explicit formula for such an infinitesimal generator that most closely upper-bounds the initial process. The adaptation of these results for a lower-bounding process is presented in Section 3. In Section 4, the method is used to derive results on the existence of a stationary distribution and on the error made in numerically solving the master equation on truncated state spaces and all these results are illustrated on an example.

**Notation** We fix a countable state space $\mathcal{S}$ and a partition $(S_\ell)_{\ell \in \mathbb{N}}$ of $\mathcal{S}$, i.e., for all $\ell \in \mathbb{N}$, $S_\ell \subset \mathcal{S}$ and

$$ \mathcal{S} = \bigcup_{\ell \in \mathbb{N}} S_\ell, \quad \text{and} \quad S_\ell \cap S_m = \emptyset, \quad \forall \ell \neq m. $$

In the following, we will assume that **for all** $\ell \in \mathbb{N}$, **the set S_\ell is finite**. No other assumption on the partition will be required.

**Definition 1.1**

For any $i \in \mathcal{S}$, we will call "class" of $i$ the unique integer $\ell \in \mathbb{N}$ such that $i \in S_\ell$. We denote it by $\text{cl}(i)$.

**Definition 1.2**

Let $X = (X_t)_{t \ge 0}$ and $Y = (Y_t)_{t \ge 0}$ be two Markov chains on $\mathcal{S}$ and $\mathbb{N}$, respectively. We say that $Y$ is
an upper-bounding process (of X) if

$$Y_t \ge \mathrm{cl}(X_t), \quad \forall t \ge 0, \text{ almost surely.}$$

We say that Y is a lower-bounding process (of X) if

$$Y_t \le \mathrm{cl}(X_t), \quad \forall t \ge 0, \text{ almost surely.}$$

As we will consider Markov chains on $S$, $N$, and $S \times N$, we will make use of the following convention:

the letters $i$ and $j$ denote states in $S$,

the letters $k$, $\ell$, and $m$ denote states in $N$.

In the article, we will abusively call *matrix* any real-valued mapping with 2 or more input variables, taking values in a countable set. For any countable set $\Omega$, we denote as $Q(\Omega)$ the set of matrices $Q = (Q_{i,j})_{i,j \in \Omega}$ such that

* for all $i \in \Omega$, for all $j \in \Omega \setminus \{i\}$, $Q_{i,j} \ge 0$,

* for all $i \in \Omega$, $\sum_{j \in \Omega \setminus \{i\}} Q_{i,j} < \infty$,

* for all $i \in \Omega$, $\sum_{j \in \Omega} Q_{i,j} = 0$.

We will call such matrices *transition-rate matrices*. Note that we allow the matrix $Q$ to have an infinite number of non-zero values. In the case of reaction networks, this generality allows us to consider an infinite number of reactions. Some special types of reaction networks are composed of an infinite number of reactions (such as polymerization reactions [15]). However, most reaction networks are composed of a finite number of reactions but may be composed of non-standard reactions that allow for transitions to infinitely many different states (such as the production of proteins in bursts with stochastic burst sizes) [34].

In the following, the transition-rate matrix of the Markov chain of interest will be denoted by $Q$. The aim of this article is to construct the transition-rate matrix of the upper- and lower-bounding processes that will be denoted by $\tilde{Q} = (\tilde{Q}_{\ell,m})_{\ell,m \in N} \in Q(N)$. We will make use of the following notations for $Q$ and $\tilde{Q}$:

$$Q_{i,\ell} = \sum_{j \in S_{\ell}} Q_{i,j}, \quad \forall i \in S, \forall \ell \in N,$$

$$\tilde{Q}_{\ell,0:m} = \sum_{k=0}^{m} \tilde{Q}_{\ell,k}, \quad \forall \ell, m \in N, \qquad (1)$$

$$\tilde{Q}_{\ell,m:\infty} = \sum_{k=m}^{\infty} \tilde{Q}_{\ell,k}, \quad \forall \ell, m \in N.$$

These notations will also be used for general matrices (which are not necessarily transition-rate matrices). For example, we will use the notation

$$Q_{i,0:m} = \sum_{k=0}^{m} Q_{i,k} = \sum_{k=0}^{m} \sum_{j \in S_k} Q_{i,j}, \quad \forall i \in S, \forall m \in N.$$

Note that for all $\hat{Q} \in Q(\Sigma)$, for all $(\ell, m) \in \Sigma^2$, $\hat{Q}_{\ell,0:m} = -\hat{Q}_{\ell,m+1:\infty}$.

Finally, we recall the master equation governing the time evolution of the probability distribution of a continuous-time Markov chain $X$ on $S$ with transition-rate matrix $Q$:

$$\frac{dp}{dt}(t,j) = \sum_{i \in S} Q_{i,j} p(t,j),$$

where $p(t,i) = \mathbb{P}(X_t = i)$. This equation can be formulated in a vector form:

$$\frac{dp}{dt}(t) = p(t)Q,$$

where $p(t) = (p(t,i))_{i \in S}$ is a row vector.
# 2 Construction of an upper-bounding process

## 2.1 Main result

We state here our first main result, which gives two simple conditions on the matrix $\tilde{Q}$ of the process Y for which Y is an upper-bounding process of X.

**Theorem 2.1**

Let $Q \in Q(S)$ and let $\tilde{Q} \in Q(\mathbb{N})$. Let $p^0$ be a distribution on $S$. Assume the two following properties:

$$ (\text{A1}) \quad \forall \ell, m \in \mathbb{N}, \forall i \in S_{\ell}, \tilde{Q}_{\ell,0:m} \le Q_{i,0:m}, $$

$$ (\text{A2}) \quad \forall \ell, m \in \mathbb{N} \text{ such that } \ell \neq m, \tilde{Q}_{\ell+1,0:m} \le \tilde{Q}_{\ell,0:m}. $$

Then there exist two continuous-time Markov chains $X = (X_t)_{t \ge 0}$ and $Y = (Y_t)_{t \ge 0}$, defined on $S$ and $\mathbb{N}$ respectively, with transition-rate matrices $Q$ and $\tilde{Q}$, such that $X_0$ is distributed according to $p^0$, and such that Y is an upper-bounding process of X.

Subsections 2.2 and 2.3 are dedicated to the proof of the theorem. The main idea of our analysis consists in constructing $(X_t, Y_t)_{t \ge 0}$ as a Markov chain on $S \times \mathbb{N}$ with a transition-rate matrix $R$ satisfying certain properties (see Lemma 2.3). In particular, $R$ will be such that no transition towards a state $(i, \ell)$ such that $cl(i) > \ell$ is possible. We will provide the reader with an explicit construction of $R$, based on a technique from transport theory (see Proposition 2.1). Subsection 2.4 will be dedicated to the construction of an explicit transition-rate matrix $\tilde{Q}$ satisfying Assumptions **A1** and **A2**.

**Remark 2.1**

Theorem 2.1 shows that for constructing an upper-bounding process, it suffices to find a transition-rate matrix $\tilde{Q}$ such that Assumptions **A1** and **A2** are satisfied. For a given $Q$ and given partition $(S_\ell)_{\ell \in \mathbb{N}}$, we provide in Proposition 2.2 a construction of a transition-rate matrix $\tilde{Q}$ such that these assumptions hold true. Let us stress that this construction will not require any additional assumption on $\tilde{Q}$ and $(S_\ell)_{\ell \in \mathbb{N}}$, and so the two theoretical results (Theorem 2.1 and Proposition 2.2) building our approach do not induce any restriction on the nature of the reaction network.

Note that different partitions lead to different upper-bounding processes through the assumptions **A1** and **A2**. So a suitable choice of a partition may lead to a tighter upper-bounding process, as discussed in Section 4.2 in a concrete example.

We give an interpretation of the two standing assumptions in the next remarks. Note that when $m \ge \ell$, the term $\tilde{Q}_{\ell,0:m}$ has no straightforward interpretation, as the underlying sum contains the diagonal negative term $\tilde{Q}_{\ell,\ell}$. So in this case, we reformulate our assumptions with the help of the equality $\tilde{Q}_{\ell,m+1:\infty} = -\tilde{Q}_{\ell,0:m}$.

**Remark 2.2**

Assumption **A1** essentially says that Y decreases at a slower rate than X moves from states to other states that belong to a lower class. Vice-versa, Y increases at a faster rate than X moves from states to other states that belong to a higher class. This can be formalized as follows.

* When $m < \ell$, the assumption is equivalent to:

$$ \tilde{Q}_{\ell,0:m} \le \min_{i \in S_{\ell}} Q_{i,0:m}. $$

The quantity $\tilde{Q}_{\ell,0:m}$ is the rate at which Y moves from $\ell$ to classes that are smaller or equal to m and the quantity $\min_{i \in S_{\ell}} Q_{i,0:m}$ is the minimal rate at which X transitions from a state in class $\ell$ to another state in a class lower or equal to m.

* When $m \ge \ell$, the assumption is equivalent to:

$$ \tilde{Q}_{\ell,m+1:\infty} \ge \max_{i \in S_{\ell}} Q_{i,m+1:\infty}. $$

The quantity $\tilde{Q}_{\ell,m+1:\infty}$ is the rate at which Y moves from $\ell$ to classes that are strictly greater than m and the quantity $\max_{i \in S_{\ell}} Q_{i,m+1:\infty}$ is the maximal rate at which X transitions from a state in class $\ell$ to another state in a class strictly higher than m.
**Remark 2.3**

Assumption A2 essentially says that Y moves to low classes at a faster rate when the initial state is small, and vice-versa: Y moves to larger classes at a faster rate when the initial state is large. This can be formalized as follows.

* When $m < \ell$, the assumption is equivalent to:

$$ \tilde{Q}_{\ell+1,0:m} \le \tilde{Q}_{\ell,0:m}. $$

The quantities $\tilde{Q}_{\ell+1,0:m}$ and $\tilde{Q}_{\ell,0:m}$ describe the rate at which Y moves from class $\ell+1$, respectively from class $\ell$, to classes that are lower or equal to $m$.

* When $m > \ell$, the assumption is equivalent to:

$$ \tilde{Q}_{\ell+1,m+1:\infty} \ge \tilde{Q}_{\ell,m+1:\infty}. $$

The quantities $\tilde{Q}_{\ell+1,m+1:\infty}$ and $\tilde{Q}_{\ell,m+1:\infty}$ describe the rates at which Y moves from class $\ell+1$, respectively from class $\ell$, to classes that are strictly greater than $m$.

## 2.2 A preliminary result

We denote by $\ell^1(\mathbb{N}, \mathbb{R}_+)$ the set of sequences $u = (u_k)_{k \in \mathbb{N}}$ of non-negative real numbers whose sum is finite. Following the convention introduced in the introduction, we denote $u_{0:\infty} = \sum_{k=0}^{\infty} u_k$ for a given $u \in \ell^1(\mathbb{N}, \mathbb{R}_+)$. Similarly, we define $\ell^1(\mathbb{N}^2, \mathbb{R}_+)$ as the set of doubly-indexed sequences $\pi = (\pi_{k,\ell})_{(k,\ell) \in \mathbb{N}^2}$ of non-negative numbers whose sum is finite. Given $\pi \in \ell^1(\mathbb{N}^2, \mathbb{R}_+)$, $k \in \mathbb{N}$, and $\ell \in \mathbb{N}$, we denote $\pi_{k,0:\infty} = \sum_{\ell=0}^{\infty} \pi_{k,\ell}$ and $\pi_{0:\infty,\ell} = \sum_{k=0}^{\infty} \pi_{k,\ell}$.

**Proposition 2.1**

Let $a$ and $b$ be in $\ell^1(\mathbb{N}, \mathbb{R}_+)$. Assume that $a_{0:k} \ge b_{0:k}$, for all $k \in \mathbb{N}$. Assume further that $a_{0:\infty} = b_{0:\infty}$. Then there exists $\pi \in \ell^1(\mathbb{N}^2, \mathbb{R}_+)$ satisfying the following properties:

$$ \pi_{k,0:\infty} = a_k, \quad \pi_{0:\infty,\ell} = b_\ell, \quad k > \ell \Rightarrow \pi_{k,\ell} = 0, \quad \forall (k, \ell) \in \mathbb{N}^2. \quad (2) $$

The rest of the subsection is dedicated to the construction of such a sequence $\pi$. One can understand Proposition 2.1 from the point of view of optimal transportation theory: a commodity is available in sources, labeled by $k=0, 1, \dots$, and must be shipped to some factories, labeled by $\ell=0, 1, \dots$. Then $a_k$ represents the amount of commodity available at source $k$, $b_\ell$ the capacity of factory $\ell$, and $\pi_{k,\ell}$ is the amount of commodity that is shipped from the source $k$ to the factory $\ell$. The first two conditions in (2) simply say that there is no loss of commodity. The third condition says that the commodity should be shipped from a given source to a factory with a greater or equal label.

In order to prove Proposition 2.1, we first consider the problem of shipping an amount $x$ of commodity (from a single source) to factories of capacities $u_0, u_1, \dots$ exploiting them at their maximal capacity in increasing order. Concretely, if $x \le u_0$, we ship the whole amount to factory 0. Otherwise, we ship the amount $u_0$ to the factory 0. Then the remaining amount $x - u_0$ is entirely shipped to the factory 1 if $x - u_0 \le u_1$. Otherwise, we ship the amount $u_1$, and so on. This leads to the following mathematical definition. Given $u \in \ell^1(\mathbb{N}, \mathbb{R}_+)$ and $x \in [0, u_{0:\infty}]$, we define the sequence $(\Pi[x, u]_k)_{k \in \mathbb{N}}$ by

$$ \Pi[x, u]_k = \min(\max(x - u_{0:k-1}, 0), u_k), $$

with the convention that $u_{0:-1} = 0$. Two cases must be considered in order to interpret the formula. If $x \le u_{0:k-1}$, then the total capacity of the factories $0, \dots, k-1$ is larger than the amount of commodity, which does not need to be shipped to factory $k$. In this case we indeed have $\Pi[x, u]_k = \min(0, u_k) = 0$. Otherwise, when $x > u_{0:k-1}$, the factories $0, \dots, k-1$ should be used at their full capacity. The above formula simplifies to $\Pi[x, u]_k = \min(x - u_{0:k-1}, u_k)$. The quantity $x - u_{0:k-1} > 0$ represents the remaining amount to be allocated. If this quantity is less than $u_k$, it can be fully allocated to factory $k$, otherwise, we restrict it to $u_k$.

We have the following technical result.

**Lemma 2.1**

Let $u \in \ell^1(\mathbb{N}, \mathbb{R}_+)$ and let $x \in [0, u_{0:\infty}]$. The five following statements hold true.
1. For all $k \in \mathbb{N}$, $\Pi[x, u]_k \le u_k$.

2. For any $k \in \mathbb{N}$, if $u_{0:k} \le x$, then $\Pi[x, u]_k = u_k$.

3. If $x = u_{0:\infty}$, then $\Pi[x, u] = u$.

4. We have $\Pi[x, u] \in \ell^1(\mathbb{N}, \mathbb{R}_+)$ and in particular, $\Pi[x, u]_{0:\infty} = x$.

5. For any $x' \in [0, x]$, for all $k \in \mathbb{N}$, we have

$$0 \le \Pi[x, u]_k - \Pi[x', u]_k \le x - x'$$

**Proof:** 1. Straightforward.

2. Let $k \in \mathbb{N}$ be such that $u_{0:k} \le x$. Then we have $x - u_{0:k-1} = x - u_{0:k} + u_k \ge u_k \ge 0$. Therefore $\Pi[x, u]_k = \min(x - u_{0:k-1}, u_k) = u_k$ as was to be shown.

3. If $x = u_{0:\infty}$, then we have $u_{0:k} \le x$ for all $k \in \mathbb{N}$, and therefore $\Pi[x, u]_k = u_k$ for any $k \in \mathbb{N}$ (by the previous point).

4. From the definition of $\Pi[x, u]$, we have $\Pi[x, u]_k \ge 0$ for any $k \in \mathbb{N}$. We only have to prove that $\Pi[x, u]_{0:\infty} = x$.

If $x = u_{0:\infty}$, we have $\Pi[x, u] = u$ (by the previous point) and the claim holds true. Otherwise, if $x < u_{0:\infty}$, there exists $K \in \mathbb{N}$ such that

$$u_{0:k-1} < x \le u_{0:k}.$$

It is easy to verify that for $k = 0, \dots, K-1$, $\Pi[x, u]_k = u_k$, that for $k = K$, $\Pi[x, u]_K = x - u_{0:K-1}$ and that for $k > K$, $\Pi[x, u]_k = 0$. The claim follows.

5. Let $k \in \mathbb{N}$. Then $\Pi[x, u]_k$ results from the composition of two functions, that is, $\Pi[x, u]_k = \alpha_k \circ \beta_k(x)$, where

$$\beta_k(x) = \max(x - u_{0:k-1}, 0) \quad \text{and} \quad \alpha_k(y) = \min(y, u_k).$$

Both functions $\alpha_k$ and $\beta_k$ are non-decreasing and Lipschitz continuous with modulus 1. Therefore, the composition $\alpha_k \circ \beta_k$ is also non-decreasing and Lipschitz continuous with modulus 1. The claim immediately follows.

Given $a$ and $b$ in $\ell^1(\mathbb{N}, \mathbb{R}_+)$ satisfying the assumptions of Proposition 2.1, we define $\bar{\Pi}[a, b]: (k, \ell) \in \mathbb{N}^2 \to \bar{\Pi}[a, b]_{k, \ell} \in \mathbb{R}$ as

$$\bar{\Pi}[a, b]_{k, \ell} = \Pi[a_{0:k}, b]_{\ell} - \Pi[a_{0:k-1}, b]_{\ell}, \quad \forall (k, \ell) \in \mathbb{N}^2,$$

using again the convention $a_{0:-1} = 0$. Since the sequences $a$ and $b$ satisfy the assumptions of Proposition 2.1, for any $k \in \mathbb{N}$, $a_{0:k} \le b_{0:\infty}$. Thus for any $k$, the two terms $\Pi[a_{0:k}, b]$ and $\Pi[a_{0:k-1}, b]$ are well-defined.

**Lemma 2.2**

Let $a$ and $b$ satisfy the assumptions of Proposition 2.1. Then $\bar{\Pi}[a, b]$ lies in $\ell^1(\mathbb{N}^2, \mathbb{R}_+)$ and satisfies (2).

**Proof:** Let us define $\Lambda: (k, \ell) \in (\{-1\} \cup \mathbb{N}) \times \mathbb{N} \mapsto \Lambda_{k, \ell} \in \mathbb{R}$ as

$$\Lambda_{k, \ell} = \Pi[a_{0:k}, b]_{\ell}.$$

For convenience, we set $\pi = \bar{\Pi}[a, b]$. We have $\pi_{k, \ell} = \Lambda_{k, \ell} - \Lambda_{k-1, \ell}$ for all $(k, \ell) \in \mathbb{N}^2$. First we verify that $\pi \in \ell^1(\mathbb{N}^2, \mathbb{R}_+)$. As $a_{0:k} \ge a_{0:k-1}$, for any $k$, we have that $\Lambda_{k, \ell} \ge \Lambda_{k-1, \ell}$, by Lemma 2.1.5. Thus $\pi_{k, \ell} \ge 0$. Note that if $\pi$ satisfies (2) then $\pi$ has finite sum, so it only remains to verify that (2) holds true.

* Let $k \in \mathbb{N}$. Using Lemma 2.1.4, we have that for any $\ell \in \mathbb{N}$,

$$\pi_{k,0:\ell} = \Lambda_{k,0:\ell} - \Lambda_{k-1,0:\ell} \xrightarrow{\ell \to +\infty} a_{0:k} - a_{0:k-1} = a_k.$$

* For $k \in \mathbb{N}$,

$$\pi_{0:k, \ell} = \Lambda_{k, \ell} - \Lambda_{-1, \ell} = \Lambda_{k, \ell} = \Pi[a_{0:k}, b]_{\ell}.$$

Using Lemma 2.1.5, the function $x \mapsto \Pi[x, b]_{\ell}$ is 1-Lipschitz so it is continuous. Therefore,

$$\lim_{k \to +\infty} \pi_{0:k, \ell} = \Pi[a_{0:\infty}, b]_{\ell} = \Pi[b_{0:\infty}, b]_{\ell} = b_{\ell}.$$

The last equality follows from Lemma 2.1.3.

* Let $k > \ell$. We have

$$b_{0:\ell} \le a_{0:\ell} \le a_{0:k-1} \le a_{0:k}.$$

Therefore, by Lemma 2.1.2, we have $\Lambda_{k, \ell} = b_{\ell}$ and $\Lambda_{k-1, \ell} = b_{\ell}$ and finally $\pi_{k, \ell} = b_{\ell} - b_{\ell} = 0$.
## 2.3 Proof of the main result

We prove in this subsection Theorem 2.1. We assume that a transition-rate matrix $\tilde{Q}$ satisfying Assumptions A1 and A2 is given (its construction will be investigated in Section 2.4).

**Lemma 2.3**

Assume that there exists a transition-rate matrix $R = (R_{(i,\ell),(j,m)})_{(i,\ell),(j,m) \in S \times \mathbb{N}}$ satisfying the following properties:

$$ 
\begin{array}{ll}
\text{a) } R_{(i,\ell),(j,0:\infty)} = Q_{i,j}, & \forall i \in S, \forall \ell \in \mathbb{N}, \forall j \in S, \\
\text{b) } R_{(i,\ell),(0:\infty,m)} = \tilde{Q}_{\ell,m}, & \forall i \in S, \forall \ell \in \mathbb{N}, \forall m \in \mathbb{N}, \quad (3) \\
\text{c) } \text{cl}(i) \le \ell \text{ and } \text{cl}(j) > m \Rightarrow R_{(i,\ell),(j,m)} = 0, & \forall (i, \ell) \in S \times \mathbb{N}, \forall (j, m) \in S \times \mathbb{N}.
\end{array}
$$

Then for any $(X_0, Y_0) \in S \times \mathbb{N}$ such that $\text{cl}(X_0) \le Y_0$, the continuous-time Markov chain $Y$ with transition-rate matrix $\tilde{Q}$ starting in $Y_0$ is an upper-bounding process of $X$.

**Proof:** The transition-rate matrix allows us to define a coupling between the processes $X$ and $Y$.

The conditions (3.a) and (3.b) ensure that the marginals of this coupling are well-defined and have respectively the same law than the processes $X$ and $Y$.

Moreover the condition (3.c) ensures that the jumps to state $(j, m) \in S \times \mathbb{N}$ such that $\text{cl}(j) > m$ are not possible.

Thus since the initial condition verifies $\text{cl}(X_0) \le Y_0$, the order $\text{cl}(X) \le Y$ is preserved over time. ■

The proof of Theorem 2.1 amounts to finding a transition matrix $R$ satisfying (3). The search for a transition matrix $R$ will be successively simplified in the following lemmas, so as to boil down to an application of Proposition 2.1.

The first transformation of the problem that we realize simply consists of a translation of $R$, $Q$ and $\tilde{Q}$. For any pair $(i, \ell) \in S \times \mathbb{N}$, we define $M_{i,\ell}$ as

$$ M_{i,\ell} = -Q_{i,i} - \tilde{Q}_{\ell,\ell} \ge 0. $$

For any pair $(i, j) \in S^2$ and for any pair $(\ell, m) \in \mathbb{N}^2$, we define $Q'_{(i,\ell),j}$ and $\tilde{Q}'_{(i,\ell),m}$ as

$$ Q'_{(i,\ell),j} = Q_{i,j} + M_{i,\ell} \delta_{i=j}, \quad \forall (i,j) \in S^2, $$

$$ \tilde{Q}'_{(i,\ell),m} = \tilde{Q}_{\ell,m} + M_{i,\ell} \delta_{\ell=m}, \quad \forall (\ell,m) \in \mathbb{N}^2. $$

We use here the Kronecker symbol, i.e., $\delta_{i=j} = 1$ if $i = j$, otherwise $\delta_{i=j} = 0$. We directly see from the above definitions that $Q'_{(i,\ell),j} \ge 0$ and $\tilde{Q}'_{(i,\ell),m} \ge 0$.
**Lemma 2.4**

Let $R' = (R'_{(i,\ell),(j,m)}(i,\ell),(j,m))_{(i,\ell),(j,m) \in S \times \mathbb{N}}$ satisfy the following properties:

$$
\begin{align}
a) \quad & R'_{(i,\ell),(j,0:\infty)} = Q'_{(i,\ell),j}, && \forall i \in S, \forall \ell \in \mathbb{N}, \forall j \in S, \nonumber \\
b) \quad & R'_{(i,\ell),(0:\infty,m)} = \tilde{Q}'_{(i,\ell),m}, && \forall i \in S, \forall \ell \in \mathbb{N}, \forall m \in \mathbb{N}, \tag{4} \\
c) \quad & \text{cl}(i) \le \ell \text{ and } \text{cl}(j) > m \Rightarrow R'_{(i,\ell),(j,m)} = 0, && \forall (i, \ell) \in S \times \mathbb{N}, \forall (j, m) \in S \times \mathbb{N}, \nonumber \\
d) \quad & R'_{(i,\ell),(j,m)} \ge 0, && \forall (i, \ell) \in S \times \mathbb{N}, \forall (j, m) \in S \times \mathbb{N}. \nonumber
\end{align}
$$

$$
c) \quad \text{cl}(i) \le \ell \text{ and cl}(j) > m \Rightarrow R'_{(i,\ell),(j,m)} = 0, \quad \forall (i, \ell) \in S \times \mathbb{N}, \forall (j, m) \in S \times \mathbb{N},
$$

$$
d) \quad R'_{(i,\ell),(j,m)} \ge 0, \quad \forall (i, \ell) \in S \times \mathbb{N}, \forall (j, m) \in S \times \mathbb{N}.
$$

Then the matrix $R = (R_{(i,l),(j,m)}_{(i,l),(j,m)})_{(i,l),(j,m) \in S \times \mathbb{N}}$ defined by

$$
R_{(i,\ell),(j,m)} = R'_{(i,\ell),(j,m)} - M_{i,\ell} \delta_{(i,\ell)}(j,m),
$$

is a transition-rate matrix satisfying (3).

**Proof:** Let us verify that $R$ is a transition-rate matrix. The condition (4.d) and the definition of $R$ imply the nonnegativity of the non-diagonal terms of $R$. Moreover, we have

$$
\sum_{(j,m) \in S \times \mathbb{N}} R_{(i,\ell),(j,m)} = \sum_{(j,m) \in S \times \mathbb{N}} R'_{(i,\ell),(j,m)} - M_{i,\ell} = \sum_{j \in S} Q'_{(i,\ell),j} - M_{i,\ell} = \sum_{j \in S} Q_{i,j} + M_{i,\ell} - M_{i,\ell} = 0,
$$

and

$$
\sum_{\substack{(j,m) \in S \times \mathbb{N} \\ (j,m) \neq (i,\ell)}} R_{(i,\ell),(j,m)} = -R'_{(i,\ell),(i,\ell)} + M_{i,\ell} < \infty,
$$

which proves that $R$ is a transition-rate matrix.

Next, using condition (4.a) and the definition of $Q'$, we have that

$$
\sum_{m \in \mathbb{N}} R_{(i,\ell),(j,m)} = \sum_{m \in \mathbb{N}} R'_{(i,\ell),(j,m)} - M_{i,\ell} \delta_{i=j} = Q'_{(i,\ell),j} - M_{i,\ell} \delta_{i=j} = Q_{i,j},
$$

which implies condition (3.a). Condition (3.b) can be verified in a similar fashion. Finally, consider $i$, $j$, $\ell$, and $m$ such that $\text{cl}(i) \le \ell$ and $\text{cl}(j) > m$. Then $(i, \ell) \ne (j, m)$, which implies that $R_{(i,\ell),(j,m)} = R'_{(i,\ell),(j,m)} = 0$, which proves (3.c) and concludes the proof. ■

The following lemma shows that a matrix $R'$ satisfying (4) can be deduced from a matrix $R'' = (R''_{(i,\ell),(k,m)}(i,\ell))_{(i,\ell) \in S \times \mathbb{N}, (k,m) \in \mathbb{N}^2}$ that is of simpler form, in so far as the third index of $R''$ lies in $\mathbb{N}$ (while the third index of $R'$ lies in $S$). Note that in the statement of condition (5.a) below, we make use of the convention (1), so we have $Q'_{(i,\ell),k} = \sum_{j \in S_k} Q_{(i,\ell),j}$.

**Lemma 2.5**

Let $R'' = (R''_{(i,\ell),(k,m)}(i,\ell))_{(i,\ell) \in S \times \mathbb{N}, (k,m) \in \mathbb{N}^2}$ be such that

$$
\begin{align*}
a) \quad & R''_{(i,\ell),(k,0:\infty)} = Q'_{(i,\ell),k}, && \forall i \in S, \forall \ell \in \mathbb{N}, \forall k \in \mathbb{N}, \\
b) \quad & R''_{(i,\ell),(0:\infty,m)} = \tilde{Q}'_{(i,\ell),m}, && \forall i \in S, \forall \ell \in \mathbb{N}, \forall m \in \mathbb{N}, \tag{5} \\
c) \quad & \text{cl}(i) \le \ell \text{ and } k > m \Rightarrow R''_{(i,\ell),(k,m)} = 0, && \forall (i, \ell) \in S \times \mathbb{N}, \forall (k, m) \in \mathbb{N}^2, \\
d) \quad & R''_{(i,\ell),(k,m)} \ge 0, && \forall (i, \ell) \in S \times \mathbb{N}, \forall (k, m) \in \mathbb{N}^2.
\end{align*}
$$

$$
c) \quad \text{cl}(i) \le \ell \text{ and } k > m \Rightarrow R''_{(i,\ell),(k,m)} = 0, \quad \forall (i, \ell) \in S \times \mathbb{N}, \forall (k, m) \in \mathbb{N}^2,
$$

$$
d) \quad R''_{(i,\ell),(k,m)} \ge 0, \quad \forall (i, \ell) \in S \times \mathbb{N}, \forall (k, m) \in \mathbb{N}^2.
$$

Then, the matrix $R'$ defined for all $((i, \ell), (j, m)) \in (S \times \mathbb{N})^2$ by

$$
R'_{(i,\ell),(j,m)} = \begin{cases} \frac{Q'_{(i,\ell),j}}{Q'_{(i,\ell),\text{cl}(j)}} R''_{(i,\ell),(\text{cl}(j),m)} & \text{if } Q'_{(i,\ell),\text{cl}(j)} > 0, \\ 0 & \text{otherwise,} \end{cases}
$$

satisfies (4).
**Proof:** Let us verify (4.a). If $Q'_{(i,\ell),cl(j)} > 0$, we have with the definition of $R'$ and condition (5.a) that

$$ \sum_{m \in \mathbb{N}} R'_{(i,\ell),(j,m)} = \frac{Q'_{(i,\ell),j}}{Q'_{(i,\ell),cl(j)}} \sum_{m \in \mathbb{N}} R''_{(i,\ell),(cl(j),m)} = \frac{Q'_{i,j}}{Q'_{i,cl(j)}} Q'_{(i,\ell),cl(j)} = Q'_{(i,\ell),j}, $$

as was to be proved. Otherwise, if $Q'_{(i,\ell),cl(j)} = 0$, then necessarily $Q'_{(i,\ell),j} = 0$, and by definition, $R'_{(i,\ell),(j,m)} = 0$, which shows that (4.a) also holds in that case.

Let us verify (4.b). We have

$$ R'_{(i,\ell),(0:\infty,m)} = \sum_{k \in \mathbb{N}} \sum_{j \in S_k} R'_{(i,\ell),(j,m)}. $$

If $Q'_{(i,\ell),k} \neq 0$, we have

$$ \sum_{j \in S_k} R'_{(i,\ell),(j,m)} = \sum_{j \in S_k} \left( \frac{Q'_{(i,\ell),j}}{Q'_{(i,\ell),cl(j)}} R''_{(i,\ell),(cl(j),m)} \right) = \frac{R''_{(i,\ell),(k,m)}}{Q'_{(i,\ell),k}} \sum_{j \in S_k} Q'_{(i,\ell),j} = R''_{(i,\ell),(k,m)}. $$

If $Q'_{(i,\ell),k} = 0$, then $\sum_{j \in S_k} R'_{(i,\ell),(j,m)} = 0$ and $R''_{(i,\ell),(k,m)} = 0$ (by conditions a and d). Thus again, we have

$$ \sum_{j \in S_k} R'_{(i,\ell),(j,m)} = R''_{(i,\ell),(k,m)}. $$

Combining the three above equalities and (5.b), we deduce that

$$ R'_{(i,\ell),(0:\infty,m)} = \sum_{k \in \mathbb{N}} R''_{(i,\ell),(k,m)} = \tilde{Q}'_{(i,\ell),m}, $$

as was to be shown.

Let $i, \ell, j$, and $m$ be such that $cl(i) \le \ell$ and $cl(j) > m$. Then, by condition (5.c), $R''_{(i,\ell),(cl(j),m)} = 0$ and thus $R'_{(i,\ell),(j,m)} = 0$, which proves that (4.c) is satisfied. Condition (4.d) is also satisfied, as a direct consequence of (5.d).

## Lemma 2.6

Let Assumptions A1 and A2 hold true. Then there exists $R''$ satisfying (5).

**Proof:** The construction of $R''$ is realized for each pair $(i, \ell)$, separately from one another. So we fix $(i, \ell)$ and consider two cases.

* Case 1: $cl(i) \le \ell$. The construction of $R''_{(i,\ell),(\cdot,\cdot)}$ follows from Proposition 2.1, applied with the sequences $a := (Q'_{(i,\ell),j})_{j \in \mathbb{N}}$ and $b := (\tilde{Q}'_{(i,\ell),m})_{m \in \mathbb{N}}$. The result of Proposition 2.1 indeed implies all statements of (5). It remains to verify the applicability of the proposition. First we note that $a$ and $b$ are nonnegative sequences. Moreover,

$$ a_{0:\infty} = Q'_{(i,\ell),0:\infty} = Q_{i,0:\infty} + M_{i,\ell} = M_{i,\ell}, $$

and similarly, $b_{0:\infty} = M_{i,\ell}$; so $a, b \in \ell^1(\mathbb{N}, \mathbb{R}_+)$. It remains to verify that for any $k \in \mathbb{N}$, $a_{0:k} \ge b_{0:k}$. Let $k \in \mathbb{N}$. Note first that

$$ \begin{aligned} a_{0:k} &= Q'_{(i,\ell),0:k} = Q_{i,0:k} + M_{i,\ell} \delta_{k \ge cl(i)}, \\ b_{0:k} &= \tilde{Q}'_{(i,\ell),0:k} = \tilde{Q}_{\ell,0:k} + M_{i,\ell} \delta_{k \ge \ell}. \end{aligned} $$

Next we consider two cases.

- Case (a): $k \notin \{cl(i), \dots, \ell-1\}$. Using Assumption A1 once and Assumption A2 repeatedly, we obtain that

$$ Q_{i,0:k} \ge \tilde{Q}_{cl(i),0:k} \ge \tilde{Q}_{cl(i)+1,0:k} \ge \dots \ge \tilde{Q}_{\ell,0:k}. $$

Next, as $cl(i) \le \ell$, we have

$$ \delta_{k \ge cl(i)} \ge \delta_{k \ge \ell}. $$

Combining the two derived inequalities with the expressions of $a_{0:k}$ and $b_{0:k}$, we deduce that $a_{0:k} \ge b_{0:k}$ as was to be shown.

- Case (b): $k \in \{cl(i), \dots, \ell-1\}$. Note first that the above proof does not work anymore, as it would require us to utilize the inequality $\tilde{Q}_{k,0:k} \ge \tilde{Q}_{k+1,0:k}$, which we do not assume because it is wrong in general (indeed: $\tilde{Q}_{k,0:k} \le 0$ and $\tilde{Q}_{k+1,0:k} \ge 0$). So we proceed differently. We have the inequalities:

$$ Q_{i,0:k} \ge Q_{i,i} \quad \text{and} \quad \tilde{Q}_{\ell,0:k} \le -\tilde{Q}_{\ell,\ell}, $$
which are actually true for any $k$ in $\mathbb{N}$. As we consider here the case with $k \ge \text{cl}(i)$ and $k < \ell$, we have:

$$ 
\begin{aligned}
a_{0:k} &\ge Q_{i,i} + M_{i,\ell} \delta_{k \ge \text{cl}(i)} = Q_{i,i} + M_{i,\ell} = -\tilde{Q}_{\ell,\ell}, \\
b_{0:k} &\le -\tilde{Q}_{\ell,\ell} + M_{i,\ell} \delta_{k \ge \ell} = -\tilde{Q}_{\ell,\ell}.
\end{aligned}
$$

So $a_{0:k} \ge b_{0:k}$, as was to be shown.

* Case 2: $\text{cl}(i) > \ell$. This case is easier to treat, as condition (5.c) is trivially satisfied. Define

$$ R''_{(i,\ell)(k,m)} = \begin{cases} \frac{1}{M_{i,\ell}} Q'_{(i,\ell),k} \tilde{Q}'_{(i,\ell),m}, & \text{if } M_{i,\ell} > 0, \\ 0, & \text{otherwise.} \end{cases} $$

The nonnegativity of $R''_{(i,\ell)(k,m)}$ is clear in both cases. From the definition of $\tilde{Q}_{(i,\ell),j}$, we deduce that

$$ \tilde{Q}'_{(i,\ell),0:\infty} = \tilde{Q}_{\ell,0:\infty} + M_{i,\ell} \left( \sum_{m=0}^{\infty} \delta_{\ell=m} \right) = M_{i,\ell}. $$

If $M_{i,\ell} > 0$, then

$$ R''_{(i,\ell),(k,0:\infty)} = \frac{1}{M_{i,\ell}} Q'_{(i,\ell),k} \tilde{Q}'_{(i,\ell),0:\infty} = \frac{1}{M_{i,\ell}} Q'_{(i,\ell),k} M_{i,\ell} = Q'_{(i,\ell),k}, $$

so condition (5.a) holds true. If $M_{i,\ell} = 0$, then $Q_{i,i} = 0$ and $Q_{i,j} = 0$, for any $j$. This further implies that $Q'_{(i,\ell),k} = 0$, and so condition (5.a) also holds true in that case. Condition (5.b) can be verified with similar arguments. ■

**Proof of Theorem 2.1:** Lemma 2.6 gives the existence of a matrix $R''$ satisfying the assumptions of Lemma 2.5, which in turn gives the existence of a matrix $R'$ satisfying the assumptions of Lemma 2.4, which gives the existence of a matrix $R$ satisfying the assumptions of Lemma 2.3, which finally implies Theorem 2.1. ■

## 2.4 Construction of an upper-bounding transition-rate matrix

We address the issue of finding a transition-rate matrix $\tilde{Q}$ that satisfies Assumptions A1 and A2, for a given transition-rate matrix $Q \in Q(S)$, which we suppose is fixed for the whole section. Our analysis relies on a transformation of the transition-rate matrices $\tilde{Q} \in Q(\mathbb{N})$, which we define first. Consider the sets:

$$ \Omega^{-} = \{ (\ell, m) \in \mathbb{N}^{2} | m < \ell \} \quad \text{and} \quad \Omega^{+} = \{ (\ell, m) \in \mathbb{N}^{2} | m > \ell \}. $$

Consider the map $\Phi: \tilde{Q} \in Q(\mathbb{N}) \mapsto \Phi(\tilde{Q})$, where $\Phi(\tilde{Q}) = (U^{-}, U^{+})$ and where

$$ U^{-} = (\tilde{Q}_{\ell,0:m})_{(\ell,m) \in \Omega^{-}}, $$

$$ U^{+} = (\tilde{Q}_{\ell,m:\infty})_{(\ell,m) \in \Omega^{+}}. $$

Note that the function $\Phi$ takes as an argument a matrix in $Q(\mathbb{N})$ and returns a pair of double-indexed sequences of non-negative numbers (the first is indexed by $\Sigma^{-}$ and the second by $\Sigma^{+}$).

We introduce next the set $U$, defined as the set of pairs $(U^{-}, U^{+})$ where $U^{-} = (U_{\ell,m}^{-})_{(\ell,m) \in \Omega^{-}}$ and $U^{+} = (U_{\ell,m}^{+})_{(\ell,m) \in \Omega^{+}}$ satisfy the following:

$$ a) \quad U_{\ell,m-1}^{-} \le U_{\ell,m}^{-} \quad \forall (\ell, m) \in \Omega^{-} \text{ such that } m > 0, \qquad (6) $$

$$ b) \quad U_{\ell,m}^{+} \ge U_{\ell,m+1}^{+}, \quad \forall (\ell, m) \in \Omega^{+}, $$

and

$$ a) \quad U_{\ell,0}^{-} \ge 0, \quad \forall \ell \in \mathbb{N} \setminus \{0\}, \qquad (7) $$

$$ b) \quad \lim_{m \to +\infty} U_{\ell,m}^{+} = 0, \quad \forall \ell \in \mathbb{N}. $$
**Lemma 2.7**

The mapping $\Phi$ is a bijection from $\mathcal{Q}(\mathbb{N})$ to $\mathbf{U}$. Its inverse is given by

$$ \Phi^{-1}(U^{-}, U^{+})_{\ell,m} = \begin{cases} U_{\ell,m}^{-} - U_{\ell,m-1}^{-} & \text{if } (\ell, m) \in \Omega^{-}, \\ -U_{\ell,\ell-1}^{-} - U_{\ell,\ell+1}^{+} & \text{if } \ell = m, \\ U_{\ell,m}^{+} - U_{\ell,m+1}^{+} & \text{if } (\ell, m) \in \Omega^{+}. \end{cases} $$

with the convention for all $\ell \in \mathbb{N}$, $U_{\ell, -1}^{-} = 0$.

**Proof:** Let us call $\Psi$ the function from $\mathbf{U}$ to $\mathcal{Q}(\mathbb{N})$ defined in the lemma. As a first step, we will prove that the mappings $\Phi$ and $\Psi$ are well-defined.

Let $\tilde{Q} \in \mathcal{Q}(\mathbb{N})$ and $(U^{-}, U^{+}) = \Phi(\tilde{Q})$. We need to check that $(U^{-}, U^{+})$ is in $\mathbf{U}$, that is it verifies (6) and (7).

Since for any $\ell \neq m \in \mathbb{N}$, $\tilde{Q}_{\ell,m} \geq 0$, by definition of $U_{\ell,m}^{-}$ and $U_{\ell,m}^{+}$, the equation (6) is verified by $(U^{-}, U^{+})$.

Moreover, for all $\ell \geq 1$, $U_{\ell,0}^{-} = \tilde{Q}_{\ell,0:0} = \tilde{Q}_{\ell,0} \geq 0$. Finally, for any $\ell \in \mathbb{N}$, by definition of $\mathcal{Q}(\mathbb{N})$, the sum of the sequence $(\tilde{Q}_{\ell,m})_{m>\ell}$ is finite, so the limit of the rest of this sum goes to 0 and we have proved the condition (7)b.

Let $(U^{-}, U^{+}) \in \mathbf{U}$ and $\tilde{Q} = \Psi(U^{-}, U^{+})$. We need to check that $\tilde{Q} \in \mathcal{Q}(\mathbb{N})$.

With (6) and (7)a, the positivity of $\tilde{Q}_{\ell,m}$ for any $\ell \neq m$ is straightforward. Let us compute the sum according to the second index, we have on one hand for any $(\ell, m) \in \Omega^{-}$,

$$ \tilde{Q}_{\ell,0:\ell-1} = U_{\ell,0}^{-} + \sum_{k=1}^{\ell-1} (U_{\ell,k}^{-} - U_{\ell,k-1}^{-}) = U_{\ell,\ell-1}^{-}, $$

and on the other hand, with (7)b, for any $(\ell, m) \in \Omega^{+}$:

$$ \sum_{k=\ell+1}^{m} \tilde{Q}_{\ell,k} = \sum_{k=\ell+1}^{m} (U_{\ell,k}^{+} - U_{\ell,k+1}^{+}) = U_{\ell,\ell+1}^{+} - U_{\ell,m+1}^{+} \underset{m \to +\infty}{\longrightarrow} U_{\ell,\ell+1}^{+}. $$

$$ \text{So } \sum_{m \neq \ell} \tilde{Q}_{\ell,m} = U_{\ell,\ell-1}^{-} + U_{\ell,\ell+1}^{+} < \infty \text{ and } \sum_{m \in \mathbb{N}} \tilde{Q}_{\ell,m} = 0. $$

It only remains to prove that for any $\tilde{Q} \in \mathcal{Q}(\mathbb{N})$, $\Psi(\Phi(\tilde{Q})) = \tilde{Q}$. For convenience, we set $V_{\ell,m} = \Psi(\Phi(\tilde{Q}))_{\ell,m}$ for all $(\ell, m) \in \mathbb{N}^2$ and $\Phi(\tilde{Q}) = (U^{-}, U^{+})$.

Let $(\ell, m) \in \mathbb{N}^2$, three cases must be considered.

* If $(\ell, m) \in \Omega^{-}$, $V_{\ell,m} = U_{\ell,m}^{-} - U_{\ell,m-1}^{-} = \tilde{Q}_{\ell,0:m} - \tilde{Q}_{\ell,0:m-1} = \tilde{Q}_{\ell,m}$.

* If $m = \ell$, $V_{\ell,\ell} = -U_{\ell,\ell-1}^{-} - U_{\ell,\ell+1}^{+} = -\tilde{Q}_{\ell,0:\ell-1} - \tilde{Q}_{\ell,\ell+1:\infty} = \tilde{Q}_{\ell,\ell}$ where the last equality comes from the fact that the sum according to the second coordinate of $\tilde{Q}$ is equal to 0.

* If $(\ell, m) \in \Omega^{+}$, $V_{\ell,m} = U_{\ell,m}^{+} - U_{\ell,m+1}^{+} = \tilde{Q}_{\ell,m:\infty} - \tilde{Q}_{\ell,m+1:\infty} = \tilde{Q}_{\ell,m}$.

This concludes the proof.

<figure><img src="image_3.png" /></figure>

Recall that we have fixed a transition-rate matrix $Q$. Next we define $f^{-} = (f_{\ell,m}^{-})_{(\ell,m) \in \Omega^{-}}$ and $f^{+} = (f_{\ell,m}^{+})_{(\ell,m) \in \Omega^{+}}$ as

$$ f_{\ell,m}^{-} = \min_{i \in S_{\ell}} Q_{i,0:m} \quad \text{and} \quad f_{\ell,m}^{+} = \max_{i \in S_{\ell}} Q_{i,m:\infty}, \qquad (8) $$

which are well-defined since each class is finite. For simplicity, we do not mention explicitly the dependence of $f^{-}$ and $f^{+}$ with respect to $Q$ in the notation.

**Lemma 2.8**

Let $\tilde{Q} \in \mathcal{Q}(\mathbb{N})$ and let $(U^{-}, U^{+}) = \Phi(\tilde{Q})$. Then $\tilde{Q}$ satisfies Assumption A1 if and only if

$$ a) \quad U_{\ell,m}^{-} \leq f_{\ell,m}^{-}, \quad \forall (\ell, m) \in \Omega^{-}, \qquad (9) $$

$$ b) \quad U_{\ell,m}^{+} \geq f_{\ell,m}^{+}, \quad \forall (\ell, m) \in \Omega^{+}. $$
Moreover, $\tilde{Q}$ satisfies Assumption A2 if and only if

$$ a) \quad U_{\ell+1,m}^{-} \le U_{\ell,m}^{-} \quad \forall (\ell, m) \in \Omega^{-}, \qquad (10) $$

$$ b) \quad U_{\ell+1,m+1}^{+} \ge U_{\ell,m+1}^{+}, \quad \forall (\ell, m) \in \Omega^{+}. $$

**Proof:** The result directly follows from the case disjunction realized in Remarks 2.2 and 2.3. ■

A consequence of the three above lemmas is that the problem of finding a transition-rate matrix $\tilde{Q}$ satisfying Assumptions A1 and A2 is equivalent to finding $U^{-}$ and $U^{+}$ satisfying (6), (7), (9), and (10). Note that the conditions to be satisfied by $U^{-}$ are independent from those to be satisfied by $U^{+}$. We define next two candidates $\bar{U}^{-}$ and $\bar{U}^{+}$ as follows:

$$ \bar{U}_{\ell,m}^{-} = \min_{\ell' \in \{m+1, \dots, \ell\}} f_{\ell',m}^{-} \quad \forall (\ell, m) \in \Omega^{-}, \qquad (11) $$

$$ \bar{U}_{\ell,m}^{+} = \max_{\ell' \in \{0, \dots, \ell\}} f_{\ell',m}^{+} \quad \forall (\ell, m) \in \Omega^{+}. $$

Finally, we provide in the next proposition an explicit expression for a transition-rate matrix $\tilde{Q}$ satisfying Assumptions A1 and A2. Let us highlight the fact that the provided construction does not require any assumption on the given partition (except the finiteness of the classes).

**Proposition 2.2**

Firstly, the quantities $\bar{U}^{-}$ and $\bar{U}^{+}$ are well-defined and satisfy (6), (7), (9), and (10). Therefore $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})$ is a transition-rate matrix satisfying Assumptions A1 and A2.

Moreover, the transition-rate matrix $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})$ is optimal in the following sense: any other transition-rate matrix $\tilde{Q}$ satisfying Assumptions A1 and A2 is necessarily such that

$$ \tilde{Q}_{\ell,0:m} \le \Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m}, \quad \forall (\ell, m) \in \mathbb{N}^{2}. $$

**Proof:** Since $\bar{U}^{-}$ (respectively $\bar{U}^{+}$) is a minimum (resp. a maximum) on a finite non empty set, it is well-defined.

It is direct to see that $\bar{U}^{-}$ and $\bar{U}^{+}$ satisfy (6), (9) and (10). Condition (7.a) is straightforward since $f_{\ell,m}^{-}$ is positive for all $(\ell, m) \in \Omega^{-}$.

Finally, for any $\ell \in \mathbb{N}$, the sequence $(f_{\ell,m}^{+})_{m>\ell}$ is decreasing, so we have for all $(\ell, m) \in \Omega^{+}$,

$$ \bar{U}_{\ell,m}^{+} = \max_{\ell' \in \{0, \dots, \ell\}} f_{\ell',m}^{+} = \max_{\ell' \le \ell} \max_{i \in S_{\ell'}} Q_{i,m:\infty}. $$

Then since the sequence $(Q_{i,m:\infty})_{m}$ goes to 0 and the two above maximum are taken on finite sets, the condition (7.b) is fulfilled.

By Lemma 2.8, $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})$ is a transition-rate matrix satisfying Assumptions A1 and A2.

It remains to prove the optimality of the matrix $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})$. Let $\tilde{Q}$ be a transition-rate matrix satisfying Assumptions A1 and A2 and suppose that there exists $(\ell, m) \in \mathbb{N}^{2}$ such that

$$ \tilde{Q}_{\ell,0:m} > \Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m}. \qquad (12) $$

Two cases must be considered.

* If $m < \ell$, with Lemma 2.7, $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m} = \bar{U}_{\ell,m}^{-} = \min_{\ell' \in \{m+1, \dots, \ell\}} f_{\ell',m}^{-}$.

Thus there exists an integer $\ell_0 \in [m+1, \ell]$ such that $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m} = f_{\ell_0,m}^{-}$. Moreover, with Assumptions A1 and A2, and since $\ell_0 > m$ we have

$$ \tilde{Q}_{\ell,0:m} \le \tilde{Q}_{\ell_0,0:m} \le f_{\ell_0,m}^{-} = \Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m}. \qquad (13) $$

There is a contradiction between (12) and (13).

* If $m \ge \ell$, with Lemma 2.7, $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m} = -\bar{U}_{\ell,m+1}^{+} = - \max_{\ell' \in \{0, \dots, \ell\}} f_{\ell',m+1}^{+}$.

Thus there exists an integer $\ell_0 \in [0, \ell]$ such that $\Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m} = -f_{\ell_0,m+1}^{+}$. Moreover, with Assumptions A1 and A2, and since $\ell_0 \le \ell \le m$ we have

$$ \tilde{Q}_{\ell,0:m} \le \tilde{Q}_{\ell_0,0:m} = -\tilde{Q}_{\ell_0,m+1:\infty} \le -f_{\ell_0,m+1}^{+}. \qquad (14) $$

Again, there is a contradiction between (12) and (14).
So we have proven by contradiction that

$$ \tilde{Q}_{\ell,0:m} \le \Phi^{-1}(\bar{U}^{-}, \bar{U}^{+})_{\ell,0:m}, \quad \forall (\ell, m) \in \mathbb{N}^{2}. $$

This concludes the proof.

<figure><img src="image_4.png" /></figure>

## 2.5 Uniform bounds of solutions of the master equation

In this section, we derive from our main result (Theorem 2.1) a comparison on the probability distribution of processes in a family parameterized by a parameter $\theta$.

Let $(X^\theta)_{\theta \in \Theta}$ be a family of processes parameterized by a set $\Theta$ of suitable parameters. For any $\theta \in \Theta$, the master equation determines the probability distribution of $X^\theta$,

$$ \forall t \ge 0, \quad \frac{dp^\theta}{dt}(t) = p^\theta(t)Q^\theta, $$

where $p^\theta(t) = (p_i^\theta(t))_{i \in S}$ is a row vector such that $p_i^\theta(t) = \mathbb{P}(X_t^\theta = i | X_0^\theta)$.

In the following, the solution of the master equation linked to the upper-bounding process $Y$ (i.e. for the matrix $\tilde{Q}$) will be denoted by $\tilde{p}_t$. Then we can state an inequality for the probability distributions of the family $(X^\theta)_{\theta \in \Theta}$ and of $Y$.

**Theorem 2.2**

Assume that for any $\theta \in \Theta$, the pair $(Q^\theta, \tilde{Q})$ satisfies the Assumptions A1 and A2.

Assume further that

$$ \forall \ell \in \mathbb{N}, \quad \tilde{p}_{0:\ell}(0) \le \inf_{\theta \in \Theta} p_{0:\ell}^{\theta}(0). $$

Then for any $t \ge 0$ and $\ell \in \mathbb{N}$,

$$ \tilde{p}_{0:\ell}(t) \le \inf_{\theta \in \Theta} p_{0:\ell}^{\theta}(t). \quad (15) $$

The above inequality means that the process $\tilde{p}$ keeps less mass in the $\ell$ first classes than $p^\theta$ that is the mass go faster to infinity for the process $\tilde{p}$.

**Remarks:**

* If the set $\Theta$ is a singleton, Theorem 2.2 is a weaker result than Theorem 2.1 since it provides an inequality between the probability distributions $p$ and $\tilde{p}$.

However, Theorem 2.1 cannot be generalised to a family of processes $(X^\theta)_{\theta \in \Theta}$.

* Since $\tilde{p}$ and $p^\theta$ are probability distributions, the inequality (15) can be rewritten in terms of the distribution tails:

$$ \forall t \ge 0, \forall \ell \in \mathbb{N}, \inf_{\theta \in \Theta} p_{\ell:\infty}^{\theta}(t) \le \tilde{p}_{\ell:\infty}(t). $$

**Proof:** Let $\theta \in \Theta$, with Theorem 2.1, for any $t \ge 0$,

$$ \mathrm{cl}(X_t^\theta) \le Y_t \quad \text{a.s.} $$

So for any $\ell \in \mathbb{N}$, we have

$$ \{Y_t \le \ell\} \subseteq \{\mathrm{cl}(X_t^\theta) \le \ell\}. $$

By definition of the probability distributions,

$$ \forall \theta \in \Theta, \forall t \ge 0, \forall \ell \in \mathbb{N}, \tilde{p}_{0:\ell}(t) \le p_{0:\ell}^{\theta}(t). $$

By taking the infimum on the set $\Theta$, we obtain the result.

<figure><img src="image_5.png" /></figure>

## 3 Construction of a lower-bounding process

In this section, by adapting proofs and formulas, we give simple conditions on the transition-rate matrix of a lower-bounding process and then construct such a matrix.

For the sake of simplicity, we use the same notation than in the previous section on the upper-bounding process. Let us denote by $\tilde{Q}$ the transition-rate matrix of the lower-bounding process that we construct.
To be able to construct the lower-bounding process, the following hypothesis on the transition-rate matrix $Q$ is necessary:

$$ \forall m \in \mathbb{N}, \exists B_m \ge 0, \forall i \in S, \mathrm{cl}(i) > m \Rightarrow Q_{i,m} \le B_m. \qquad (16) $$

This hypothesis says that the rate of all reactions leading to a decrease of classes and ending in a class $m$ are uniformly bounded with respect to the starting class. It is directly verified when the number of reactions is finite.

## 3.1 Main result

### Theorem 3.1

Let $Q \in \mathcal{Q}(S)$ satisfying Equation (16) and let $\tilde{Q} \in \mathcal{Q}(\mathbb{N})$. Let $p^0$ be a distribution on $S$. Assume the two following properties:

(B1) $\forall \ell, m \in \mathbb{N}, \forall i \in S_{\ell}, \tilde{Q}_{\ell,0:m} \ge Q_{i,0:m},$

(B2) $\forall \ell, m \in \mathbb{N}$ such that $\ell \ne m, \tilde{Q}_{\ell+1,0:m} \le \tilde{Q}_{\ell,0:m}.$

Then there exist two continuous-time Markov chains $X = (X_t)_{t \ge 0}$ and $Y = (Y_t)_{t \ge 0}$, defined on $S$ and $\mathbb{N}$ respectively, with transition-rate matrices $Q$ and $\tilde{Q}$, such that $X_0$ is distributed according to $p^0$, and such that $Y$ is a lower-bounding process of $X$.

As in the previous section, we can give an interpretation of the two assumptions. Firstly note that the Assumption B2 is the same as Assumption A2. Conversely, the inequality is inverted between Assumptions A1 and B1.

### Remark 3.1

Assumption B1 essentially says that $Y$ decreases at a faster rate than $X$ moves from states to other states that belong to a lower class. Vice-versa, $Y$ increases at a smaller rate than $X$ moves from states to other states that belong to a higher class. This can be formalized as follows.

* When $m < \ell$, the assumption is equivalent to:

$$ \tilde{Q}_{\ell,0:m} \ge \min_{i \in S_{\ell}} Q_{i,0:m}. $$

The quantity $\tilde{Q}_{\ell,0:m}$ is the rate at which $Y$ moves from $\ell$ to classes that are smaller or equal to $m$ and the quantity $\min_{i \in S_{\ell}} Q_{i,0:m}$ is the minimal rate at which $X$ transitions from a state in class $\ell$ to another state in a class lower or equal to $m$.

* When $m \ge \ell$, the assumption is equivalent to:

$$ \tilde{Q}_{\ell,m+1:\infty} \le \max_{i \in S_{\ell}} Q_{i,m+1:\infty}. $$

The quantity $\tilde{Q}_{\ell,m+1:\infty}$ is the rate at which $Y$ moves from $\ell$ to classes that are strictly greater than $m$ and the quantity $\max_{i \in S_{\ell}} Q_{i,m+1:\infty}$ is the maximal rate at which $X$ transitions from a state in class $\ell$ to another state in a class strictly higher than $m$.

The interpretation of the Assumption B2 has already been stated in Remark 2.3 since this assumption is the same as Assumption A2.

## 3.2 Construction of a lower-bounding transition rate-matrix

By adapting the approach and formulas of Subsection 2.4, we find a transition-rate matrix $\tilde{Q}$ that satisfies Assumptions B1 and B2.

With the same notations as in Subsection 2.4, we define the two doubly-indexed sequences $\hat{f}^- = (\hat{f}_{\ell,m}^-)_{(\ell,m) \in \Omega^-}$ and $\hat{f}^+ = (\hat{f}_{\ell,m}^+)_{(\ell,m) \in \Omega^+}$ as

$$ \hat{f}_{\ell,m}^- = \max_{i \in S_{\ell}} Q_{i,0:m} \quad \text{and} \quad \hat{f}_{\ell,m}^+ = \min_{i \in S_{\ell}} Q_{i,m:\infty}. \qquad (17) $$
We define next two candidates $\hat{U}^-$ and $\hat{U}^+$ for the image of the matrix $\tilde{Q}$ by the function $\Phi$ (see Lemma 2.7) as follows:

$$ \hat{U}_{\ell,m}^{-} = \sup_{\ell' \in \{\ell, \ell+1, \dots\}} \hat{f}_{\ell',m}^{-} \quad \forall (\ell, m) \in \Omega^{-}, \qquad (18) $$

$$ \hat{U}_{\ell,m}^{+} = \min_{\ell' \in \{\ell, \dots, m-1\}} \hat{f}_{\ell',m}^{+} \quad \forall (\ell, m) \in \Omega^{+}. $$

**Proposition 3.1**

The quantities $\hat{U}^-$ and $\hat{U}^+$ are well-defined and are in the set $\mathbf{U}$.

Then $\Phi^{-1}(\hat{U}^-, \hat{U}^+)$ is a transition-rate matrix satisfying Assumptions B1 and B2.

Moreover, this transition-rate matrix $\Phi^{-1}(\hat{U}^-, \hat{U}^+)$ is optimal in the following sense: any other transition-rate matrix $\tilde{Q}$ satisfying Equation (16) and Assumptions B1 and B2 is necessarily such that

$$ \Phi^{-1}(\hat{U}^{-}, \hat{U}^{+})_{\ell,0:m} \leq \tilde{Q}_{\ell,0:m}, \quad \forall (\ell, m) \in \mathbb{N}^{2}. $$

**Proof:** Since $\hat{U}^+$ is a minimum on a finite set, it is well-defined. With Equation (16),

$$ \forall (\ell, m) \in \Omega^{-}, \quad \hat{f}_{\ell,m}^{-} \leq B_m, $$

where $B_m$ is independent of $\ell$. Thus $\hat{U}^-$ is well-defined.

The rest of the proof is an adaptation of the proof of Proposition 2.2. ■

<figure><img src="image_6.png" /></figure>

# 4 Applications

In this section, we will detail how different properties of the original process can be derived from the lower and upper-bounding processes. As an example, we will use the following chemical reaction network composed of three species $X_1$, $X_2$ and $X_3$

$$ \emptyset \xrightarrow{b_1(x_2+x_3)+b_2} X_1 \qquad 2X_1 \xrightarrow{\alpha x_1(x_1-1)} 3X_2 \qquad X_i \xrightarrow{d_i x_i} \emptyset \quad \text{for } 1 \le i \le 3 \qquad (19) $$
$$ X_1 + X_2 \xrightarrow{\beta x_1 x_2} X_3 $$

with all parameters ($b_1, b_2, \alpha, \beta, (d_i)_{1 \le i \le 3}$) strictly positive. In the following, we will call $\theta$ the parameter vector.

This network is composed of six reactions: one reaction creating the species $X_1$, two second-order reactions for the creation of species $X_2$ and $X_3$, and three reactions for the degradation of each species. The dynamics of each species depend on those of the other two due to the reaction rates chosen. In particular, the creation of the species $X_1$ is positively regulated by the sum of the number of molecules $X_2$ and $X_3$.

The motivation for choosing this reaction network as an example is that it poses significant challenges: the non-linearity of the reaction propensities and the dimensionality of the network make it difficult to analyze the dynamics and to establish properties of the CTMC. Notably, the creation of molecules of $X_2$ at a quadratic rate implies that molecule numbers could potentially diverge to infinity in finite time for some choices of model parameters.

Let $X$ be the continuous-time Markov chain on $\mathbb{N}^3$ associated to this network. In the following, we will apply our method with different sets of classes in order to derive results on the original process. For the upper-bounding processes, we will consider the two following sets of classes

$$ \forall \ell \in \mathbb{N}, \quad S_{\ell}^{+} = \{2x_1 + x_2 + x_3 = \ell\} \quad \text{and} \quad \hat{S}_{\ell}^{+} = \{x_1 + x_2 + x_3 = \ell\}. \qquad (20) $$

The partition $(\hat{S}_{\ell}^{+})_{\ell \in \mathbb{N}}$ is the naive choice of classes that groups states according to the total number of molecules. Contrary to that, the partition $(S_{\ell}^{+})_{\ell \in \mathbb{N}}$ has a linear form too but gives more weight to the molecules of species $X_1$. The main motivation for weighting more the species $X_1$ is to take into account
the structure of the reaction network by ensuring that the reaction $2X_1 \rightarrow 3X_2$ decreases classes. For the lower-bounding process, we will consider only one set of classes

$$ \forall \ell \in \mathbb{N}, \quad S_{\ell}^{-} = \{2x_1 + 2x_2 + 5x_3 = \ell\}. \qquad (21) $$

In the same spirit as for the two previous partitions, the weights in the partition $(S_{\ell}^{-})_{\ell \in \mathbb{N}}$ are chosen such that the two reactions $2X_1 \rightarrow 3X_2$ and $X_1 + X_2 \rightarrow X_3$ increase classes.

In a first step, we will illustrate the usefulness of our method to obtain conditions on the parameters for non-explosion and ergodicity of the multi-dimensional CTMC $X$. Then, in a second step, we will compare the optimal upper-bounding processes obtained with the two partitions $(S_{\ell}^{+})_{\ell \in \mathbb{N}}$ and $(\hat{S}_{\ell}^{+})_{\ell \in \mathbb{N}}$ introduced before which then allows us to discuss the importance of the choice of the classes for the coupling and for the study of the stationary distribution. Finally, we will bound the error made in numerically solving the master equation on a truncated state space.

# 4.1 Asymptotic behavior

## 4.1.1 Theoretical results

Let us recall some links between the asymptotic behavior of a CTMC $X$ and hitting times. Let $X$ be an irreducible continuous-time Markov chain on a state space $S$.

The sequence $J = (J_n)_{n \in \mathbb{N}}$ of jump times of $X$ are defined by

$$ J_0 = 0 \quad \text{and for any } n \ge 1, \quad J_n = \inf\{t \ge J_{n-1}, X_t \ne X_{J_{n-1}}\}, $$

where $\inf \emptyset = \infty$ by convention.

The process $X$ is explosive if and only if $\mathbb{P}(\sup_n J_n < \infty) > 0$.

For any state $i \in S$, the first hitting time of $i$ is the first time at which the chain reaches state $i$ (excluding time 0), that is

$$ T_i^X = \inf(t > 0, X_t = i), $$

with the convention that $T_i^X = \infty$ if the chain never reaches state $i$.

The asymptotic behavior of the irreducible continuous-time Markov chain $X$ can be described with the hitting times:

1. $X$ is transient $\iff$ there exists $i \in S$ such that $\mathbb{P}(T_i^X = \infty) > 0$.

2. $X$ is recurrent $\iff$ there exists $i \in S$ such that $\mathbb{P}(T_i^X = \infty) = 0$.

We distinguish between two types of recurrent processes:

* $X$ is null-recurrent $\iff$ there exists $i \in S$ such that $\mathbb{E}(T_i^X) = \infty$.

* $X$ is positive recurrent $\iff$ there exists $i \in S$ such that $\mathbb{E}(T_i^X) < \infty$.

Since $X$ is irreducible, the existence of a stationary distribution is equivalent to positive recurrence. Note that only transient process can be explosive, so a non-explosive process can be transient, null-recurrent or positive recurrent.

In the following, the coupling presented in the previous sections is used in order to determine the asymptotic behavior of $X$. Let $Z$ be a lower-bounding process for classes $(S_{\ell}^{-})_{\ell \in \mathbb{N}}$ and $Y$ a upper-bounding process for classes $(S_{\ell}^{+})_{\ell \in \mathbb{N}}$ (see Sections 2 and 3), that is for any $t \ge 0$,

$$ \left\{ \begin{array}{ll} Z_t & \le \text{cl}^{+}(X_t), \quad \text{for } (S_{\ell}^{-})_{\ell \in \mathbb{N}}, \\ \\ \text{cl}^{-}(X_t) & \le Y_t, \quad \text{for } (S_{\ell}^{+})_{\ell \in \mathbb{N}} \end{array} \right. $$

where the functions $\text{cl}^{+}$ (resp. $\text{cl}^{-}$) is the function defined in Definition 1.1 for classes $(S_{\ell}^{+})_{\ell \in \mathbb{N}}$ (resp. for $(S_{\ell}^{-})_{\ell \in \mathbb{N}}$).
**Theorem 4.1**

Suppose that processes Y and Z are irreducible. Then according to the asymptotic behavior of the lower and upper-bounding processes Z and Y, we can deduce properties on the asymptotic behavior of the process X presented in the Table 1:

<table><caption>Table 1: Asymptotic behavior of X according to that of the lower and upper-bounding processes Z and Y. Grey regions represent impossible combinations. In the blue regions, the asymptotic behavior of X can be fully deduced from the bounding processes. In the white regions, we only have access to partial information.</caption><thead><tr><th colspan="2">Upper process Y<sub>t</sub></th><th colspan="2">Transient</th><th colspan="2">Recurrent</th></tr><tr><th colspan="2">Lower process Z<sub>t</sub></th><th>Explosive</th><th>Non-explosive</th><th>Null-recurrent</th><th>Positive recurrent</th></tr></thead><tbody><tr><th rowspan="2">Transient</th><th>Explosive</th><td>Explosive</td><td></td><td></td><td></td></tr><tr><th>Non-explosive</th><td>Transient (Explosive or not)</td><td>Transient and non-explosive</td><td></td><td></td></tr><tr><th rowspan="2">Recurrent</th><th>Null-recurrent</th><td>Transient or null-recurrent</td><td>Non-explosive</td><td>Null-recurrent</td><td></td></tr><tr><th>Positive recurrent</th><td>No information</td><td>Non-explosive</td><td>Recurrent</td><td>Positive recurrent</td></tr></tbody></table>

**Proof:** The results on the asymptotic behavior of the process X presented in the Table 1 are equivalent to the following set of implications:

$$ \left\{ \begin{array}{lclclcl} Y & \text{non-explosive} & \Rightarrow & X & \text{non-explosive}, & & Z & \text{explosive} & \Rightarrow & X & \text{explosive}, \\ Y & \text{recurrent} & \Rightarrow & X & \text{recurrent}, & & Z & \text{transient} & \Rightarrow & X & \text{transient}, \\ Y & \text{positive recurrent} & \Rightarrow & X & \text{positive recurrent}. & & Z & \text{null-recurrent} & \Rightarrow & X & \text{non positive recurrent}. \end{array} \right. $$

Since Y is an upper-bounding process of X, for any $i \in S_0^+$, $\{T_i^X = \infty\} \Rightarrow \{T_0^Y = \infty\}$, so

$$ \mathbb{P}(T_i^X = \infty) \leq \mathbb{P}(T_0^Y = \infty). \qquad (22) $$

The same inequality holds for the expectation,

$$ \begin{aligned} \mathbb{E}(T_i^X) &= \int_0^\infty \mathbb{P}(T_i^X \ge t) dt \\ &\le \int_0^\infty \mathbb{P}(T_0^Y \ge t) dt \\ &\le \mathbb{E}(T_0^Y). \end{aligned} \qquad (23) $$

With the same calculations, we obtain corresponding inequalities for the lower-bounding process Z:

$$ \forall i' \in S_0^-, \quad \left\{ \begin{array}{lclcl} \mathbb{P}(T_0^{Z} = \infty) & \le & \mathbb{P}(T_{i'}^{X} = \infty) & & \\ \mathbb{E}(T_0^{Z}) & \le & \mathbb{E}(T_{i'}^{X}) & & \end{array} \right. . \qquad (24) $$

Then the implications on the left follow from Equations (22) and (23), since we suppose that processes X and Y are irreducible. The implications on the right follow from Equation (24) since we suppose that processes X and Z are irreducible. ■

**Remark 4.1**

When the bounding processes are constructed specifically to study the asymptotic behavior of the initial Markov chain X (by applying Theorem 4.1), the relevant property is the asymptotic behavior of the bounding processes. Then, it is often much simpler to compute the quantities $(\bar{U}_{\ell,m})_{(\ell,m) \in \Omega^-}$
and $(\bar{U}_{\ell,m}^{+})_{(\ell,m)\in\Omega^+}$ (see Equation (11)), $(\hat{U}_{\ell,m}^{-})_{(\ell,m)\in\Omega^-}$ and $(\hat{U}_{\ell,m}^{+})_{(\ell,m)\in\Omega^+}$ (see Equation (18)) only for large $\ell$, that is for any $\ell$ bigger than a constant $L \in \mathbb{N}$.

### 4.1.2 Application of the method to the toy model

In this subsection, we compute respectively the transition-rate matrices $\tilde{Q}^-$ and $\tilde{Q}^+$ of the optimal lower and upper-bounding processes Z and Y for two different sets of classes and then deduce results on the asymptotic behavior of the process X.

Let us start with the lower-bounding process Z and recall the set of classes defined by Equation (21),

$$ \forall \ell \in \mathbb{N}, \quad S_{\ell}^{-} = \{x \in \mathbb{N}^3, 2x_1 + 2x_2 + 5x_3 = \ell\}. $$

Using the results of Section 3, the transition-rate matrix $Q^-$ of the lower-bounding process Z can be expressed as

$$ \tilde{Q}^{-} = \Phi(\hat{U}^{-}, \hat{U}^{+}), $$

where $(\hat{U}_{\ell,m}^{-})_{(\ell,m)\in\Omega^-}$ and $(\hat{U}_{\ell,m}^{+})_{(\ell,m)\in\Omega^+}$ are defined by Equation (18). Since we focus on the long-term behavior of the process Z and to keep formulas simple, the only values that matter are $(\hat{U}_{\ell,m}^{-})_{(\ell,m)\in\Omega^-}$ and $(\hat{U}_{\ell,m}^{+})_{(\ell,m)\in\Omega^+}$ for $\ell$ large enough (see Remark 4.1).

**Proposition 4.1**

Consider the reaction network defined in Equation (19) and the partition defined in Equation (21). There exists $L \in \mathbb{N}$ such that for any $\ell \ge L$, the quantities $\hat{U}^-$ and $\hat{U}^+$ are explicitly given by

$$ \hat{U}_{\ell,m}^{-} = \begin{cases} 0 & \text{if } m < \ell - 5, \\ \frac{d_3}{5}(m+5) & \text{if } \ell - 5 \le m \le l - 3, \\ B(\theta)(m+2) & \text{if } m \in \{\ell - 2, \ell - 1\}. \end{cases} $$

$$ \hat{U}_{\ell,m}^{+} = \begin{cases} \frac{b_1}{5}\ell + A(\theta) & \text{if } \ell + 1 \le m \le \ell + 2, \\ 0 & \text{if } \ell + 3 \le m. \end{cases} $$

$$ \text{with } A(\theta) = b_2 - \frac{(\alpha + \frac{2}{5}b_1)^2}{4\alpha} \quad \text{and} \quad B(\theta) = \max\left(\frac{d_1}{2}, \frac{d_2}{2}, \frac{d_3}{5}\right). $$

**Proof:** The proof of Proposition 4.1 is detailed in Subsection A.1. ■

Note that the previous proposition gives the expression of the lower-bounding process only for large integers $\ell$ corresponding to large classes of the initial process. The computation of the exact expression of the transition rates is possible for any integer $\ell$ but leads to many case disjunctions depending on the parameter values. Since transition rates for small classes are not relevant for the study of the asymptotic behavior of the process, we refrain from listing them here.

Now, we construct the transition-rate matrix of the upper-bounding process Y. Let us recall the following set of classes defined by Equation (20),

$$ \forall \ell \in \mathbb{N}, \quad S_{\ell}^{+} = \{2x_1 + x_2 + x_3 = \ell\}. $$

**Proposition 4.2**

Consider the reaction network defined in Equation (19) and the partition $(S_{\ell}^{+})_{\ell \in \mathbb{N}}$ defined in Equation (20). There exists $L \in \mathbb{N}$ such that for any $\ell \ge L$, the matrix $\tilde{Q}^+$ of the optimal upper-bounding process Y for the set of classes $(S_{\ell}^{+})_{\ell \in \mathbb{N}}$ verifies,

$$ \begin{align*} \tilde{Q}_{\ell, \ell-1}^{+} &= \min(d_2, d_3)\ell - C(\theta), & \tilde{Q}_{\ell, \ell}^{+} &= -\tilde{Q}_{\ell, \ell-1}^{+} - \tilde{Q}_{\ell, \ell+2}^{+}, \\ \tilde{Q}_{\ell, \ell+2}^{+} &= b_1\ell + b_2, & \tilde{Q}_{\ell, m}^{+} &= 0 \quad \text{if } m \notin \{\ell-1, \ell+2\}, \end{align*} $$

$$ \text{where } C(\theta) = \frac{(\alpha + 2 \min(d_2, d_3) - d_1)^2}{4\alpha} \mathbf{1}_{\alpha + 2 \min(d_2, d_3) - d_1 \ge 0}(\theta). $$
**Proof:** The proof of Proposition 4.2 is detailed in Subsection A.2.

<figure><img src="image_7.png" /></figure>

In order to be able to apply Theorem 4.1 to obtain results on the asymptotic behavior of the process X, it only remains to study the asymptotic behavior of the two one-dimensional CTMCs Y and Z. Based on [36], we defined the following quantities:

• for the process Y,

$$B'_{1} = 2b_{1} - \min(d_{2}, d_{3}) \quad \text{and} \quad B'_{3} = 2b_{2} + C(\theta) - \frac{3}{2} \min(d_{2}, d_{3}),$$

• and for the process Z,

$$B_{1} := \lim_{\ell \to \infty} \frac{\hat{U}_{\ell, \ell+1:\infty}^{+} - \hat{U}_{\ell, 0:\ell-1}^{-}}{\ell} = \frac{2}{5}b_{1} - \frac{3}{5}d_{3} - 2B(\theta),$$

$$B_2 := \lim_{\ell \to \infty} (\hat{U}_{\ell, \ell+1:\infty}^{+} - \hat{U}_{\ell, 0:\ell-1}^{-} - B_1 \ell) = 2A(\theta) - \frac{3}{5}d_3 - B(\theta),$$

$$B_3 := B_2 - \frac{1}{2} \lim_{\ell \to \infty} \frac{\sum_{k=1}^{\infty} (2k-1) \hat{U}_{\ell, \ell+k}^{+} + \sum_{k=1}^{\ell} (2k-1) \hat{U}_{\ell, \ell-k}^{-}}{\ell} = 2A(\theta) - \frac{2}{5}b_1 - \frac{27}{10}d_3 - 3B(\theta).$$

Parameters $B_1$ and $B'_1$ represent the average jump size of the processes. $B_2$ is the average jump size of the process if the linear contributions are removed. The expression of the average jump size of the process Y ($B'_2$) is not detailed here since this parameter is always strictly positive, meaning that it has no impact on the asymptotic behavior of the process. The sign of these five parameters determine the asymptotic behavior of the processes Y and Z.

<table><caption>Table 2: Asymptotic behavior of X according to the parameter values. The grey regions represent impossible parameter combinations. In the blue regions, the asymptotic behavior of X can be fully deduced from the bounding processes. In the white regions, we only have access to partial information.</caption><thead><tr><th></th><th>Upper process Y<sub>t</sub></th><th>Transient and non-explosive</th><th>Null-recurrent</th><th>Recurrent</th></tr></thead><tbody><tr><th>Lower process Z<sub>t</sub></th><td></td><td>B'<sub>1</sub> > 0 or { B'<sub>1</sub> = 0<br>B'<sub>3</sub> > 0</td><td>{ B'<sub>1</sub> = 0<br>B'<sub>3</sub> ≤ 0</td><td>B'<sub>1</sub> < 0</td></tr><tr><th>Transient</th><td>B<sub>1</sub> > 0 or { B<sub>1</sub> = 0<br>B<sub>3</sub> > 0</td><td>Transient and non-explosive</td><td></td><td></td></tr><tr><th>Recurrent</th><td>{ B<sub>1</sub> = 0<br>B<sub>2</sub> ≥ 0<br>B<sub>3</sub> ≤ 0</td><td>Non-explosive</td><td>Null-recurrent</td><td></td></tr><tr><th></th><td>{ B<sub>1</sub> = 0<br>B<sub>2</sub> < 0</td><td>or B<sub>1</sub> < 0</td><td>Non-explosive</td><td>Recurrent</td><td>Positive recurrent</td></tr></tbody></table>

Table 2 illustrates the kind of condition that we can obtain with the construction of the lower and upper-bounding processes. In this case, the processes Y and Z are non-explosive for all parameter values so the corresponding column and row in Table 1 have been removed.

## 4.2 Importance of the choice of classes

Let us recall the following partition defined by Equation (20),

$$\forall \ell \in \mathbb{N}, \hat{S}_{\ell}^{+} = \{x_1 + x_2 + x_3 = \ell\}.$$

In this partition $\hat{S}^+$, the states are ordered by the total number of molecules. This partition is the simplest partition but does not take into account the structure of the reaction network.
**Proposition 4.3**

Consider the reaction network defined in Equation (19) and the partition $(\hat{S}_{\ell}^{+})_{\ell \in \mathbb{N}}$ defined in Equation (20). Then the matrix $\tilde{Q}^{+}$ of the optimal upper-bounding process $\hat{Y}$ for the set of classes $(\hat{S}_{\ell}^{+})_{\ell \in \mathbb{N}}$ verifies, for any $\ell \ge 0$,

$$ \tilde{Q}_{\ell,\ell-1}^{+} = d_2 \ell, \quad \tilde{Q}_{\ell,\ell+1}^{+} = \alpha \ell^2 - \alpha \ell + b_2, \\ \tilde{Q}_{\ell,\ell}^{+} = -\tilde{Q}_{\ell,\ell-1}^{+} - \tilde{Q}_{\ell,\ell+1}^{+}. $$

For any set of parameters $\theta$, the process $\hat{Y}$ is explosive.

**Proof:** The calculations of the matrix $\tilde{Q}^{+}$ is detailed in Subsection A.3.

With the expressions of the birth and death rates $\tilde{Q}_{\ell,l+1}^{+}$ and $\tilde{Q}_{\ell,l-1}^{+}$ and Theorems (3.1), (3.3) and (3.7) of [36], we obtain directly the same the process $\hat{Y}$ is explosive. ■

One can see that using the classes $\hat{S}^{+}$ leads to a different optimal upper-bounding process compared to the process in Proposition 4.2 that was obtained with the classes $S^{+}$. Concretely, we obtain a process with a quadratic birth rate and linear death rate as opposed to the linear birth and death rates in Proposition 4.2.

According to Table 1, when the upper-bounding process is explosive (first column of Table 1), the original process $X$ can have any asymptotic behavior: explosive, transient, null-recurrent or positive recurrent. Thus the naive choice of classes gives no information about the asymptotic behavior of the process $X$. On the contrary, the upper-bounding process $Y$ used in Subsection 4.1.2 gives relevant conditions on the asymptotic behavior of the process $X$ (see Table 2). Thus, the choice of the set of classes influences the quality of the results obtained on the stationary distribution.

### 4.2.1 Coupling with two sets of classes

In this subsection, numerical simulations are used to illustrate the coupling used in the proof of Theorem 2.1 and the difference between the two upper-bounding processes $Y$ and $\hat{Y}$ (see Propositions 4.2 and 4.3) related respectively to the partition $S^{+}$ and to the naive partition $\hat{S}^{+}$.

For the numerical simulations, we used the set of parameters $\theta$ detailed in Table 3 for which the initial process $X$ and the upper-bounding process $Y$ have a stationary distribution (see Table 2) as opposed to the explosive process $\hat{Y}$.

<table><caption>Table 3: Parameter values used for the numerical illustrations.</caption><thead><tr><th></th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>α</th><th>β</th><th>d<sub>1</sub></th><th>d<sub>2</sub></th><th>d<sub>3</sub></th></tr></thead><tbody><tr><th>θ</th><td>1</td><td>5/2</td><td>5/2</td><td>2</td><td>5/2</td><td>5/2</td><td>3</td></tr></tbody></table>

The initial condition is fixed to $X_0 = (15, 5, 5)$ that corresponds to $Y_0 = 40$ and $\hat{Y}_0 = 25$.

Figure 1 illustrates the difference in behavior at infinity between the two upper-bounding processes: the existence of a stationary distribution for $Y$ compared to the explosion of process $\hat{Y}$. The three trajectories of $\hat{Y}$ gradually move away from the initial process $X$ (see right-hand panel) and can only be stochastically simulated for very small time ($T_f = 0.013$ here) due to the explosion of the process. This is due to the quadratic birth rate of $\hat{Y}$ (see Proposition 4.3). On the contrary, the trajectories of $Y$ stabilise and remain close to $X$ (see left-hand panel).

## 4.3 Truncation error

### 4.3.1 Theoretical result

In principle, the time evolution of the probability distribution of a CTMC can be computed by solving the corresponding master equation. However, for most reaction networks the state space is countably infinite and solving the master equation is impossible without first truncating the state space [32]. Truncating the state space gives rise to a linear system of ordinary differential equations whose solution approximates the solution of the full master equation. In the case of a multi-dimensional CTMC, even after truncation
<figure><img src="image_8.png" /><figcaption>Figure 1: Three trajectories of the coupling process for the two sets of classes $S^+$ (left panel) and $\hat{S}^+$ (right panel). On each panel, we plot the trajectory of the upper-bounding process (second component of the coupling obtained with Theorem 2.1) in dashed lines and the class of the original process in solid lines with respect to time.</figcaption></figure>

the state space is often still large leading to a time-consuming numerical solving of this linear system. It is therefore impossible to test different truncations when we want to guarantee a given approximation error due to the truncation.

In this section, we show how the bounding processes can be used to guide the choice of truncation and to quantify the approximation error. Since the bounding processes are one-dimensional, they are easy to analyse and simulate so our approach allows us to deal with high-dimensional initial CTMC without any numerical problems.

Let $\Omega$ be the truncated state space and $T_f$ the final time. The truncation error is defined as the probability that the process of interest $X$ will exit $\Omega$ over the whole time interval $[0, T_f]$,

$$e_{trunc}(\Omega) := \mathbb{P}(\exists t \le T_f, X_t \notin \Omega).$$

The goal of this section is to use the upper-bounding process constructed in Section 2 to bound this truncation error. For the sake of simplicity, we will only consider truncated state spaces $\Omega$ defined through the classes, i.e for any $N \in \mathbb{N}$,

$$\Omega_N = \bigcup_{\ell=0}^{N} S_{\ell}.$$

With Theorem 2.1, we can construct a one-dimensional continuous-time Markov chain $Y_t$ on $\mathbb{N}$ such that

$$\forall t \ge 0, \text{cl}(X_t) \le Y_t.$$

We have a first obvious bound given by

$$e_{trunc}(\Omega_N) \le \mathbb{P}(\exists t \le T_f, Y_t > N).$$

It is possible to compute an upper bound of $\mathbb{P}(\exists t \le T_f, Y_t > N)$ in function of the solution to the CME associated with $Y$. In practice, only a truncated version $\tilde{p}^M$ of the CME can be computed. For a given $M \in \mathbb{N}$ with $M \ge N$, we define it as follows:

$$\forall t \ge 0, \forall m \in [0, M], \quad \frac{d\tilde{p}_m^M}{dt}(t) = \sum_{\ell=0}^{M} \tilde{p}_\ell^M(t) \tilde{Q}_{\ell,m}, \quad \tilde{p}_m^M(0) = \mathbb{P}[X_0 \in S_m].$$
Finally, we introduce the two following quantities,

• the exit rate $F_{\Omega_N}^{\text{out}}(t)$ at time $t$ of the set $\Omega_N$,

$$F_{\Omega_N}^{\text{out}}(t) = \sum_{\ell=0}^{N} \sum_{m=N+1}^{M} \tilde{p}_{\ell}^{M}(t) \tilde{Q}_{\ell,m},$$

• and the integral of $F_{\Omega_N}^{\text{out}}$ on $[0, T_f]$,

$$F_{\Omega_N} = \int_0^{T_f} F_{\Omega_N}^{\text{out}}(s) \, ds.$$

Our next theorem gives an upper bound of $e_{\text{trunc}}(\Omega_N)$ in function of $\tilde{p}^M$.

**Theorem 4.2**

Let $(M, N) \in \mathbb{N}^2$ such that $M \ge N$. Then the truncated error is upper bounded by

$$e_{\text{trunc}}(\Omega_N) \le \underbrace{(1 - \tilde{p}_{0:M}^M(T_f)) + \tilde{p}_{N+1:M}^M(0)}_{=:C(N)} + F_{\Omega_N}. \quad (25)$$

So for a fixed error $\epsilon$, the minimal state space that guarantees a truncation error smaller than $\epsilon$ is $\Omega_{\tilde{N}}$ with

$$\hat{N} = \min E(\epsilon), \quad (26)$$

where $E(\epsilon) = \{N, 0 \le N \le M \text{ s.t } C(N) + F_{\Omega_N} < \epsilon\}$.

**Proof:** Let us call $\tilde{p}^N$ the solution to the CME associated with $Y$ on the truncated state space $[0, N]$, i.e.,

$$\forall t \ge 0, \forall m \in [0, N], \quad \frac{d\tilde{p}_m^N}{dt}(t) = \sum_{\ell=0}^{N} \tilde{p}_\ell^N(t) \tilde{Q}_{\ell,m}, \quad \tilde{p}_m^N(0) = \mathbb{P}[X_0 \in S_m].$$

We have a probabilistic interpretation of $\tilde{p}_m^N$, given by

$$\tilde{p}_m^N(t) = \mathbb{P}[Y_t = m \text{ and } Y_s \le N, \forall s \in [0, t]]. \quad (27)$$

Therefore we have

$$e_{\text{trunc}}(\Omega_N) \le 1 - \tilde{p}_{0:N}^N(T_f).$$

We need now to find a lower bound of the quantity $\tilde{p}_{0:N}^N(T_f)$. We have

$$\tilde{p}_{0:N}^N(T_f) - \tilde{p}_{0:M}^M(T_f) = \tilde{p}_{0:N}^N(0) - \tilde{p}_{0:M}^M(0) + \int_0^{T_f} \frac{d}{dt} (\tilde{p}_{0:N}^N(t) - \tilde{p}_{0:M}^M(t)) dt \\ = - \tilde{p}_{N+1:M}^M(0) + \int_0^{T_f} \left[ \underbrace{\sum_{\ell=0}^N (\tilde{p}_\ell^N(t) \sum_{m=0}^N \tilde{Q}_{\ell,m})}_{=:\phi(t)} - \sum_{\ell=0}^M (\tilde{p}_\ell^M(t) \sum_{m=0}^M \tilde{Q}_{\ell,m}) \right] dt. \quad (28)$$

We have the following inequality, for any $m \in [0, N]$:

$$\tilde{p}_m^N(t) = \mathbb{P}[Y_t = m \text{ and } Y_s \le N, \forall s \in [0, t]] \le \mathbb{P}[Y_t = m \text{ and } Y_s \le M, \forall s \in [0, t]] = \tilde{p}_m^M(t).$$

As $\sum_{m=0}^N \tilde{Q}_{\ell,m} \le 0$, we deduce that

$$\phi(t) \ge \sum_{\ell=0}^{N} (\tilde{p}_{\ell}^{M}(t) \sum_{m=0}^{N} \tilde{Q}_{\ell,m}) - \sum_{\ell=0}^{M} (\tilde{p}_{\ell}^{M}(t) \sum_{m=0}^{M} \tilde{Q}_{\ell,m}).$$

Re-arranging the right-hand side, we obtain that $\phi(t) \ge \phi_1(t) + \phi_2(t)$, where

$$\phi_1(t) = \sum_{\ell=0}^{N} \left[ \tilde{p}_{\ell}^{M}(t) \left( \sum_{m=0}^{N} \tilde{Q}_{\ell,m} - \sum_{m=0}^{M} \tilde{Q}_{\ell,m} \right) \right] \quad \text{and} \quad \phi_2(t) = \sum_{\ell=N+1}^{M} \left( \tilde{p}_{\ell}^{M}(t) \sum_{m=0}^{M} \tilde{Q}_{\ell,m} \right).$$

Obviously, we have

$$\phi_1(t) = - \sum_{\ell=0}^{N} \tilde{p}_{\ell}^{M}(t) \left( \sum_{m=N+1}^{M} \tilde{Q}_{\ell,m} \right) \quad \text{and} \quad \phi_2(t) \le 0. \quad (29)$$
Combining (28) with (29), we obtain that

$$ \tilde{p}_{0:N}^{N}(T_f) - \tilde{p}_{0:M}^{M}(T_f) = -\tilde{p}_{N+1:M}^{M}(0) - \sum_{\ell=0}^{N} \sum_{m=N+1}^{M} \int_{0}^{T_f} \tilde{p}_{\ell}^{M}(t) \tilde{Q}_{\ell,m} dt, $$

from which (25) follows, thanks to (27).

<figure><img src="image_9.png" /></figure>

Theorem 4.2 can be used to devise a method to determine for a fixed error $\epsilon$ the smallest truncated state space $\Omega_{\hat{N}}$ by only simulating the process $Y$ as follows:

* Calculate the transition-rate matrix $\tilde{Q}$ of the upper-bounding process $Y$ using Proposition 2.2,

* Solve the master equation for the process $Y$ on $[0, T_f]$ on the finite space $[0, M]$ for $M$ such that

$$ 1 - \tilde{p}_{0:M}^{M}(T_f) \le \epsilon. $$

* Compute the size $\hat{N}$ of the truncated state space using Equation (26).

Finding a value for $M$ such that the quantity $1 - \tilde{p}_{0:M}^{M}(T_f)$ is smaller than $\epsilon$ is easy because $Y$ is a one-dimensional CTMC so a large value for $M$ can be picked directly. Just directly choosing a very large truncation for $X$ would not be feasible since the high-dimensional state space would make solving the master equation too costly in terms of calculation time. Moreover, the same value of $M$ can be used to compute the minimal truncated state space for different values of $\epsilon$ which avoids having to solve a new master equation.

### 4.3.2 Application of the method to the toy model

In this subsection, we illustrate the method derived from Theorem 4.2 on the estimation of the error coming from the truncation of the state space with the upper-bounding process used in the previous section.

Let us recall the following set of classes defined by Equation (20),

$$ \forall \ell \ge 0, S_{\ell}^{+} = \{2x_1 + x_2 + x_3 = \ell\}, $$

for which the optimal upper-bounding process $Y$ is defined through its transition-rate matrix $Q^+$ by Proposition 4.2.

For the numerical simulations, we used the set of parameters $\theta$ detailed in Table 3, $M = 2000$ and $T_f = 4$. With the same notation as in Theorem 4.2, we denote by $\tilde{p}^M$ the solution of the master equation of $Y$ truncated on $[0, M]$ until time $T_f$. The initial condition $X_0$ of the original process is set to (30, 10, 10) which corresponds to $\tilde{p}^M(0) = \delta_{80}$.

The upper bound $E^T$ of the truncation error, defined by Equation (25), is expressed for any $0 \le N \le M$ as

$$ E^T(N) = (1 - \tilde{p}_{0:M}^M(T_f)) + \delta_{N \le 80} + \int_0^{T_f} \sum_{x \in \Omega_N} \sum_{y \notin \Omega_N} \tilde{p}_x^M(t) Q^+(x, y), $$

$$ \text{with } \Omega_N = \bigcup_{\ell=0}^{N} S_{\ell}^{+}. $$

The quantity $E^T(N)$ is plotted in Figure 2 for different final times $T$ (x-axis) and for different values of $N$ (y-axis). Note that the $E^T(N)$ is equal to 1 for $N$ smaller than 80 which is consistent with the initial condition being equal to 80. The orange line (respectively the red line) represents the smallest $N$ for which the upper bound of the error is smaller than a fixed error $\epsilon$ equal to 0.1 (resp. 0.01). Note that the same value of $M$ was used to compute the orange and red lines meaning that we solve the master equation links to $Y$ once. The quantity $E^T(N)$ preserves the stationarity of the process since the orange and red lines does not go to infinity with time.

If we had chosen the classes $\tilde{S}^+$ (see Equation (20)) and the corresponding optimal upper-bounding process $\tilde{Y}$, we would not have obtained a converging truncation bound since $\tilde{Y}$ is explosive.
<figure><img src="image_10.png" /><figcaption>Figure 2: Upper-bound $E^T(N)$ of the truncation error as a function of the final time $T$ and truncation parameter $N$ for the well-selected classes $S^+$. The y-axis represents the parameter $N$ chosen for truncation. The orange line (respectively the red line) represents the smallest $N$ obtained with Equation (26) that leads to an error smaller than 0.1 (resp. 0.01).</figcaption></figure>

Note that since $E^T(N)$ is an upper-bound of a probability, it could take values greater than 1 even if we know that the estimated probability is lower than 1. Knowing that values of $E^T(N)$ greater than 1 are set to 1 (yellow part of the graph).

# 5 An example of enzyme kinetics

In this section, we study the long-time behaviour of a classical reaction network describing the Michaelis-Menton enzyme kinetics [16]. In this network, a substrate is converted into a product by the enzyme. Let us consider the following reaction network composed of four species, $S$ (substrate), $E$ (enzyme), $SE$ (substrate-enzyme complex) and $P$ (product):

$$S + E \underset{c_2 n_{SE}}{\overset{c_1 n_S n_E}{\rightleftharpoons}} SE \underset{a_1 n_P n_E}{\overset{a_2 n_{SE}}{\rightleftharpoons}} P + E \qquad \emptyset \underset{d n_S}{\overset{b}{\rightleftharpoons}} S \qquad P \underset{\emptyset}{\overset{d n_P}{\rightleftharpoons}} \emptyset$$

where all parameters $a_1, a_2, b, c_1, c_2$, and $d$ are strictly positive and where $n_S, n_E, n_{SE}$ and $n_P$ denote the number of molecules of each of the species.

We denote by $X$ the continuous-time Markov chain on $\mathbb{N}^4$ describing the evolution of the vector $(n_S, n_E, n_{SE}, n_P)$. The state space of the process $X$ can be reduced by using the conservation law present in the reaction network:

$$n_E + n_{SE} = E_{tot}$$

where $E_{tot} = n_E(0) + n_{SE}(0)$ is the initial quantity of enzyme and enzyme-substrate complex. Then the process $X$ takes values in the following countable state space:

$$S = \{ (n_S, n_E, n_{SE}, n_P) \in \mathbb{N}^4, n_E + n_{SE} = E_{tot} \}.$$

In the following, we apply our method to construct an upper-bounding process for the set of classes defined by

$$S_{\ell} = \{ (n_S, n_E, n_{SE}, n_P) \in S | n_S + n_P = \ell \}, \quad \forall \ell \in \mathbb{N}.$$

Note that the sets $(S_\ell)_{\ell \in \mathbb{N}}$ build a partition of $S$. The quantities $n_E$ and $n_{SE}$ are bounded by definition of $S$, which implies that all sets $S_\ell$ are finite.

## Proposition 5.1

The matrix $\tilde{Q}$ of the optimal upper-bounding process $Y$ for the set of classes $(S_\ell)_{\ell \in \mathbb{N}}$ verifies, for any
$$ \begin{array}{l} \ell \ge 0, \\ \tilde{Q}_{\ell, \ell+1} = (a_2 + c_2)E_{tot} + b, \quad \tilde{Q}_{\ell, \ell-1} = d\ell, \\ \text{and } \tilde{Q}_{\ell, m} = 0 \text{ if } m \notin \{\ell-1, \ell, \ell+1\}. \text{ Moreover, whatever the choice of parameter values, the process } Y \text{ is positive recurrent, so the process } X \text{ associated with the enzyme kinetics network is also positive recurrent.} \end{array} $$

**Proof:** Starting from a state $(n_S, n_E, n_{SE}, n_P)$ in class $\ell \in \mathbb{N} \setminus \{0\}$,

* after the reactions $S + E \rightarrow SE, P + E \rightarrow SE, S \rightarrow \emptyset$ or $P \rightarrow \emptyset$, the process ends in the class $\ell-1$,

* after the reactions $SE \rightarrow S + E, SE \rightarrow P + E$ or $\emptyset \rightarrow S$, the process ends in the class $\ell+1$.

Thus we have, with Equation (8),

$$ f_{\ell,m}^{-} = \begin{cases} 0 & \text{if } m \le \ell - 2, \\ \min_{x \in S_{\ell}} (c_1 n_S n_E + a_1 n_P n_E + d n_S + d n_P) = d\ell & \text{if } m = \ell - 1, \end{cases} $$

$$ f_{\ell,m}^{+} = \begin{cases} \max_{x \in S_{\ell}} (c_2 n_{SE} + a_2 n_{SE} + b) = (a_2 + c_2)E_{tot} + b & \text{if } m = \ell + 1, \\ 0 & \text{if } m \ge \ell + 2. \end{cases} $$

Using Equation (11), Proposition 2.2 and Lemma 2.7, the process defined on $\mathbb{N}$ by the transition-rate matrix $\tilde{Q}$

$$ \tilde{Q}_{\ell, \ell+1} = (a_2 + c_2)E_{tot} + b, \quad \tilde{Q}_{\ell, \ell-1} = d\ell, $$

and $\tilde{Q}_{\ell,m} = 0$ for $m \notin \{\ell-1, \ell, \ell+1\}$, is an upper-bounding process of the initial process $X$.

The process $Y$ is a birth-death one-dimensional process with constant birth rate and linear death rate so it is always positive recurrent. The result on the ergodicity of the process $X$ follows from Theorem 4.1. ■

# 6 Discussion

Stochastic reaction networks are tremendously flexible and can serve as models to represent stochastic processes in various different scientific areas. However, in most applications, reaction networks comprise a multitude of interacting species, which implies that their dynamics follow a multi-dimensional continuous-time Markov chain that is notoriously difficult to analyze. In most cases, analytical solutions to the governing master equation are not available and even just stochastic simulation may be computationally too expensive to be practical.

In this paper, we have shown that despite these difficulties it may still be possible to establish properties of the CTMC such as its long-term behavior and the existence of a stationary distribution. To achieve this, we defined a partition of the state space that we referred to as classes throughout the manuscript, and showed that the original CMTC transferred to the classes can be upper and lower bounded by one-dimensional birth-death processes that are easy to analyze analytically and numerically. We proved that properties of the bounding processes can be used to gain insights on the original CTMC. For instance, for irreducible CTMCs, the existence of a stationary distribution for the upper-bounding process guarantees that the original CTMC also admits a stationary distribution. These results allow us to gain some level of understanding about complex intractable CTMCs from the study of easily tractable one-dimensional bounding processes.

The quality and accuracy of so-obtained results will generally depend on how tightly the bounding processes bound the original CTMC. We have shown how, for a given partition of the state space of the original CTMC, one can define the rates of the bounding processes that are guaranteed to lead to optimal bounds (Subsections 2.4 and 3.2). However, whether or not these optimal bounds are practically useful still depends on how exactly the classes that partition the state space are chosen. Maybe the most natural choice is to define classes by the total number of individuals/molecules of all interacting species that are present in the system. We have shown, however, that this is a choice that is agnostic to the structure of the reaction network and may lead to bad results. For instance, for the case study in Section 4, this naive choice of classes leads to an (optimal) upper-bounding process that is generally explosive despite the fact that the original CTMC, depending on the parameter values that determine the reaction rates, may admit a stationary distribution. On the other hand, as we have shown in Subsection 4.1, a better choice of the classes that respects the structure of the reaction network allows one to use the lower and upper-bounding processes to determine the asymptotic behavior of the original process in different regions
of parameter space.

In addition to using bounding processes to establish properties such as recurrence or existence of a
stationary distribution, we have shown how the upper-bounding process can be used to guide the choice
of finite state projections for the numerical analysis of CTMCs with countably infinite state spaces. Trun-
cation of the state space is often the only possible approach for numerically approximating the solution
of the master equation of the process. Thereby, approximation errors depend on how the truncation is
chosen and increase with time such that it is often difficult to say a priori how large the truncation needs
to be to achieve a given error. We have shown in Subsection 4.3 that the error that a given truncation
will create in the approximation of the solution of the master equation of the original CTMC can be
bounded based on the solution of the master equation of the upper-bounding process, which is easy to
approximate numerically since the upper-bounding process is one-dimensional. As we have shown in our
case study, this allows us to straightforwardly calculate necessary truncation sizes to achieve given error
guarantees for any desired time horizon (Figure 2).

In sum, we believe that our scheme for constructing one-dimensional bounding processes for multi-
dimensional CTMCs, together with our proofs on how properties of these bounding processes carry over
to the original CTMC, will help researchers to gain a better understanding of the properties of the com-
plex stochastic reaction networks that emerge from applications in biology and other fields.

**Acknowledgments:** This work was supported by the ANR (Opt-MC project, ANR-22-CE45-0019) and the European Research Council (ERC-2022-STG, BridgingScales, grant agreement 101075989).

References

[1] G. Albi, R. Chignola, and F. Ferrarese. Efficient ensemble stochastic algorithms for agent-based models with spatial predator–prey dynamics. *Math. Comput. Simul.*, 199:317–340, 2022.

[2] L. J. Allen. *An introduction to stochastic processes with applications to biology*. CRC press, 2010.

[3] D. F. Anderson, D. Cappelletti, and T. G. Koyama, M. and Kurtz. Non-explosivity of Stochastically Modeled Reaction Networks that are Complex Balanced. *Bull. Math. Biol.*, 80:2561–2579, 2018.

[4] D. F. Anderson and S. Cotter. Product-Form Stationary Distributions for Deficiency Zero Networks with Non-mass Action Kinetics. *Bull. Math. Biol.*, 78:2390–2407, 2016.

[5] D. F. Anderson and T. G. Kurtz. Continuous time Markov chain models for chemical reaction networks. In *Design and Analysis of Biomolecular Circuits*, pages 3–42. Springer New York, 2011.

[6] D. F. Anderson and T. D. Nguyen. Results on stochastic reaction networks with non-mass action kinetics. *Math. Biosci. Eng.*, 16(4):2118–2140, 2019.

[7] R. Anderson and R. May. *Infectious diseases of humans: dynamics and control*. Oxford University Press, 1992.

[8] W. J. Anderson. *Continuous-time Markov chains*. Springer Series in Statistics. Springer New York, 1991.

[9] A. Andreychenko, L. Bortolussi, R. Grima, P. Thomas, and V. Wolf. *Distribution approximations for the chemical master equation: comparison of the method of moments and the system size expansion*, page 39–66. Springer International Publishing, 2017.

[10] J. Bai and P. Wang. Conditional Markov chain and its application in economic time series analysis. *J. Appl. Econ.*, 26(5):715–734, 2011.

[11] A. M. Barreto and M. D. Fragoso. Lumping the states of a finite Markov chain through stochastic factorization. *IFAC Proc.*, 44(1):4206–4211, 2011.

[12] A. Briggs and M. Sculpher. An introduction to Markov modelling for economic evaluation. *Phar-macoeconomics*, 13(4):397–409, 1998.

[13] F. Campos, S. Bruno, Y. Fu, D. Del Vecchio, and R. Williams. Comparison theorems for stochastic chemical reaction networks. *Bull. Math. Biol.*, 85(39), 2023.
[14] N. Champagnat and D. Villemonais. General criteria for the study of quasi-stationarity. *Electron. J. Probab.*, 28:1–84, 2023.

[15] B. D. Coleman and T. G. Fox. A multistate mechanism for homogeneous ionic polymerization. ii. the molecular weight distribution. *J. Am. Chem. Soc.*, 85(9):1241–1244, 1963.

[16] A. Cornish-Bowden. *Fundamentals of Enzyme Kinetics*. Portland Press, London, 2004.

[17] E. B. Davies. *Linear operators and their spectra*. Cambridge University Press, 2007.

[18] L. Einkemmer, J. Mangott, and M. Prugger. A hierarchical dynamical low-rank algorithm for the stochastic description of large reaction networks. *arXiv preprint*, 2024.

[19] M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain. Stochastic gene expression in a single cell. *Science*, 297(5584):1183–1186, 2002.

[20] M. Függer, T. Nowak, and J. Rybicki. Majority consensus thresholds in competitive Lotka-Volterra populations. In *Proceedings of the 43rd ACM Symposium on Principles of Distributed Computing*, PODC '24, pages 76–86. Association for Computing Machinery, 2024.

[21] A. Ganguly, T. Petrov, and H. Koepp. Markov chain aggregation and its applications to combinatorial reaction networks. *J. Math. Biol.*, 69:767–797, 2014.

[22] D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. *J. Phys. Chem.*, 81(25):2340–2361, 1977.

[23] D. T. Gillespie. Approximate accelerated stochastic simulation of chemically reacting systems. *J. Chem. Phys.*, 115(4):1716–1733, 2001.

[24] A. Gupta and M. Khammash. Computational identification of irreducible state-spaces for stochastic reaction networks. *SIAM J. Appl. Dyn. Syst.*, 17(2):1213–1266, 2018.

[25] A. Gupta, J. Mikelson, and M. Khammash. A finite state projection algorithm for the stationary solution of the chemical master equation. *J. Chem. Phys.*, 147(15), 2017.

[26] M. Hansen, C. Wiuf, and C. Xu. Stationary measures of continuous-time markov chains with applications to stochastic reaction networks. *Adv. Appl. Probab.*, 57(3):1031–1067, 2025.

[27] L. Hoessly. Stationary distributions via decomposition of stochastic reaction networks. *J. Math. Biol.*, 82:67, 2021.

[28] J. Kuntz, P. Thomas, G.-B. Stan, and M. Barahona. Stationary distributions of continuous-time Markov chains: a review of theory and truncation-based approximations. *SIAM Review*, 63(1):3–64, 2021.

[29] T. Lindvall. *Lectures on the coupling method*. Dover Publications, Inc, 2002.

[30] Y. Lu, C. Gao, and D. Dochain. Chemical reaction network decomposition technique for stability analysis. *Automatica*, 146:110650, 2022.

[31] S. P. Meyn and R. L. Tweedie. Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes. *Adv. Appl. Probab.*, 25(3):518–548, 1993.

[32] B. Munsky and M. Khammash. The finite state projection algorithm for the solution of the chemical master equation. *J. Chem. Phys.*, 124(4):044104, 2006.

[33] G. Rubino and B. Sericola. A finite characterization of weak lumpable Markov processes. Part II: The continuous time case. *Stoch. Process. Their Appl.*, 45(1):115–125, 1993.

[34] V. Shahrezaei and P. S. Swain. Analytical distributions for stochastic gene expression. *Proc. Natl. Acad. Sci. U.S.A.*, 105(45):17256–17261, 2008.

[35] V. Wolf, R. Goel, M. Mateescu, and T. A. Henzinger. Solving the chemical master equation using sliding windows. *BMC Syst. Biol.*, 4:1–19, 2010.

[36] C. Xu, M. C. Hansen, and C. Wiuf. Full classification of dynamics for one-dimensional continuous-time Markov chains with polynomial transition rates. *Adv. Appl. Probab.*, 55(1):321–355, 2023.
# A Technical calculations

## Lemma A.1

Let $(a, b, c_1, c_2) \in (\mathbb{R}^+)^2 \times \mathbb{R}^2$ and for $x \in \mathbb{R}^2$,

$$F(x) = ax_1^2 + bx_1x_2 + c_1x_1 + c_2x_2.$$

The goal of this lemma is to compute the minimum of the function $F$ on a subset of the upper right quadrant of the plane. To this purpose, we define for any $(\ell, \alpha_1, \alpha_2) \in \mathbb{N}_*^3$,

$$M_{\ell} = \min_{x \in D_{\ell}} F(x),$$

with $D_{\ell} = \{x \in \mathbb{R}_+^2, \alpha_1 x_1 + \alpha_2 x_2 \le \ell\}$.

Then the elements of $D_{\ell}$ where the minimum of $F$ can be reached (infinite values are allowed if the denominator is null) are

$$p_1 = \left(-\frac{c_2}{b}, \frac{2ac_2-bc_1}{b^2}\right), \quad p_2 = (0, 0), \quad p_3 = \left(0, \frac{\ell}{\alpha_2}\right),$$

$$p_4 = \left(-\frac{c_1}{2a}, 0\right), \quad p_5 = \left(\frac{\ell}{\alpha_1}, 0\right), \quad p_6 = \left(-d_1\ell + d_2, \frac{1+d_1\alpha_1}{\alpha_2}\ell - \frac{\alpha_1}{\alpha_2}d_2\right),$$

with $d_1 = \frac{b}{2(a\alpha_2 - b\alpha_1)}$ and $d_2 = \frac{\alpha_1 c_2 - \alpha_2 c_1}{2(a\alpha_2 - b\alpha_1)}$.

So we have, with $\mathcal{P} = \{p_i, 1 \le i \le 6\}$,

$$M_{\ell} = \min_{p \in \mathcal{P} \cap D_{\ell}} F(p).$$

**Proof:** Let $\bar{D}_{\ell}$ be the interior of $D_{\ell}$ and $\mathring{D}_{\ell}$ the boundary of $D_{\ell}$, that is

$$\bar{D}_{\ell} = \{x \in \mathbb{R}_{+}^{2}, \alpha_{1}x_{1} + \alpha_{2}x_{2} < \ell\},$$

$$\mathring{D}_{\ell} = \left\{ (0, x_2), \ 0 \le x_2 \le \frac{\ell}{\alpha_2} \right\} \cup \left\{ (x_1, 0), \ 0 \le x_1 \le \frac{\ell}{\alpha_1} \right\} \cup \{x \in \mathbb{R}_{+}^{2}, \ \alpha_1 x_1 + \alpha_2 x_2 = \ell\}.$$

Two cases must be considered.

1. For $x \in \bar{D}_{\ell}$, the minimum $M_{\ell}$ is reached on $x$ if $x$ is a critical point that is

$$\frac{\partial F}{\partial x_1}(x) = 2ax_1 + bx_2 + c_1 = 0,$$

$$\frac{\partial F}{\partial x_2}(x) = bx_1 + c_2 = 0.$$

The only solution of this system is $p_1 = \left(-\frac{c_2}{b}, \frac{2ac_2-bc_1}{b^2}\right)$.

2. For $\mathring{D}_{\ell}$, we will consider separately each set in the definition of $\mathring{D}_{\ell}$.

* For $x \in \{(0, x_2), 0 \le x_2 \le \frac{\ell}{\alpha_2}\}$, $F(x) = c_2x_2$.

The minimum is reached for $x_2 = 0$ or $x_2 = \frac{\ell}{\alpha_2}$ according to the sign of $c_2$. Thus, the corresponding points in the plane are $p_2 = (0, 0)$ and $p_3 = (0, \frac{\ell}{\alpha_2})$.

* For $x \in \{(x_1, 0), 0 \le x_1 \le \frac{\ell}{\alpha_1}\}$, $F(x) = ax_1^2 + c_1x_1 = x_1(ax_1 + c_1)$.

If $a > 0$, the minimum is reached for $x_1 = -\frac{c_1}{2a}$, which corresponds to $p_4 = (-\frac{c_1}{2a}, 0)$.

If $a = 0$, the minimum is reached for $x_1 = 0$ or $x_1 = \frac{\ell}{\alpha_1}$, which corresponds to $p_2$ and $p_5$.

* For $x \in \{x \in \mathbb{R}_+^2, \alpha_1 x_1 + \alpha_2 x_2 = \ell\}$,

$$F(x) = ax_1^2 + b \frac{x_1}{\alpha_2} (\ell - \alpha_1 x_1) + c_1 x_1 + \frac{c_2}{\alpha_2} (\ell - \alpha_1 x_1) \\ = \left(a - b \frac{\alpha_1}{\alpha_2}\right) x_1^2 + \left(\frac{b}{\alpha_2} \ell + c_1 - c_2 \frac{\alpha_1}{\alpha_2}\right) x_1 + \frac{c_2}{\alpha_2} \ell.$$

If $(a - b \frac{\alpha_1}{\alpha_2}) \le 0$, the minimum is reached for $x_1 = 0$ or $x_1 = \frac{\ell}{\alpha_1}$, which corresponds to $p_3$ and $p_5$.

If $(a - b \frac{\alpha_1}{\alpha_2}) > 0$, an easy computation shows that the minimum is reached for $x_1 = -d_1\ell + d_2$ which corresponds to $p_6$.
## A.1 Proof of Proposition 4.1

The proof is composed of two steps: first the calculation of the quantities $\hat{f}^-$ and $\hat{f}^+$ defined by Equation (17), then with Equation (18), we obtain the expressions of $\hat{U}^-$ and $\hat{U}^+$.

In this proposition, the set of classes used is

$$ \forall \ell \in \mathbb{N}, \quad S_{\ell}^{-} = \{x \in \mathbb{N}^3, 2x_1 + 2x_2 + 5x_3 = \ell\}. $$

Starting from a state $x$ in class $\ell \ge 5$,

* after the reaction $X_3 \to \emptyset$, the process ends in the class $\ell - 5$,

* after the reaction $X_1 \to \emptyset$ or $X_2 \to \emptyset$, the process ends in the class $\ell - 2$,

* after the reaction $X_1 + X_2 \to X_3$, the process ends in the class $\ell + 1$,

* after the reactions $\emptyset \to X_1$ and $2X_1 \to 3X_2$, the process ends in the class $\ell + 2$.

Thus we have, with Equation (17),

$$ \begin{aligned} \hat{f}_{\ell,m}^{-} &= \begin{cases} 0 & \text{if } m < \ell - 5, \\ \max_{x \in S_{\ell}^{-}} d_3 x_3 = \frac{d_3}{5} \ell & \text{if } \ell - 5 \le m \le \ell - 3, \\ \max_{x \in S_{\ell}^{-}} (d_1 x_1 + d_2 x_2 + d_3 x_3) = \max \left( \frac{d_1}{2}, \frac{d_2}{2}, \frac{d_3}{5} \right) \ell & \text{if } \ell - 2 \le m < \ell, \end{cases} \\ \hat{f}_{\ell,m}^{+} &= \begin{cases} \min_{x \in S_{\ell}^{-}} (b_1(x_2+x_3) + b_2 + \alpha x_1(x_1-1) + \beta x_1 x_2) & \text{if } m = \ell+1, \\ \min_{x \in S_{\ell}^{-}} (b_1(x_2+x_3) + b_2 + \alpha x_1(x_1-1)) & \text{if } m = \ell+2, \\ 0 & \text{if } \ell+3 \le m. \end{cases} \end{aligned} $$

Moreover,

$$ \begin{aligned} \hat{f}_{\ell, \ell+2}^{+} &= \min_{x \in S_{\ell}^{-}} (b_1(x_2 + x_3) + b_2 + \alpha x_1(x_1 - 1)) \\ &= \min_{x \in D_{\ell}} \left( \alpha x_1^2 - \left(\alpha + \frac{2}{5}b_1\right)x_1 + \frac{3}{5}b_1x_2 \right) + b_2 + \frac{b_1}{5}\ell \quad \text{where } D_{\ell} = \{x \in \mathbb{R}_+^2, 2x_1 + 2x_2 \le \ell\} \end{aligned} $$

Using Lemma A.1, for $\ell$ large enough the minimum is reached in $p_4 = \left(\frac{\alpha + \frac{2}{5}b_1}{2\alpha}, 0\right)$ and

$$ \hat{f}_{\ell, \ell+2}^{+} = \frac{b_1}{5}\ell + A(\theta), $$

$$ \text{with } A(\theta) = b_2 - \frac{(\alpha + \frac{2}{5}b_1)^2}{4\alpha}. $$

With the same computation, $\hat{f}_{\ell, \ell+1}^{+} = \frac{b_1}{5}\ell + A(\theta)$.

The result follows from Equation (18).

## A.2 Proof of Proposition 4.2

The proof is composed of two steps: first the calculation of the quantities $f^-$ and $f^+$ defined by Equation (8), then with Equation (11), we obtain the expressions of $U^-$ and $U^+$ from which we deduce the expression of the matrix $Q^+$.

In this proposition, the set of classes used is

$$ \forall \ell \in \mathbb{N}, \quad S_{\ell}^{+} = \{x \in \mathbb{N}^3, 2x_1 + x_2 + x_3 = \ell\}. $$

Starting from a state $x$ in class $\ell \ge 2$,

* after the reactions $X_1 \to \emptyset$ or $X_1 + X_2 \to X_3$, the process ends in the class $\ell - 2$,

* after the reactions $2X_1 \to 3X_2$, $X_2 \to \emptyset$ or $X_3 \to \emptyset$, the process ends in the class $\ell - 1$,
• after the reaction $\emptyset \rightarrow X_1$, the process ends in the class $\ell + 2$.

Thus we have, with Equation (8),

$$f_{\ell,m}^{-} = \begin{cases} 0 & \text{if } m < \ell - 2, \\ \min_{x \in S_{\ell}^{+}} (d_1 x_1 + \beta x_1 x_2) = 0 & \text{if } m = \ell - 2, \\ \min_{x \in S_{\ell}^{+}} \left( \alpha x_1 (x_1 - 1) + \sum_{i=1}^{3} d_i x_i \right) & \text{if } m = \ell - 1, \end{cases}$$

$$f_{\ell,m}^{+} = \begin{cases} \max_{x \in S_{\ell}^{+}} (b_{1}(x_{2} + x_{3}) + b_{2}) = b_{1}\ell + b_{2} & \text{if } \ell + 1 \le m \le \ell + 2, \\ 0 & \text{if } \ell + 3 \le m. \end{cases}$$

Moreover,

$$\begin{aligned} f_{\ell,\ell-1}^{-} &= \min_{x \in S_{\ell}^{+}} \left( \alpha x_{1}(x_{1}-1) + \sum_{i=1}^{3} d_{i}x_{i} \right) \\ &= \min_{x \in D_{\ell}} \left( \alpha x_{1}^{2} + (d_{1}-\alpha-2d_{3})x_{1} + (d_{2}-d_{3})x_{2} \right) + d_{3}\ell, \quad \text{where } D_{\ell} = \{x \in \mathbb{R}_{+}^{2}, 2x_{1} + x_{2} \le \ell\}. \end{aligned}$$

Using Lemma A.1, for $\ell$ large enough the minimum is reached in

$$p_4 = \left( \frac{\alpha + 2d_3 - d_1}{2\alpha}, 0 \right) \quad \text{or} \quad p_6 = \left( \ell - \frac{\alpha + 2d_2 - d_1}{\alpha}, 0 \right),$$

and

$$f_{\ell,\ell-1}^{-} = \min(d_2, d_3)\ell - C(\theta),$$

$$\text{with } C(\theta) = \frac{(\alpha + 2 \min(d_2, d_3) - d_1)^2}{4\alpha} \mathbf{1}_{\alpha + 2 \min(d_2, d_3) - d_1 \ge 0}(\theta).$$

The result follows directly from Equation (11).

## A.3 Proof of Proposition 4.3

The proof is composed of two steps: first the calculation of the quantities $f^-$ and $f^+$ defined by Equation (8), then with Equation (11), we obtain the expressions of $U^-$ and $U^+$ from which we deduce the expression of the matrix $\tilde{Q}^+$.

In this proposition, the set of classes used is

$$\forall \ell \in \mathbb{N}, \quad \tilde{S}_{\ell}^{+} = \{x \in \mathbb{N}^{3}, x_{1} + x_{2} + x_{3} = \ell\}.$$

Starting from a state $x$ in class $\ell \ge 1$,

• after the reactions $X_1 + X_2 \rightarrow X_3$, $X_1 \rightarrow \emptyset$, $X_2 \rightarrow \emptyset$ or $X_3 \rightarrow \emptyset$, the process ends in the class $\ell - 1$,

• after the reactions $\emptyset \rightarrow X_1$ or $2X_1 \rightarrow 3X_2$, the process ends in the class $\ell + 1$.

Thus we have, with Equation (8),

$$f_{\ell,m}^{-} = \begin{cases} 0 & \text{if } m < \ell - 1, \\ \min_{x \in \tilde{S}_{\ell}^{+}} \left( \beta x_{1}x_{2} + \sum_{i=1}^{3} d_{i}x_{i} \right) = \min(d_{1}, d_{2}, d_{3})\ell & \text{if } m = \ell - 1, \end{cases}$$

$$f_{\ell,m}^{+} = \begin{cases} \max_{x \in \tilde{S}_{\ell}^{+}} (\alpha x_{1}(x_{1}-1) + b_{1}(x_{2} + x_{3}) + b_{2}) = \alpha\ell^{2} - \alpha\ell + b_{2} & \text{if } m = \ell + 1, \\ 0 & \text{if } \ell + 2 \le m. \end{cases}$$

The result follows directly from Equation (11).