Fast reactions with non-interacting species in stochastic reaction networks

We consider stochastic reaction networks modeled by continuous-time Markov chains. Such reaction networks often contain many reactions, potentially occurring at different time scales, and have unknown parameters (kinetic rates, total amounts). This makes their analysis complex. We examine stochastic reaction networks with non-interacting species that often appear in examples of interest (e.g. in the two-substrate Michaelis Menten mechanism). Non-interacting species typically appear as intermediate (or transient) chemical complexes that are depleted at a fast rate. We embed the Markov process of the reaction network into a one-parameter family under a two time-scale approach, such that molecules of non-interacting species are degraded fast. We derive simplified reaction networks where the non-interacting species are eliminated and that approximate the scaled Markov process in the limit as the parameter becomes small. Then, we derive sufficient conditions for such reductions based on the reaction network structure for both homogeneous and time-varying stochastic settings, and study examples and properties of the reduction.


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 [25].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 Z n ≥0 , 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 [6].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) [32] and singular perturbation theory in the sense of Tikhonov-Fenichel [35,9] have been the main means to derive lower dimensional reduced models.(See [13,11] and references therein for the relationship between the two approaches and when the QSSA is valid.)A motivation for our work is the situation for which reactions involving so-called non-interacting species [8,30] in the reactant has reaction rate constants scaled by 1/ , where > 0 is a small number [7].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 [31,19,12].One stream of research eliminates species using heuristic projection arguments [16,17].Rigorous simplifications and reductions often follow scaling limits of Markov processes in a multi-scale setting [23].These might be applied to concrete examples, if certain conditions are satisfied and the different scaling parameters are balanced (in a specific sense) [3,20,28].Scaling laws for a special class of reaction networks with intermediate species (a special type of noninteracting species) and their explicit reduction are given in [4].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 λ y→y (x) := 1 λ y→y (x).
• Slow reactions with unscaled transition intensity λ y→y (x) := λ y→y (x).The fast reactions are determined by non-interacting species.Our work is inspired by previous work on non-interacting species in deterministic systems [8,30,7] and intermediate species in stochastic models [4].The technical part is based on singularly perturbed Markov chains [37,38] and watched Markov chains [10].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 κ1κ2z S 1 (z S 2 +1) κ2(z S 2 +1)+κ3 for the first reaction and κ1κ3z S 1 κ2(z S 2 +1)+κ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 −−→ S 3 , or without consumption of the S 2 molecule, resulting in the net reaction S 1 −−→ 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 § 2 preliminaries on graph theory and reaction networks are covered.In § 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 § 4, we study transient approximations for stochastic reaction networks with non-interacting species via the previously introduced reduction.In § 5, we give realistic examples, study sufficient conditions for reductions and compare stationary properties of the reduction with the original reaction network.Finally in § 6, we discuss the results, approach, and elaborate on the relation to the literature.In the Appendices § A, § B, § C, we give proofs as well as brief introductions to the theory on singularly perturbed CTMCs and watched CTMCs.

Preliminaries
2.1.Notation.Let R p be the real p-dimensional space, and R p >0 (R p ≥0 ) the subset of elements of R p with strictly positive (non-negative) entries in all components.A vector y ∈ R p is written as (y 1 , • • • , y p ), where y i is the i-th component.For vectors denotes the component-wise maximum, and y 1 ≥ y 2 denotes component-wise inequality.The inner product between y 1 and y 2 is denoted y 1 , y 2 .The cardinality of a set A is denoted |A|.

potentially listed as the corresponding sequence of edges).
A multi-digraph 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.

Reaction networks (RNs
where S is a finite set of species S = {S Reactions are directed edges between complexes, written as r i : y i → y i or, generically, as r : y → y , potentially omitting r i , r.A reaction is said to consume the reactant y and create the product y .An RN is said to be finite if R is finite, and otherwise it is infinite. We diverge in two ways from the standard definition of RNs: trivial reactions r : y → 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 [34,33].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.(i) Two species interact if they both appear in a complex of C.
(ii) A subset U ⊆ S is non-interacting if it contains no pair of interacting species, and the stoichiometric coefficients of the species in U in all complexes are either 0 or 1.The species of U are said to be non-interacting.
(iii) If U is non-interacting, the species in S \ U are said to be core species.
Example 2.2.Consider a two-substrate Michaelis Menten mechanism [15, Section 3.1.2]:In this paper, it is convenient to work with the directed stoichiometric subspace of N = (C, R), defined as (note, this definition allows the RN to be infinite).For v ∈ R n , the set (v +T )∩R n ≥0 defines a directed stoichiometric compatibility class of N .The RN is said to be conservative, respectively, sub-conservative, if there exists a positive vector c ∈ R n >0 , such that c, y − y = 0, respectively, c, y − y ≤ 0 for any reaction r : y → y ∈ 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.

Stochastic reaction networks (SNRs
).An SRN is an RN together with a CTMC X(t), t ≥ 0, on Z n ≥0 , modeling the number of molecules of each species over time.A reaction r : y → y fires with transition intensity λ r (x), in which case the state jumps from X(t) = x to x + y − y [2].The Markov process with transition intensities λ r : For (stochastic) mass-action kinetics, the transition intensity for r : y → y is , where x! := n i=1 x i !, and κ r is a positive reaction rate constant [2].If there are infinitely many reactions, we assume (1) such that the corresponding CTMC is well-defined in the sense that it has no instantaneous jumps [26].
Assumption 1.For all reactions, r : y → y ∈ R, the transition intensity λ r : Z n ≥0 → R ≥0 satisfies the following In particular, (X(t)) t≥0 , is confined to the directed stoichiometric compatibility class (x + T ) ∩ Z n ≥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 U ⊆ S be a noninteracting subset of species.For simplicity, we let be the projections onto the first p coordinates and last m coordinates of R n , respectively.Consequently, we denote a state by x = (z, u) ∈ Z p ≥0 × Z m ≥0 = Z n ≥0 .Define the following sets of reactions such that R U ∪ R U are the reactions that involve species in U.
We consider a subset of fast (a terminology to be motivated below) reactions F ⊆ R U with the following structural property.Definition 3.1.The reactions in F are proper w.r.t.U, that is, any non-interacting species is part of a sequence of reactions in (R U \ R U ) ∪ F of the form 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 U 1 = {EA, EAB} is a non-interacting set and F = R U is a set of fast proper reactions.For example, the following is a fast chain: Also, U 2 = {EA, P } is a non-interacting set of species, but F = R U 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 [30].We later use this graph to define a reduced RN and a reduced SRN by elimination of non-interacting species.Definition 3.3.Let N = (C, R) be an RN on S, U ⊆ S a set of non-interacting species, and F ⊆ R U a set of fast proper reactions.Let G U ,F = (V U ,F , E U ,F ) be the labeled multi-digraph with vertex set and edge set Furthermore, let L : E U ,F → {1, . ..} be the function that maps γ ∈ E U ,F to the corresponding index of the reaction label, that is, if γ ∈ E U ,F has label r i , then The construction of the graph G U ,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.
. There are infinitely many walks from * in to * 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 G U ,F with start vertex * in and end vertex * out is a fast chain.Define the set of such walks by Note that W U ,F might be infinite, even for finite RNs, as in Example 3.4 above.
We will consider the situation in which the reactions of F occur fast compared to the remaining reactions R \ F, in the sense that the transition intensities of the reactions in F are scaled by 1/ for small .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 U ,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.
Then the following holds: • r(Γ) ≥ 0, and p(Γ) := r(Γ) + Proof.We give complete proofs here for convenience, but note that the results also follow from [14].The second item follows by definition.As we take the maximum coordinate-wise and w 1 ≥ 0, then r(Γ) has non-negative coordinate in each species.To see that p(Γ) has non-negative coordinate in each species, we note that by definition y to both sides, we get 0 ≤ y L(γ l ) ≤ p(Γ), and the result follows.By Lemma 3.5, we might consider r(Γ) and p(Γ) as elements in Z p ≥0 , and r(Γ) → p(Γ) as the reduced reaction obtained by contraction along Γ. 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 Γ 0 of G U ,F consisting of the two edges associated to r 1 , r 2 .Then a trivial reaction.Likewise, the walk Γ 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 Γ B 1 consisting of the three edges associated to r 1 , r 3 and consists either of a sequence of edges with labels r 1 , r 3 , r 4 , . . ., r 3 , r 4 , r 3 , r 5 or a sequence of edges with labels r 1 , r 3 , r 4 , . . ., 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 Γ A 1 , respectively, Γ 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 Γ The lemma implies that the reactions corresponding to a walk in W U ,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 r(Γ).A similar statement appears in [14,Corollary 3.2].
Lemma 3.7 allows us to define transition intensities of the reduced reactions in a natural way.For Γ = (γ 1 , γ 2 , . . ., γ l ) ∈ W U ,F , define the function Λ Γ : , , and out(s(γ j )) denotes the set of outgoing edges of s(γ j ) in G U ,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 ( 2) is well-defined due to (1).
Definition 3.8.Let N = (C, R) be an RN on S, let U ⊆ S be a non-interacting set of species, and and C U ,F being the set of vertices of R U ,F .
For r ∈ R U ,F , define Definition 3.9.Let N = (C, R) be an SRN on S with transition intensities λ r (•), r ∈ R, satisfying Assumption 1.Let U ⊆ S be a set of non-interacting species, and F ⊆ R U a set of proper fast reactions.Then, the possibly infinite SRN is the reduced SRN, obtained by elimination of (U, F) from N .
As a reduced SRN has species set O = S \ U, the associated CTMC lives in Z p ≥0 .Even if the reduced SRN has infinitely many reactions, the associated CTMC is always well-defined, that is, it has no instantaneous states [26].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, We note that the last inequality is not necessarily an equality, 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 → y ∈ R U ,F of the reduced SRN, it holds that τ r (z) > 0 if and only if z ≥ y.
Only one direction is not obvious.Assume z ≥ r(Γ).Using the definition of Λ Γ (•), it is sufficient to show that the numerators in the fraction are all non-zero.This holds by Lemma 3.7, hence Λ Γ (z) > 0.
Example 3.12.We continue with Examples 3.4, 3.6 assuming stochastic massaction kinetics.There are three reduced reactions with the first two being trivial, with infinitely many walks underlying the second 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 Analogously, we compute the sums over walks for the trivial reactions.We have τ 1 (z) = Λ Γ0 (z) and .
Furthermore, a simple calculation gives the following identity for z ∈ Z p ≥0 , 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.
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 § 5.1.

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 N = (C, R) be an SRN on a species set S (of size n) with non-interacting species U ⊆ S, fast proper reactions F ⊆ R U , and transition intensities λ r , r ∈ R. From this SRN, we construct a family of SNRs, indexed by a parameter > 0, with corresponding Q-matrix Q .In particular, the transition intensities of the reactions in F consuming non-interacting species are scaled in as follows: • λ r (x) := 1 λ r (x), r ∈ F, that is, for small the reactions are fast compared to the remaining reactions, • λ r (x) := λ r (x), r ∈ R\F, that is, the transition intensities are independent of , and the reactions are slow.Thus, we might write the corresponding Q-matrix as as sum of two terms, a fast part Q, scaled by 1 , and a slow part Q, It follows from Definition 3.9 that the reduced SRN, obtained from the SRN with transition intensities λ r , r ∈ R, has transition intensities τ r , r ∈ R U ,F , independent of .Let (X (t)) t≥0 be Markov chains on Z n ≥0 with generators Q , > 0, respectively.We note that all Q , > 0, are dynamically equivalent in the sense that the connectivity of the state space and the decomposition of the state space Z n ≥0 into communicating classes are independent of .Furthermore, let (X 0 (t)) t≥0 be a Markov chain on Z p ≥0 (excluding the m = n − p coordinates for non-interacting species) with generator Q 0 , obtained from the transition intensities τ r , r ∈ R U ,F .In the following, we provide conditions that guarantee that the dynamics of (X (t)) t≥0 is similar (in a sense to be made precise) to the dynamics of (X 0 (t)) t≥0 , whenever is small.Technically, we will consider the limit as → 0.
We restrict attention to a particular closed set E of Q , such that E ∩ (Z p ≥0 × {0}) = ∅.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) ∈ E 0 := ρ O (E∩(Z p ≥0 ×{0})), then by construction of the reduced reactions, X 0 (t) ∈ 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 ∈ E 0 satisfy the following: (3) 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 : 0 → U, r 2 : S + U → 0; if there 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 : 0 → U , r 2 : 2S + U → 3S + U , r 3 : U → 0; there is positive probability of an infinite walk r 1 , r 2 , r 2 , . . .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 λ r , r ∈ R, but only on the structure of the underlying RN (together with U and F).Furthermore, as (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 ⊆ E ∩ (Z p ≥0 × {0}), we define B 0 = ρ O (B).Theorem 4.1.Let N = (C, R) be an SRN on S with transition intensities λ r , r ∈ R, U ⊆ S a set of non-interacting species, and F ⊆ R U a set of fast proper reactions.Suppose Assumptions 1 and 2 hold.Let π be a probability distribution on E ∩ (Z p ≥0 × {0}), and π 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 ⊆ E ∩ (Z p ≥0 × {0}) and any 0 < T : The proof is in Appendix § B.3.As a consequence of the theorem we have: Then, for all t ≥ 0 it holds that: lim If the initial state X (0) has more than one molecule of the non-interacting species, then these will in general be depleted quickly for small .This is however not always the case, see Example 5.2.

The case of time-dependent transition intensities.
In the previous section, we considered the scaling limit of time-homogeneous SRNs.While timeconstant 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 [1].
The setup is essentially the same as in Section 4.1, except we now allow timedependent transition intensities, λ r (t, x), r ∈ R, for t in a compact time interval t ∈ [0, T ].We make the following assumption.
Assumption 3. The time-dependent transition intensities, λ r (t, x), r ∈ R, are such that λ r (t, •) satisfies Assumption 1 for all t ∈ [0, T ] and such that for all 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 ∈ [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.

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.

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 κ i as labels of the reactions, and omit the reaction names r i .A key source book for realistic examples is [15], that contains reaction networks used in systems biology.Other useful references are [5], which focuses on enzyme kinetics, and [25], which has examples from a range of areas in biology.
Example 5.1.We first consider the two-substrate Michaelis-Menten mechanism [15, 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) U = {EA, EAB} and fast reactions F = R U = {r 2 , r 3 , r 4 , r 5 }.Assuming the same in the stochastic setting gives the reduced RN in Example 3.6 with one non-trivial reaction, Assumption 2 is fulfilled for all z ∈ Z 5 ≥0 (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 Z 5 ≥0 × {0} non-trivially.
Example 5.2.Consider the two-substrate Michaelis-Menten mechanism with U = {EA, EAB} as above, but fast reactions F = {r 3 , r 4 , r 5 } ⊆ R U , that is, 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 Note, that this example does not have mass-action-like kinetics, as τ 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 [15, 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 U = {ES, EI, ESI}.In our setting, taking U to be non-interacting species with F = R U , then the reduced SRN becomes For any closed set E in a directed stoichiometric compatibility class that intersects Z 4 ≥0 × {0} non-trivially, Assumption 2 is satisfied.Example 5.4.Consider an example of allosteric activation [15,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) U = {ER, ERS}.Choosing analogously F = R U in the stochastic case gives a proper set of fast reactions.Then, the reduced RN has infinitely many reactions, given by: All z ∈ Z 4 ≥0 (p = 4) satisfy (3).As 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 Z 4 ≥0 × {0} non-trivially, Assumption 2 is satisfied.
Example 5.5.As a final example, consider a mechanism-based inhibitor system [25, 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 [25,Section 6.4].For this, often the QSSA is applied to the short-lived species U = {X, Y }.Analogously, we consider the stochastic system with F = R U being fast and proper, and study the corresponding stochastic reduction.The reduced SRN has reactions with transition intensities: .
As U = {X, Y } are intermediate species, all z ∈ Z 4 ≥0 (p = 4) satisfy (3).As 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 F being proper are both easy to check, Assumption 2 is nontrivial 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 § A for proof).
Proposition 5.6.Assume N is a sub-conservative SRN on a species set S, U ⊆ S is a set of non-interacting species, and F ⊆ R U is a proper set of reactions.Suppose Assumption 1 holds.Furthermore, assume that one of the following holds: 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 → 0. First one might ask whether states x = (z, 0) ∈ Z p ≥0 × {0} and z are of the same type (positive recurrent/null recurrent/transient) for N and N U ,F , respectively.Furthermore, the stationary distribution (if it exists) for the scaled SRN N as → 0 might match the stationary distribution of N U ,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 N with mass-action kinetics, Example 5.9.Consider the SRN N with mass-action kinetics, The next result follows from [14,Theorem 5.6].We recall that irreducible components of sub-conservative RNs are finite [2].For an irreducible component ). Theorem 5.10.Let N = (C, R) be a sub-conservative SRN on S, U ⊆ S a set of intermediate species, and F = R U a set of proper reactions.Suppose Assumption 1 holds.Then x = (z, 0) ∈ Z p ≥0 × {0} is transient for N if and only if z is transient for N U ,F , and the same holds for positive recurrence.In particular, if a subset E is an irreducible component for N then E 0 is an irreducible component for N U ,F .
From Theorem 5.10, Lemma 5.7 and [18, Theorem 3], we get convergence of the stationary distribution.
Corollary 5.11.Let N = (C, R) be a sub-conservative SRN on S, U ⊆ S a set of intermediate species, and F = R 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 π the unique stationary distribution of Q on E, and let π 0 be the unique stationary distribution of the reduced SRN on E 0 .Then, π ((z, 0)) → π 0 (z), for z ∈ E 0 , and → 0.

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 on [37,38], 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 [4].In contrast to our setting, [4] 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 [22,23,3,20,31].This is not required in our setting (nor in [4] 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 → 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 [4,Example 5.3] or [20, § 6.5], because even for small , 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.
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 [37,38].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 [36], 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 [18] 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 can be written as two reactions with transition intensities κ 1 z S1 1 {0} (z S3 ) for the first reaction and 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 F and the example falls outside our setting.As a matter of fact, the scaling parameter 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 ).This distribution is not concentrated on the part of the state space without molecules of non-interacting species.Rescaling time by would retrieve the original chain.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 [37,38] for homogeneous CTMCs on a finite state space, using their notation.For more on the two-scale setting for CTMCs, we refer to [27,37,38].
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 [38, § 4.3, § 4.4, § 4.5].
Let a sequence of homogeneous generators of a CTMC be given by on a finite state space I.By rearrangement, the states can be divided into blocks of states I A , I W , I T according to their role in Q, which are absorbing, irreducible, and transient, respectively.The matrices Q, Q are partitioned accordingly into block-matrices and take the form The states in I A with Q AA = 0, Q AW = 0, Q AT = 0 correspond to absorbing states of Q.The states in I T correspond to transient states of Q (hence Q T T is a Hurwitz matrix, furthermore stable and non-singular see Lemma B.5).Finally, I W consists of the states that are part of (non-trivial) closed irreducible components, so Q W W in (4) can further be decomposed as ( 5) , where Q i is an m i × m i matrix.We denote by ν i the stationary distribution of Q i for i = 1, . . ., l (which exist by assumptions on the state space).
Remark B.1.Our presentation differs slightly from [38].In that work, so-called weakly irreducible classes is used [38,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 Q via I T , this is not a restriction.We note that our setting is included in [38, § 4.5].
For > 0, the forward equation of the CTMC gives an ODE where π 0 is the initial probability distribution.
In [38], they construct sequences of functions that approximate p (t) in the limit for → 0 uniformly on [0, T ] for different powers of .One contribution comes from the so-called outer expansion, that approximates p (t) for t > 0. Another part comes from the so-called initial-layer correction, that approximates p (t) in a neighborhood of t = 0.This might be ignored as it is zero in our case [38].Hence, we will focus on the zeroth order outer expansion as this is enough to obtain an approximation to order O( ) (see Proposition B.3).
Following [38,Theorem 4.45], the distribution of the CTMC converges to the zeroth order outer expansion ϕ 0 (t) for t > 0. In the case we consider, it is given as where ϑ i (t) are scalar functions and ν i are the stationary distributions from matrix (5), and 0 T is the vector of zero-entries.Then, let (with m 0 = 0) where and define the matrix (6) ).
Note that Q * 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. diag(1, Then the zeroth order outer expansion ϕ 0 (t) (corresponding to [38, (4.86)]) with initial distribution π 0 = (p A , p W , p T ) is determined by the following ODE, which defines the generator of a CTMC that is given by Q * .This CTMC has reduced state space I A ∪ {s 1 , • • • , s l }, where each state s i corresponds to the assembled irreducible component, which has stationary distribution ν i from Q W W , For any probability distribution π on E and any state w ∈ E, P π (X(t) = w) = x∈E π x P x (X(t) = w).Therefore it is enough to prove the statement for: -a closed set of the form Ẽ := N (x) with x ∈ N 0 ∩ E. This follows as E = ∪ x ∈E∩N0 N (x) by definition.-an initial distribution with support on M(x ) ∩ N 0 for x ∈ Ẽ ∩ N 0 .This follows as by definition Ẽ ∩ N 0 = (∪ x ∈ Ẽ∩N0 M(x )) ∩ N 0 , see Proposition B.7.Hence, we might assume that N is finite, that E := N (x) for some x ∈ N 0 and that π has support on M(x ) ∩ N 0 ⊆ E for x ∈ E. We denote M := M(x ).B.3.2.The proof.With the preparations above, the matrix Q from § B.1 is comprised of the transition intensities from reactions of F with a non-interacting species of U in the reactant and Q is comprised of the other reactions R \ F, where we divide the states of E as follows: where by definition E = I S1 ∪ I S2 ∪ I F1 ∪ I F2 .We next illustrate the possible slow transitions (that is, transitions from Q) outgoing from I S1 in blue on the left and all possible fast transitions (that is, transitions from Q) displayed in red on the right by their transition state diagrams.
These transition diagrams 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 § B.1) divided into blocks: Proof.Recall that the slow transitions come from R \ F and the fast ones from F by assumption.We go through the slow transitions and then the fast ones.
-By contradiction there are no slow transitions from I S1 to I S2 or I F2 .
-By Assumption 1, Q has zero entries on E 0 = I S1 ∪ I S2 (the coordinates of non-interacting species are zero).By the definition of the sets I F1 , I S2 and I F2 there are no fast transitions from I F1 to I S2 or I F2 .
By assumption, the initial distribution π 0 has support contained in I S1 .We next go through the roles of the states in I S1 , I S2 , I F1 , I F2 .
States in Next, we make the following general observation.
Remark B.10.Let (X(t)) t≥0 be a CTMC with generator Q on a state space S. Let π 0 be a probability distribution with supp (π) ⊆ Z ⊆ S, where |Z| < ∞.Assume that there are no x ∈ Z, y ∈ S \ Z such that x → y.Then, Q| Z×Z is a Q-matrix, and denoting by X restr the corresponding CTMC, we have P π0 (X(t) ∈ B) = P π0 (X restr (t) ∈ B), for B ⊆ S, t ≥ 0.
Then, it is enough to consider the restricted Q-matrix Q| A×A and X restr for computations of such probabilities.Let S be the reduced state space of Q * (cf., § B.1). Clearly A ⊆ S , as these are the absorbing states in Q, and at least the transient states have been eliminated from S .Hence I S1 ⊆ S .We next focus on Q * | I S 1 ×S , which contains the possible outgoing transitions from I S1 .It is enough to show that Q * | I S 1 ×S \I S 1 has only zero entries.For this, it suffices to consider the matrices in the middle of ( 6), which we recall here, (10) We next check the right summand of (10), which we denote Then, consider the restriction to R| I S 1 ×S .Keeping the decomposition of A, T, W in (9) in mind, the following hold for the involved matrices Q AT , Q −1 T T , Q T A , Q T W by Lemma B.9: and we get that Hence, the property holds for R, that is, the right hand side of (10) has the required property.
Next we look at the left hand side of (10), and again restrict to I S1 × S \ T .By Lemma B.9, the slow part has only non-zero transitions to I S1 or I F1 , hence also this matrix has the required property.As this holds for both summands of (10), we have shown that the property holds for Q * .
Finally, noting the form of equation ( 8) and equation (11) we see that only the slow and fast parts outgoing from I F1 and I S1 contribute to Q * | I S 1 ×I S 1 .Then we can compute Q * for the state space I F1 ∪ I S1 (i.e. with restricted Q , Q, Q) or for the state space I S1 ∪ I S2 ∪ I F1 ∪ I F2 .However, in both cases the expression we get for Q * | I S 1 ×I S 1 is the same, hence the following holds.
Lemma B.12.If the initial distribution has support on I S1 , for the computation of ϕ 0 (t), it is enough to restrict the state space of Q to I F1 ∪ I S1 .
We treat the restriction of the state space to I F1 ∪I S1 in the following.Restricting to I F1 ∪ I S1 , we have where the states I S1 are absorbing in Q, and the states of I F1 are transient states.
Note that this setting corresponds to the situation of § B.1 where we removed the irreducible part (that is, I W in the notation of § B.1), I F1 corresponds to T and I S1 corresponds to A of § B.1.The scaled CTMC converges to the zeroth order outer expansion ϕ 0 (t) for t > 0, giving ϕ 0 (t) := (ϕ S1 (t), ϕ F1 (t)) = (ϕ S1 (t), 0 F1 ) which satisfies the ODE from § B.1 with (12) φS1 (t) = ϕ S1 (t)( with initial distribution π 0 = (p S1 , p F1 ).As the initial distribution has support on I S1 , we get the following from Proposition B.3, where again X (t) is the CTMC

r 4 :
EAB −−→ EA + B, r 5 : EAB −−→ E + P + Q, 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 S = {E, A, B, EA, EAB, P, Q} and complex set C = {E + A, EA, EA + B, EAB, E + P + Q}.The sets U 1 = {EA, EAB}, U 2 = {EA, P } and U 3 = {EAB} are non-interacting.

3 .
weakly reversible, and for any reaction r : y → y in ∈ R U \ R U , the reverse reaction r : y → y is in F. (c) N is weakly reversible and F = R 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 N is an SRN on a species set S, and that U ⊆ S is a set of intermediate species.Suppose Assumption 1 holds.Then, F ⊆ R U is proper if and only if all z ∈ Z p ≥0 satisfy (3) of Assumption 2. We summarize the properties of the examples from § 5.1.In examples 5.1, 5.3 and 5.4 we can directly conclude by Proposition 5.6 (b) that Assumption 2 holds.Example N sub-cons (3) fulfilled N U ,F finite 3Comparison of the original RN and the reduced RN.
with transition intensities satisfying Assumption 1 by Corollary 3.11.Any (a, b, c, 0) with a + b ≥ 1 is transient for N , but recurrent for N U ,F .Reaction r 5 in 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 N U ,F .

Appendix B .
Proof of Theorem 4.1 F1 are transient in Q by Assumption 2 and Lemma B.8, and Q| I F 1 ×I F 1 is non-singular (see Lemma B.5). States in I S1 , I S2 are absorbing in Q. States in I F2 can be transient, absorbing or part of a closed communicating class in Q. Correspondingly, ϕ 0 (t) is zero on I F1 and on the transient part of I F2 .

Lemma B. 11 .
The CTMC of the generator Q * (see (6) of § B.1) in our setting has the property of Remark B.10 with Z := I S1 = M ∩ N 0 .Proof.First, we decompose I F2 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 Q.Then, writing (9)A = I S1 ∪ I S2 ∪ F 2,a , T = I F1 ∪ F 2,t , W = F 2,w ,corresponds to the situation of § B.1. 2