1407/1407.8316.md

Validity conditions for moment closure approximations in stochastic chemical kinetics

David Schnoerr ^1,2, Guido Sanguinetti ², Ramon Grima ¹

¹ School of Biological Sciences, University of Edinburgh, UK

² School of Informatics, University of Edinburgh, UK

Abstract

Approximations based on moment-closure (MA) are commonly used to obtain estimates of the mean molecule numbers and of the variance of fluctuations in the number of molecules of chemical systems. The advantage of this approach is that it can be far less computationally expensive than exact stochastic simulations of the chemical master equation. Here we numerically study the conditions under which the MA equations yield results reflecting the true stochastic dynamics of the system. We show that for bistable and oscillatory chemical systems with deterministic initial conditions, the solution of the MA equations can be interpreted as a valid approximation to the true moments of the CME, only when the steady-state mean molecule numbers obtained from the chemical master equation fall within a certain finite range. The same validity criterion for monostable systems implies that the steady-state mean molecule numbers obtained from the chemical master equation must be above a certain threshold. For mean molecule numbers outside of this range of validity, the MA equations lead to either qualitatively wrong oscillatory dynamics or to unphysical predictions such as negative variances in the molecule numbers or multiple steady-state moments of the stationary distribution as the initial conditions are varied. Our results clarify the range of validity of the MA approach and show that pitfalls in the interpretation of the results can only be overcome through the systematic comparison of the solutions of the MA equations of a certain order with those of higher orders.## I. INTRODUCTION

The chemical master equation (CME) is the well accepted mesoscopic description of chemical systems in well-mixed and dilute conditions [1]. However, for most systems, analytic solutions are unknown. The stochastic simulation algorithm (SSA [2]) is a popular Monte Carlo method for sampling from the probability distribution of the CME, but the SSA is computationally expensive and is tractable only for chemical reaction systems with a small number of reactions and for parameters such that the number of molecules is not too large.

For chemical systems solely composed of unimolecular reactions, the problems of the SSA are immaterial since the equations for the moments derived from the CME can be solved exactly [3]. However this is not the case for chemical systems with at least one bimolecular reaction which are the norm rather than the exception in nature. In this case the equations for the moments of the CME constitute an infinite hierarchy of coupled equations and hence they cannot generally be solved exactly. The roughest approximation to the problem involves solving the deterministic rate equations for the chemical system which leads to accurate estimates for the mean concentrations in the limit of large molecule numbers [4]. Such a formalism is however inadequate when one is interested in estimating the size of the fluctuations, i.e., the variance in the fluctuations about the mean concentrations or when the goal consists in obtaining an approximate closed form solution for the probability distribution of the CME. These estimates are particularly important when the dynamics are strongly affected by the inherent stochasticity in the timing of chemical reaction events (intrinsic noise) such as the case when one or more chemical species are present in low abundances.

Over the past few decades, several methods have been developed which provide a more accurate and complete picture than that obtained from deterministic rate equations. These methods involve solving a set of deterministic ordinary differential equations whose solution provides an approximation to the moments of the probability distribution of the CME. Two popular methods of this type are the linear-noise approximation [5, 6] and moment-closure approximations (MAs) [7–11]. The linear-noise approximation corresponds to the leading order term of the system-size expansion of the CME and hence the accuracy of its predictions and the range of its applicability has been deduced by considering the next to leading order terms of the expansion [12–14]. In contrast MAs follow from an ad-hoc truncation of theinfinite hierarchy of coupled moment equations of the CME and hence little is known about their range of validity and the accuracy of the moment estimates that they provide. Error estimates for the moments of the distribution of the CME of monostable chemical systems in the limit of large molecule numbers have been derived by Grima [7]. A more fundamental question which has not been studied to-date is: when can we trust MAs to lead to physically meaningful estimates? i.e., positive real mean concentrations and positive real even central moments of the fluctuations in molecule numbers (note that by the latter we mean all moments of the type $\langle \prod_i (x_i - \langle x_i \rangle)^{k_i} \rangle$ such that $k_i$ is even for all $i$ ; this convention will be used throughout the paper).

In this paper we report on a numerical study of a class of MAs which seeks to answer the aforementioned question. The paper is organised as follows. In Section II we provide an introduction to the mathematical framework of MAs, specify the class of MA methods which we will be concerned with and define a set of criteria which guarantee physical admissibility of the solution of the MA equations. In Section III, we numerically investigate the properties of the MA equations for four chemical reaction systems which are representative of deterministically monostable, bistable and oscillatory systems; we show that the MA equations lead to physically meaningful solutions only above a certain critical molecule number for deterministic monostable systems and only in a finite range of molecule numbers for deterministic bistable and oscillatory systems. We conclude in Section IV.

II. MOMENT CLOSURE APPROXIMATIONS: DEFINITIONS AND CRITERIA FOR PHYSICAL VALIDITY

A. Background

Consider a chemical system involving species $X_i$ ( $i = 1, \dots, N$ ) interacting via $R$ chemical reactions:

$$\sum _{i=1}^N{s}_{ij}{X}_i\stackrel{c_j}{\to }\sum _{i=1}^N{r}_{ij}{X}_i,j=1,\dots ,R\mathrm{.}(1)\backslash sum\_\{i=1\}^N\; s\_\{ij\}\; X\_i\; \backslash xrightarrow\{c\_j\}\; \backslash sum\_\{i=1\}^N\; r\_\{ij\}\; X\_i,\; \backslash quad\; j\; =\; 1,\; \backslash dots,\; R.\; \backslash quad\; (1)$$

Here, $c_j$ is the rate constant of reaction $j$ . The stoichiometric matrix is defined as $S_{ij} = r_{ij} - s_{ij}$ . Under well-mixed and dilute conditions the system can be described by the joint probability distribution at time $t$ , $P(\mathbf{n}, t)$ , where $\mathbf{n} = (n_1, \dots, n_N)$ is the state vector of the system and$n_i$ is the number of molecules of species $X_i$ . Its time evolution is governed by the CME [1]:

$${\mathrm{\partial }}_{t}P(\mathbf{n},t)=\sum _{r=1}^R{f}_r(\mathbf{n}-{\mathbf{S}}_r)P(\mathbf{n}-{\mathbf{S}}_r,t)-\sum _{r=1}^R{f}_r(\mathbf{n})P(\mathbf{n},t),(2)\backslash partial\_t\; P(\backslash mathbf\{n\},\; t)\; =\; \backslash sum\_\{r=1\}^R\; f\_r(\backslash mathbf\{n\}\; -\; \backslash mathbf\{S\}\_r)\; P(\backslash mathbf\{n\}\; -\; \backslash mathbf\{S\}\_r,\; t)\; -\; \backslash sum\_\{r=1\}^R\; f\_r(\backslash mathbf\{n\})\; P(\backslash mathbf\{n\},\; t),\; \backslash quad\; (2)$$

where $\mathbf{S}_r$ is the $r$ th column vector of $S$ and $f_r(\mathbf{n})$ is the propensity function of reaction $r$ defined as [15]:

$$f_r(\mathbf{n})=c_r\mathrm{\Omega }\prod _{k=1}^N\frac{n_k!}{(n_k-s_{kj})!{\mathrm{\Omega }}^{s_{kj}}},(3)f\_r(\backslash mathbf\{n\})\; =\; c\_r\; \backslash Omega\; \backslash prod\_\{k=1\}^N\; \backslash frac\{n\_k!\}\{(n\_k\; -\; s\_\{kj\})!\; \backslash Omega^\{s\_\{kj\}\}\},\; \backslash quad\; (3)$$

and $\Omega$ is the volume of the system. To obtain the time evolution equation for the moment $\langle n_i \dots n_l \rangle$ we multiply Eq. (2) by $n_i \dots n_l$ and sum over all molecule numbers:

$${\mathrm{\partial }}_t⟨n_i\dots n_l⟩=\sum _{n_1,\dots ,n_N=0}^{\mathrm{\infty }}{n}_i\dots n_l[\sum _{r=1}^R{f}_r(\mathbf{n}-{\mathbf{S}}_r)P(\mathbf{n}-{\mathbf{S}}_r,t)-\sum _{r=1}^R{f}_r(\mathbf{n})P(\mathbf{n},t)],(4)\backslash partial\_t\; \backslash langle\; n\_i\; \backslash dots\; n\_l\; \backslash rangle\; =\; \backslash sum\_\{n\_1,\; \backslash dots,\; n\_N=0\}^\{\backslash infty\}\; n\_i\; \backslash dots\; n\_l\; \backslash left[\; \backslash sum\_\{r=1\}^R\; f\_r(\backslash mathbf\{n\}\; -\; \backslash mathbf\{S\}\_r)\; P(\backslash mathbf\{n\}\; -\; \backslash mathbf\{S\}\_r,\; t)\; -\; \backslash sum\_\{r=1\}^R\; f\_r(\backslash mathbf\{n\})\; P(\backslash mathbf\{n\},\; t)\; \backslash right],\; \backslash quad\; (4)$$

$$=\sum _{r=1}^R⟨(n_i+S_{ir})\dots (n_l+S_{lr})f_r(\mathbf{n})⟩-\sum _{r=1}^R⟨n_i\dots n_l{f}_r(\mathbf{n})⟩\mathrm{.}(5)=\; \backslash sum\_\{r=1\}^R\; \backslash langle\; (n\_i\; +\; S\_\{ir\})\; \backslash dots\; (n\_l\; +\; S\_\{lr\})\; f\_r(\backslash mathbf\{n\})\; \backslash rangle\; -\; \backslash sum\_\{r=1\}^R\; \backslash langle\; n\_i\; \backslash dots\; n\_l\; f\_r(\backslash mathbf\{n\})\; \backslash rangle.\; \backslash quad\; (5)$$

For the first two moments one obtains:

$${\mathrm{\partial }}_t⟨n_i⟩=\sum _{r=1}^R{S}_{ir}⟨f_r(\mathbf{n})⟩,(6)\backslash partial\_t\; \backslash langle\; n\_i\; \backslash rangle\; =\; \backslash sum\_\{r=1\}^R\; S\_\{ir\}\; \backslash langle\; f\_r(\backslash mathbf\{n\})\; \backslash rangle,\; \backslash quad\; (6)$$

$${\mathrm{\partial }}_t⟨n_i{n}_j⟩=\sum _{r=1}^R[S_{jr}⟨n_i{f}_r(\mathbf{n})⟩+S_{ir}⟨f_r(\mathbf{n})n_j⟩+S_{ir}{S}_{jr}⟨f_r(\mathbf{n})⟩]\mathrm{.}(7)\backslash partial\_t\; \backslash langle\; n\_i\; n\_j\; \backslash rangle\; =\; \backslash sum\_\{r=1\}^R\; \backslash left[\; S\_\{jr\}\; \backslash langle\; n\_i\; f\_r(\backslash mathbf\{n\})\; \backslash rangle\; +\; S\_\{ir\}\; \backslash langle\; f\_r(\backslash mathbf\{n\})\; n\_j\; \backslash rangle\; +\; S\_\{ir\}\; S\_\{jr\}\; \backslash langle\; f\_r(\backslash mathbf\{n\})\; \backslash rangle\; \backslash right].\; \backslash quad\; (7)$$

For chemical systems composed of only unimolecular reactions, the moment equations are closed and can be solved explicitly. However for systems with at least one bimolecular reaction, this is not the case: the equation for a certain moment will depend on higher-order moments thus leading to an infinite hierarchy of equations which cannot be solved. The idea behind moment closure approximation involves the artificial truncation of this hierarchy at a certain order to obtain a finite set of equations that can be solved numerically. The truncation involves replacing all moments above a certain order by a function of the lower order moments. The latter function is ad-hoc and hence the accuracy of such approximations is not clear. A popular means of imposing the truncation involves setting the $(N + 1)$ th and higher-order cumulants to zero which leads to a closed set of equations for the first $N$ moments (see for example [8–11]). We will refer to this method as the $N$ -moment approximation ( $N$ -MA) and exclusively focus on this class of moment closure approximations for the rest of this article. The conventional deterministic rate equations correspond to the 1MA, i.e, setting the variance to zero, and ignoring any factors of the form $\Omega^{-\alpha}$ .## B. Criteria for physical admissibility of the MA equations

Here we formulate a set of criteria which guarantee physically meaningful predictions of MA approximations and which we will repeatedly use through the rest of this article. Given deterministic initial conditions, i.e., variance and all higher-order central moments are initially equal to zero, and provided the CME has a stationary solution for the probability distribution function, the MA equations should converge to a single steady-state in the limit of long time, and the trajectories should preserve a positive mean and even central moments in the molecule numbers for all times and for all initial conditions.

Note that a unique steady-state solution of the MA equations is generally expected for all systems independent of whether they are deterministically monostable, bistable or oscillatory. What we here mean by “unique” (and throughout the rest of this article) is that given a fixed set of parameters, the time-dependent solution of the MA equations should converge in the limit of long times to the same fixed point for all possible initial conditions. Clearly this has to be the case since the stationary probability distribution of the CME is independent of initial conditions and hence the same steady-state moments must be reachable from all initial conditions.

Note also that the long time solution of the MA equations should not show sustained oscillations for systems with time-independent rate constants. This is since the moments of the CME in the limit of long time are always non-oscillatory (though the approach to the steady-state can be oscillatory) independent of whether the deterministic rate equations exhibit sustained oscillations or not. The explanation behind this phenomenon is that even if single trajectories of the SSA display sustained oscillations, independent trajectories get out of phase as time progresses and hence the ensemble-average over all the trajectories can only lead to non-oscillatory moments in the limit of long time.

III. NUMERICAL ANALYSIS OF THE MA EQUATIONS

In this section, we show that the criteria set forth in Section II are not met by the MA equations for a number of chemical systems including some which are of biological relevance. We consider three types of systems: those whose rate equations have a single steady-state (deterministic monostable systems), those whose rate equations have two steady-states (de-terministic bistable systems) and those whose rate equations predict sustained oscillations (deterministic oscillatory systems). The chosen systems were selected since they are simple enough to study in depth while at the same time their behavior is representative of a large class of systems encountered in chemistry and biochemistry. Stochastic simulations for all systems except that in Section III A.1 were done using the software package iNA [13].

A. Deterministic monostable systems

1. Bursty gene expression with no feedback

We consider a model of bursty gene expression followed by a post-translational protein dimerisation reaction:

$$\mathrm{\varnothing }\stackrel{c_1{q}_m}{\to }mP,m\in \mathbb{N}\setminus \{0\},P+P\stackrel{c_2}{\to }\mathrm{\varnothing },(8)\backslash emptyset\; \backslash xrightarrow\{c\_1\; q\_m\}\; mP,\; \backslash quad\; m\; \backslash in\; \backslash mathbb\{N\}\; \backslash setminus\; \backslash \{0\backslash \},\; \backslash quad\; P\; +\; P\; \backslash xrightarrow\{c\_2\}\; \backslash emptyset,\; \backslash quad\; (8)$$

where $P$ is the protein species, $c_2$ is the dimerisation rate constant, and $m$ is the burst size. Experimental [16] and theoretical [17, 18] evidence indicates a geometric distribution $q_m = p(1-p)^{m-1}$ with constant parameter $p$ ( $p < 1$ ) as an appropriate model for bursting; $c_1 q_m$ are then the rates at which bursts of size $m$ are created. The production step can be viewed as either an infinite number of input reactions, or equivalently as a single input reaction with input size $m$ being a random variable. Note that $1/p$ is equal to the mean burst size. Note also that in the limit of $p \approx 1$ , the expression is non-bursty and the set of protein production reactions in (8) reduces to the single reaction $\emptyset \xrightarrow{c_1} P$ .

We rescale time as $\tau = tc_2/\Omega$ and define the dimensionless constant $k = \Omega^2 c_1/c_2$ , where $\Omega$ is the volume of the system. It is easy to show that the rate equations for this system have a unique positive fixed point which is globally attractive for all $k$ ; exact stochastic simulations using the stochastic simulation algorithm (SSA [1]) also show that the CME has a stationary solution for all values of $k$ .

Starting from the CME for reaction scheme (8), we can derive the time evolution equations for the first moment $\langle n \rangle$ and the second moment $\langle n^2 \rangle$ as described in Section II:

$${\mathrm{\partial }}_{\tau }⟨n⟩=-2⟨n^2⟩+2⟨n⟩+k⟨m⟩,(9)\backslash partial\_\backslash tau\; \backslash langle\; n\; \backslash rangle\; =\; -2\backslash langle\; n^2\; \backslash rangle\; +\; 2\backslash langle\; n\; \backslash rangle\; +\; k\backslash langle\; m\; \backslash rangle,\; \backslash quad\; (9)$$

$${\mathrm{\partial }}_{\tau }⟨n^2⟩=-4⟨n^3⟩+8⟨n^2⟩-4⟨n⟩+2⟨n⟩k⟨m⟩+k⟨m^2⟩,(10)\backslash partial\_\backslash tau\; \backslash langle\; n^2\; \backslash rangle\; =\; -4\backslash langle\; n^3\; \backslash rangle\; +\; 8\backslash langle\; n^2\; \backslash rangle\; -\; 4\backslash langle\; n\; \backslash rangle\; +\; 2\backslash langle\; n\; \backslash rangle\; k\backslash langle\; m\; \backslash rangle\; +\; k\backslash langle\; m^2\; \backslash rangle,\; \backslash quad\; (10)$$Figure 1: Stability of the 2MA equations describing bursty gene expression (reaction scheme (8)). The three plots in the upper panel show the real and imaginary parts of the mean protein number $\mu$ and the variance $\Sigma$ in protein number fluctuations of the three fixed points $S_i$ for the 2MA Eqs. (11)-(12) as a function of $k$ for a mean burst size of 20 (i.e., $p = 1/20$ ). The three plots in the lower panel show the associated real and imaginary parts of the eigenvalues of the Jacobian, $\lambda_1$ and $\lambda_2$ , respectively. We find that $S_2$ and $S_3$ are always unstable, while $S_1$ is stable and the mean and variance are positive, if and only if $k > k_1$ .

where $\langle m \rangle = 1/p$ and $\langle m^2 \rangle = (2 - p)/p^2$ (these follow from the definition of $q_m$ ). These equations are then closed by using the 2MA (setting the third cumulant of $n$ to zero), leading to:

$${\mathrm{\partial }}_{\tau }\mu =⟨m⟩k+2\mu -2{\mu }^2-2\mathrm{\Sigma },(11)\backslash partial\_\backslash tau\; \backslash mu\; =\; \backslash langle\; m\; \backslash rangle\; k\; +\; 2\backslash mu\; -\; 2\backslash mu^2\; -\; 2\backslash Sigma,\; \backslash quad\; (11)$$

$${\mathrm{\partial }}_{\tau }\mathrm{\Sigma }=⟨m^2⟩k+4\mu (\mu -1)+8(1-\mu )\mathrm{\Sigma },(12)\backslash partial\_\backslash tau\; \backslash Sigma\; =\; \backslash langle\; m^2\; \backslash rangle\; k\; +\; 4\backslash mu(\backslash mu\; -\; 1)\; +\; 8(1\; -\; \backslash mu)\backslash Sigma,\; \backslash quad\; (12)$$

where $\mu = \langle n \rangle$ and $\Sigma = \langle n^2 \rangle - \langle n \rangle^2$ are the mean and variance in protein numbers.

Setting the left hand side of Eqs. (11)-(12) to zero, and solving simultaneously, one finds that there are three possible solutions which we call $S_i = (\mu_i, \Sigma_i)$ , $i = 1, \dots, 3$ . Figure 1 shows the real and imaginary parts of the mean and variance of these three fixed points, as well as the real and imaginary parts of the associated eigenvalues of the Jacobian of Eqs. (11)-(12) for the case $p = 1/20$ (this corresponds to a mean burst size of 20 which has been measured experimentally for gene expression [16]). By inspection of Figure 1, we see that of the three possible steady-state solutions only $S_1$ is physically admissible; this is since it is the onlysteady-state solution which displays a positive mean and variance of molecule numbers and which is locally stable (negative real part of the eigenvalues of the Jacobian). However note that these properties only manifest for $k$ larger than a certain critical value $k = k_1 \approx 70$ . This would lead one to surmise that the CME has a stationary solution only for $k$ greater than this critical value. However, as noted earlier this is not the case: the CME has a stationary solution for all values of $k$ . These results taken together imply that the 2MA does not give a physically meaningful steady-state solution for all values of $k$ .

Next we study the time-evolution leading to the steady-state. Figure 2 shows the numerically integrated time trajectories for several initial conditions for three different $k$ values and for $p = 1/20$ . We find that for $k < k_2 \approx 1.2k_1 \approx 85$ , some of the trajectories diverge as time goes to infinity which is unphysical since the CME has a stable fixed point. This instability manifests for all deterministic initial conditions for $k < k_1$ and for initial conditions characterised by a small initial mean number of protein molecules for $k_1 \leq k < k_2$ . For $k \geq k_2$ , however, the trajectories converge to the stable fixed point and are non-negative at all times. We verified this numerically for initial conditions up to $\mu(\tau = 0) = 10^6$ . Hence in coincidence with the steady-state analysis above, we find that the 2MA only gives physically admissible solutions for a certain range of $k$ . Note that the requirement that the time trajectories of the moments are physically admissible at all times is harder to satisfy than the requirement that there exists a single physically admissible steady-state; this is since the critical value of $k$ above which the former is satisfied ( $k_2$ ) is larger than the critical value of $k$ above which the steady-state criteria are satisfied ( $k_1$ ). For the rest of this article we shall refer to the first requirement noted above as the time-dependent criterion and the second requirement as the steady-state criterion.

The analysis described above was specifically for the case of $p = 1/20$ . Qualitatively similar results for the 2MA equations are found for all values of $p$ , i.e, there exist a critical $p$ -dependent values $k_1$ and $k_2$ such that for $k > k_1$ the steady-state criterion is satisfied and for $k > k_2$ the time-dependent criterion is satisfied; $k_2 \geq k_1$ such that the time-dependent criterion is generally the more difficult of the two criteria to satisfy. Figure 3 shows the $p$ -dependence of the critical value $k_2$ as well as the $p$ -dependence of the corresponding mean particle number $\langle n \rangle_{\text{CME}}$ . Note that both $k_2$ and $\langle n \rangle_{\text{CME}}$ increase with $1/p$ implying that the larger the burstiness in protein expression, the larger is the critical molecule number above which the 2MA gives physically meaningful results. In particular for the case $p = 1/20$ , weFigure 2: Time trajectories of the 2MA equations and of the ensemble-averaged SSA in the $\mu$ - $\Sigma$ plane for different values of $k$ characterising bursty gene expression (reaction scheme (8)). The initial conditions are deterministic. The blue and magenta dots show the stable and unstable fixed points (FPs) of the 2MA equations, respectively. While for $k = 65$ there is no fixed point in the positive orthant and the time trajectories diverge for all initial conditions, for $k = 75$ most initial conditions lead to trajectories converging to a unique stable fixed point. However, a small initial mean number of molecules still leads to divergence. Finally, for $k = 85$ all initial conditions lead to trajectories converging to a unique fixed point with positive mean and variance. In contrast as SSA simulations show, the CME has a stationary solution for all values of $k$ . The mean burst size is 20 ( $p = 1/20$ ) as in Figure 1.

had earlier found that $k_2 = 85$ which corresponds to a mean steady-state protein number of $\langle n \rangle_{CME} \approx 25$ , i.e., the 2MA equations for a mean protein burst size of 20 give physically meaningful results for the time-evolution of the system only when the number of protein molecules in steady-state exceeds 25. It is well known that protein numbers per cell can be very small, even of the order of a few molecules and hence our results show that one must be careful in the use of the moment-closure approximation to understand cell level phenomena.

Similar results to the ones found for the 2MA are found for the higher-order moment closure approximations. In Figure 4 we show the 3MA analog of the 2MA time-evolution analysis shown in Figure 2. The qualitative similarity between the two figures is evident: the time-evolution criteria are only satisfied for $k$ greater than a certain critical value ( $\approx 100$ ) since for smaller values of $k$ , we have a non-positive variance of molecule numbers as the steady-state is approached. We have also verified the same qualitative behaviour for theFigure 3: Critical values ( $k_2$ , $\langle n \rangle_{\text{CME}}$ ) of the 2MA as a function of mean burst size $1/p$ for reaction scheme (8). Both curves are monotonically increasing with increasing $1/p$ implying that the critical molecule number above which the 2MA equations lead to physically meaningful results (both criteria in Section II B are satisfied) increases with burstiness in protein expression.

4MA and 5MA (in addition to the mean and variance, we require the fourth central moment to be positive for these MA equations) which suggests that this behaviour exists for any order of the MA method.

Adopting the terminology that $k_1$ and $k_2$ are the values of $k$ above which the MA equations of any order satisfy the steady-state and time-dependent criteria, we show in Figure 5 the dependence of these two values of $k$ on the mean burst size $1/p$ for all four MAs. Two observations can be made: (i) for all orders of the MA, $k_2 > k_1$ implying that the time-dependent criteria are the more difficult of the two to satisfy; this also implies that the mean molecule number associated with $k_2$ is the critical molecule number above which the MA equations give physically meaningful results. (ii) for a given $1/p$ , the value of $k_2$ is monotonically increasing with increasing order of the MA, i.e., the range of molecule numbers over which one observes physically meaningful results decreases with increasing order of the MA. For $1/p = 20$ , for example, we obtain the values $k_2 = 85, 116, 420, 643$ for the 2-5 MA equations respectively. Simulating the system using the SSA we find that these values correspond to mean protein numbers of 26, 31, 61 and 75, respectively. Finally in Figure 6 we show a plot of the relative error in the mean number of molecules predicted by the MA as a function of $k$ and of the order of the MA; as expected accuracy increases with the order but the range of $k$ over which the approximation is valid decreases with increasing order.Figure 4: Time trajectories of the 3MA equations and the ensemble-averaged SSA in the $\mu - \Sigma$ plane for different values of $k$ characterising bursty gene expression (reaction scheme (8)). For $k = 20$ there is an unphysical stable fixed point of the MA equations with negative variance. For $k = 60$ the stable fixed point has positive mean and variance but the variance becomes negative as the fixed point is approached. Finally, for $k = 150$ all trajectories converge to a physically admissible steady-state and are physically meaningful (positive mean and variance at all times). In contrast the CME has a stationary solution for all values of $k$ . The mean burst size is 20 ( $p = 1/20$ ).

Summarising our analysis in this section shows that moment closure approximations for a gene circuit involving a bimolecular reaction give physically meaningful results (satisfy both criteria set forth in Section II B) only when the protein molecule numbers are above a certain critical threshold. This threshold increases with the mean burst size and closure order.

2. Gene expression with negative feedback

We next consider the following gene regulatory network:

$$D_u\stackrel{c_1}{\to }{D}_u+X,D_u+X\underset{c_3}{\overset{c_2}{\rightleftharpoons }}{D}_b,X\stackrel{c_4}{\to }\mathrm{\varnothing }\mathrm{.}(13)D\_u\; \backslash xrightarrow\{c\_1\}\; D\_u\; +\; X,\; \backslash quad\; D\_u\; +\; X\; \backslash xrightleftharpoons[c\_3]\{c\_2\}\; D\_b,\; \backslash quad\; X\; \backslash xrightarrow\{c\_4\}\; \backslash emptyset.\; \backslash quad\; (13)$$

A single gene in the unbound state $D_u$ expresses a protein $X$ which subsequently binds to the same gene and forms a non-expressive complex $D_b$ . This is a negative feedback loop since the protein suppresses its own expression. Note that the mRNA is here not modelled explicitly for simplicity. It has been analytically shown that the rate equations possessFigure 5: Critical values $k_1$ and $k_2$ as a function of the mean burst size $1/p$ for the 2-5 MA equations describing bursty gene expression (reaction scheme (8)). In panel (a) we show that $k_2 > k_1$ for all orders of the MA and for all $p$ ; hence the time-evolution criterion is satisfied over a smaller range of $k$ than the steady-state criterion. In panel (b) we show that $k_2$ increases monotonically with burstiness ( $1/p$ ) and with the order of the MA. Hence the higher the order of the MA, and the larger the noise in gene expression, the smaller is the range of $k$ over which the MA equations give physically meaningful results for the time-evolution of the system.

a single steady-state solution and the CME has a stationary solution for the probability distribution function for all parameter values [19]. Hence as for the previous example, we can check if the steady-state and time-dependent criteria are satisfied by the MA equations for this circuit. Note that this example unlike the previous one is two dimensional.

We fix the rate constants to $c_1/c_4 = 10$ , $c_2/c_4 = 1$ , $c_3/c_4 = 0.5$ and study the nature of the solutions of the MA equations as a function of the cell volume $\Omega$ . Solving the 2-5MA equations in steady-state we find that they possess a unique stable fixed point with positive means and even central moments only above a certain critical volume $\Omega = \Omega_1$ ; in particular we find the values $\Omega_1 = 5.8, 16, 28, 39$ , for the 2-5MA, respectively. Solving the 2-5MA equations as a function of time we find that for volumes less than the critical volumes stated above, the trajectories of the moments either diverge or converge to an unphysical fixed point characterised by a negative mean. In contrast for volumes larger than the critical volumes $\Omega_1$ all initial conditions tested give rise to convergent and physically meaningful time trajectories. These results imply that for this set of constants, the critical volumes above which the time-dependent criterion is satisfied are the same as the critical volumesFigure 6: Relative error in the mean number of molecules of the MA ( $\Delta\mu$ ) as a function of the non-dimensional parameter $k$ and of the order of the MA. The error is computed as the absolute difference between the steady-state prediction of the first moment of the MA equations and the ensemble averaged SSA result for the mean number of molecules divided by the latter. The parameter $p$ is fixed to $1/20$ . The accuracy of the MA increases with the order of the approximation while the range of validity of the MA equations decreases with the order (since the minimum value of $k$ at which the MA gives a physically meaningful result increases with order).

above which the steady-state criterion is satisfied.

From the exact solution of the CME for this circuit [19], we find that the mean protein numbers corresponding to the critical volumes of the 2-5MA are 3.3, 5.2, 6.3 and 6.9, respectively. Hence as for the previous example of bursty gene expression with no feedback, (i) the MA equations possess physically admissible solutions only when the steady-state molecule numbers predicted by the CME are above a certain minimum; (ii) these critical molecule numbers increase with the order of the MA equations.

B. A deterministic bistable system

Next, we consider a bistable reaction system that has been extensively studied in the literature, the Schlögl system [20]

$$2X\underset{c_2}{\overset{c_1}{\rightleftharpoons }}3X,\mathrm{\varnothing }\underset{c_4}{\overset{c_3}{\rightleftharpoons }}X\mathrm{.}(14)2X\; \backslash xrightleftharpoons[c\_2]\{c\_1\}\; 3X,\; \backslash quad\; \backslash emptyset\; \backslash xrightleftharpoons[c\_4]\{c\_3\}\; X.\; \backslash quad\; (14)$$Figure 7: Monostable and bistable regimes of the MA solution as a function of the volume $\Omega$ for the Schlögl reaction scheme in (14) with parameter set S1. The crosses denote the critical volumes above which all time trajectories for the mean and the even central moments are positive and convergent to a stable steady-state.

The deterministic rate equations of the system (14) have two steady-state solutions (bistability) in certain parameter regimes and one steady-state solution (monostability) in the rest of parameter space [21]. The CME can also be solved exactly in steady state since the system is in detailed balance [21, 22]. In what follows, we study the nature of the solutions of the MA equations for a parameter set in which the system is deterministically bistable and a parameter set for which it is monostable.

First, we consider the set of rate constants $c_1/c_4 = 2.6, c_2/c_4 = 1.7, c_3/c_4 = 0.1$ which we shall refer to as “S1”. The deterministic rate equations have two positive stable fixed points for approximate concentrations, 0.16 and 1.0, respectively. In what follows we study the nature of the solutions of the MA equations for S1 as a function of the volume $\Omega$ . Solving the time-dependent MA equations, we find that the requirement of convergent and physically meaningful time trajectories is only fulfilled above a critical volume. We find that the value of this critical volume increases with increasing closure order, ranging from 0.32 for the 2MA to 5.6 for the 5MA which corresponds to a range of mean molecule numbers (as calculated from the CME) equal to $\langle n \rangle \approx 0.02 - 1$ . Hence as for the monostable circuits studied earlier, it is clear that the MA equations give valid meaningful results only above a certain critical molecule number.

However there is a major difference between the solution of the MA equations for this system and the MA equations for the previous two monostable systems: the time-dependentFigure 8: Time trajectories of the SSA and 2-4MA equations for the Schlögl reaction scheme in (14) with parameter set S1 and volume $\Omega = 235$ . The inset for the 2MA shows how the time trajectories for small initial mean molecule numbers converge to one of two stable fixed points; time trajectories for large initial mean molecule numbers converge to the other stable fixed point. This bistable behaviour evident in the solution of the 2MA is “washed away” as the order of the MA is increased and the time trajectories become similar to those of the SSA. This shows that the bistability in the 2MA solution at this volume should be considered as an artifact of the moment-closure approximation method.

MA equations converge to either one or two physically admissible steady-state values of the moments, depending on the volume. This bifurcation behaviour together with the critical molecule numbers described above are schematically summarised in Figure 7; the time-evolution of the 2-4MA equations are compared to the SSA in Figure 8. The non-uniqueness of a steady-state with positive mean and even central moments breaks the physical admissibility criteria set forth in Section II B since the moments of any probability distribution (independent of whether the probability distribution of the CME is unimodal or multimodal) should be single-valued. Hence on this basis one would conclude that the MA equations givephysically admissible results only for volumes above the critical one and for which all trajectories converge to a single-valued solution. An inspection of Figure 7 shows that this criterion implies that the 2-5MA equations lead to physically admissible solutions in the volume ranges 3.3 – 147, 2.1 – 346, 2.8 – 530 and 5.6 – 707 respectively. Thus for this deterministically bistable system, there are two critical volumes and not one as for monostable systems. It is also the case that the size of the volume range increases with the order of the MA, this being mostly due to the fact that ceiling of this range (the volume at which the MA equations switch from one to two physically admissible steady-state solutions) increases with the order of the MA. The above results have also been verified for 8 other parameter sets for which the deterministic rate equations of the Schlögl model are bistable (see Table I). Note that in a single case (first parameter set in Table I), there is no regime where the steady-state solution of the 2MA equations is physically meaningful since the solution is bistable for all volumes.

The fact that the upper critical volume increases rapidly with the order of the MA is indeed an indirect verification that the regime in which there are two physically admissible steady-state MA solutions is an artifact of the approximation. The origin of the unphysical bifurcation is currently unclear although we have confirmed (see Appendix A) that it is not related to a sudden breakdown of the cumulant neglect assumption at the heart of the MA equations.

We note that another interpretation of the phenomenon may at first appear possible. When the deterministic rate equations possess two steady-state solutions the probability distribution of the CME for large volumes is expected to be bimodal; similarly when the rate equations possess one steady-state then the distribution is unimodal for large volumes. Now if one sees the MA equations as a refinement of the conventional rate equations, in the sense that they are valid not only for large volumes but over a wider range of volumes, then one may postulate that the number of physically admissible steady-state solutions of the MA equations reflects the number of maxima of the probability distribution of the CME. As we now show this interpretation is not valid for the Schlögl model. For the parameter set $S1$ , the probability distribution of the CME is unimodal for volumes below $\Omega = 30$ and bimodal for larger volumes. In contrast the 2MA equations show two transitions from bimodal at low volumes to unimodal at intermediate volumes to bimodal at large volumes (see Figure (7)); furthermore the volumes at which the transition from unimodal to bimodal

parameters			2MA		3MA		4MA		5MA
$c\_1/c\_4$	$c\_2/c\_4$	$c\_3/c\_4$	$\backslash Omega\_1$	$\backslash Omega\_2$	$\backslash Omega\_1$	$\backslash Omega\_2$	$\backslash Omega\_1$	$\backslash Omega\_2$	$\backslash Omega\_1$	$\backslash Omega\_2$
1	$\backslash frac\{1\}\{8\}$	$\backslash frac\{1\}\{8\}$	0.0094	0	1.1	8.7	1.5	17.	1.4	29000
1	$\backslash frac\{1\}\{8\}$	$\backslash frac\{1\}\{4\}$	1.4	180	0.83	410	1.1	610	1.1	19000
1	$\backslash frac\{1\}\{4\}$	$\backslash frac\{1\}\{8\}$	39.	76.	1.1	140	1.5	190	1.6	$1.4\; \backslash times\; 10^6$
1	$\backslash frac\{1\}\{4\}$	$\backslash frac\{1\}\{4\}$	1.4	37.	0.8	91.	1.1	140	1.6	190
2	$\backslash frac\{1\}\{8\}$	$\backslash frac\{1\}\{8\}$	2.6	4300	1.7	8600	2.1	72000	1.9	44000
2	$\backslash frac\{1\}\{4\}$	$\backslash frac\{1\}\{8\}$	2.6	1300	1.7	2800	2.1	4100	2.5	23000
2	$\backslash frac\{1\}\{2\}$	$\backslash frac\{1\}\{8\}$	2.6	370	1.7	810	2.1	1200	2.2	38000
2	1	$\backslash frac\{1\}\{8\}$	2.8	74.	1.6	180	2.1	280	3.3	370

Table I: The table shows the two volumes ( $\Omega_1$ and $\Omega_2$ ) between which the MA equations have a physically meaningful solution for the Schlögl reaction scheme in Eq. (14). The data is generated for 8 distinct parameter sets. In the volume range $\Omega_1$ – $\Omega_2$ the time trajectories of the MA are convergent to a unique steady-state and exhibit physically meaningful moments at all times (positive mean and even central moments). Above $\Omega_2$ , the moments are not unique and below $\Omega_1$ the time trajectories do not converge or else there exist negative mean and even central moments at a point in time. We consider all combinations for the parameters $c_1/c_4$ , $c_2/c_4$ and $c_3/c_4$ drawn from the set ${2^{-3}, 2^{-2}, \dots, 2^3}$ for which the rate equations are bistable.

occurs ( $\Omega = 147, 346, 530, 707$ for the 2-5MA respectively) bear no relationship to the actual transition volume $\Omega = 30$ obtained from the CME. Hence we arrive to the conclusion as before, namely that generally the bifurcations in the solutions of the MA equations are artificial in the sense that they are not indicative of any real transition in the number of maxima of the probability distribution of the CME. Preliminary analysis (see Appendix B) does however suggest that the two steady-state solutions of the MA equations contain information (position and width) on the two peaks of the bimodal distribution of the CME but not on the height of the two peaks; this information is partial in the sense that it cannot be used to reconstruct the probability distribution or to calculate any of its moments and furthermore the information about the individual peaks is only valid over a subset of the volumes over which the probability distribution of the CME is bimodal. Hence in line withFigure 9: Monostable and bistable regimes of the MA solution as a function of the volume $\Omega$ for the Schlögl reaction scheme in (14) with parameter set S2. The crosses denote the critical volumes above which all time trajectories for the mean and the even central moments are positive and convergent to a stable steady-state.

our previous analysis, it can be concluded that it is only safe to trust information from the MA equations when the time-evolution of the trajectories is physically meaningful (positive mean and even central moments) and when they converge to a single steady-state.

Next, we consider the set of rate constants $c_1/c_4 = 10, c_2/c_4 = 0.1, c_3/c_4 = 0.1$ which we shall refer to as “S2”. The deterministic rate equations have now one positive stable fixed point and the probability distribution of the CME is unimodal for all volumes that we checked. As for the monostable circuits studied in Section III A, the MA equations only give physically meaningful results (time-dependent criterion is satisfied) above a certain critical volume. These together with the number of steady-state solutions of the MA equations are summarised in Fig. 9. The critical volumes range from $\Omega \approx 3.3$ for the 2MA to $\Omega \approx 5$ for the 5MA and the associated mean molecule numbers (calculated from the CME) range from $\langle n \rangle \approx 330$ to $\langle n \rangle \approx 510$ . The bifurcation from single to two steady-states in the MA equations is completely unphysical since (i) moments can only be single-valued, (ii) there is no corresponding transition in the number of maxima of the probability distribution of the CME, and (iii) the two solutions of the MA equations do not provide any meaningful local information on the individual modes of the distribution since the distribution is unimodal (unlike the case for bistable parameter set S1).### C. A deterministic oscillatory system

Consider next the Brusselator, a well known oscillating chemical system [23, 24]

$$2X+Y\stackrel{c_1}{\to }3X,X\stackrel{c_2}{\to }Y,\mathrm{\varnothing }\underset{c_4}{\overset{c_3}{\rightleftharpoons }}X\mathrm{.}(15)2X\; +\; Y\; \backslash xrightarrow\{c\_1\}\; 3X,\; \backslash quad\; X\; \backslash xrightarrow\{c\_2\}\; Y,\; \backslash quad\; \backslash emptyset\; \backslash xrightleftharpoons[c\_4]\{c\_3\}\; X.\; \backslash quad\; (15)$$

The deterministic rate equations predict sustained oscillations in certain parameter regimes, and damped or no oscillations in other regimes. In contrast, stochastic simulations show that the mean molecule numbers either exhibit damped oscillations or no oscillations [25]; sustained oscillations are only seen in individual trajectories of the SSA but due to dephasing between independent trajectories, the mean molecule numbers (calculated over an ensemble of trajectories) only show damped oscillations. We now study the MA predictions for this system.

We fix the rate constants to $c_1/c_4 = 1, c_2/c_4 = 3, c_3/c_4 = 1$ and study the MA equation properties as a function of the volume $\Omega$ . The deterministic rate equations exhibit sustained oscillations while the mean molecule numbers computed from the SSA show damped oscillations (for all tested volumes) before settling to a non-oscillatory steady-state (see top panel of four graphs in Figure 10). In the same figure, we show the solutions of the 2-5MA as a function of the volume. Note that the MA solutions show a transition from sustained oscillations to damped oscillations at a certain critical volume; the MA solution is an approximation for the moments, i.e., ensemble-averaged behaviour of the CME, however no such transition is seen in the CME solution and hence the sustained oscillation solution of the MA equations is to be treated as an artifact. A further proof of the artifactual nature of the transition is that the critical volume at which it occurs increases rapidly with the order of the MA; these critical volumes are $\Omega = 22, 89, 320, 4900$ for the 2-5MA equations (denoted collectively as $\Omega_2$ ). This picture is consistent with the critical volume going to infinity in the limit of infinite order of the MA which would imply agreement with the ensemble-average results of the CME in this limit. As for the unphysical bifurcation in the Schlögl system, the origin of the artificial transition from sustained to damped oscillations is unclear since it does not appear to be related to a sudden breakdown of the cumulant neglect assumption at the heart of the MA equations.

As for previous systems, another critical volume is found below which the trajectories of the MA equations do not converge or lead to unphysical moments. These critical volumes are found to be $\Omega = 1.4, 0.41, 0.57$ and $0.79$ for the 2-5MAs (denoted collectively as $\Omega_1$ ).Figure 10: Comparison of the time evolution of the mean molecule numbers according to the CME and to the MA equations for the Brusselator system (15) as a function of the volume $\Omega$ . The rate parameters are $c_1/c_4 = 1, c_2/c_4 = 3, c_3/c_4 = 1$ . Blue and red lines denote the solutions for species $X$ and $Y$ respectively. Note that the MA equations predict a transition from sustained to damped oscillations as the volume decreases but the CME predicts no such transition; the transition is thus artifactual also since the critical volume at which it occurs increases with the order of the MA.

These correspond to a range of mean molecule numbers of species $X$ between 0.4 and 1.4 and for species $Y$ between 8 to 9. A schematic summary of our numerical analysis, including the transition behaviour earlier discussed, is shown in Fig. 11.

The time trajectory and transition analysis put together lead us to the conclusion that the MA equations for the Brusselator lead to physically meaningful predictions only in the finite range of volumes $(\Omega_1, \Omega_2)$ (and the associated finite range of molecule numbers); this conclusion is similar to that for the bistable parameter set S1 of the Schlögl system in Section III B. The analysis here was specifically done for the rate constants $c_1/c_4 = 1, c_2/c_4 = 3, c_3/c_4 =$Figure 11: Damped and sustained oscillation regimes of the MA solution as a function of the volume $\Omega$ for the Brusselator reaction scheme (15) with rate constants $c_1/c_4 = 1, c_2/c_4 = 3, c_3/c_4 = 1$ . The crosses denote the critical volumes above which all time trajectories for the mean and even central moments are positive and convergent to a steady-state.

1; we have verified that the same picture and conclusions emerge from studying twelve other sets of rate constants (see Table II). Note that in 8 cases, there is no regime where the steady-state solution of the 2MA equations is physically meaningful since the solution is oscillatory for all volumes; in contrast there is just one case where the 3MA suffers a similar problem and no such cases are found using the 4 and 5MA. These results suggests that our analysis broadly holds for all parameters such that the deterministic rate equations predict sustained oscillations.

IV. SUMMARY AND CONCLUSION

In this paper we have elucidated, by means of several exemplary reaction systems, the range of validity of a popular class of MA equations. In particular our numerical results suggest that the solutions of these equations are only physically meaningful when the steady-state mean molecule numbers obtained from the CME are above a certain threshold for deterministic monostable systems and when these molecule numbers fall within a certain finite range of mean molecule numbers for deterministic bistable and oscillatory systems. Our results have important implications for the use of MA approaches to either predict the

parameters			2MA		3MA		4MA		5MA
$c\_1/c\_4$	$c\_2/c\_4$	$c\_3/c\_4$	$\backslash Omega\_1$	$\backslash Omega\_2$	$\backslash Omega\_1$	$\backslash Omega\_2$	$\backslash Omega\_1$	$\backslash Omega\_2$	$\backslash Omega\_1$	$\backslash Omega\_2$
$\backslash frac\{1\}\{2\}$	2	$\backslash frac\{1\}\{2\}$	1.8	1.8	0.67	5.7	0.52	17	0.52	44
$\backslash frac\{1\}\{2\}$	2	1	0.9	14	0.37	350	0.44	4900	0.67	$>\; 10^6$
$\backslash frac\{1\}\{2\}$	4	$\backslash frac\{1\}\{2\}$	3	3	1.1	1.1	0.4	4.6	0.38	8
$\backslash frac\{1\}\{2\}$	4	1	1.5	1.5	0.56	5.9	0.4	12	0.61	21
$\backslash frac\{1\}\{2\}$	4	2	0.77	28	0.23	88	0.34	210	0.38	410
1	2	$\backslash frac\{1\}\{2\}$	1.8	1.8	0.67	23	0.67	103	0.65	$>\; 10^6$
1	4	$\backslash frac\{1\}\{2\}$	1.8	1.8	1.1	3.9	0.57	10	0.9	17
1	4	1	1.5	1.5	0.56	20	0.52	40	0.79	$>\; 10^6$
2	2	$\backslash frac\{1\}\{2\}$	1.8	27	0.72	710	0.88	9700	1.3	$>\; 10^6$
2	4	$\backslash frac\{1\}\{2\}$	1.5	1.5	1.1	11	0.79	25	1.2	42
2	4	1	1.5	56	0.59	170	0.54	420	0.77	830
4	4	$\backslash frac\{1\}\{2\}$	3	3	1.1	41	1	80	1.6	130

Table II: The table shows the two volumes ( $\Omega_1$ and $\Omega_2$ ) between which the MA equations have a physically meaningful solution for the Brusselator reaction scheme in Eq. (15). The data is generated for 12 distinct parameter sets. In the volume range $\Omega_1 - \Omega_2$ the time trajectories of the MA are convergent to a unique steady-state and exhibit physically meaningful moments at all times (positive mean and even central moments). Above $\Omega_2$ , the moment solutions of the MA exhibit artificial sustained oscillations (the SSA always has damped oscillations) and below $\Omega_1$ the time trajectories do not converge or else there exist negative mean and even central moments at a point in time. We consider all combinations for the parameters $c_1/c_4$ , $c_2/c_4$ and $c_3/c_4$ drawn from the set ${1/2, 1, 2, 4}$ for which the deterministic rate equations predict sustained oscillations. The initial condition is zero mean and higher-order moments.

stochastic dynamics of chemical systems [9, 10] or for parameter inference [26, 27].

By physically meaningful solutions we specifically mean that (i) the mean molecule numbers and the even central moments of the fluctuations in molecule numbers predicted by the MA equations are positive real numbers at all times and they converge to steady-state values whenever the CME has a stationary solution, (ii) the moments are unique in thesense that the same steady-state moments can be reached from all initial conditions, and (iii) the moments do not exhibit sustained oscillations in the limit of long times. The first two properties are self-evident while the third may not be immediately obvious – this stems from the fact that though individual trajectories of the SSA can exhibit sustained oscillations, all these noisy trajectories are not in phase since they are independent (provided the rate constants are time independent) and hence an ensemble average leads to damped oscillations through destructive interference between the trajectories. We found that deterministic monostable systems suffer from a breakdown of property (i) for mean molecule numbers below a certain value, deterministic bistable systems suffer from the same and as well from a breakdown of property (ii) for mean molecule numbers above a certain value, while deterministic oscillatory systems have the same property as monostable systems and as well suffer from a breakdown of property (iii) for mean molecule numbers above a certain value.

We have shown that the transition in the number of steady-state solutions (from 1 to 2) which lead to a breakdown of property (ii) and the transition from damped to oscillatory long time behaviour which lead to the breakdown of property (iii) occur at increasing higher mean molecule numbers as the order of the MA equations is increased – this strongly suggests that these transitions are unphysical. This conclusion was further supported by showing that the transitions are uncorrelated with sudden changes in the number of maxima of the probability distribution of the CME or in the moment dynamics of the CME.

We note that above the mean molecule numbers for which properties (ii) and (iii) breakdown, the solution of the moment equations closely resembles that of the deterministic rate equations rather than the moments of the probability distribution of the CME. However as noted above, the range of mean molecule numbers for which all three properties are satisfied is found to increase rapidly as the order of the MA equations is increased – this means that the solution of the MA resembles less the deterministic rate equations and more the CME as the order is increased. This result is consistent with the fact that the first-order MA is either the same or approximately equal to the deterministic rate equations (since the 1MA involves setting the variance of fluctuations to zero) while the infinite order MA is equivalent to the exact moments of the CME.

We emphasise that the results here reported are specifically for bursty gene expression and a genetic negative feedback loop (monostable systems), the Schlögl model (a bistablesystem) and the Brusselator (an oscillatory system); it remains to be seen whether the conclusions obtained herein are common to all deterministic monostable, bistable and oscillatory systems, albeit its unlikely that the latter can be conclusively proved since the MA framework is typically only amenable to numerical analysis.

In conclusion, we have here elucidated the conditions necessary for the validity of the solution of the MA equations. Our results suggest that though this approach presents an efficient computational means of obtaining approximate solutions to the moments of the CME, there are several pitfalls in the interpretation of the results particularly for deterministic bistable and oscillatory systems; these difficulties can be only be overcome through the systematic comparison of the solutions of the MA equations of a certain order with those of higher orders.

Acknowledgments

G.S. acknowledges support from the European Research Council under grant MLCS 306999.

---

• [1] D. T. Gillespie, Annu. Rev. Phys. Chem. 58, 35 (2007)
• [2] D. T. Gillespie, J. Phys. Chem. 81, 2340 (1977)
• [3] D. A. McQuarrie, J. App. Prob. 4, 413 (1967)
• [4] T. G. Kurtz, J. Chem. Phys. 57, 2976 (1972)
• [5] N. van Kampen, Can. J. Phys. 39, 551 (1961)
• [6] J. Elf and M. Ehrenberg, Gen. Res. 13, 2475 (2003)
• [7] R. Grima, J. Chem. Phys. 136, 154105 (2012)
• [8] L. Ferm, P. Lötstedt and A. Hellander, J. Sci. Comput. 34, 127 (2008)
• [9] M. Ullah and O. Wolkenhauer, J. Theor. Biol. 260, 340 (2009)
• [10] C. A. Gomez-Urbe and G. C. Verghese, J. Chem. Phys. 126, 024109 (2007)
• [11] A. Ale, P. Kirk and M. P. H. Stumpf, J. Chem. Phys. 138, 174101 (2013)
• [12] R. Grima, J. Chem. Phys. 133, 035101 (2010)
• [13] P. Thomas, H. Matuschek and R. Grima, PLoS ONE 7(6): e38518 (2012)- [14] P. Thomas, H. Matuschek and R. Grima, BMC Genomics 14 (Suppl 4):S5 (2013)
• [15] N. G. van Kampen, Stochastic Processes in Physics and Chemistry (Elsevier, 2007)
• [16] L. Cai, N. Friedman, and X.S. Xie, Nature 440, 358 (2006)
• [17] M. Thattai and A. van Oudenaarden, Proc. Natl. Acad. Sci. USA 98, 8614 (2001)
• [18] V. Shahrezaei and P.S. Swain, Proc. Natl. Acad. Sci. USA 105, 17256 (2008)
• [19] R. Grima, D. R. Schmidt and T. J. Newman, J. Chem. Phys. 137, 035104 (2012)
• [20] F. Schlögl; Z. Physik 253, 147 (1972)
• [21] I. Matheson, D. F. Walls and C. W. Gardiner, J. Stat. Phys. 12, 21 (1975)
• [22] W. Ebeling and L. Schimansky-Geier, Physica A 98, 587 (1979)
• [23] I. Prigogine and R. Lefever, J. Chem. Phys. 48, 1695 (1968).
• [24] R. Lefever, G. Nicolis, and P. Borckmans, J. Chem. Soc., Faraday Trans. I 84, 1013 (1988).
• [25] D. L. K. Toner and R. Grima, J. Chem. Phys. 138, 055101 (2013)
• [26] P. Milner, C. S. Gillespie and D. J. Wilkinson, Statistics and Computing 23, 287 (2013)
• [27] C. Zechner et al., Proc. Natl. Acad. Sci. USA 109, 8340 (2012)

Appendix A: Origin of the unphysical bifurcation in the MA equations

One is naturally led to the question if there is a way to predict the unphysical bifurcation in the MA equations for the Schlögl model studied in Section III B. Recall that the derivation of the MA equations requires that some of the cumulants are zero. This is of course an assumption and hence one could surmise that perhaps this assumption breaks down dramatically when the MA equations experience the unphysical bifurcation, e.g., the mentioned cumulants of the probability distribution solution of the CME may suddenly take large values as the volume crosses a certain threshold. We now test this hypothesis. Since the reaction scheme (14) includes bimolecular and trimolecular reactions, to get a closed set of equations for the first $N$ moments we have to set the $(N + 1)$ th and $(N + 2)$ th cumulants to zero. We denote the $N$ th cumulant by $\kappa_N$ in the following. We compute the exact values of these two cumulants from the stationary solution of the CME and in particular plot their dependence with the volume. The results are shown in Figure 12. We find that the behaviour of the cumulants gives no apparent hint to the bifurcation of the MAs from one to two positive stable fixed points. Thus our present analysis is inconclusive as to origin ofFigure 12: Cumulants of the concentrations for the Schlogl reaction scheme in (14) with parameter set S1, as a function of the volume $\Omega$ . $\kappa_n$ denotes the cumulant of order $n$ which is computed from the exact solution of the CME. The solid vertical lines show the critical volumes ( $\Omega_2$ ) at which the MA equations change from having a single to two stable fixed points with positive mean and even central moments. Note that there is no sudden change in the cumulants as a function of volume which anticipates the transition observed in the solution of the MA equations.

the unphysical bifurcation in the MA equation solutions.

Appendix B: Information contained in the bistable solution of the 2MA equations

Since the bistable solution of the MA equations cannot be interpreted as an approximation to the moment of the probability distribution of the CME, the question arises if they can be meaningfully interpreted otherwise. In the bistable regime which appears when solving the 2MA equations in small volumes, the solution of the CME is unimodal, and there is no interpretation of the bistability. For large volumes, however, it turns out that the moments of the two 2MA fixed points provide a good approximation to the moments of the individual peaks of the bimodal probability distribution of the CME. This can be seen by constructing probability distribution functions from the MAs in the bistable regime in the following way.First, we construct two Gaussians with the respective mean and variance obtained from the two solutions of the 2MA. Secondly a linear superposition of these two Gaussians is created; the two superposition weights are calculated from the CME probability distribution by numerically calculating the area under each peak.

Figure 13 shows the exact CME solution (magenta dots) and the probability distribution constructed as detailed above from the 2MA (blue curve) for four different volumes ( $\Omega = 150, 200, 250, 300$ ). There is reasonably good agreement between the two probability distribution functions. However, note that to our knowledge there is no method to obtain the weights self-consistently from the MA equations (we estimated these from the CME itself) and hence all one can conclude is that *in the large volume bistable regime, the 2MA provides information on the individual peaks of the bimodal CME distribution but not on the whole distribution itself. It is also the case that there are volumes for which the CME distribution is bimodal but for which the 2MA equations are monostable ( $30 < \Omega < 147$ ), in which case no reconstruction of the probability distribution function by the aforementioned method is possible.*Figure 13: A comparison of the exact probability distribution function (pdf) obtained from CME (magenta dots) and the reconstructed pdf from the 2MA equations (solid blue line) for the Schlögl reaction scheme in (14) with parameter set S1. The volumes are $\Omega = 150, 200, 250, 300$ for (a) - (d) respectively. The 2MA pdf is constructed from the two steady-state solution of the 2MA as detailed in Appendix B.