Stochastic approximations of higher-molecular by bi-molecular reactions

Reactions involving three or more reactants, called higher-molecular reactions, play an important role in mathematical modelling in systems and synthetic biology. In particular, such reactions underpin a variety of important bio-dynamical phenomena, such as multi-stability/multi-modality, oscillations, bifurcations, and noise-induced effects. However, as opposed to reactions involving at most two reactants, called bi-molecular reactions, higher-molecular reactions are biochemically improbable. To bridge the gap, in this paper we put forward an algorithm for systematically approximating arbitrary higher-molecular reactions with bi-molecular ones, while preserving the underlying stochastic dynamics. Properties of the algorithm and convergence are established via singular perturbation theory. The algorithm is applied to a variety of higher-molecular biochemical networks, and is shown to play an important role in synthetic biology.

(1973), Cook et al. (1989), and Wilhelm (2000) in context of reactions of arbitrary order at the stochastic level.
The paper is organized as follows. In Sect. 2, we prove that any one-species thirdorder reaction can be approximated with a suitable family of second-order networks, and we apply the results in Sect. 3 on the Schlögl system (Schlögl 1972). In Sect. 4, we generalize the results from Sect. 2 to arbitrary multi-species higher-order reactions under mass-action kinetics. In Sect. 5, we apply the generalized results to higher-order reaction networks displaying noise-induced phenomena. Finally, we conclude with a summary and discussion in Sect. 6. The notation and background theory used in the paper are introduced as needed, and are summarized in Appendix 1. More detailed analyses underlying Sect. 4 are provided in Appendices 2-3.

Special case: one-species third-order reactions
Let us consider an arbitrary input reaction network R 0 = R 0 (X ), under mass-action kinetics (Feinberg 1979), involving a single biochemical species X , given by where ν j ∈ Z ≥ andν j ∈ Z ≥ are the reactant and product stoichiometric coefficients of the j-reaction, respectively, while k j ∈ R > is the corresponding dimensionless rate coefficient; Z ≥ and R > are the sets of nonnegative integers and positive real numbers, respectively. See also Appendix 1 for notation and reaction network theory. Here, (3X k 1 − →ν 1 X ) is an arbitrary one-species third-order (tri-molecular) reaction, which we wish to approximate, while R ρ = R ρ (X ), called the residual network, contains the remaining reactions from R 0 that we do not wish to approximate. Let us now consider the output mass-action network R ε = R ε (X , Y ), containing an auxiliary species Y , given by where we denote the irreversible forward and backward reactions (2X κ 1
Remark. In the original coordinates (x, y), operator L 0 corresponds to the reaction Y 1 − → 2X , and equation (6) has infinitely many solutions. This degeneracy arises from the fact that the process induced by L 0 satisfies a local linear conservation law x + 2y =x, wherex is time-independent. Using this conservation law as a coordinate change leads to equation (6), wherex is only a parameter, with the solution p 0 (y) = δ y,0 which is unique up to a multiplicative (x, τ )-dependent constant. Remark. Let l K 0 = {p : Z ≥ → R | p(y) = 0 ∀y ≥ K } be the space of sequences whose elements beyond K ≥ 0 are zero. If L 0 is as given in (4) and if f ∈ l K 0 , then any solution of L 0 p = f satisfies p ∈ l K 0 . Therefore, since p 0 (x, ·, τ ) ∈ l 1 0 , it follows that (7) and (8) are both finite-dimensional systems of linear equations with respect to the variable y for any parameter choice (x, τ ). Equation (7). Using (4) and (9), it follows that Considering the form of L 0 and (10), we seek a solution of (7) in a separable form, Substitution into (7) leads to the general solution Equation (8). As remarked previously, (8) is a finite-dimensional system of linear equations; in particular, p 2 (x, ·, τ ) ∈ l max{2,γ }+1 0 ; therefore, the Fredholm alternative theorem holds (Kreyszig 1989). The adjoint (backward) operator corresponding to L 0 is given by L * 0 = y E −1 y − 1 , and its null-space N is given by N (L * 0 ) = {1}; see also Appendix 1. Therefore, it follows from the Fredholm alternative theorem that (8) has a solution if and only if the solvability condition 1, RHS y = 0 holds, where 1, · y = y∈Z ≥ ·, and RHS denotes the right-hand side of (8); substituting (4), (9), and (11) into the solvability condition, one obtains the effective CME where L ρ is the forward operator induced by the residual network R ρ . Remark. In the original coordinate t, the Fredholm alternative theorem applied to (7) enforces a trivial effective CME d/dt p 0 (x, t) = 0; to capture non-trivial dynamics, we have rescaled time to a longer scale. Remark. Equation (12) describes a time-evolution of the PMF for the stochastic process and not the original copy-number X (t). However, (9) implies that process Y (t) spends most of the time at y = 0 as ε → 0, so that the PMFs for X (t) and X (t) match as ε → 0.

Kinetic and stoichiometric conditions
To ensure that the dynamics of the output network (2) matches that of the input network (1), coefficients κ 1 and κ 2 have to suitably scale with ε. In particular, the CME for (1) reads In order for (12) and (13) to match, we impose the kinetic condition, given by and the stoichiometric condition, given bỹ where o(·) is the "little-o" asymptotic symbol, see also Appendix 1. (14) respectively ensure that the operators L 1 and L 2 remain slower than ε −1 L 0 in (3); otherwise, in the degenerate case when is the "big-O" symbol, one obtains families of perturbation problems distinct from the one considered in this section.

Convergence
The formal perturbation analysis from Sect. 2.1 has been performed under the assump- , and that κ 1 and κ 2 are independent of ε, which is inconsistent with the kinetic condition (14). We stress that an objective of the analysis in Sect. 2.1 was precisely to uncover admissible ε-scalings of κ 1 and κ 2 which ensure that (1) and (2) match. Having formally obtained such candidates, we now perform a convergence analysis, without the aforementioned assumptions, under a particular scaling which satisfies (14), given by whereκ 1 ,κ 2 are suitable ε-independent parameters. The CME for network (2) under (16) reads Substituting into (17) the fractional-power perturbation series Since systems (6)- (8) and (19) have the same form, the same is true for their solutions.
Remark. Constant c = c(K, Sx × S y , T ) appearing in (22) increases linearly with T (see (26)), i.e. condition ε 1/T 2 is sufficient for achieving accuracy for all t ∈ [0, T ]. In Sect. 3, we demonstrate with numerical simulations that, while sufficient, the condition ε 1/T 2 is not necessary for the error to satisfy a bound of the form (22). In particular, we present an example biochemical network for which lim t→∞ p ε (·, ·, t) − p 0 (·, ·, t) l 1 (Sx ×S y ) = O(ε 1/2 ), showing that even stationary PMFs can obey an error bound of the form (22). Remark. Result (22) remains valid if for all t ∈ [0, T ] the rate coefficientsκ =κ(t) and k ρ = k ρ (t) are nonnegative-valued functions of time with continuous first derivatives. Remark. One can derive error bounds analogous to (22) for scalings other than (16), which are consistent with (14); however, note that such scalings lead to systems of perturbation equations whose form is not the same as (19), see also Appendix 2. In particular, one can show that, under the fractional-power scaling κ 1 =κ 1 ε −n/d and κ 2 =κ 2 ε −(1−n/d) , with n, d ∈ Z > and n/d < 1, the error can be asymptotically bounded by cε 1/d for some c > 0; scaling (16) is a special case obtained by taking n = 1 and d = 2.
The joint-PMF error estimate (22) holds for every choice of the rate coefficients κ 1 andκ 2 , and for every choice of the stoichiometric coefficientsν andγ . Under the particular choices (14) and (15), Proposition 2.1 implies the following marginal-PMF error estimate. In what follows, we let p x ε (x, t) ≡ 1, p ε (x, ·, t) be the x-marginal PMF of network (2).

Multi-reaction approximations
The results from Sects. (2.1) and (2.2) have been achieved when a single input trimolecular reaction (3X k 1 − →ν 1 X ) is replaced by the output bi-molecular network R ε 1 (X ; Y , κ 1 ) ∪ R 2 (X ; Y , κ 2 ,ν,γ ) involving one auxiliary species Y . Assume now that we wish to approximate multiple input tri-molecular reactions, M j=1 (3X Then, performing analogous perturbation analysis as in Sects. (2.1)-(2.2), one can readily prove that analogous convergence results hold if the M input reactions are replaced by 3M output reactions M j=1 More efficiently, one can readily prove that the same convergence results also hold if the M input reactions are replaced by

Example 2.1 Consider the third-order input network
By applying the algorithm (2) independently to each of the two reactions from (29), one obtains the second-order output network The kinetic and stoichiometric conditions (14) and (15) are respectively given by: Output network (30) contains 6 reactions and 2 auxiliary species Y 1 and Y 2 . More efficiently, an alternative output network is given by and contains 4 reactions and 1 auxiliary species Y . The kinetic and stoichiometric conditions read

Example: The Schlögl network
In this section, we apply the results developed in Sect. 2 to the one-species third-order Schlögl network (Schlögl 1972), given by where ∅ represents species or processes that are not explicitly modelled. In Fig. 1a, we display as a black curve the stationary PMF p 0 = p 0 (x; k 1 , (k 2 , k 3 , k 4 )) for the input network (34) under a particular choice of the rate coefficients, which has been obtained by numerically solving the underlying stationary CME; one T. Plesa

Fig. 2
Panel a displays the stationary PMF of the input network (34) as a black curve, and the x-marginal PMF for the output network (35), under (37) with ε = 10 −2 , as a blue histogram. Panel b displays a log-log plot of the l 1 -distance between the PMFs for networks (34) and (35) as a function of ε. Panel c displays the error as a function of the scaling factor s for three different values of ε, and with (κ 1 ,κ 2 ) and (k 1 , k 2 , k 3 , k 4 ) fixed as in Fig. 1 can notice that p 0 displays two maxima (bi-modality). Approximating the third- (2), while preserving the residual network The stoichiometric condition (15) demands thatν = (2 − 2γ ), and there are two choices: takingγ = 0 implies thatν = 2, takingγ = 1 implies thatν = 0, while takingγ ≥ 2 implies thatν < 0, which is biochemically infeasible. In what follows, we take (ν,γ ) = (0, 1), and consider different scaling factors s to satisfy the kinetic condition (14). Let us first satisfy (14) by setting In Fig. 1a, we display the stationary x-marginal PMF of the output network (35) under (36) for different values of the parameter ε. In particular, when ε = 10 −2 , the PMF is shown as a dashed purple curve; while bi-modal, this intermediate PMF is inaccurately distributed. On the other hand, when ε = 10 −6 , the PMF p x ε is shown as a blue histogram, and is in an excellent match with target PMF p 0 . In Fig. 1b, we show a log-log plot of a numerically approximated error p x ε − p 0 l 1 as a function of ε. Also shown, as a dashed blue line, is the reference curve p x ε − p 0 1 = ε 1/2 ; one can notice an excellent match in the slopes of the two curves, in accordance with the finite-time result (27) from Corollary 2.1. In Fig. 2c, we display the stationary y-marginal PMF for network (35) when ε = 10 −6 , which is shown in Fig. 2d to converge to the Kronecker-delta function centered at zero in accordance with the finite-time result (28).
Corollary 2.1 provides information about the error p x ε − p 0 l 1 in the limit ε → 0. Let us now discuss how one may decrease this error for a fixed ε by choosing an appropriate scaling factor s. To this end, note that network (2) consists of an ordered chain of reactions: in order for R 2 to fire, and mimic (1), one requires that the forward reaction from R ε 1 fires first. The reactant of the forward reaction from R ε 1 , forming the start of the chain, is given by 2X , and the propensity function is given by α 1 (x) = κ 1 x(x − 1). On the other hand, R 2 involves as a reactant the auxiliary species Y , with the propensity function α 2 (x, y 1 ) = κ 2 x y. Since Y (t) spends most of the time at y = 0 for sufficiently small ε, it follows that the underlying joint-PMF is concentrated in the neighborhood of the x-axis, and α 2 (x, y)/κ 2 < α 1 (x)/κ 1 . This observation suggests that, for a fixed smaller ε, there is an optimal ratio κ 1 /κ 2 , sufficiently small to speed up reaction R 2 , and sufficiently large to ensure that R ε 1 is triggered often enough. To this end, let us now satisfy (14) by setting In Fig. 2a, we display the stationary x-marginal PMF of the output network (35) under (37) (36), with the former occurring at a rate ε 1/3 ; see also a remark below Proposition 2.1. In Fig. 2c, we display the l 1 -distance between the input and output PMFs as a function of the scaling factor s for three different values of the parameter ε. One can notice that the error appears to be minimized approximately at s = 3/10 when ε = 10 −2 , and that, for larger values of ε, an overall better performance is achieved by taking s < 1/2. Figure 2c also suggests that the error does not converge to zero at the degenerate points s = 0 and s = 1.

General case: Multi-species higher-order reactions
Let us consider an arbitrary input mass-action reaction network R 0 = R 0 (X ) involving N biochemical species X = {X 1 , X 2 , . . . , X N }, given by In particular, we assume that the first reaction involves m ≥ 1 distinct reactants {X 1 , X 2 , . . . , X m } ⊆ X , which we enforce by demanding that the underlying reactant stoichiometric coefficients are nonzero; we also assume that the first reaction is of order n ≥ 3. Furthermore, for convenience, we assume that the reactant species in the first reaction are ordered according to nondecreasing stoichiometric coefficients, i.e. ν 1,i ≤ ν 1, j if i < j, for all i, j ∈ {1, 2, . . . , m}. We wish to approximate the n-th order target (first) reaction from (38) with a suitable second-order network, while leaving the residual network R ρ unchanged, analogously as in Sect. 2. To this end, consider the output mass-action reaction network R ε = R ε (X , Y) given by where ∅ is the empty set, and with the sub-networks . . , Y n−2 }, and consists of (2n − 3) reactions, (n − 2) of which are firstorder, and (n − 1) of second-order. Reaction network (2) from Sect. 2 is a special case of (39)-(40) with m = N = 1 and n = 3.

Kinetic and stoichiometric conditions
In Appendix 2, we generalize the formal perturbation analysis from Sect. 2.1 and derive kinetic and stoichiometric conditions, analogous to (14) and (15), respectively, which ensure that the CMEs for the input network (38) and the output network (39)-(40) match. In particular, the generalized kinetic condition is given by Requirement (41) states that the product of the rate coefficients of the slower reactions from (40), n−1 i=1 κ i , divided by the product of the rate coefficients of the faster reactions, 1/ε n−2 , must be equal to the rate coefficient of the target reaction from (38), k 1 . On the other hand, when the reaction R n−1 (X j ) from (40) does not contain the auxiliary species {Y 1 , Y 2 , . . . , Y n−3 }, then the generalized stoichiometric conditions are given bỹ The stoichiometric conditions valid when (γ 1 ,γ 2 , . . . ,γ n−3 ) = 0 take a more complicated form, and can be obtained algebraically as explained in Appendix 2. One can also readily obtain the stoichiometric conditions graphically, as we now outline via an example.
= (2X 1 + X 2 )) to the products in (47) as many times as desired, as long as the resulting complex contains nonnegative stoichiometric coefficients. For example, by adding the complex On the other hand, adding ∅ . = (Y 2 − 2X 1 − X 2 ) three times to the products in (47) leads to which is not a biochemical reaction, as the product complex is not nonnegative.
Remark. The graphical approach taken in Example 4.1 applies generally: one can extract the formal equalities, such as ∅ .
, directly from the fastest reactions in R ε . Writing the final reaction from R ε with the same product complex as in the first reaction from R 0 , one can then add the formal equalities as many times as desired to the products of the final reaction from R ε , provided the resulting complex remains nonnegative.
Remark. Theorem 4.1 also holds under more general choices of {ν i } m i=1 and {γ i } n−2 i=1 , which have been outlined in Example 4.1; see also Appendices 2 and 3. In particular, the fact that a given higher-order input network may be approximated by a variety of different second-order output networks is favorable, since there is a greater flexibility to meet various biochemical constraints which may be necessary for successful experimental implementations; see also Example 4.2. Remark. Algorithm (39)-(40) and Theorem 4.1 extend naturally to the case when multiple higher-order input reactions are approximated by second-order ones; see Sections 2.3 and 5.2. Remark. Analogous remarks as those under Proposition 2.1 also apply to Theorem 4.1. In particular, constant c = c(K, S x × S y , T ) from (53) increases linearly with time T , so that a sufficient condition for achieving accuracy for all t ∈ [0, T ] is that ε 1/T n−1 ; hence, to meet this sufficient condition for a fixed T , the higher the order of the target reaction, n, the smaller the asymptotic parameter, ε, must be chosen. Furthermore, note that (53) remains valid if for all t ∈ [0, T ] the rate coefficientsκ = κ(t) and k ρ = k ρ (t) are continuously differentiable nonnegative-valued functions of time.
To formulate Theorem 4.1, we have assumed a fixed ordering of the reactants and reactions in (39)-(40). One can readily prove analogous results for other suitable orderings.

Example 4.2 Consider the third-order input reaction
Output network (39)-(40) is given by in particular, the forward reaction from R ε 1 is a second-order hetero-reaction, involving two distinct reactants X 1 and X 2 . One can readily show that the results presented in this section also hold for the output network R ε 1 (X 2 , X 2 ) ∪ R 2 (X 1 ), given by for which the forward reaction from R ε 1 is a second-order homo-reaction, involving only X 2 as reactants. Another admissible output network is R ε for which the forward reaction from R ε 1 is of first-order. Using analogous perturbation analysis to that underpinning Theorem 4.1, one can show that the (x 1 , x 2 )-marginal PMF of (56) also converges to the PMF of (54) at a rate ε 1/2 as ε → 0; the same convergence occurs for the output network (57) at a rate ε 1/3 , since two, and not only one, auxiliary species are introduced. Depending on the application area, a particular output network might be more desirable than others; for example, given a set of molecular species {X 1 , X 2 , Y 1 } with predefined biophysical properties, it may be easier to experimentally realize a reaction of the form 2X 2 → Y 1 than X 1 + X 2 → Y 1 .

Examples: Noise-induced phenomena
In this section, we apply the results from Sect. 4 to two test networks arising from theoretical synthetic biology and displaying noise-induced phenomena that are absent at the deterministic level. In particular, the first network, given by (58), plays an important role in the stochastic morpher controller (Plesa et al. 2021) that can globally morph the PMF of a suitable reaction network into any desired form. The second network, given by (63), is part of the noise-control algorithm (Plesa et al. 2018) that can redesign a given reaction network to locally reshape the underlying PMF in a mean-preserving manner.

Biochemical Kronecker-delta function
Let us consider the fourth-order mass-action input reaction network Long-time PMF of (58), under a particular choice of the rate coefficients k 1 < k 2 , is shown in Fig. 3a as black dots interpolated with solid lines. The PMF is close to the Kronecker-delta function centered at x = 3. In particular, when there are less than four molecules of X present, x < 4, only the first reaction from (58) fires and X experiences a constant positive drift until four molecules are present. When x ≥ 4, both reactions from (58) fire, with the second one, having a larger propensity function, overpowering the first one and generating a net-negative drift. The combined effect of the two reactions forces X to spend most of the time at the state x = 3.  (58) with (k 1 , k 2 ) = (10 −3 , 10 −2 ) as black dots interpolated with solid lines; the long-time x-marginal PMF of the output network (59), with (ν,γ 1 ,γ 2 ) = (0, 0, 1), under rate coefficients (62), is shown as a blue histogram when ε = 10 −3 . Panel b displays a log-log plot of the l 1 -distance between the long-time PMFs for networks (58) and (59). Panel c displays the long-time PMF of the input network (63)-(64) with (k 1 , k 2 , k 2,5 ,k 4,2 ) = (1, 1, 1, 1), and the x-marginal PMF of the output network (65)-(66), under rate coefficients (68) with β = 1/12, for two different values of ε. The conservation constant for (63) is fixed to c = 7, and (65) is initialized with zero copy-numbers of the species {Y i } 5 i=1 andỸ 4 . The plots have been generated using the Gillespie algorithm (Gillespie 1977) Applying the algorithm (39)-(40) on the fourth-order input network (58), one obtains a suitable second-order output network given by The stoichiometric condition (42) for network (59) withγ 1 = 0 is given bỹ we fixγ 2 = 1, so thatν = 0. On the other hand, the kinetic condition (41) for (59) reads which, using (52) with e.g. (κ 1 ,κ 2 ,κ 3 ) = (k 2 , 1, 1), is satisfied with In Fig. 3a, we display the long-time x-marginal PMF of the output network (59) with (ν,γ 1 ,γ 2 ) = (0, 0, 1) and the rate coefficients (62) with ε = 10 −3 , which is in an excellent agreement with the input PMF. In Fig. 3b, we show that the output PMF converges to the input one at a rate ε 1/3 , consistent with Theorem 4.1.

Noise-induced tri-modality
In this section, we demonstrate how the approach from Sect. 4, developed to approximate a single higher-order target reaction from an input network, generalizes to the case of multiple target input reactions; analogous ideas have been discussed in Sect. 2.3 for third-order reactions. In particular, one can apply the algorithm (39)-(40) to each of the target input reactions independently; however, such an approach may lead to output networks with larger number of reactions and auxiliary species Y, which may be biochemically expensive to engineer. More efficiently, if some of the higher-order input reactions involve common reactant sub-complexes, then some of the intermediate species Y can be re-used to simultaneously approximate multiple input reactions.
To illustrate these ideas, consider the seventh-order input network R 0 = R 0 (X 1 , X 2 ) given by where Species X 1 and X 2 from (63) are conserved, x 1 + x 2 = c; in what follows, we fix the conservation constant to c = 7. Network (63) has been obtained by applying the noise-control algorithm (Plesa et al. 2018) on the residual network R ρ ; in particular, sub-networks R 2,5 andR 4,2 , called zero-drift networks, introduce a state-dependent noise and decrease the PMF of R ρ at x 1 = 2 and x 1 ∈ {4, 5}, respectively, while preserving the underlying mean. In Fig. 3c, we display the long-time PMF of (63) as black dots interpolated with solid lines. One can notice that the network displays noise-induced tri-modality, with the modes x 1 ∈ {1, 3, 6}. Applying algorithm (39)-(40) to each of the four reactions from R 2,5 ∪R 4,2 independently requires 18 auxiliary species and 40 reactions in total, i.e. 5 auxiliary species and 11 reactions for each of the reaction from R 2,5 , and 4 auxiliary species and 9 reactions for each of the reaction fromR 4,2 . However, since both reactions from R 2,5 involve the same reactants, one can reduce their order simultaneously by using 5 auxiliary species; similarly, 4 auxiliary species suffice to reduce the order of both reactions fromR 4,2 . Furthermore, all of the reactions from R 2,5 ∪R 4,2 involve a common reactant sub-complex, namely (2X 1 +2X 2 ), so that, instead of using 9 auxiliary species for R 2,5 ∪R 4,2 , one can use only 6. These considerations give rise to the output network where The sub-network R 2,5 from (64) 3 encodes the common sub-complex (2X 1 +2X 2 ). Instead of applying (39)-(40) independently to each of the reaction from R 2,5 ∪R 4,2 , which requires 18 auxiliary species and 40 reactions, we have achieved the same goal in (65) with 6 auxiliary species and 16 reactions.

Discussion
In this paper, we have shown that, by introducing auxiliary species (dimension expansion) and suitable time-scaled-separated reactions, any higher-order input reaction network can be mapped to a second-order output one, with the underlying stochastic dynamics being preserved. This order-reduction algorithm has been previously formally established at the deterministic level for third-and fourth-order reactions (Tyson 1973;Cook et al. 1989;Wilhelm 2000). In this paper, we have generalized this algorithm to reactions with arbitrary number and composition of reactants, and we have augmented our formal results with rigorous convergence analyses at the stochastic level. In particular, we have shown that the time-dependent probability distributions of the input and output networks are arbitrarily close over suitable bounded domains in an appropriate asymptotic limit of some of the underlying rate coefficients. In Sect. 2, we have shown that an arbitrary one-species input reaction of order n = 3, given in (1), can be approximated by a family of second-order output networks, given in (2), provided the kinetic and stoichiometric conditions (14) and (15), respectively, are satisfied. Convergence for a family of the output networks has been proved in Corollary 2.1. In Sect. 3, we have numerically verified the results from Sect. 2 on the Schlögl network (34). In Appendices 2-3, the results from Sect. 2 have been generalized to arbitrary multi-species reactions of order n ≥ 3, and these results have been presented in Sect. 4. In particular, we have shown that an arbitrary multi-species input reaction of order n ≥ 3, given in (38), can be approximated with a family of second-order output networks, given in (39)-(40), provided the kinetic and stoichiometric conditions (41) and (42), respectively, hold. Convergence for a particular family of output networks has been established in Theorem 4.1, where we have shown that, for an input reaction of order n ≥ 3, the order of convergence is given by (n − 1) −1 ; hence, the higher the order of the input reaction, the slower the convergence. In Sect. 5, we have applied the results from Sect. 4 to the fourth-and seventh-order input networks (58) and (63), respectively, arising from theoretical synthetic biology (Plesa et al. 2018(Plesa et al. , 2021, and displaying noise-induced phenomena. The results established in this paper may play an important role in synthetic biology, and particularly in nucleic-acid-based synthetic biology, also known as DNA computing (Zhang and Winfree 2009). In this context, it has been proved that, assuming one can experimentally vary reaction rate coefficients over a sufficiently large range, any abstract reaction network of up to second-order, under mass-action kinetics, can be experimentally compiled into a physical second-order network with DNA molecules, with the underlying deterministic dynamics being preserved over compact time-intervals (Soloveichik et al. 2010). This molecular compiler has been proved to also preserve the underlying stochastic dynamics (Plesa et al. 2018). In this context, results from Sect. 4, and Theorem 4.1 in particular, imply the following corollary.
Corollary 6.1 (Universal molecular compiler) Assume that the reaction rate coefficients in the DNA compiler from Soloveichik et al. (2010) can be varied over arbitrarily large range. Then, any mass-action input reaction network, of arbitrary order, can be compiled into a second-order DNA-based output network with the compiler from Soloveichik et al. (2010), in such a way that the probability distributions for the input and output networks are arbitrarily close over sufficiently large compact state-spaces and any compact time-interval.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: Background
Notation. Union and intersection of sets A 1 and A 2 are denoted by A 1 ∪ A 2 and A 1 ∩ A 2 , respectively. The empty set is denoted by ∅. Set R is the space of real numbers, R ≥ the space of nonnegative real numbers, and R > the space of positive real numbers. Similarly, Z is the space of integer numbers, Z ≥ the space of nonnegative integer numbers, and Z > the space of positive integer numbers. When convenient, Euclidean vector x : {1, 2, . . . , N } → R is displayed in the column-form , where x i = x(i) and · denotes the transpose operator; the zero vector is given by 0 = (0, 0, . . . , 0) ∈ R N . Given any two sequences p, q : Z N ≥ → R, we define a bilinear form p, q ≡ x∈Z N ≥ p(x)q(x); the l 1 -norm of p is given by

Biochemical reaction networks
We consider reaction networks R = R(X ) firing in well-mixed unit-volume reactors under mass-action kinetics (Feinberg 1979), involving N biochemical species X = {X 1 , X 2 , . . . , X N } interacting via M reactions, given by Here, k j ∈ R > is the rate coefficient of the j-reaction, and we let k ≡ (k 1 , k 2 , . . . , k M ) ∈ R M > . Integers ν j,l ,ν j,l ∈ Z ≥ are the reactant and product stoichiometric coefficients of the species X l in the j-reaction, respectively, and we let ν j ≡ (ν j,1 , ν j,2 , . . . , ν j,N ) ∈ Z N ≥ andν j ≡ (ν j,1 ,ν j,2 , . . . ,ν j,N ) ∈ Z N ≥ . If ν j = 0 (respectively,ν j = 0), then the reactant (respectively, product) of the j-reaction is the zero species, denoted by ∅, representing species that are not explicitly modelled or non-biochemical processes. When convenient, we denote two irreversible reactions The order of j-reaction from network R is given by 1, ν j ∈ Z ≥ . The order of reaction network R is given by the order of its highest-order reaction; R of order higher than two is said to be a higher-order network.

Stochastic model of reaction networks
Under suitable conditions, copy-numbers of the biochemical species from (69), denoted by X(t) = (X 1 (t), X 2 (t), . . . , X N (t)) ∈ Z N ≥ , where t ∈ R ≥ is the timevariable, can be modelled as a continuous-time discrete-space Markov chain . The probability-mass function (PMF) p(x, t), i.e. the probability that the copynumber vector X(t) ∈ Z N ≥ at time t > 0 is given by x ∈ Z N ≥ , satisfies a partial difference-differential equation, called the chemical master equation (CME) (Erban et al. 2019;Van Kampen 2007), given by Here, the step operator E Linear operator L : D(L) ⊆ l 1 → l 1 , defined in (70)
Let us write the solution of (75) as a perturbation series Substituting (79) into (75), and equating terms of equal powers in ε, and defining p −1 (x, y, τ ) ≡ 0, one obtains the following system of n equations: Order 1/ε n−1 equation. Operator L 0 , given in (76), acts and depends only on the variable y, and each summand L i 0 is multiplied on the right by a factor y i ; therefore Order 1/ε n−2 equation. Substituting (81) into the right-hand side, one obtains where we use the fact that each of the operators {L i } n−1 i=2 from (77) is multiplied on the right by a nonconstant factor y i−1 , and we use the identity y i δ y i ,0 ≡ 0.
Since the right-hand side of (94) is a linear combination of a particular solution of the inhomogeneous linear equation (91) and the general solution of the underlying homogeneous equation, it follows that (94) is the general solution of (91).

Kinetic and stoichiometric conditions
In order for the effective CME (103) to match the CME of the input network (38), the kinetic condition (41) must hold. Furthermore, we require that x = x 1 = (ν 1 −ν 1 ).