Fast reactions with non-interacting species in stochastic reaction networks

1.   Introduction

Reaction network theory offers a quantitative framework for biochemistry, systems biology, and cellular biology, by enabling the modeling of biological systems. Deterministic models have been the main focus area with contributions going back more than hundred years ^[1]. However, the increasing interest in living systems at the cellular level motivates the use of stochastic models to describe the variation and noise found in systems with low molecule numbers. Typically, stochastic models are continuous-time Markov chains (CTMCs) on the state space $\mathbb{Z}^n_{\ge0}$, where a state represents the vector of molecule numbers of the $n$ species in the system.

We are here concerned with the transient behavior of such CTMCs, in contrast to the stationary behavior (the existence of stationary distributions). In practice, variations in molecule numbers and reaction rates might yield phenomena that evolve on different time-scales, enabling simplifications ^[2]. In particular, we are interested in systems with two time-scales, also known as slow-fast systems, where a set of reactions are fast (in a relative sense) compared to the remaining (slow) reactions. The objective is to approximate, in a mathematically rigorous way, the dynamics of the original CTMC, with another CTMC in smaller dimension (with fewer species). This other CTMC should ideally be interpreted as a reaction network, obtained by reduction of the original reaction network. We will give conditions for when this can be done.

In the deterministic setting, the heuristic quasi-steady state approximation (QSSA) ^[3] and singular perturbation theory in the sense of Tikhonov-Fenichel ^[4,5] have been the main means to derive lower dimensional reduced models, where ^[5] is the classic reference of Tikhonov, while ^[4] is from Fenichel. (See ^[6] or ^[7] and references therein for the relationship between the two approaches and when the QSSA is valid, and ^[8] for comparison with stochastic QSSA.) A motivation for our work is the situation for which reactions involving so-called non-interacting species ^[9] in the reactant has reaction rate constants scaled by $1/\epsilon$, where $\epsilon > 0$ is a small number ^[10]. We will consider a similar situation for stochastic reaction networks and point out differences between the two settings.

Many stochastic studies consider physical or heuristic-based derivations to extract reduced reaction networks, for example, SDEs or hybrid models, or ad-hoc reductions ^[11,12,13]. ^[11] comments on general deterministic approaches, ^[12] studies heuristics for stochastic single scale reductions, while ^[13] reviews practical aspects of deterministic and stochastic reductions. One stream of research eliminates species using heuristic projection arguments ^[14,15], where ^[14] treats some reversible cases and ^[15] irreversible examples. Rigorous simplifications and reductions often follow scaling limits of Markov processes in a multi-scale setting ^[16]. These might be applied to concrete examples, if certain conditions are satisfied and the different scaling parameters are balanced (in a specific sense), as done for specific examples of reaction networks ^[17], more generally ^[18], or for spatial reaction networks ^[19]. Scaling laws for a special class of reaction networks with intermediate species (a special type of non-interacting species) and their explicit reduction are given in ^[20]. We extends this work to reaction networks with non-interacting species on two time-scales.

These two time-scales separate the transition intensities of the CTMC into a fast and a slow component such that the $Q$-matrix has the following form

The transition intensities of the reactions in the network are divided into two kinds:

● Fast reactions with scaled transition intensity $\lambda^ \epsilon_{y \to y{^{\prime}}}(x): = \frac{1}{ \epsilon} \lambda_{y \to y{^{\prime}}}(x)$.

● Slow reactions with unscaled transition intensity $\lambda^ \epsilon_{y \to y{^{\prime}}}(x): = \lambda_{y \to y{^{\prime}}}(x)$.

The fast reactions are determined by non-interacting species. Our work is inspired by previous work on non-interacting species in deterministic systems that focussed on algebraic reductions ^[9], applications to biological examples ^[21], the connection to singular perturbation ^[10], and intermediate species in stochastic models ^[20]. The technical part is based on singularly perturbed Markov chains ^[22] and watched Markov chains ^[23]. In order to derive the limiting dynamics systematically, we associate a graph to the reaction network that captures the dynamics of the fast reactions. This enables the definition of the reduced reaction network and its dynamics. We give limit theorems for the approximation of the original CTMC to the reduced CTMC on compact time intervals. Furthermore, we study the case of time-heterogeneous CTMCs, and show that the same reductions work.

As an example, consider a mass-action reaction network as follows

with two fast reactions determined by the presence of the non-interacting species $U_1$, that degrades fast. The reduced reaction network is given as

with transition intensities $\frac{ \kappa_1 \kappa_2z_{S_1}(z_{S_2}+1)}{ \kappa_2(z_{S_2}+1)+ \kappa_3}$ for the first reaction and $\frac{ \kappa_1 \kappa_3z_{S_1}}{ \kappa_2(z_{S_2}+1)+ \kappa_3}$ for the second. After creation of a $U_1$ molecule (and a $S_2$ molecule), then it might be degraded either by consumption of the $S_2$ molecule, resulting in the net reaction $S_1\xrightarrow{ }S_3$, or without consumption of the $S_2$ molecule, resulting in the net reaction $S_1\xrightarrow{ }S_2+S_4$. In both cases, the transition intensities reflect that the presence of $S_1$ is required for the reactions to take place.

We next outline the content, where in $\S$ 2 preliminaries on graph theory and reaction networks are covered. In $\S$ 3, we introduce non-interacting species and introduce a graph that is used to define the reduced reaction network by elimination of the non-interacting species. In $\S$ 4, we study transient approximations for stochastic reaction networks with non-interacting species via the previously introduced reduction. In $\S$ 5, we give realistic examples, study sufficient conditions for reductions and compare stationary properties of the reduction with the original reaction network. Finally in $\S$ 6, we discuss the results, approach, and elaborate on the relation to the literature. In the Appendices $\S$ A, $\S$ B, $\S$ C, we give proofs as well as brief introductions to the theory on singularly perturbed CTMCs and watched CTMCs.

2.   Preliminaries

2.1.   Notation

Let $\mathbb{R}^p$ be the real $p$-dimensional space, and $\mathbb{R}^p_{ > 0}$ ($\mathbb{R}^p_{\geq0}$) the subset of elements of $\mathbb{R}^p$ with strictly positive (non-negative) entries in all components. A vector $y\in \mathbb{R}^p$ is written as $(y^1, \cdots, y^p)$, where $y^i$ is the $i$-th component. For vectors $y_1, \cdots, y_k \in \mathbb{R}^p$, $\rm{max}(y_1, \cdots, y_k)$ denotes the component-wise maximum, and $y_1\ge y_2$ denotes component-wise inequality. The inner product between $y_1$ and $y_2$ is denoted $\langle y_1, y_2 \rangle$. The cardinality of a set $A$ is denoted $|A|$.

2.2.   Graph theory

$\mathcal{G} = (\mathcal{V}, \mathcal{E})$ a directed graph consists of the set of vertices $\mathcal{V}$ and edge set $\mathcal{E}$. A directed subgraph of $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ is a directed graph $\mathcal{G}{^{\prime}} = (\mathcal{V}{^{\prime}}, \mathcal{E}{^{\prime}})$ with $\mathcal{V}{^{\prime}}\subseteq \mathcal{V}, \mathcal{E}{^{\prime}}\subseteq \mathcal{E}$ with $\mathcal{E}{^{\prime}}$ on $\mathcal{V}{^{\prime}}$. A walk is a directed path $V_{i_1}\rightarrow V_{i_2}\rightarrow \cdots \rightarrow V_{i_{l-1}}\rightarrow V_{i_l}$ (potentially listed as the corresponding sequence of edges).

A multi-digraph $\mathcal{G}$ is a directed graph where multiple edges between the same vertices are allowed. In particular, a multi-digraph comes with two functions

where the functions $s, t$ are the source and target function, respectively. Both self-edges (edges $e$ with $t(e) = s(e)$) and parallel edges are possible.

2.3.   Reaction networks (RNs)

A RN on a finite set $\mathcal{S}$ is a digraph $\mathcal{N} = (\mathcal{C}, \mathcal{R})$, where $\mathcal{S}$ is a finite set of species $\mathcal{S} = \{S_1, \cdots, S_n\}$, $\mathcal{C}$ a potentially infinite set of complexes and $\mathcal{R}$ a potentially infinite set of reactions $\mathcal{R} = \{r_1, r_2\cdots\}$. Complexes are non-negative linear combinations of species, $y = \sum_{i = 1}^ny^i S_i$, identified with vectors $y = (y^1, \ldots, y^n)$ in $\mathbb{Z}_{\geq 0}^n$. Reactions are directed edges between complexes, written as $r_i\colon y_i\to y_i{^{\prime}}$ or, generically, as $r\colon y\to y{^{\prime}}$, potentially omitting $r_i, r$. A reaction is said to consume the reactant $y$ and create the product $y{^{\prime}}$. An RN is said to be finite if $\mathcal{R}$ is finite, and otherwise it is infinite.

We diverge in two ways from the standard definition of RNs: trivial reactions $r\colon y\to y$ (self-loops) are allowed, and the numbers of complexes and reactions are allowed to be infinite. Both extensions are useful when dealing with reduced RNs. From a dynamical point of view, trivial reactions might always be ignored as the dynamics is the same with and without them. Realistic model of bursty gene expression with infinitely many reactions have been proposed in the literature ^[24,25], where ^[24] focusses on corresponding theory for extensions of density-dependent Markov processes and ^[25] heuristically investigate such models. However, our motivation is not to accommodate such examples, but rather to ensure that an RN obtained by reduction of a finite RN is also an RN, finite or infinite. The construction of the reduced RN also holds even if the original RN is infinite.

Definition 2.1. (ⅰ) Two species interact if they both appear in a complex of $\mathcal{C}$.

(ⅱ) A subset $\mathcal{U}\subseteq \mathcal{S}$ is non-interacting if it contains no pair of interacting species, and the stoichiometric coefficients of the species in $\mathcal{U}$ in all complexes are either $0$ or $1$. The species of $\mathcal{U}$ are said to be non-interacting, and the species in $\mathcal{S}\setminus \mathcal{U}$ are said to be core species.

Example 2.2. Consider a two-substrate Michaelis Menten mechanism [26,Section 3.1.2]:

where $E$ is an enzyme catalyzing the conversion of two substrates $A, B$ into two other substrates $P, Q$ by means of transient (or intermediate) steps; here $EA$ and $EAB$ are known as intermediate complexes formed by binding of the molecules in the reactants.

This RN has species set $\mathcal{S} = \{E, A, B, EA, EAB, P, Q\}$ and complex set $\mathcal{C} = \{E+A, EA, EA+B, EAB, E+P+Q\}$.

The sets $\mathcal{U}_1 = \{EA, EAB\}$, $\mathcal{U}_2 = \{EA, P\}$ and $\mathcal{U}_3 = \{EAB\}$ are non-interacting, with corresponding core species $\mathcal{O}_1 = \mathcal{S}\setminus \mathcal{U}_1 = \{E, A, B, P, Q\}$, $\mathcal{O}_2 = \mathcal{S}\setminus \mathcal{U}_2 = \{E, A, B, EAB, Q\}$, and $\mathcal{O}_3 = \mathcal{S}\setminus \mathcal{U}_3 = \{E, A, EA, B, P, Q\}$. In particular, a species might be a core or a non-interacting species, depending on the choice of $\mathcal{U}$.

In this paper, it is convenient to work with the directed stoichiometric subspace of $\mathcal{N} = (\mathcal{C}, \mathcal{R})$, defined as

(note, this definition allows the RN to be infinite). For $v\in \mathbb{R}^n$, the set $(v+\mathcal{T})\cap \mathbb{R}^n_{\geq 0}$ defines a directed stoichiometric compatibility class of $\mathcal{N}$. The RN is said to be conservative, respectively, sub-conservative, if there exists a positive vector $c\in \mathbb{R}^n_{ > 0}$, such that $\langle c, y'-y\rangle = 0$, respectively, $\langle c, y'-y\rangle\le 0$ for any reaction $r\colon y \to y{^{\prime}}\in \mathcal{R}$. A sub-conservative RN has compact directed stoichiometric compatibility classes (but not stoichiometric compatibility classes). For a weakly reversible RN, the directed and the undirected stoichiometric subspaces are the same.

2.4.   Stochastic reaction networks (SNRs)

An SRN is an RN together with a CTMC $X(t), t\ge 0$, on $\mathbb{Z}_{\geq 0}^n$, modeling the number of molecules of each species over time. A reaction $r\colon y\to y{^{\prime}}$ fires with transition intensity $\lambda_r(x)$, in which case the state jumps from $X(t) = x$ to $x+y{^{\prime}}-y$ ^[27]. The Markov process with transition intensities $\lambda_r\colon \mathbb{Z}_{\geq 0}^n\to \mathbb{R}_{\geq 0}$, $r\in \mathcal{R}$, has $Q$-matrix

For (stochastic) mass-action kinetics, the transition intensity for $r\colon y\to y{^{\prime}}$ is

where $x!: = \prod_{i = 1}^nx_i!$, and $\kappa_r$ is a positive reaction rate constant ^[27]. If there are infinitely many reactions, we assume

such that the corresponding CTMC is well-defined in the sense that it has no instantaneous jumps ^[28].

When the reactions are indexed, $\mathcal{R} = \{r_1, r_2, \ldots\}$, we occasionally write $\lambda_i$ and $\kappa_i$, $i = 1, 2, \ldots$, for convenience. The following assumption holds in particular for stochastic mass-action kinetics.

Assumption 1. For all reactions, $r\colon y\to y{^{\prime}}\in \mathcal{R}$, the transition intensity $\lambda_r\colon \mathbb{Z}_{\geq 0}^n\to \mathbb{R}_{\geq 0}$ satisfies the following

Under Assumption 1, $(X(t))_{t\ge 0}$, stays in $\mathbb{Z}^n_{\geq 0}$, if $X(0) = x\in \mathbb{Z}^n_{\geq 0}$. In particular, $(X(t))_{t\ge 0}$, is confined to the directed stoichiometric compatibility class $(x+\mathcal{T})\cap \mathbb{Z}^n_{\geq 0}$. Assumption 1 is fundamental and enforces the reactions to be compatible with the transition intensities of the CTMC.

3.   Elimination of non-interacting species through fast reactions

3.1.   Notation for non-interacting species

In the following, we focus on RNs and SRNs with non-interacting species. To fix notation, let $\mathcal{U}\subseteq \mathcal{S}$ be a non-interacting subset of species. For simplicity, we let $\mathcal{U} = \{U_1, \dots, U_m\}$ and $\mathcal{O} = \mathcal{S}\setminus \mathcal{U} = \{S_1, \dots, S_p\}$ ($p = n-m$), such that $\mathcal{S} = \mathcal{O}\cup \mathcal{U} = \{S_1, \dots, S_p, U_1, \dots, U_m\}$. We let

be the projections onto the first $p$ coordinates and last $m$ coordinates of $\mathbb{R}^n$, respectively. Consequently, we denote a state by $x = (z, u)\in \mathbb{Z}^p_{\ge0}\times \mathbb{Z}^m_{\ge0} = \mathbb{Z}^n_{\ge0}$.

Define the following sets of reactions

such that $\mathcal{R}_{\mathcal{U}}\cup \mathcal{R}_{\mathcal{U}}'$ are the reactions that involve species in $\mathcal{U}$.

We consider a subset of fast (a terminology to be motivated below) reactions $\mathcal{F}\subseteq \mathcal{R}_{\mathcal{U}}$ with the following structural property.

Definition 3.1. The reactions in $\mathcal{F}$ are proper w.r.t. $\mathcal{U}$, that is, any non-interacting species is part of a sequence of reactions in $(\mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}})\cup \mathcal{F}$ of the form

with $i_0, \ldots, i_l\in\{1, \ldots, k\}$, $\rho_ \mathcal{U}(y_{i_0}) = \rho_ \mathcal{U}(y_{i_l}{^{\prime}}) = 0$ and $\rho_ \mathcal{U}(y_{i_j}) = \rho_ \mathcal{U}(y_{i_{j-1}}{^{\prime}})\neq 0$ for $j = 1, \cdots, l$. Such a sequence of reactions is called a fast chain.

Hence, any molecule of a non-interacting species can be degraded through a sequence of fast reactions (provided sufficient molecule numbers of core species).

Example 3.2. Recall Example 2.2. The set $\mathcal{U}_1 = \{EA, EAB\}$ is a non-interacting set and $\mathcal{F} = \mathcal{R}_{ \mathcal{U}_1}$ is a set of fast proper reactions. For example, the following is a fast chain:

Also, $\mathcal{U}_2 = \{EA, P\}$ is a non-interacting set of species, but $\mathcal{F} = \mathcal{R}_{ \mathcal{U}_2}$ is not a set of fast proper reactions, as there are no fast reactions with $P$ in the reactant.

3.2.   The reduced reaction network

We introduce a labeled multi-digraph to capture the conversion and creation of non-interacting species through fast chains. This graph is similar to the one introduced in ^[21]. We later use this graph to define a reduced RN and a reduced SRN by elimination of non-interacting species.

Definition 3.3. Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be an RN on $\mathcal{S}$, $\mathcal{U}\subseteq \mathcal{S}$ a set of non-interacting species, and $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ a set of fast proper reactions. Let $\mathcal{G}_{\mathcal{U}, \mathcal{F}} = (\mathcal{V}_{\mathcal{U}, \mathcal{F}}, \mathcal{E}_{\mathcal{U}, \mathcal{F}})$ be the labeled multi-digraph with vertex set

and edge set

Furthermore, let $L\colon \mathcal{E}_{\mathcal{U}, \mathcal{F}}\to \{1, \ldots\}$ be the function that maps $\gamma\in \mathcal{E}_{\mathcal{U}, \mathcal{F}}$ to the corresponding index of the reaction label, that is, if $\gamma\in \mathcal{E}_{\mathcal{U}, \mathcal{F}}$ has label $r_i$, then $L(\gamma) = i$.

The construction of the graph $\mathcal{G}_{\mathcal{U}, \mathcal{F}}$ also makes sense in the case of an infinite RN. The vertex set is finite, but the number of edges between any two vertices might be infinite.

Example 3.4. Consider Example 2.2 with $\mathcal{U} = \mathcal{U}_1 = \{EA, EAB\}$ and $\mathcal{F} = \mathcal{R}_{\mathcal{U}}$. Then $\mathcal{G}_{\mathcal{U}, \mathcal{F}}$ is

There are infinitely many walks from $\ast_{in}$ to $\ast_{out}$, as one might take arbitrary many 'rounds' in the loop before exiting.

By definition, any finite sequence of reactions that corresponds to the reaction labels of a walk in $\mathcal{G}_{\mathcal{U}, \mathcal{F}}$ with start vertex $\ast_{in}$ and end vertex $\ast_{out}$ is a fast chain. Define the set of such walks by

Note that $W_{\mathcal{U}, \mathcal{F}}$ might be infinite, even for finite RNs, as in Example 3.4 above.

We will consider the situation in which the reactions of $\mathcal{F}$ occur fast compared to the remaining reactions $\mathcal{R}\setminus \mathcal{F}$, in the sense that the transition intensities of the reactions in $\mathcal{F}$ are scaled by $1/ \epsilon$ for small $\epsilon$. Thus, we consider a fast-slow dynamic regime. In that case, it is natural to expect that whenever a non-interacting molecule is created, then it will be degraded almost instantaneously through fast reactions, before any other non-fast reaction occurs. Such sequences of reactions (creation and degradation) are encoded in the walks of $W_{\mathcal{U}, \mathcal{F}}$. To understand the fast dynamics, it is therefore important to understand the net gain of core species in the walks and their probabilities of occurring.

In preparation for this, consider a walk $\Gamma = (\gamma_1, \gamma_2, \ldots, \gamma_l)\in W_{\mathcal{U}, \mathcal{F}}$, and denote

where $L(\gamma_j)$ is the reaction index of the edge $\gamma_j$. Define

Note that $\mathrm{r}(\Gamma)$ depends on the order of the elements of $\Gamma$.

Lemma 3.5. Let $\Gamma = (\gamma_1, \gamma_2, \ldots, \gamma_l)\in W_{\mathcal{U}, \mathcal{F}}.$ Then the following holds:

● $\mathrm{r}(\Gamma)\ge 0$, and $\mathrm{p}(\Gamma): = \mathrm{r}(\Gamma)+\sum_{i = 1}^ly'_{L(\gamma_i)}-y_{L(\gamma_i)}\ge 0$

● $p_ \mathcal{U}(\mathrm{r}(\Gamma)) = 0$, $p_ \mathcal{U}(\mathrm{p}(\Gamma)) = 0$.

Proof. We give complete proofs here for convenience, but note that the results also follow from ^[29].The second item follows by definition. As we take the maximum coordinate-wise and $w_1\ge 0$, then $\mathrm{r}(\Gamma)$ has non-negative coordinate in each species. To see that $\mathrm{p}(\Gamma)$ has non-negative coordinate in each species, we note that by definition $y_{L(\gamma_l)}+\sum_{i = 1}^{l-1}y_{L(\gamma_i)}-y'_{L(\gamma_i)}\leq \mathrm{r}(\Gamma)$. Adding $\sum_{i = 1}^ly_{L(\gamma_i)}'-y_{L(\gamma_i)}$ to both sides, we get $0\le y_{L(\gamma_l)}'\leq \mathrm{p}(\Gamma)$, and the result follows.

By Lemma 3.5, we might consider $\mathrm{r}(\Gamma)$ and $\mathrm{p}(\Gamma)$ as elements of $\mathbb{Z}_{\geq 0}^p$, and as reactant and product, respectively, of the reduced reaction $\mathrm{r}(\Gamma)\to \mathrm{p}(\Gamma)$, obtained by contraction along $\Gamma$. Lemma 3.7 below justifies this view, and relates the compatibility requirement in Assumption 1 to the reduced reactions.

Example 3.6. We continue with Example 3.4. Consider the walk $\Gamma_0$ of $\mathcal{G}_{ \mathcal{U}, \mathcal{F}}$ consisting of the two edges associated to $r_1, r_2$. Then

a trivial reaction. Likewise, the walk $\Gamma^A_1$ consisting of the edges associated to $r_1, r_3, r_4, r_2$, also results in a trivial reaction,

In contrast, the walk $\Gamma^B_1$ consisting of the three edges associated to $r_1, r_3$ and $r_5$ yields

Any other walk of $W_{\mathcal{U}, \mathcal{F}}$ consists either of a sequence of edges with labels $r_1, r_3, r_4, \ldots,$ $r_3, r_4, r_3, r_5$ or a sequence of edges with labels $r_1, r_3, r_4, \ldots, r_3, r_4, r_2$. As the net gain of core species in the contraction of the reactions $r_3, r_4$ is zero, the corresponding reduced reactions are the same as that of $\Gamma^A_1$, respectively, $\Gamma^B_1$. Hence, there is only one non-trivial reduced reaction. For future reference, denote the corresponding walks with $n$ instances of $r_3$ by $\Gamma^A_n$, respectively, $\Gamma^B_n$, for $n\ge 1$.

Lemma 3.7. Suppose Assumption 1 holds. Let$\Gamma = (\gamma_1, \gamma_1, \ldots, \gamma_l)\in W_{\mathcal{U}, \mathcal{F}}$. Furthermore, let $z\in \mathbb{Z}_{\geq 0}^p$ and $x = (z, 0)\in \mathbb{Z}^n_{\geq0}$. Then

if and only if $z\geq \mathrm{r}(\Gamma)$.

The lemma implies that the reactions corresponding to a walk in $W_{\mathcal{U}, \mathcal{F}}$ can fire in succession of each other (without other non-fast reactions firing in between), if and only if the present molecule counts of core species is larger or equal to $\mathrm{r}(\Gamma)$. A similar statement appears in [29,Corollary 3.2].

Lemma 3.7 allows us to define transition intensities of the reduced reactions in a natural way. For $\Gamma = (\gamma_1, \gamma_2, \ldots, \gamma_l)\in W_{\mathcal{U}, \mathcal{F}}$, define the function $\Lambda_{ \Gamma}\colon \mathbb{Z}_{\geq 0}^p\to \mathbb{R}_{\geq0}$,

where $\xi_{i}: = y_{L(\gamma_i)}'-y_{L(\gamma_i)}$, and ${\rm{out}}(s(\gamma_{j}))$ denotes the set of outgoing edges of $s(\gamma_{j})$ in $\mathcal{G}_{\mathcal{U}, \mathcal{F}}$. Each term in the product is the probability that the desired reaction is chosen out of all possible fast reactions with the same non-interacting species in the reactant. The convention $0/0 = 0$ is used. Even if the RN is infinite, then (3.1) is well-defined due to (2.1).

Definition 3.8. Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be an RN on $\mathcal{S}$, let $\mathcal{U}\subseteq \mathcal{S}$ be a non-interacting set of species, and $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ a set of proper fast reactions. The reduced RN $\mathcal{N}_{\mathcal{U}, \mathcal{F}} = (\mathcal{C}_{\mathcal{U}, \mathcal{F}}, \mathcal{R}_{\mathcal{U}, \mathcal{F}})$ on $\mathcal{O} = \mathcal{S}\setminus \mathcal{U}$, obtained by elimination of $(\mathcal{U}, \mathcal{F})$ from $\mathcal{N}$, is the possibly infinite RN defined by

and $\mathcal{C}_{ \mathcal{U}, \mathcal{F}}$ being the set of vertices of $\mathcal{R}_{ \mathcal{U}, \mathcal{F}}$.

For $r\in \mathcal{R}_{ \mathcal{U}, \mathcal{F}}$, define

Definition 3.9. Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be an SRN on $\mathcal{S}$ with transition intensities $\lambda_r(\cdot), r\in \mathcal{R}$, satisfying Assumption 1. Let $\mathcal{U}\subseteq \mathcal{S}$ be a set of non-interacting species, and $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ a set of proper fast reactions. Then, the possibly infinite SRN $\mathcal{N}_{\mathcal{U}, \mathcal{F}} = (\mathcal{C}_{\mathcal{U}, \mathcal{F}}, \mathcal{R}_{\mathcal{U}, \mathcal{F}})$ on $\mathcal{O} = \mathcal{S}\setminus \mathcal{U}$ with transition intensities $\tau_r(\cdot), r\in \mathcal{R}_{\mathcal{U}, \mathcal{F}}$, given by

is the reduced SRN, obtained by elimination of $(\mathcal{U}, \mathcal{F})$ from $\mathcal{N}$.

As a reduced SRN has species set $\mathcal{O} = \mathcal{S}\setminus \mathcal{U}$, the associated CTMC lives in $\mathbb{Z}_{\geq 0}^p$. Even if the reduced SRN has infinitely many reactions, the associated CTMC is always well-defined, that is, it has no instantaneous states ^[28]. If $\mathcal{U}$ consists of intermediate species, then the reduced reaction network and the corresponding transition intensities are easily found [20,Theorem 3.1]. For intermediate species, the number of molecules of core species does not change along a walk, in contrast to for instance Example 3.6. In particular, if the original SRN has mass-action kinetics, then the reduced SRN has mass-action kinetics as well.

As remarked earlier, one might discard any trivial reaction. However, the transition intensities of the trivial reactions play a crucial role in Assumption 2 in Section 4.1.

Lemma 3.10. Suppose Assumption 1 holds. Then, the reduced SRN has Q-matrix with all off-diagonal row sums finite.

Proof. By definition, we need to show that

For this, we note that

Hence, it is enough to show that the second summand in the second term is finite as the first is finite by assumption. By construction of the reduced transition intensities,

The left hand side is the sum of all transition intensities of walks (reduced reactions) with initial reaction in $\mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}}$. This sum is at most equal to the right hand side: Let $\Gamma = (\gamma_1, \ldots, \gamma_l)$ be a walk and $y_{L(\gamma_1)}$ be the reactant of the first reaction in the walk. If $y_{L(\gamma_1)}\le z\le \mathrm{r}(\Gamma)$ and $\mathrm{r}(\Gamma)\not = z$, then the first reaction can fire, but not all reactions in the walk can fire in succession, hence the corresponding $\Lambda_{ \Gamma}(z)$ is zero. For a concrete illustration, see Example 5.2.

As a consequence of Lemma 3.7, Assumption 1 holds for a reduced SRN, provided it holds for the original SRN.

Corollary 3.11. Suppose Assumption 1 holds. Then, for any reaction $r:y\to y{^{\prime}}\in \mathcal{R}_{\mathcal{U}, \mathcal{F}}$ of the reduced SRN, it holds that $\tau_{r}(z) > 0$ if and only if $z\geq y$.

Proof. It is enough to prove it for $\Lambda_{ \Gamma}(z),$ where $\Gamma\in W_{\mathcal{U}, \mathcal{F}}$ and $z\in \mathbb{Z}_{\geq 0}^p$, that is,

Only one direction is not obvious. Assume $z\geq\mathrm{r}(\Gamma)$. Using the definition of $\Lambda_{ \Gamma}(\cdot)$, it is sufficient to show that the numerators in the fraction are all non-zero. This holds by Lemma 3.7, hence $\Lambda_{ \Gamma}(z) > 0$.

Example 3.12. We continue with Examples 3.4, 3.6, assuming stochastic mass-action kinetics. There are three reduced reactions, the first reaction being non-trivial and the last two trivial,

with infinitely many walks underlying the first and third reduced reaction. Here, $s_i$ is used to denote the reactions of the reduced SRN, rather than $r_i$, to distinguish the reactions of the reduced RN from those of the original RN. We find

such that

Analogously, we compute the sums over walks for the trivial reactions. We have $\tau_2(z) = \Lambda_{ \Gamma_0}(z)$ and

Furthermore, a simple calculation gives the following identity for $z\in \mathbb{Z}_{\geq 0}^p$,

An interesting observation is that the kinetics is mass-action-like in the sense that the numerator alone is of mass-action form, while the denominator is positive (for any state). This is, however, not true in general, see Example 5.2.

Example 3.13. Consider the RN,

with $\mathcal{U} = \{U_1\}$ and $\mathcal{F} = \mathcal{R}_{\mathcal{U}}$ proper. Taken with mass-action kinetics with corresponding rate constants, $\kappa_1, \kappa_2, \kappa_3$, the reduced SRN has infinitely many reactions, $s_n\colon S_1+nS_2\to S_3$, $n\ge0$, with transition intensities,

where by convention $\prod_{i = 0}^{-1} = 1$. On any particular state, at most finitely many reactions can be active, that is, have non-zero transition intensity.

For examples from the biochemical literature, we refer to $\S$ 5.1.

4.   Two-scale SRNs

We study transient approximability of the CTMC via the reduced SRNs under appropriate assumptions. We cover two settings, the standard time-homogeneous CTMC setting for SRNs, and afterwards the setting where the transition intensities of the SRN are allowed to be time-dependent.

4.1.   Transient approximation

Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be an SRN on a species set $\mathcal{S}$ (of size $n$) with non-interacting species $\mathcal{U}\subseteq \mathcal{S}$, fast proper reactions $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$, and transition intensities $\lambda_r, r\in \mathcal{R}$. From this SRN, we construct a family of SNRs, indexed by a parameter $\epsilon > 0$, with corresponding $Q$-matrix $Q^\epsilon$. In particular, the transition intensities of the reactions in $\mathcal{F}$ consuming non-interacting species are scaled in $\epsilon$ as follows:

● $\lambda^ \epsilon_{r}(x): = \frac{1}{ \epsilon} \lambda_{r}(x)$, $r\in \mathcal{F}$, that is, for small $\epsilon$ the reactions are fast compared to the remaining reactions,

● $\lambda^ \epsilon_{r}(x): = \lambda_{r}(x)$, $r\in \mathcal{R}\setminus \mathcal{F}$, that is, the transition intensities are independent of $\epsilon$, and the reactions are slow.

Thus, we might write the corresponding $Q$-matrix as as sum of two terms, a fast part $\widetilde{Q}$, scaled by $\frac{1}{ \epsilon}$, and a slow part $\widehat{Q}$,

It follows from Definition 3.9 that the reduced SRN, obtained from the SRN with transition intensities $\lambda^ \epsilon_{r}$, $r\in \mathcal{R}$, has transition intensities $\tau_r$, $r\in \mathcal{R}_{ \mathcal{U}, \mathcal{F}}$, independent of $\epsilon$.

Let $(X_ \epsilon(t))_{t\ge 0}$ be Markov chains on $\mathbb{Z}_{\ge0}^n$ with generators $Q^ \epsilon$, $\epsilon > 0$, respectively. We note that all $Q^ \epsilon$, $\epsilon > 0$, are dynamically equivalent in the sense that the connectivity of the state space and the decomposition of the state space $\mathbb{Z}^n_{\ge0}$ into communicating classes are independent of $\epsilon$. Furthermore, let $(X_0(t))_{t\ge0}$ be a Markov chain on $\mathbb{Z}_{\ge0}^p$ (excluding the $m = n-p$ coordinates for non-interacting species) with generator $Q^0$, obtained from the transition intensities $\tau_r$, $r\in \mathcal{R}_{ \mathcal{U}, \mathcal{F}}$. In the following, we provide conditions that guarantee that the dynamics of $(X_ \epsilon(t))_{t\ge 0}$ is similar (in a sense to be made precise) to the dynamics of $(X_0(t))_{t\ge0}$, whenever $\epsilon$ is small. Technically, we will consider the limit as $\epsilon\to0$.

We restrict attention to a particular closed set $E$ of $(X_ \epsilon(t))_{t\ge 0}$, such that $E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\})\not = \emptyset$. That is, there are states in $E$ with no molecules of non-interaction species being present. If the reduced SRN starts in $X_0(0)\in E_0: = \rho_ \mathcal{O}(E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\}))$, then by construction of the reduced reactions, $X_0(t)\in E_0$ for all $t > 0$. Furthermore, $E_0$ is a closed set of $Q^0$.

We require the following.

Assumption 2. The closed set $E$ is finite, and all $z\in E_0$ satisfy the following:

The sum on the left hand side also includes walks giving rise to trivial reactions, as in Example 3.12. The condition might fail in two ways. Either a walk is blocked because of lack of molecules of core species (for example $r_1\colon0\to U,$ $r_2\colon S+U\to 0$; if there is no molecules of $S$, then the fast chain $r_1, r_2$ is blocked and a molecule of $U$ cannot be degraded), or one might be trapped in infinite walks (for example, $r_1\colon0\to U$, $r_2\colon2S+U\to 3S+U$, $r_3\colon U\to 0$; there is positive probability of an infinite walk $r_1, r_2, r_2, \ldots$ without degradation of the $U$ molecule). The latter cannot occur if $E$, and hence $E_0$, are finite as this implies that at most finitely many (reduced) reactions can be active on any state of $E$ ($E_0$).

For a sub-conservative RN, any closed set of states $E$ (not necessarily a communicating class) in a directed stoichiometric compatibility class is finite. If Assumption 1 is satisfied, then the conditions in Assumption 2 do not depend on the transition intensities $\lambda_r$, $r\in \mathcal{R}$, but only on the structure of the underlying RN (together with $\mathcal{U}$ and $\mathcal{F}$). Furthermore, as (2.1) holds by assumption, then Assumption 2 is meaningful even for infinite RNs. In particular, Theorem 4.1 also holds for infinite RNs as $E$ is finite.

For a subset $B\subseteq E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\})$, we define $B_0 = \rho_ \mathcal{O}(B)$.

Theorem 4.1. Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be an SRN on $\mathcal{S}$ with transition intensities $\lambda_r$, $r\in \mathcal{R}$, $\mathcal{U}\subseteq \mathcal{S}$ a set of non-interacting species, and $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ a set of fast proper reactions. Suppose Assumptions 1 and 2 hold. Let $\pi$ be a probability distribution on $E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\})$, and $\pi_0$ the induced probability distribution on $E_0$ by omitting the last $n-p$ coordinates of the states of $E$. Then, the following holds for any $0 < T$:

In particular, for any $B\subseteq E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\})$ and any $0 < T$:

The proof is in Appendix $\S$ B.3. As a consequence of the theorem we have:

Corollary 4.2. Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be an SRN on $\mathcal{S}$ with transition intensities $\lambda_r$, $r\in \mathcal{R}$, $\mathcal{U}\subseteq \mathcal{S}$ a set of non-interacting species, and $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ a set of fast proper reactions. Suppose Assumptions 1 and 2 hold. Let $x = (z, 0)\in E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\})$ and $B\subseteq E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\})$. Then, for all $t\geq0$ it holds that:

If the initial state $X_ \epsilon(0)$ has more than one molecule of the non-interacting species, then these will in general be depleted quickly for small $\epsilon$. This is however not always the case, see Example 5.2.

4.2.   The case of time-dependent transition intensities

In the previous section, we considered the scaling limit of time-homogeneous SRNs. While time-constant transition intensities are reasonable in many situations, time-dependence becomes relevant, for example, in connection with variation in experimental setups, cycle-dependent mechanisms or temperature changes ^[30].

The setup is essentially the same as in Section 4.1, except we now allow time-dependent transition intensities, $\lambda_r(t, x)$, $r\in \mathcal{R}$, for $t$ in a compact time interval $t\in [0, T]$. We make the following assumption.

Assumption 3. The time-dependent transition intensities, $\lambda_r(t, x)$, $r\in \mathcal{R}$, are such that $\lambda_r(t, \cdot)$ satisfies Assumption 1 for all $t\in [0, T]$ and such that for all $x\in \mathbb{Z}^n_{\geq0}$, $\lambda_r(\cdot, x)$, $r\in \mathcal{R}$, are $C^1$-functions from $[0, T]\to \mathbb{R}_{\geq 0}$ with Lipschitz derivatives.

The transition intensities are scaled as in the homogeneous case. Under Assumption 3, the reachability properties of the CTMC are invariant over time and match the reachability properties of the homogeneous SRNs under Assumption 1. Assumption 2 holds for all $t\in[0, T]$, provided it holds for one $t$, as the assumption is structural.

With these changes and amendments, then Theorem 4.1 and Corollary 4.2 hold as well for the time-dependent case. A proof is sketched in Appendix C.

5.   Examples and properties

We next illustrate the reduction procedure with various examples. Furthermore, we derive sufficient conditions for Assumption 2 to hold, and compare the long-term behavior of the CTMC of the reduced SRN with the one of the original SRN.

5.1.   Zoo of examples

In this section, we give realistic examples to elaborate on the reduction. All examples are taken with stochastic mass-action kinetics, hence they satisfy Assumption 1. To indicate mass-action kinetics, we put the reaction rate constants $\kappa_i$ as labels of the reactions, and omit the reaction names $r_i$. A key source book for realistic examples is ^[26], that contains reaction networks used in systems biology. Other useful references are ^[31], which focuses on enzyme kinetics, and ^[1], which has examples from a range of areas in biology.

Example 5.1. We first consider the two-substrate Michaelis-Menten mechanism [26,Section 3.1.2]; also discussed in Examples 2.2, 3.4, 3.6, and 3.12, and repeated here for convenience with reaction names replaced by their reaction constants:

The mechanism is commonly studied as a deterministic mass-action system, using Tikhonov-Fenichel singular perturbation theory or the QSSA, with short-lived intermediate complexes (non-interacting species) $\mathcal{U} = \{EA, EAB\}$ and fast reactions $\mathcal{F} = \mathcal{R}_ \mathcal{U} = \{r_2, r_3, r_4, r_5\}$. Assuming the same in the stochastic setting gives the reduced RN in Example 3.6{or Example 3.12 with} one non-trivial reaction,

with transition intensity

We omit the trivial reactions; further details are given in example 3.12. Assumption 2 is fulfilled for all $z\in \mathbb{Z}_{\geq 0}^5$ ($p = 5$, see Example 3.12), in particular the original RN is conservative. Then, Theorem 4.1 applies to any closed set $E$ in a directed stoichiometric compatibility class that intersects $\mathbb{Z}_{\geq 0}^5\times \{0\}$ non-trivially.

Example 5.2. Consider the two-substrate Michaelis-Menten mechanism with $\mathcal{U} = \{EA, EAB\}$ as above, but with set of fast reactions $\mathcal{F} = \{r_3, r_4, r_5\}\subseteq \mathcal{R}_ \mathcal{U}$. We remove $r_2$ from the fast reactions in the previous example. Then there is only one reduced reaction in contrast to Example 5.1 that additionally had two trivial reactions, given by

with transition intensity

Assumption 2 holds for all $z = (z_A, z_B, z_E, z_P, z_Q)\in \mathbb{Z}_{\ge 0}^5$ with $z_B > 0$:

Also note that for $z_A\geq1, z_E\geq 1$ and $z_B = 0$, we have $0 = \tau_{1}(z) < \kappa_1z_Ez_A = \lambda_1(z, 0)\neq 0$. This is an example of the inequality in the proof of Lemma 3.10.

This example does not have mass-action-like kinetics, as $\tau_1(z)$ cannot be expressed as a fraction with positive denominator such that the numerator is of mass-action form. Alternatively, we might extend by multiplication and division by $z_B$.

Example 5.3. Consider a non-competitive inhibition network [26,Section 3.2.2]:

where $E$ is an enzyme that catalyze the conversion of a substrate $S$ into a substrate $P$. The conversion is delayed by an inhibitor $I$ that binds to the enzyme and to the substrate through the intermediate complex $EI$. It is non-competitive in the sense that $I$ does not compete with $E$ for substrate binding, but rather acts on $E$ directly; a phenomenon known as allosteric regulation.

This example is typically analyzed by means of the QSSA with short-lived species; here $\mathcal{U} = \{ES, EI, ESI\}$. In our setting, taking $\mathcal{U}$ to be non-interacting species with $\mathcal{F} = \mathcal{R}_{\mathcal{U}}$, then the reduced SRN becomes

with transition intensities

\medskip For any closed set $E$ in a directed stoichiometric compatibility class that intersects $\mathbb{Z}_{\geq 0}^4\times \{0\}$ non-trivially, Assumption 2 is satisfied.

Example 5.4. Consider an example of allosteric activation [26,Problem 3.7.8] that models the regulation of an enzyme by binding of an allosteric activator before the enzyme can bind a substrate:

Here, $R$ is the allosteric activator, $S$ the substrate, $E$ the enzyme, and $P$ the product. While the RN is reminiscent of the the two-substrate example 5.1, the enzyme-activator complex $ER$ stays intact after the product $P$ dissociates in reaction $r_5$.

To analyze the behavior of the system, often the QSSA is applied to a mass-action ODE system with short-lived intermediate complexes (non-interacting species) $\mathcal{U} = \{ER, ERS\}$. Choosing analogously $\mathcal{F} = \mathcal{R}_{\mathcal{U}}$ in the stochastic case gives a proper set of fast reactions. Then, the reduced RN has infinitely many reactions, given by:

for $k\ge 1$, with transition intensities

All $z\in \mathbb{Z}_{\geq 0}^4$ ($p = 4$) satisfy (4.1). As $\mathcal{N}$ is sub-conservative, any closed set $E$ of a directed stoichiometric compatibility class is finite. Hence, for any closed set $E$ in a directed stoichiometric compatibility class that intersects $\mathbb{Z}_{\geq 0}^4\times \{0\}$ non-trivially, Assumption 2 is satisfied.

Example 5.5. As a final example, consider a mechanism-based inhibitor system [1,Section 6.4],

where $S, P$ are substrate and product, respectively, $X, Y$ intermediate complexes, $E$ an enzyme in active form (implying it might bind to the substrate) and $E_i$, the enzyme in its inactivated form. A main interest is to know the final ratio of the product to the inactivated enzyme [1,Section 6.4]. For this, often the QSSA is applied to the short-lived species $\mathcal{U} = \{X, Y\}$. Analogously, we consider the stochastic system with $\mathcal{F} = \mathcal{R}_{\mathcal{U}}$ being fast and proper, and study the corresponding stochastic reduction. The reduced SRN has reactions

with transition intensities:

As $\mathcal{U} = \{X, Y\}$ are intermediate species, all $z\in \mathbb{Z}_{\geq 0}^4$ ($p = 4$) satisfy (4.1). As $\mathcal{N}$ is sub-conservative, for any closed sets $E$ of a directed stoichiometric compatibility class Assumption 2 is satisfied.

5.2.   Properties enabling sufficient conditions for simplification

While Assumption 1 and $\mathcal{F}$ being proper are both easy to check, Assumption 2 is non-trivial in general. Assumption 2 requires a finite closed set of the CTMC and that molecules of non-interacting species can be degraded through chains of fast reactions. The following gives sufficient conditions for this to hold (see Appendix $\S$ A for proof).

Proposition 5.6. Assume $\mathcal{N}$ is a sub-conservative SRN on a species set $\mathcal{S}$, $\mathcal{U}\subseteq \mathcal{S}$ is a set of non-interacting species, and $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ is a proper set of reactions. Suppose Assumption 1 holds. Furthermore, assume that one of the following holds:

$(a)$ $\mathcal{F}\cup (\mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}})$ is weakly reversible.

$(b)$ $\mathcal{F}\cap \mathcal{R}_{\mathcal{U}}\cap \mathcal{R}_{\mathcal{U}}'$ is weakly reversible, and for any reaction $r\colon y\to y'$ in $\in \mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}}$, the reverse reaction $r'\colon y'\to y$ is in $\mathcal{F}$.

$(c)$ $\mathcal{N}$ is weakly reversible and $\mathcal{F} = \mathcal{R}_{\mathcal{U}}$.

Then, for an arbitrary closed set in a (directed) stoichiometric compatibility class, Assumption 2 is satisfied.

None of the conditions in the theorem are necessary for Assumption 2 to hold. For intermediate species a stronger result holds.

Lemma 5.7. Assume $\mathcal{N}$ is an SRN on a species set $\mathcal{S}$, and that $\mathcal{U}\subseteq \mathcal{S}$ is a set of intermediate species. Suppose Assumption 1 holds. Then, $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ is proper if and only if all $z\in \mathbb{Z}_{\geq 0}^p$ satisfy (4.1) of Assumption 2.

We summarize the properties of the examples from $\S$ 5.1. In examples 5.1, 5.3 and 5.4 we can directly conclude by Proposition 5.6 (b) that Assumption 2 holds.

5.3.   Comparison of the original RN and the reduced RN

A natural question following the transient approximability of SRNs with non-interacting species is whether the reduced SRN approximates the stationary behavior in the limit as $\epsilon\to 0$. First one might ask whether states $x = (z, 0)\in \mathbb{Z}_{\geq 0}^p\times \{0\}$ and $z$ are of the same type (positive recurrent/null recurrent/transient) for $\mathcal{N}$ and $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$, respectively. Furthermore, the stationary distribution (if it exists) for the scaled SRN $\mathcal{N}$ as $\epsilon\to 0$ might match the stationary distribution of $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$. The following examples show that such correspondences do not hold in general, even for reduction by intermediates.

Example 5.8. Consider the following SRN $\mathcal{N}$ with mass-action kinetics,

Choosing $\mathcal{U} = \{U_1\}$, $\mathcal{F} = \{r_3\}$, then the reduced SRN $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$ is $S_1\underset{ }{\overset{ }{\mathop{\rightleftharpoons }}}\, S_2$ with transition intensities satisfying Assumption 1 by Corollary 3.11. Any $(a, b, c, 0)$ with $a+b\geq 1$ is transient for $\mathcal{N}$, but recurrent for $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$. Reaction $r_5$ in $\mathcal{N}$ causes absorption into $(0, 0, a+b+c+d, 0)$. As $r_5$ is not present in the reduced SRN, this cannot happen in $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$.

Example 5.9. Consider the SRN $\mathcal{N}$ with mass-action kinetics,

Choosing $\mathcal{U} = \{U_1\}$, $\mathcal{F} = \{r_3\}$, then the reduced SRN $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$ becomes $S_1\to S_2$. The state $(a, b, 0)$ with $a\geq 1$ is recurrent for $\mathcal{N}$, but transient for $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$.

The next result follows from [29,Theorem 5.6]. We recall that irreducible components of sub-conservative RNs are finite ^[27]. For an irreducible component $E\subseteq \mathbb{Z}^n_{\ge 0}$, define $E_0 = \rho_{ \mathcal{O}}(E\cap (\mathbb{Z}_{\geq 0}^p\times \{0\}))$.

Theorem 5.10. Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be a sub-conservative SRN on $\mathcal{S}$, $\mathcal{U}\subseteq \mathcal{S}$ a set of intermediate species, and $\mathcal{F} = \mathcal{R}_ \mathcal{U}$ a set of proper reactions. Suppose Assumption 1 holds. Then $x = (z, 0)\in \mathbb{Z}_{\geq 0}^p\times \{0\}$ is transient for $\mathcal{N}$ if and only if $z$ is transient for $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$, and the same holds for positive recurrence. In particular, if a subset $E$ is an irreducible component for $\mathcal{N}$ then $E_0$ is an irreducible component for $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$.

From Theorem 5.10, Lemma 5.7 and [32,Theorem 3], we get convergence of the stationary distribution.

Corollary 5.11. Let $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ be a sub-conservative SRN on $\mathcal{S}$, $\mathcal{U}\subseteq \mathcal{S}$ a set of intermediate species, and $\mathcal{F} = \mathcal{R}_ \mathcal{U}$ a set of fast proper reactions. Suppose Assumption 1 holds. Let $E$ be an irreducible component, and let $E_0$ be the corresponding projected set. Furthermore, denote by $\pi_ \epsilon$ the unique stationary distribution of $Q_ \epsilon$ on $E$, and let $\pi_0$ be the unique stationary distribution of the reduced SRN on $E_0$. Then, $\pi_ \epsilon((z, 0))\to \pi_0(z)$, for $z\in E_0,$ and $\epsilon\to 0.$

6.   Relation to previous work and discussion

We compare our setting for reduction of SRNs with non-interacting species to other approaches of model simplification. Then we discuss assumptions, extensions and examples that are not covered by our approach.

In our treatment, we consider SRNs with fast reactions consuming non-interacting species, generalizing earlier settings in three ways: 1) extending from intermediate species to non-interacting species, 2) allowing arbitrary transition intensities satisfying a compatibility condition, and 3) allowing non-homogeneous Markov dynamics, that is, time-dependent transition intensities. As the set-up is based mostly on singular perturbation of Markov chains with transient states ^[33], which is also available in the book ^[22], we loose the possibility to study multiple time scales, but are restricted to two scales.

Reductions of SRNs with non-interacting species resembles reduction by intermediate species as intermediate species form a special class of non-interacting species ^[20]. In contrast to our setting, ^[20] allows a multi-scale setting. Other multi-scale reductions require the parameters to fulfill a balancing equation between scales of species concentrations and scales of transition intensities ^{[13,16,17,18,36]}, where ^[36] treats the classical scaling for convergence to the ODE. This is not required in our setting (nor in ^[20] for the reactions involving intermediate species), in fact such equation will not hold. Balancing is violated because fast reactions do not have 'limits' themselves as $\epsilon\to 0$, but only jointly through contraction of sequences of reactions. Furthermore, at the process-level, we do not have convergence in the Skorohod topology, see [20,Example 5.3] or [18,$\S$ 6.5], because even for small $\epsilon$, molecules of non-interacting species are created and exist for small amounts of time, while this is not possible in the reduced SRN. Convergence in Skorohod topology typically require that species abundance is measured in concentrations, rather than in molecule numbers. In contrast, we might have convergence in the Jakubowski topology from ^[34], also see [18,Remark A.14] for more in the context of reaction networks. Whenever, a molecule of a non-interacting species is produced, the Markov chain $(X_\epsilon(t))_{t\ge 0}$ will jump a finite number of times (stochastically bounded in $\epsilon$) in a small time interval before the molecule is degraded again. This implies the (so-called) up- and downcrossing numbers of ^[34] is bounded, a prerequisite for convergence.

The results on the transient dynamics in our two-scale setting could potentially be extended by allowing general infinite state space (that is, weaken Assumption 2). In so, the transition rates of the original CTMC become unbounded for mass-action kinetics, and the theory on singular perturbations for CTMCs are not applicable ^[22]. Furthermore, the original process might have sample paths that diverge to infinity in a finite amount of time. Hence, theory applicable to unbounded cases as well as conditions ensuring non-explositivity would be of interest to develop. We note that even 1-dimensional CTMCs arising from RNs can have both positive recurrent and explosive irreducible components ^[35], implying that different parts of the state space might have to be analyzed separately.

Another generalization would be to allow many scales rather than two scales. This could be done using ^[32] in the finite state space case. However, it becomes difficult to find the reduced reactions, as these depend on the concrete form of the transition intensities and their scalings. Also, a reduction might not be interpreted as SRNs fulfilling Assumption 1. As an example, consider

with stochastic mass-action kinetics. The reduction for small $\epsilon$ can be written as two reactions

with transition intensities $\kappa_1z_{S_1}{\mathbb{1}}_{\{0\}}(z_{S_3})$ for the first reaction and $\kappa_1z_{S_1}{\mathbb{1}}_{ \mathbb{N}\setminus\{0\}}(z_{S_3})$ for the second. In particular, the first transition intensity invalidates Assumption 1.

In our approach, it is important that non-interacting species are produced and degraded. Consider the SRN with reactions

Neither $U_1$ nor $U_2$ are produced nor degraded, hence there does not exist a proper set of reactions $\mathcal{F}$ and the example falls outside our setting. As a matter of fact, the scaling parameter $\epsilon$ would apply uniformly to all reactions as they all transform one non-interacting species into another. Hence, the distribution of the CTMC approaches the stationary distribution within a short time span (for small $\epsilon$). This distribution is not concentrated on the part of the state space without molecules of non-interacting species. Rescaling time by $\epsilon$ would retrieve the original chain.

A.   Complements for sufficient conditions for Assumption 2

The following lemma proves Proposition 5.6 on sufficient conditions for Assumption 2.

Lemma A.1. Assume $\mathcal{N} = (\mathcal{C}, \mathcal{R})$ is a sub-conservative SRN on a species set $\mathcal{S}$ with transition intensities $\lambda_r$, $r\in \mathcal{R}$, that $\mathcal{U}\subseteq \mathcal{S}$ is a set of non-interacting species, and $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$ is a proper set of fast reactions. Suppose Assumption 1 holds.Then either of the following are sufficient for Assumption 2 to be satisfied for an arbitrary closed set in a (directed) stoichiometric compatibility class:

$(a)$ $\mathcal{F}\cup (\mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}})$ is weakly reversible.

$(b)$ $\mathcal{F}\cap \mathcal{R}_{\mathcal{U}}\cap \mathcal{R}_{\mathcal{U}}'$ is weakly reversible, and for any reaction $r\in \mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}}$, the inverse reaction is in $\mathcal{F}$.

$(c)$ $\mathcal{N}$ is weakly reversible and $\mathcal{F} = \mathcal{R}_{\mathcal{U}}$.

Proof. As $\mathcal{N}$ is sub-conservative, any closed set of a directed stoichiometric compatibility class is finite. Hence, it is enough to show that (4.1) holds for arbitrary $x = (z, 0)\in \mathbb{Z}_{\geq0}^p\times \{0\}$. As (2.1) holds by assumption, it is enough to show that for arbitrary $r\in \mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}}$ with $\lambda_{r}(z, 0)\neq 0$ the following holds:

Consider the SRN with species set $\mathcal{S}$ and reactions $\mathcal{F}$, denoted $\mathcal{N}|_ \mathcal{F}$, that we take with the same transition intensities as $\mathcal{N}$. Consider the associated discrete-time Markov chain (DTMC) from the jump chain of $\mathcal{N}|_ \mathcal{F}$ (cf. [28,p. 82]), which is the sequence of states taken by the CTMC. This will be denoted by $(Y_n)_{n\in \mathbb{N}}$. By definition of $\Lambda_{ \Gamma}(\cdot)$, $\frac{\Lambda_{ \Gamma}(z)}{ \lambda_{r}(z, 0)}$ corresponds to the probability that the transitions of the jump chain $(Y_n)_{n\in \mathbb{N}}$ when starting from $(z, 0)+y'-y$ comes from the sequence of reactions in $\Gamma$.

Denote the set of reachable states from $(z, 0)+y'-y$ via $\mathcal{N}|_ \mathcal{F}$ by $\mathcal{N}|_ \mathcal{F}((z, 0)+y'-y)$. As $\mathcal{N}$ is sub-conservative, $\mathcal{N}|_ \mathcal{F}$ is sub-conservative and $\mathcal{N}|_ \mathcal{F}((z, 0)+y'-y)$ is finite. In the following, let $e_i\in \mathbb{Z}^n_{\geq 0}$ denote the vector with $U_i$-coordinate entry equal to one and all other entries equal to zero.

(a) Assume $\mathcal{F}\cup (\mathcal{R}_{\mathcal{U}}'\setminus \mathcal{R}_{\mathcal{U}})$ is weakly reversible. Then, if $(z{^{\prime}}, 0)+e_i\in \mathcal{N}|_ \mathcal{F}((z, 0)+y'-y)$, then also $(z, 0)\in \mathcal{N}|_ \mathcal{F}((z{^{\prime}}, 0)+e_i)$. Hence, any such $(z{^{\prime}}, 0)+e_i$ is transient in the jump chain $(Y_n)_{n\in \mathbb{N}}$ of $\mathcal{N}|_ \mathcal{F}$. On the other hand by definition, any state of the form $(z', 0)$ is absorbing for the jump chain of $\mathcal{N}|_ \mathcal{F}$. As the set of reachable states is finite, the time until absorption is a.s. finite. So overall the jump chain ends in one of the absorbing states after a finite time, and

For (b) and (c), the proofs are similar to (a).

B.   Proof of Theorem 4.1

The proof of Theorem 4.1 is based on singularly perturbed and watched Markov chains. We therefore first introduce these concepts for the convenience of the reader.

B.1.   Singularly perturbed CTMCs

Singular perturbation theory for Markov chains is a well-developed topic in probability theory. For the convenience of the reader we restate and summarize results of ^[22,33] for homogeneous CTMCs on a finite state space, using their notation. For more on the two-scale setting for CTMCs, we refer to ^[22,37], where ^[37] is a book more generally treating multiscale methods in mathematics.

The setting we introduce has a Markov chain with fast and slow components subject to so-called weak and strong interactions, where the fast components consist of absorbing, weakly irreducible or transient states of the fast dynamics. These cases are treated in detail in [22,$\S$ 4.3,$\S$ 4.4,$\S$ 4.5].

Let a sequence of homogeneous generators of a CTMC be given by

on a finite state space $\mathcal{I}$. By rearrangement, the states can be divided into blocks of states $\mathcal{I}_A, \mathcal{I}_W, \mathcal{I}_T$ according to their role in the Markov chain associated to $\widetilde{Q}$, which are absorbing, irreducible, and transient, respectively. The matrices $\widetilde{Q}, \widehat{Q}$ are partitioned accordingly into block-matrices and take the form

The states in $\mathcal{I}_A$ with $\widetilde{Q}_{AA} = 0, \widetilde{Q}_{AW} = 0, \widetilde{Q}_{AT} = 0$ correspond to absorbing states of $\widetilde{Q}$. The states in $\mathcal{I}_T$ correspond to transient states of $\widetilde{Q}$ (hence $\widetilde{Q}_{TT}$ is a Hurwitz matrix, furthermore stable and non-singular see Lemma B.5). Finally, $\mathcal{I}_W$ consists of the states that are part of (non-trivial) closed irreducible components, so $\widetilde{Q}_{WW}$ in (B.1) can further be decomposed as

where $\widetilde{Q}^i$ is an $m_i\times m_i$ matrix. We denote by $\nu^i$ the stationary distribution of $\widetilde{Q}^i$ for $i = 1, \ldots, l$ (which exist by assumptions on the state space).

Remark B.1. Our presentation differs slightly from ^[22]. In that work, so-called weakly irreducible classes is used [22,Definition 2.7a]; classes containing exactly one closed irreducible component and possibly other states that lead to this component (which are then transient) in the case of time-homogeneous CTMCs. As we include transient states of $\widetilde{Q}$ via $\mathcal{I}_T$, this is not a restriction. We note that our setting is included in [22,$\S$ 4.5].

For $\epsilon > 0$, the forward equation of the CTMC gives an ODE

where $\pi_0$ is the initial probability distribution.

In ^[22], they construct sequences of functions that approximate $p^ \epsilon(t)$ in the limit for $\epsilon\to 0$ uniformly on $[0, T]$ for different powers of $\epsilon$. One contribution comes from the so-called outer expansion, that approximates $p^ \epsilon(t)$ for $t > 0$. Another part comes from the so-called initial-layer correction, that approximates $p^ \epsilon(t)$ in a neighborhood of $t = 0$. This might be ignored as it is zero in our case ^[22]. Hence, we will focus on the zeroth order outer expansion as this is enough to obtain an approximation to order $O(\epsilon)$ (see Proposition B.3).

Following [22,Theorem 4.45], the distribution of the CTMC converges to the zeroth order outer expansion $\varphi_0(t)$ for $t > 0$. In the case we consider, it is given as

where $\vartheta^i(t)$ are scalar functions and $\nu^i$ are the stationary distributions from matrix (B.2), and $0_T$ is the vector of zero-entries. Then, let (with $m_0 = 0$)

where

and define the matrix

Note that $\overline{Q}_\ast$ is not the same as $Q^0$ (which was previously defined for SRNs).

Remark B.2. The matrices in the product above have different dimensions, i.e. ${\rm{diag}}(1, \cdots, 1, \nu^1, \cdots, \nu^l)$ is a $(| \mathcal{I}_A|+l)\times (| \mathcal{I}_A|+| \mathcal{I}_W|)$ matrix, the middle matrix is a $(| \mathcal{I}_A|+| \mathcal{I}_W|)\times (| \mathcal{I}_A|+| \mathcal{I}_W|)$ matrix whereas ${\rm{diag}}(1, \cdots, 1, \widetilde{{\mathbb{1}}})$ is a $(| \mathcal{I}_A|+| \mathcal{I}_W|)\times (| \mathcal{I}_A|+l)$ matrix.

Then the zeroth order outer expansion $\varphi_0(t)$ (corresponding to [22,(4.86)]) with initial distribution $\pi_0 = (p_A, p_W, p_T)$ is determined by the following ODE,

which defines the generator of a CTMC that is given by $\overline{Q}_*$. This CTMC has reduced state space $\mathcal{I}_A\cup\{s_{1}, \cdots, s_{l}\}$, where each state $s_i$ corresponds to the assembled irreducible component, which has stationary distribution $\nu^i$ from $\widetilde{Q}_{WW}$, see (B.2). Then, denoting the projection onto the reduced state space of the initial value of (B.4) by ${\rm{Pr}}(\pi_0): = (\varphi_A(0), \vartheta^1(0), \cdots, \vartheta^l(0))$, we get a solution

Denote the sequence of CTMCs associated to $Q^ \epsilon$ by $X_ \epsilon(t)$ and the CTMC from $\overline{Q}_\ast$ by $X_0(t)$. Under these conditions the following holds [22,Theorem 4.45].

Proposition B.3. Let $\pi_0$ be a probability distribution on $\mathcal{I}_A\cup \mathcal{I}_W\cup \mathcal{I}_T$ with support on $\mathcal{I}_A$. Then, for all $T > 0$ and all $B\subseteq \mathcal{I}_A$:

The contribution of the zeroth order outer expansion is sufficient to approximate the above events to order $O(\epsilon)$ as the initial-layer corrections are zero ^[22]. Hence, the CTMC generated by $\overline{Q}_*$ can be seen as the limit CTMC (in the sense of Proposition B.3). The above will be the main reference of this section for the proof of Theorem 4.1.

Remark B.4. In the inhomogeneous case, the generator of the scaled CTMC of ^[22] has the form $Q^ \epsilon_t = \frac{1}{ \epsilon}\widetilde{Q}_t+\widehat{Q}_t$, where $\widetilde{Q}_t, \widehat{Q}_t$ satisfy some regularity conditions on $t\in [0, T]$. The conclusions we have stated are essentially the same in this case. In particular, the error $O(\epsilon)$ for the zeroth order outer expansion is as in Proposition B.3.

B.2.   Watched CTMCs

We introduce watched CTMCs on a finite state space. Watched Markov chains appear as a restriction of a bigger Markov chain, which we observe only if the chain is in a specific subset of the state space ^[38,39]. For a textbook treatment see ^[38], and for a review and connections to linear algebra we refer to ^[39].

In the following, we consider a CTMC $(X(t))_{t\ge 0}$ with generator $Q$ on a finite state space $S$ with $E_1\subseteq S$ a subset of the state space. Writing the $Q$-matrix in block-form according to the sets $E_1$ and $E_2 = S\setminus E_1$, we decompose $Q$ as

Lemma B.5. If under

all states in $E_1$ are transient, then $Q_{E_1E_1}$ is stable and non-singular.

Proof. Note that no diagonal entries of $Q_{E_1E_1}$ are zero. Since all states are transient, then all communicating component of the matrix have a "leak" (the row sum of $Q_{E_1E_1}$for a least one state in a component is non-zero) and the result follows from ^[40].

The censored Markov chain on $E_2$ is defined as the process that has sample paths following the process $(X(t))_{t\ge 0}$ as long as it is in $E_2$, and ignoring the parts where $X(t)\in E_1$. This so-defined process is then again a Markov process, where under the assumption above, the transition rates of the watched Markov chain are given by

Remark B.6. In the above formula, the second summand catches the contribution from when the original CTMC enters $E_1$ until it returns to $E_2$. By definition of the watched CTMC, we have

where $\tau = \inf\{t\geq0|X(t)\in E_2\}$.

B.3.   Proof of Theorem 4.1

The proof consists of several parts. We first show without loss of generality, the following might be assumed:

$-$ the closed set $E$ has a particular form,

$-$ the initial distribution has support on a special part of the state space,

$-$ the RN $\mathcal{N}$ is finite.

Then we apply singularly perturbed CTMCs of $\S$ B.1 to the SRN. As a third step, we simplify by showing, the dynamics can be restricted to a smaller part of the state space. Finally, we connect the limit CTMC to watched CTMCs of $\S$ B.2, which shows how the reduced reactions are derived from the limit CTMC of the scaled SRNs.

B.3.1.   Preparations

In preparation for the proof, we consider the following. Let $\mathcal{N}$ be an SRN with non-interacting species set $\mathcal{U}$, fast proper reactions $\mathcal{F}\subseteq \mathcal{R}_ \mathcal{U}$, and transition intensities $\lambda_r, r\in \mathcal{R}$. For $x\in \mathbb{Z}^n_{\geq0}$, denote by $\mathcal{N}(x)$ the set of reachable states from $x$ of $\mathcal{N}$, and for $z\in \mathbb{Z}^p_{\geq0 }$, denote by $\mathcal{N}_{\mathcal{U}, \mathcal{F}}(z)$ the set of reachable states from $z$ of $\mathcal{N}_{\mathcal{U}, \mathcal{F}}$.

We introduce a surrogate model $\mathcal{M}$ that reflects the behavior of the SRN of $\mathcal{N}$ when $\epsilon$ is very small. Specifically, this model is a CTMC on $\mathbb{Z}^n_{\geq 0}$ determined by the $Q$-matrix

with transition intensities,

Similarly, for $x\in \mathbb{Z}^n_{\geq0}$, we denote by $\mathcal{M}(x)$ the set of reachable states from $x$ of $\mathcal{M}$.

Let $N_i: = \{x = (z, u)\in \mathbb{Z}_{\geq 0}^n|\sum_{j = 1}^m u_j = i\}$ for $i\geq 0$. Then, we have the following.

Proposition B.7. Suppose Assumption 1 holds, that $\mathcal{F}$ is proper, and let $x = (z, 0)\in N_0$. Then, $\rho_{ \mathcal{O}}(\mathcal{M}(x)\cap N_0)) = \mathcal{N}_{\mathcal{U}, \mathcal{F}}(z)$ and $\mathcal{M}(x)\subseteq \mathcal{N}(x)$

For the proof of Theorem 4.1, the following Lemma will be useful. The proof follows by contradiction and is omitted.

Lemma B.8. Suppose Assumption 1 holds and that $\mathcal{F}$ is proper. If Assumption 2 holds for some closed set $E$, then for any $x\in E\cap N_0$, the states in $\mathcal{M}(x)\cap N_1$ are transient for the CTMC of the SRN obtained by considering the RN defined by the reaction set $\mathcal{F}$ (and taking the same transition intensities as in $\mathcal{N}$).

Consider the family of scaled CTMCs as in the setting of $\S$ 4 with notation introduced in $\S$ B.3.1. Furthermore, consider a finite closed set $E\subseteq \mathbb{Z}_{\geq 0}^{n}$. As $E$ is finite, we can assume wlog that $\mathcal{N}$ is finite.

For any probability distribution $\pi$ on $E$ and any state $w\in E$, $P_\pi(X(t) = w) = \sum_{x\in E}\pi_xP_x(X(t) = w)$. Therefore it is enough to prove the statement for:

$-$ a closed set of the form $\tilde{E}: = \mathcal{N}(x)$ with $x\in N_0\cap E$. This follows as $E = \cup_{x{^{\prime}}\in E\cap N_0} \mathcal{N}(x)$ by definition.

$-$ an initial distribution with support on $\mathcal{M}(x{^{\prime}})\cap N_0$ for $x{^{\prime}}\in \tilde{E}\cap N_0$. This follows as by definition $\tilde{E}\cap N_0 = (\cup_{x{^{\prime}}\in \tilde{E}\cap N_0}\mathcal{M}(x{^{\prime}}))\cap N_0$, see Proposition B.7.

Hence, we might assume that $\mathcal{N}$ is finite, that $E: = \mathcal{N}(x)$ for some $x\in N_0$ and that $\pi$ has support on $\mathcal{M}(x{^{\prime}})\cap N_0\subseteq E$ for $x{^{\prime}}\in E$. We denote $M: = \mathcal{M}(x{^{\prime}})$.

B.3.2.   The proof

With the preparations above, the matrix $\widetilde{Q}$ from $\S$ B.1 is comprised of the transition intensities from reactions of $\mathcal{F}$ with a non-interacting species of $\mathcal{U}$ in the reactant and $\widehat{Q}$ is comprised of the other reactions $\mathcal{R}\setminus \mathcal{F}$, where we divide the states of $E$ as follows:

where by definition $E = \mathcal{I}_{S_1}\cup \mathcal{I}_{S_2}\cup \mathcal{I}_{F_1}\cup \mathcal{I}_{F_2}$. We next illustrate the possible slow transitions (that is, transitions from $\widehat{Q}$) outgoing from $\mathcal{I}_{S_1}$ in blue on the left and all possible fast transitions (that is, transitions from $\widetilde{Q}$) displayed in red on the right by their transition state diagrams.

Note that there are other possible slow transitions (that is, transitions from $\widehat{Q}$) outgoing from $\mathcal{I}_{S_2}, \mathcal{I}_{F_1}, \mathcal{I}_{F_2}$, but the transitions from $\mathcal{I}_{S_1}$ are enough for our considerations.

The transition diagrams above determine the zero/non-zero parts of the corresponding $Q$-matrices in block-form as follows.

Lemma B.9. The following holds for the fast and slow parts of the $Q$-matrices (notation as in $\S$ B.1) divided into blocks:

Proof. Recall that the slow transitions come from $\mathcal{R}\setminus \mathcal{F}$ and the fast ones from $\mathcal{F}$ by assumption. We go through the slow transitions and then the fast ones.

$-$ By contradiction there are no slow transitions from $\mathcal{I}_{S_1}$ to $\mathcal{I}_{S_2}$ or $\mathcal{I}_{F_2}$.

$-$ By Assumption 1, $\widetilde{Q}$ has zero entries on $E_0 = \mathcal{I}_{S_1}\cup \mathcal{I}_{S_2}$ (the coordinates of non-interacting species are zero). By the definition of the sets $\mathcal{I}_{F_1}, \mathcal{I}_{S_2}$ and $\mathcal{I}_{F_2}$ there are no fast transitions from $\mathcal{I}_{F_1}$ to $\mathcal{I}_{S_2}$ or $\mathcal{I}_{F_2}$.

Following $\S$ B.1 (see also ^[22]), the distribution of the CTMC converges to the zeroth order outer expansion $\varphi_0(t)$ for $t\geq 0$, which we consider as

By assumption, the initial distribution $\pi_0$ has support contained in $\mathcal{I}_{S_1}$. We next go through the roles of the states in $\mathcal{I}_{S_1}, \mathcal{I}_{S_2}, \mathcal{I}_{F_1}, \mathcal{I}_{F_2}$.

States in $\mathcal{I}_{F_1}$ are transient in $\widetilde{Q}$ by Assumption 2 and Lemma B.8, and $\widetilde{Q}|_{ \mathcal{I}_{F_1}\times \mathcal{I}_{F_1}}$ is non-singular (see Lemma B.5). States in $\mathcal{I}_{S_1}, \mathcal{I}_{S_2}$ are absorbing in $\widetilde{Q}$. States in $\mathcal{I}_{F_2}$ can be transient, absorbing or part of a closed communicating class in $\widetilde{Q}$. Correspondingly, $\varphi_0(t)$ is zero on $\mathcal{I}_{F_1}$ and on the transient part of $\mathcal{I}_{F_2}$.

Next, we make the following general observation.

Remark B.10. Let $(X(t))_{t\ge 0}$ be a CTMC with generator $Q$ on a state space $S$. Let $\pi_0$ be a probability distribution with ${{\rm{supp\, }}}(\pi)\subseteq Z\subseteq S$, where $|Z| < \infty$. Assume that there are no $x\in Z, y\in S\setminus Z$ such that $x\to y$. Then, $Q|_{Z\times Z}$ is a $Q$-matrix, and denoting by $X^{restr}$ the corresponding CTMC, we have

Then, it is enough to consider the restricted $Q$-matrix $Q|_{A\times A}$ and $X^{restr}$ for computations of such probabilities.

Lemma B.11. The CTMC of the generator $\overline{Q}_*$ (see (B.3) of $\S$ B.1) in our setting has the property of Remark B.10 with $Z: = \mathcal{I}_{S_1} = M\cap N_0$.

Proof. First, we decompose $\mathcal{I}_{F_2}$ into three different parts $F_{2, t}, F_{2, w}$ and $F_{2, a}$ according to whether a state is transient, part of a (non-trivial) irreducible component or absorbing in $\widetilde{Q}$. Then, writing

corresponds to the situation of $\S$ B.1.

Let $S{^{\prime}}$ be the reduced state space of $\overline{Q}_*$ (cf., $\S$ B.1). Clearly $A\subseteq S{^{\prime}}$, as these are the absorbing states in $\widetilde{Q}$, and at least the transient states have been eliminated from $S{^{\prime}}$. Hence $\mathcal{I}_{S_1}\subseteq S{^{\prime}}$. We next focus on $\overline{Q}_*|_{ \mathcal{I}_{S_1}\times S{^{\prime}}}$, which contains the possible outgoing transitions from $\mathcal{I}_{S_1}$. It is enough to show that $\overline{Q}_*|_{ \mathcal{I}_{S_1}\times S{^{\prime}}\setminus \mathcal{I}_{S_1}}$ has only zero entries. For this, it suffices to consider the matrices in the middle of (B.3), which we recall here,

We next check the right summand of (B.7), which we denote

Then, consider the restriction to $R|_{ \mathcal{I}_{S_1}\times S{^{\prime}}}$. Keeping the decomposition of $A, T, W$ in (B.6) in mind, the following hold for the involved matrices $\widehat{Q}_{AT},$ $\widetilde{Q}^{-1}_{TT},$ $\widetilde{Q}_{TA},$ $\widetilde{Q}_{TW}$ by Lemma B.9:

and we get that

Hence, the property holds for $R$, that is, the right hand side of (B.7) has the required property.

Next we look at the left hand side of (B.7), and again restrict to $\mathcal{I}_{S_1}\times S\setminus T$. By Lemma B.9, the slow part has only non-zero transitions to $\mathcal{I}_{S_1}$ or $\mathcal{I}_{F_1}$, hence also this matrix has the required property. As this holds for both summands of (B.7), we have shown that the property holds for $\overline{Q}_*$.

Finally, noting the form of equation (B.5) and equation (B.8) we see that only the slow and fast parts outgoing from $\mathcal{I}_{F_1}$ and $\mathcal{I}_{S_1}$ contribute to $\overline{Q}_*|_{ \mathcal{I}_{S_1}\times \mathcal{I}_{S_1}}$. Then we can compute $\overline{Q}_*$ for the state space $\mathcal{I}_{F_1}\cup \mathcal{I}_{S_1}$ (i.e. with restricted $Q_ \epsilon, \widetilde{Q}, \widehat{Q}$) or for the state space $\mathcal{I}_{S_1}\cup \mathcal{I}_{S_2}\cup \mathcal{I}_{F_1}\cup \mathcal{I}_{F_2}$. However, in both cases the expression we get for $\overline{Q}_*|_{ \mathcal{I}_{S_1}\times \mathcal{I}_{S_1}}$ is the same, hence the following holds.

Lemma B.12. If the initial distribution has support on $\mathcal{I}_{S_1}$, for the computation of $\varphi_0(t)$, it is enough to restrict the state space of $Q_ \epsilon$ to $\mathcal{I}_{F_1}\cup \mathcal{I}_{S_1}$.

We treat the restriction of the state space to $\mathcal{I}_{F_1}\cup \mathcal{I}_{S_1}$ in the following. Restricting to $\mathcal{I}_{F_1}\cup \mathcal{I}_{S_1}$, we have

where the states $\mathcal{I}_{S_1}$ are absorbing in $\widetilde{Q}$, and the states of $\mathcal{I}_{F_1}$ are transient states. Note that this setting corresponds to the situation of $\S$ B.1 where we removed the irreducible part (that is, $\mathcal{I}_W$ in the notation of $\S$ B.1), $\mathcal{I}_{F_1}$ corresponds to $T$ and $\mathcal{I}_{S_1}$ corresponds to $A$ of $\S$ B.1.

The scaled CTMC converges to the zeroth order outer expansion $\varphi_0(t)$ for $t > 0$, giving $\varphi_0(t): = (\varphi_{S_1}(t), \varphi_{F_1}(t)) = (\varphi_{S_1}(t), 0_{F_1})$ which satisfies the ODE from $\S$ B.1 with

with initial distribution $\pi_0 = (p_{S_1}, p_{F_1})$. As the initial distribution has support on $\mathcal{I}_{S_1}$, we get the following from Proposition B.3, where again $X_ \epsilon(t)$ is the CTMC associated to $Q_ \epsilon$ while $X_0(t)$ is the CTMC from $C$ (corresponding to $\overline{Q}_\ast$ from $\S$ B.1).

Lemma B.13. Let $\pi_0$ be a probability distribution with support on $\mathcal{I}_{S_1}$, and $B\subseteq \mathcal{I}_{S_1}$. Then the following holds for all $T > 0$:

Finally, observing the form of the matrix $\widehat{Q}_{S_1S_1}+\widehat{Q}_{S_1F_1}(\widetilde{Q}_{F_1F_1})^{-1}\widetilde{Q}_{F_1S_1}$ in (B.9), we can view it as the watched Markov chain of $\S$ B.2, where it is watched when in $\mathcal{I}_{S_1}$ with $Q$-matrix

By the interpretation of the censored Markov chain available from Remark B.6, the part $\widehat{Q}_{S_1F_1}(\widetilde{Q}_{F_1F_1})^{-1}\widetilde{Q}_{F_1S_1}$ corresponds to transition intensities given by the rates entering a state in $\mathcal{I}_{F_1}$ from $\mathcal{I}_{S_1}$ times the exit probabilities to $\mathcal{I}_{S_1}$. As $\widetilde{Q}_{F_1S_1}, \widetilde{Q}_{F_1F_1}$ have the transition intensities from the reactions $\mathcal{F}$ and $\widehat{Q}_{S_1F_1}$ from $\mathcal{R}\setminus \mathcal{F}$, this corresponds exactly to the transition rates of the defined reduced RN of $\S$ 3.2. Hence Theorem 4.1 follows.

C.   Proof of Theorem 4.1 for non-homogeneous SRNs

The proof follows the same steps as the proof of Theorem 4.1, hence it is enough to check that Assumptions 2 and 3 together with $\mathcal{F}$ being proper are sufficient for the assumptions made in ^[22] to hold. Recall that since the $\lambda_{y \to y{^{\prime}}}(t, \cdot)$ satisfy Assumption 1 for each $t\in [0, T]$, the state space decomposition of the ordinary Markov process for an RN and $Q_t$ for $t$ fix agree. Hence, by Assumption 2 and Lemma B.8, any $x\in E$ with one molecule of a non-interacting species is transient in $\widehat{Q}_t$. Furthermore by Assumption 3, $\widetilde{Q}_t, \widehat{Q}_t$ are once continuous differentiable with Lipschitz derivative ([22,Assumption A4.4(respectively,A4.5)]). With these observations we can conclude.

Acknowledgement

LH acknowledges funding from the Swiss National Science Foundations Early Postdoc.Mobility grant (P2FRP2_188023). The work presented in this article is supported by Novo Nordisk Foundation, grant NNF19OC0058354.

Conflict of interest

The authors declare there are no conflicts of interest.