Stationary distributions via decomposition of stochastic reaction networks

We examine reaction networks (CRNs) through their associated continuous-time Markov processes. Studying the dynamics of such networks is in general hard, both analytically and by simulation. In particular, stationary distributions of stochastic reaction networks are only known in some cases. We analyze class properties of the underlying continuous-time Markov chain of CRNs under the operation of join and examine conditions such that the form of the stationary distributions of a CRN is derived from the parts of the decomposed CRNs. The conditions can be easily checked in examples and allow recursive application. The theory developed enables sequential decomposition of the Markov processes and calculations of stationary distributions. Since the class of processes expressible through such networks is big and only few assumptions are made, the principle also applies to other stochastic models. We give examples of interest from CRN theory to highlight the decomposition.


Introduction
Reaction networks (CRNs) form a broadly applicable paradigm to describe the interactions of different constituents through mathematical models.CRNs are vital for the prediction and analysis of data in biochemistry, systems biology and cellular biology, and have found further applications [20,27,11].Besides their relevance in applications, CRNs continue to drive the development of areas of mathematics such as dynamical systems theory, stochastic processes and applied algebraic geometry [1,10].
A CRN consists of reactions with associated reaction rates that govern the speed of the reactions.CRNs are often defined via the reaction graph, that highlights the interactions between species and their transformations.As an example consider the enzymatic Michaelis-Menten mechanism, where an enzyme E catalyzes the conversion of a substrate S into a product P through an intermediate molecule ES: (1.1) S + E ⇋ ES → P + E.
Either a deterministic or a stochastic model is chosen to represent the dynamics of CRNs.Traditionally, deterministic models have been the preferred modelling choice.However, with the emergence of systems biology, cellular biology and synthetic biology the importance of modelling systems with small molecular counts have become important.Stochastic models of CRNs are used when the molecular counts in the system are low.They typically consist of continuous time Markov chains (CTMC), which apply to many processes in living systems [8,11,21].Furthermore the efficient mathematical analysis of their stochastic properties is an invaluable tool for their application.Two realms of investigation are generally of interest for such systems.The transient behaviour describes the time-dependent dynamics, whereas the stationary behaviour describes the dynamics in the long term after the system has reached an equilibrium.
Studying the dynamics of stochastic CRNs is difficult in general, and so they are often examined via simulations [9].The stationary behaviour and its characterization are typically analysed via the master equation.In many cases, the stationary behavior of Markov chains can be described through their stationary distribution.Exact solutions for the stationary distribution (if it exists) are not known for most systems, except for some special cases.Complex balanced reaction networks are fairly well understood by now.Deterministic complex balanced CRN have their stochastic counterparts with product-form stationary distributions of Poisson-type [1].The reverse statement is essentially also true: a stochastic CRN with product-form stationary distribution of Poisson-type (on any irreducible component) is complex balanced [4].Complex balanced CRNs are in particular weakly-reversible.Apart from that, there are some results on form of stationary distributions of non-weakly reversible reaction networks, like, e.g.autocatalytic CRN [15].
Here, we study unions (or, joins) of reaction networks in the stochastic setting.Our main focus is the form and existence of stationary distributions.While [15] focussed on a particular class of interest of non-weakly reversible CRNs with applications in particle systems, life sciences and condensation, we generalise here the underlying proof principle for stationary distributions.We give tools to systematically find the stationary distributions for the joined CRN, given the stationary distributions of the smaller CRNs.To be more precise, in CRNs where the stationary distributions of the decompositions are of product form and concur in the species in common, we can derive the stationary distribution of the full CRN from its parts.These are sufficient conditions, and examples can come from any combination of CRNs as long as the stationary distributions are of product form and satisfy some condition on the state spaces.Since the class expressible through such networks is big (i.e.interacting particle systems, cf., e.g., [19]), the principle also applies to other stochastic models.As an example, consider [15] for the relation to the inclusion process.In particular, autocatalytic CRNs and more general nonweakly reversible as well as some weakly reversible (including all complex balanced) CRNs fall under the framework we consider.
One result is then that given a reaction network G that can be decomposed as a reaction-disjoint union G = G 1 ∪ G 2 , with G 1 , G 2 essentials and of product form stationary distributions such that the product-form functions agree in the species in common, the stationary distribution of G is of product form and derived from G 1 , G 2 under a summability condition.
As an illustration consider the following CRN with Mass-action kinetics.
Then, taking as G 1 the reactions between S 1 , S 2 , and G 2 the reactions between S 2 , S 3 , we can apply the above to derive the product-form stationary distribution of G = G 1 ∪ G 2 for all positive rate constants.The stationary distribution is (see Example 5.1) where the product form functions are .
As the overall CRN G = G 1 ∪ G 2 is reversible and of deficiency two, such examples show that weakly reversible non-complex balanced CRNs can have product-form stationary distributions.
As another example consider the next CRN that can be decomposed in a complex balanced (reactions between S 3 , S 5 and between S 1 , S 3 ) and a join of a non-weakly reversible and a weakly reversible non-complex balanced CRN(the rest) with product-form stationary distributions (see Example 5.4) with f 1 , f 3 , f 5 of Poisson-form, f 2 of a form from autocatalytic CRNs and f 4 as f 1 , f 2 of the previous example.
Structure.In Section 2 we introduce basic definitions and terminology for reaction networks.Then we introduce the models for CRNs in Section 3 and focus on the stochastic model by reviewing definitions, properties, and results on stationary distributions, where at the end we introduce unions of CRNs.In Section 4 we study stochastic CRNs under joins and give some results on extending the stationary distributions from smaller CRNs to their joins.Section 5 introduces some examples to outline the application of the developed theory.
Relation to existing approaches.Previous approaches for extending analytical results on stationary distributions for reaction networks have focussed on gluing one state [21] or two states [22] of finite irreducible CTMCs.

Reaction networks
A reaction network G consists of a finite set of species S = {S 1 , • • • , S n }, a finite set of complexes, and a finite set of reactions R, which is then denoted as the triple G = (S, C, R).
We represent the complexes by vectors in Z n ≥0 , and write reactions as ν → ν ′ , where we assume ν, ν ′ ∈ C and ν = ν ′ for all ν → ν ′ ∈ R. For a reaction ν → ν ′ , ν is called the reactant and ν ′ the product.Every reaction ν → ν ′ has a positive rate constant κ ν→ν ′ associated.Then, given the vector of reaction rates κ ∈ R R >0 , we denote the CRN with rates by (G, κ).
2.1.Basic terminology.We illustrate reaction networks by their reaction graph, which is the directed graph obtained by taking the vertices C and arrows R. Connected components of the reaction graph are called linkage classes.
there is a directed path in the reaction graph that begins with ν ′ and ends in ν.The molecularity of a reaction ν → ν ′ ∈ R is equal to the number of molecules in the reactant |ν| = i ν i , and correspondingly we say such reactions are unimolecular, bimolecular, three-molecular or n-molecular.The stochiometric subspace spans a subspace of R n and is given as 3. Models and kinetics for reaction networks 3.1.Stochastic model.The progression of species counts is described by a vector X(t) = x ∈ Z n ≥0 , which changes according to the 'firing' of the reactions ν → ν ′ by jumping from x to x + ν ′ − ν with transition intensity λ ν→ν ′ (x).The Markov process with intensity functions λ ν→ν ′ : Z n ≥0 → R ≥0 can then be given by with the generator A acting by The transition intensity under mass-action kinetics (more general kinetics are possible as well [1,7]) associated to the reaction ν → ν ′ is where z! := n i=1 z i ! for z ∈ Z n ≥0 , and x ≥ ν if and only if this holds for every component, i.e. x i ≥ ν i ∀S i ∈ S.
General inquiry into stochastic CRNs proceeds by inspection of the underlying CTMC.After identifying the class structure and the (so-called) stoichiometric compatibility classes where the dynamics is confined to, the state space is decomposed into different types of states ( cf., i.e., [23]).On irreducible components, positive recurrence is equivalent to non-explositivity together with existence of an invariant distribution [23].
We next introduce some terminology for stochastic reaction networks.A reaction ≥0 if there is a state x ∈ A such that the reaction is active on x.This will mostly be used for ≥0 if it can be reached from x via the underlying CTMC.We will denote this by x → G u.
A non-empty set Γ ⊂ Z n ≥0 is an irreducible component of G if for all x ∈ Γ and all u ∈ Z n ≥0 , u is accessible from x if and only if u ∈ Γ.We say G is essential if the state space is a union of irreducible components, and G is almost essential if the state space is a union of irreducible components except for a finite number of states.

3.2.
Stationary distributions of reaction networks.Let X(t) denote the underlying stochastic process associated to a reaction network on an irreducible component Γ.Then, given that the stochastic process X(t) is positive recurrent and starts in Γ, we have that the limiting distribution is the stationary distribution, i.e.
In particular, if the underlying CTMC is positive recurrent, the stationary distribution π Γ on an irreducible component Γ is unique and describes the long-term behavior ( cf., e.g.[23]).
The stationary distribution is determined by the master equation of the underlying Markov chain: for all x ∈ Γ.A popular choice as rate function is mass-action kinetics, which then gives the following master equation: Solving equation (3.2) is in general a challenging task, even when restricting to the mass-action case.
Remark 3.1.Observe that for conservative CRNs the irreducible components are finite.Therefore the CTMC dynamics are positive recurrent (e.g., by Reuters criterion, c.f., e.g.[23]) on these irreducible components and the limiting distribution is the unique stationary distribution.Recall in particular that for infinite CTMCs existence of stationary distribution does not imply positive recurrence, cf., e.g.[23,Ex 3.5.4]or [6].

Known results on stationary distributions.
Studying transient and stationary behaviour of reaction networks are formidable tasks in general, and they are often examined via simulations [9].Analytical solutions for the stationary distribution (if it exists) are not known for most systems, except for some special cases.
Some stationary distributions of weakly reversible reaction networks are wellunderstood.Complex balanced CRNs have a Poisson product-form stationary distribution [1] and can even be characterized by that.For (G, κ) a complex balanced CRN and an irreducible component Γ, the stochastic system has product-form stationary distribution of the form where c ∈ R n >0 is a point of complex balance, c x := Si∈S c xi i , and M Γ is a normalizing constant.On the other hand, by [4,Theorem 5.1] any almost essential stochastic reaction network with product-form stationary distribution of Poissontype (i.e. in the form as above) is deterministically complex balanced.Notice that since complex balanced implies weakly reversible, these results do not apply to nonweakly reversible CRNs.Results on both product-form stationary distribution and connection to the deterministic system extend to non-mass action kinetics [1,7].Hence complex balanced CRNs are fairly well-understood.
For other classes of CRNs some results are also known [15], i.e. so-called autocatalytic CRNs, a class of non-weakly-reversible CRNs also have product-form stationary distributions.Their product form functions come from an infinite family of functions, where the first one specializes to the Poisson form as above.So for a autocatalytic CRN in the sense of [15, § 3], the stochastic dynamics has the product-form stationary distribution with product-form functions on its irreducible components (λ i and β k i are determined by the autocatalytic CRN, cf.[15, § 3]) and with Z Γ the normalising constant.Some other results on the stochastic behavior of CRN beyond complex balance are in [3] or [18].
Beyond these results little is known concerning explicit stationary distributions.

Balance equations for stationary distributions of CRNs.
We start with a general definition for balance equations, and recover some classical notions in Remark 3.3.The definition below essentially states that stationary distributions factorise according to a decomposition of the reactions of the underlying CRN [15].
Definition 3.2.Consider a CRN (G, κ) with stochastic dynamics on Γ and π a stationary distribution on Γ.We say (G, κ) is generalized balanced for π on Γ if there exists such that for all i ∈ A and all x ∈ Γ we have Remark 3.3.The notion of generalized balanced includes the following: (1) reaction balanced with index given by reactions, i.e. the tuples of subsets are {(ν → ν ′ , ν ′ → ν) ν→ν ′ ∈R } (2) complex balanced with index given by complexes, i.e. the tuples of subsets are defined for (3) reaction vector balanced with index given by a ∈ Z n , i.e. the tuples of subsets are defined for but also combinations and other possibilities.
In this paper, the following will be often used.
Remark 3.4.Let the reactions of a CRN be divided into R = R 1 ∪ R 2 , then it might happen that the stationary distribution factorises through these reactions.More formally this corresponds to generalised balance with Furthermore generalized balanced distributions on irreducible components give stationary distributions for the reaction network.Proposition 3.5.[15] If (G, κ) is a CRN with stochastic dynamics on Γ that is generalized balanced for π, then π is a stationary distribution for (G, κ) on Γ.

Unions of reaction networks.
Here we look at the operation of combining two reaction networks.Such operations were already introduced and studied in the deterministic setting in [12] where they studied the effects of combining reaction networks in the ODE setting with respect to identifiability, steady-state invariants, and multistationarity.While we will use the same framework, we study stationary properties of the stochastic model under combination and focus only on the two cases of reaction-disjoint and non-reaction-disjoint union.
We next introduce the formalisation of unifying reaction networks.
Definition 3.6.The union of reaction networks The union G 1 ∪ G 2 can be built under different assumptions between the underlying reaction networks G 1 , G 2 .The following implications holds [12]: Consider now taking the union of CRNs with rates We focus on the following two cases.(1) Gluing reaction-disjoint networks: If the two networks have no reactions in common (i.e., R 1 ∩ R 2 = ∅), then the vector of reaction rates of the union of the reaction networks is equal to (2) Gluing over reactions: If the two networks have at least one reaction in common (i.e., R 1 ∩ R 2 = ∅), then the rates of the reactions of the union of the networks which are common reactions (i.e. in If the two species sets are disjoint (S 1 ∩ S 2 = ∅), then the dynamics of the reaction networks G 1 and G 2 are independent of each other, hence some properties are directly determined by the dynamics on G 1 and G 2 (cf.Remark 4.1 for more on this in the stochastic case).
Remark 3.7.It is easy to see that both detailed balanced and complex balanced reaction networks are not closed under reaction-disjoint unions.Consider, e.g., the following example: with G 1 the part with two-molecular reactions and G 2 as the three-molecular reactions.The deficiency of G = G 1 ∪ G 2 is equal to one, hence for almost all parameters it will not be complex balanced.However, it is easy to check that both G 1 and G 2 are detailed balanced and hence complex balanced by themselves.

Stochastic reaction networks under joins
Notation.Let G = G 1 ∪G 2 be a reaction network obtained from a union of networks as in Definition 3.6.We denote the projections by where p Si is the projection to the i th component.

4.1.
Properties of stochastic dynamics under joins I.We first go through the case of a join where S 1 ∩ S 2 = ∅ for the sake of exposition and to introduce the reader to the setting.For notations on CTMCs in the context of CRNs we refer to § 3.1, or, e.g., [23].
The decomposition of state space with respect to irreducible components is simple.If Γ is an irreducible component of G, then p 1 (Γ) and p 2 (Γ) are irreducible components of G 1 , G 2 such that Γ = p 1 (Γ) × p 2 (Γ).So, for Γ a positive recurrent irreducible component we have where π 1 , π 2 are the stationary distributions on p 1 (Γ) and p 2 (Γ) of G 1 , G 2 (there is no normalizing factor since the CTMC is a product).It is easy to see that the stationary distribution on the irreducible component Γ is generalized balanced with {(R i , R i ) i∈{1,2} }(cf.Remark 3.4 and Theorem 4.9 for a proof of a generalisation).Remark 4.2.Even in the simplest setting of Remark 4.1 we can not say much concerning class structure of an x ∈ Z S ≥0 given only information about the classes of p 1 (x) for G 1 and p 2 (x) for G 2 (cf., e.g., the simple symmetric random walk on Z d ).In general, x is surely transient for We next establish some simple correspondences for the decomposition of the state space where we omit the proofs.
Lemma 4.3.The following are equivalent for a CRN G: (1) G is essential.
(2) For all x ∈ Z S ≥0 either there are no active reactions on x or we have that Since for G the following part of state space {z ∈ Z 3 |z W ≥ 0, z X = 0, z Y = 2} is not part of an irreducible component, G is not almost essential.In particular (C2) does not extend to almost essential.Remark 4.6.Even if G is essential, there might be no reaction-disjoint (or nonreaction disjoint) decomposition such that G = G 1 ∪ G 2 with G 1 , G 2 essential.As an example, consider, e.g., the following CRN which can be seen as a simple model for gene-expression [25].In this example the only essential subnetworks are ∅ ⇋ S 1 and the CRN itself.Also compare Lemma 4.7 to [24], which contains a similar result (written with different notions).Furthermore we need the following Lemma which follows by the definition of irreducible component.
(B2) Assume the same (i.e. as in (B1)) for G 2 , where we denote the stationary distribution on an irreducible component of G 2 by (B3) Assume there is an α > 0 such that for all x ∈ Γ and all S i ∈ S 1 ∩ S 2 we have B2) and (B3) are satisfied, then G = G 1 ∪G 2 has a product-form stationary distribution of the form where for is satisfied with solution (4.1) for i ∈ {1, 2}, which corresponds to generalized balanced with {(R i , R i ) i∈{1,2} }.Note that it is enough to prove it for R 1 .Then we are done by the symmetry of the assumption, and (4.1) is a stationary distribution, given it is summable.
We next prove that the master equation (4.2) holds true for reaction set R i = R 1 with solution (4.1).For x ∈ Γ by assumption p 1 (x) ∈ Z S1 ≥0 is in an irreducible component of G 1 .If this irreducible component is a singleton set, then the equation is trivially true.There are no active reactions of R 1 on x and the right side of (4.2) is zero.The left side of (4.2) is zero as well since these states are transient, i.e. the stationary distribution has no support (cf.Lemma 4.8).Hence assume it is a non-trivial irreducible component of G 1 , then inserting the proposed Ansatz (4.2) (modulo normalization) gives (4.3) Since the reactions in R 1 do not change the coordinates of S 2 \ S 1 , we have for all we can factor the equation as (4.4) By assumption Si∈S2\S1 f i (x i ) is nonzero (i.e. by contradiction with the assumption on the stationary distribution), so (4.3) is satisfied if the following holds: Now we identify the left and the right hand sides of the above equation with the corresponding sides of the master equation from G 1 with the projection p 1 (x) on the irreducible component.Since the transition rates of the reactions of R 1 only depend on the coordinates of S 1 , they are the same as the transition rates of the master equation from G 1 under p 1 (x) and we get an equality by assumption (B1).
• Theorem 4.9 assumes that the stationary distributions of G 1 , G 2 are of product-form.While this is a restriction, current results on form of stationary distributions are mostly in product-form (cf.[1,15]).Nonetheless, some examples with stationary distribution of non-product form are available [18, § 4.1] or [3], but calculating it or even writing it down in small examples is demanding.• By definition, p 12 (Γ) = p 21 (Γ), and condition (B3) requires the functions f 1 i , f 2 i with S i ∈ S 1 ∩ S 2 to be proportional on p Si (Γ) ⊆ Z ≥0 .• Notice that Theorem 4.9 assumes that the union comes from reactiondisjoint networks as in Definition 3.6.By the proof of Theorem 4.9 it would also hold for unions of reaction networks where we glue over reactions with Definition 3.6 (and similarly for its consequences, i.e., Theorems 4.12, 4.16, and Corollary 4.17).However, results on gluing over reactions are only a side product of the intended scope and does not seem very practical at the moment.We refer to Remark 4.21 for issues on applicability with respect to decomposing CRNs under gluing over reactions.
Remark 4.11.[Assumptions II] Note that assumption (B3) can be stated more general and Theorem 4.9 still holds with the same proof, i.e. in the following way: Assume there are constants α i for all S i ∈ S 1 ∩ S 2 with α i > 0 such that for all x ∈ Γ we have If this more general condition together with (B1), (B2) still holds, the conclusion of Theorem 4.9 is maintained with f i := f 2 i for S i ∈ S 1 ∩ S 2 .Furthermore it is easy to see that this does not influence the summability of (4.5).The same extension then follows for Theorem 4.16.
Furthermore, by the same proof as for Theorem 4.9 we can conclude the following for a slightly generalised setting (where on the irreducible component Γ if (4.5) is summable.

4.3.
Properties of stochastic reaction networks under joins II.We want to find sufficient conditions such that the projection p 1 (Γ) is a union of irreducible components of the stochastic dynamics of G 1 , which is a part of the assumption (B1) of Theorem 4.9.
Note that this holds in particular if G 1 is weakly reversible (cf.Lemma 4.7).
Proof.Let x ∈ Γ, and let p 1 (x) be the corresponding projected element.We have to show it is part of an irreducible component of G 1 .We distinguish the following two cases: • If there are no active reactions on p 1 (x), then by Lemma 4.3 p 1 (x) is not accessible from any other z ∈ Z S1 ≥0 , hence p 1 (x) is an irreducible component.
• Assume there are active reactions on p 1 (x).Then any other z ∈ Z S1 ≥0 is accessible from p 1 (x) if and only if p 1 (x) is accessible from this z by Lemma 4.3.Therefore the communicating class of p 1 (x) is closed.
Next we further investigate the conditions of the results of § 4.2 by focussing in particular on essential reaction networks G.
Proof.For (D1) =⇒ (D2) it suffices to observe that the projection p 1 is surjective, hence as Z S ≥0 is a union of irreducible components of G, we have that is a union of irreducible components of G 1 .In particular, G 1 is essential.(D2) =⇒ (D1) follows from Lemma 4.13.
Remark 4.15.Note that this implies in particular that an essential CRN G has a decomposition into G = G 1 ∪ G 2 with state space decomposition as in Theorem 4.9 for every irreducible component if and only if there is a decomposition with both and G 1 are essential, there might still be no such decomposition(cf.the example of Remark 4.6).
Hence, if G = G 1 ∪ G 2 can be decomposed such that G 1 and G 2 are essential, we know by Lemma 4.4 that G is essential.Furthermore by Lemma 4.13 the projections of irreducible components of G are decomposed into unions of irreducible components of G i , i ∈ {1, 2}.Therefore, in this case, we can restate Theorem 4.9 in a simplified form.Theorem 4.16.Let G = G 1 ∪ G 2 be a reaction network that can be decomposed as a reaction-disjoint union such that G 1 , G 2 are essential.Let Γ be an irreducible component of G. Assume that the irreducible components of p 1 (Γ) of G 1 (i.e.
and the irreducible components of Furthermore, assume that there is an α > 0 such that for all x ∈ Γ and all S i ∈ S 1 ∩ S 2 we have where for Then consecutive applications of Theorem 4.16 along decompositions of essential CRNs gives the following.
Corollary 4.17.Let G be a reaction network that can be decomposed as a reactiondisjoint union such that G = G 1 ∪ • • • ∪ G s with all the G j essential.Denote by S only j the species that are only in S j and no other S i , i = j, and by S shared j the species in S j that are also in at least one other S i , i = j.Assume that Γ is an irreducible component of G and each G j has product-form stationary distribution of the form on its irreducible component in p j (Γ) such that, if S j ∩ S k = ∅, then there is an α > 0 such that for all S i ∈ S j ∩ S k and all x ∈ Γ we have where if S i is in S j ∩ S k we choose f i := f j i arbitrary, such that the stationary distribution on Γ is generalized balanced with Remark 4.18.By the completeness of the results for complex balanced CRN (cf. § 3.3) it is clear we can not say more about complex balanced CRN.The same holds for a similar reason for autocatalytic CRNs since we generalise the underlying proof principle of [15], cf.Example 5.5.However, we offer a framework that can combine autocatalytic, complex balanced or other CRNs, as long as the stationary distributions are of product form and agree on the species in common.Note that it is easy to find small CRNs beyond complex balance with product form stationary distribution, and we cover only some.In particular there are both reversible, weakly reversible or non-weakly reversible CRN with product-form stationary distributions which can be combined in the framework we developed (cf.§ 5).
Remark 4.19.As another example consider which is not essential, hence there is no reaction-disjoint decomposition into G = G 1 ∪G 2 such that both are essential.Hence Theorem 4.9 still applies while Theorem 4.16 does not.
Remark 4.20 (Summability).Note that for summability of (4.1), (4.5) or (4.6) it is necessary that the stationary distributions on the projections are summable.In easy cases with infinite state space summability can possibly be checked by the ratio test, see Remark 5.2 for an example.
4.4.Decomposing CRNs and applications of Theorem 4.16.While general characterisations of existence for decompositions with a view towards Theorems 4.9, 4.12 ,4.16 are not our focus, we remark on several issues.Note that two reaction-disjoint CRNs which are given by their reactions R 1 , R 2 have a union with reaction set R = R 1 ∪ R 2 if and only if R can be decomposed as R = R 1 ∪ R 2 (i.e. with the same R 1 , R 2 ) by a reaction disjoint decomposition.Hence for a given set of reactions R, there are 2 |R| − 1 such bipartitions of the reactions, which grows exponentially with |R|.Correspondingly, brute-force algorithms can be given, e.g., for CRNs consisting of a complex balanced and an autocatalytic part.Therefore if stationary distributions for more classes of CRNs (even beyond product-form but as in Theorem 4.12) are known, this can be incorporated in a similar way for essential CRNs.Similarly decompositions along Proposition 4.17 can be checked.Note that as gluing over reactions is more general, some CRNs might be decomposable in that sense into essential CRNs where it is not possible for reaction-disjoint CRNs.Feasible strategies to cope with such situations are considerably more difficult than reaction-disjoint decompositions, but might be developed at a later point (see Remark 4.21).We further note that characterisations of when such decompositions exist are mostly unknown, and even in the essential case we currently only have characterisations for CRNs with stationary distributions given by Poisson product-form functions by [4].
Remark 4.21.While the above if and only if statement still holds for CRNs where we glue over reactions, the number of possible decompositions of a CRN G (cf. Definition 3.6) where we allow gluing over reactions is uncountable, which is not very practical.We henceforth mostly focus on decomposing along reaction-disjoint unions.

Applications and examples
We will next go through some examples in order to explain and illustrate the use of the theory developed.We mostly focus on mass-action kinetics in examples 5.1, 5.3, 5.4, 5.5 and 5.6, and consider Example 5.1 with general kinetics in Example 5.7.We conclude that many such weakly reversible CRNs of arbitrary deficiency have product-form stationary distribution independent of the rate and independent of the kinetics.While we only used the theory for CRNs, it applies to other stochastic networks and CTMCs as well.Furthermore recall that irreducible components of conservative CRNs are finite, hence the limiting distribution is the unique stationary distribution (cf.Remark 3.1).

Examples with Mass-action kinetics.
Example 5.1.As a first example consider the following CRN which is reversible and of deficiency two for an application of Theorem 4.16.
We first decompose G = G 1 ∪ G 2 into two essential CRNs: Then we analyse the stationary distributions of G 1 , G 2 on their own in order to apply Theorem 4.16 at the end.
Similar to the example of Remark 3.7, G 1 is only for some values detailed balanced.It has a stationary distribution of the form (see Remark 5.2) Next consider G 2 with stationary distribution (again with d 2 > 0) .
Now we look at G = G 1 ∪G 2 in order to apply Theorem 4.16.We choose d 1 = d 2 = 1 so that the product-functions agree.Then the stationary distribution of G is as follows, , where the product form functions are We further note that the summability in this example is trivial as the irreducible components are finite.
Remark 5.2.For G 1 of Example 5.1 observe the following • On an irreducible component with a product-form stationary distribution and a conservation relation, we will mostly factor out a constant d > 0.
As an example, consider G 1 where x 1 + x 2 is constant on the irreducible components Γ 1 i , so also d x1+x2 1 is a constant along irreducible components.Then as we divide by the normalising constant the corresponding stationary distributions are all the same for different d 1 > 0.
• G 1 is reaction vector balanced independently of the rates.We can verify that (5.1) is reaction vector balance (and hence the stationary distribution) for G by checking the following • For κ 3 = ακ 1 , κ 4 = ακ 2 with α > 0, G 1 is detailed (hence complex) balanced, and we can factorize out in f 1 from (5.1) to obtain x 2 ! .
To transform this into a standard form, we can choose d 1 = κ 2 .
• We can join G 1 with the following essential CRN G 2 x2! and where c 2 = κ+ κ− is a point of complex balance.Then choosing d 1 = c 2 makes the product-form functions f 1 2 , f 2 2 equal.Therefore if the following is summable, it is the stationary distribution where we have to check that the following sum is finite Therefore it is easy to see, e.g. by the ratio test for series, that the series converges for all positive rate parameters.
Then G 1 is autacatalytic and corresponds to reactions between S 1 , S 2 and G 2 is complex balanced and corresponds to the reactions between S 1 , S 3 .Hence G 1 , G 2 are essential, and we may apply Theorem 4.16 after deriving the stationary distributions of G 1 and G 2 , giving an easy way to compute the stationary distribution of [15,Example 4.4].This shows how to systematically decompose some examples of CRNs into smaller parts where the stationary distribution is known and of product-form.As another example consider the following CRN where we glue along two species.
Example 5.4.Let G 1 be the following essential CRN We can choose the parameters to obtain an autocatalytic CRN according to [15] on S 1 , S 2 , S 3 , and join it with the CRN on S 3 , S 4 (which was Example 5.1) with stationary distribution of the form Consider as G 2 the following complex balanced (hence weakly reversible, essential) CRN: with stationary distributions of the form x 5 !with (c 1 , c 3 , c 5 ) a point of complex balance.Then, if the rates of G 1 , G 2 are such that the product-form functions f 1 1 , f 1 3 and f 2 1 , f 2 3 can be chosen to be the same, we can give the stationary distribution of G = G 1 ∪ G 2 .
Example 5.5.Now we consider autocatalytic CRNs [15].Interacting particle systems of that form are used in inclusion processes from statistical physics and the modelling of ants and swarms [13,15,17].Consider a CRN G 1 on 2 species CRN with the reactions where m ≥ 1.Note that such CRNs are essential.Then obtaining the stationary distributions for such CRNs on two species and assembling them leads to the stationary distributions of autocatalytic CRNs from [15], which can also be obtained via the decomposition into joins with Theorem 4.16.
Example 5.6.While ergodic conservative CRNs with product-form stationary distributions have a degree of freedom to choose the product-form function, that is not the case for other ergodic CRNs.As an example consider which can be decomposed into two CRNs along Example 5.1 and Remark 5.2.Then, application of Theorem 4.9 requires that the parameters match in some sense and further summability also has to be taken care of.

5.2.
Examples with more general kinetics.We recall the setting of more general intensity functions from [1] which are given as (5.3) where the κ ν→ν ′ are positive reaction rates and θ i : Z → R ≥0 are such that θ i (x) = 0 if and only if x ≤ 0 (we use the convention that −1 j=0 a j = 1 for any {a j }).Typical kinetics used in mathematical biology are, e.g., x m , called Hill-type I/II in [26], where m is an integer and k, k 1 , k 2 are positive constants.The first specialises to stochastic Michaelis-Menten kinetics for m = 1 [1].Hence from Example 5.7 (also see Example 5.1) it is easy to see that we can assemble arbitrary CRNs of this form with product form stationary distribution independently of the rates via Theorem 4.16.This then gives the following.
Corollary 5.8.Independent of the kinetics (but with θ 2 fixed), any CRN that is a disjoint union of CRNs of the form S 2 ⇋ S i , 2S i ⇋ S 2 + S i for i = 2 has product-form stationary distribution independent of the rates.Remark 5.9.[Compatibility with complex balance in S 2 ] Let G 1 be a RN obtained from Corollary 5.8 and G 2 be a weakly reversible, deficiency zero CRN that is conservative with kinetic functions as θ 2 for species S 2 such that the only species in common between G 1 , G 2 is S 2 .Then, by [1, Theorem 6.1], Corollary 5.8 and Theorem 4.16 the CRN G 1 ∪ G 2 has product form stationary distribution with the product-form function in S 2 given by f 2 (x) = c x 2 x l=1 θ2(l) , where c ∈ R S2 ≥0 is a point of complex balance for G 2 .Remark 5.10.Note that it is usually not the case that RNs with different kinetics can be joined with matching product-form functions.This comes from the fact that for example if restricting to complex balanced CRNs with Mass-action kinetics and stochastic Michaelis-Menten kinetics, the corresponding product-form functions are different(see [1,Theorem 6.1]).

Lemma 4 . 4 .
Consider G = G 1 ∪G 2 as in Definition 3.6 and let x ∈ Γ be an element of an irreducible component Γ of G. (C1) If G is a join of reaction-disjoint networks(cf.Def.3.6), then the following holds: A reaction y → y ′ ∈ R 1 is active on x if and only if it is active on p 1 (x).(C2) If both G 1 and G 2 are essential, then their union G is essential.Remark 4.5.If G 1 is almost essential and G 2 is essential, their union G is not necessarily almost essential.As an example consider the following:

Lemma 4 . 7 .
We have the following implication for a reaction network G: G reversible =⇒ G weakly reversible =⇒ 3. of Lemma 4.3 holds for G.In particular, reversible and weakly reversible reaction networks are essential.

Lemma 4 . 8 .
Let G = G 1 ∪ G 2 be as in Definition 3.6 and consider an irreducible component

4. 2 .
Stationary distributions of joins of reaction networks.Here we will generalise the setting of Remark 4.1 in a direction where we can still deduce the form of a stationary distribution of the joined network G = G 1 ∪G 2 from the combinations of the stationary distributions of the separate reaction networks G 1 , G 2 .Notice that there are no conditions on the type of kinetics.Theorem 4.9.Let G = G 1 ∪ G 2 be a reaction network obtained from a union of reaction-disjoint networks as in Definition 3.6 with S 1 ∩ S 2 = ∅.Let Γ be an irreducible component of G. Consider the following assumptions: (B1) Assume p 1 (Γ) is a union of irreducible components of the stochastic dynamics of G 1 (i.e.p 1 (Γ) = ˙ i∈I Γ 1 i ) with stationary distributions on the irreducible components of the following form

Proof.
It suffices by Definition 3.2 and Proposition 3.5 to show that for any x ∈ Γ the master equation (4.2)

Proposition 4 . 14 .
Let G = G 1 ∪ G 2 be an essential reaction network.Then the following conditions are equivalent (D1) For every irreducible component Γ of G, the projection p 1