An alternative route to the system-size expansion

The master equation is rarely exactly solvable and hence various means of approximation have been devised. A popular systematic approximation method is the system-size expansion which approximates the master equation by a generalised Fokker–Planck equation. Here we first review the use of the expansion by applying it to a simple chemical system. The example shows that the solution of the generalised Fokker–Planck equation obtained from the expansion is generally not positive definite and hence cannot be interpreted as a probability density function. Based on this observation, one may also a priori conclude that moments calculated from the solution of the generalised Fokker–Planck equation are not accurate; however calculation shows these moments to be in good agreement with those obtained from the exact solution of the master equation. We present an alternative simpler derivation which directly leads to the same moments as the system-size expansion but which bypasses the use of generalised Fokker–Planck equations, thus circumventing the problem with the probabilistic interpretation of the solution of these equations.


SSA
stochastic simulation algorithm SSE system-size expansion LNA linear-noise approximation GFPE generalised Fokker-Planck equation CME chemical master equation

Introduction
The Markovian description of stochastic systems with discrete state space is generally described by means of a master equation [1]. The solution of the latter gives the probability that the system of interest is in a given configuration at a given time. However, exact solution of the master equation is rarely possible because of the typically large (or infinite) dimensionality of the state space. The most common strategy to bypass this difficulty involves the use of stochastic simulations via the stochastic simulation algorithm (SSA) [2] or one of its several variants. However, such simulations can easily become computationally expensive when modelling realistic systems, thus making analytical approximation procedures an appealing alternative.
One of the most well-known approximation methods is van Kampen's system-size expansion (SSE) [1], a perturbative expansion of the master equation in the inverse system size, for systems which are deterministically monostable. In the limit of infinitely large system sizes, consideration of the leading order terms of the SSE shows that the mean concentrations of the master equation agree with those of the deterministic rate equations, while the distribution of fluctuations about the mean is Gaussian and given by a Fokker-Planck equation with linear drift and diffusion coefficients. The latter is often referred to as the 'linear-noise approximation' (LNA) and is widely used in the literature (see for example [1,[3][4][5]). Consideration of higher orders of the expansion lead to generalised Fokker-Planck equations (GFPE) with third or higher order derivatives [6]. These higher-order terms have recently been used [7][8][9][10][11] to compute corrections to the mean concentration solution of the rate equations and to the second moments given by the LNA (for a review of these results see [12]). These corrections are applicable when the system size is of intermediate size. Although they are computationally advantageous and typically highly accurate when compared to stochastic simulations [13][14][15], the GFPE from which they are computed has been subject to strong criticisms. Pawula's theorem [16,17] states that the solution of the GFPE is a positive quantity only if the GFPE has at most second-order derivatives or else if the GFPE has an infinite number of terms. In other words, while the solution of the Fokker-Planck equation associated with the LNA can be interpreted as a probability density function, the solution of the GFPE with third or higher derivatives can not be so interpreted. Thus the moments computed from such an equation are dubious. Based on such arguments one may dismiss any SSE computations beyond the LNA [18]. Hence the dilemma: why are the moments calculated using higher-order terms beyond the LNA quantitatively accurate despite being derived from a GFPE which has no apparent meaningful probabilistic interpretation?
In this paper, we propose a novel method, which leads to the same moment equations as obtained using the GFPE with third and higher-order derivatives derived from the SSE, but which avoids the problems outlined above. In particular, the method relies on an expansion of the moment equations directly obtained from the master equation and hence bypasses the use of a time-evolution equation for the approximate probability density function. To be specific, we shall consider a general chemical reaction system containing different species interacting via a number of chemical reactions in a well-mixed volume. In this context the master equation is known as the 'chemical master equation' (CME) and the system size corresponds to the volume of the compartment in which the reactions occur. We emphasise however that all results derived in this paper also apply to master equations that allow a SSE but are not of CME type (see for example [19,20]). This paper is structured as follows. In section 2 we introduce a simple, instructive one variable example and show that while the solution of the corresponding GFPE with up to fourthand sixth-order derivatives can become negative, the moments calculated from such a solution are highly accurate, thereby illustrating the aforementioned dilemma. In section 3, we present the general derivation of the SSE leading to a GFPE with up to fourth-order derivatives for a general multi-species chemical system. In section 4 we present our new expansion method starting directly from the moment equations of the CME and which bypass the use of any type of GFPE. Subsequently in section 5 we show that the equations from the conventional SSE derivation of section 3 match exactly the equations derived from the alternative method in section 4. Finally, we briefly summarize our results and conclude in section 6.

An illustrative one species example
In order to demonstrate the probabilistic interpretation problem of GFPEs, we consider the following simple chemical reaction system: The species X gets created at rate k 1 and irreversibly forms dimers with rate k 2 (which we do not follow and hence we do not label). The corresponding CME is given by: where P(n, t) is the probability of n particles of X being in the system at time t, E ±x are the step operators that replace n with n ± x, i.e. E ±x f (n) = f (n ± x) for a function f (n), and Ω is the volume of the compartment. The SSE makes the ansatz [1]: where φ is the solution of the deterministic rate equations. We thus divide the concentration n/Ω into a deterministic contribution φ and a contribution ε which provides corrections to the latter. The corresponding deterministic rate equations are given by: The ansatz in (3) effectively leads to a change of variables from n to ε. We thus have to transform from the distribution P(n, t) over n to a new distribution Π( , t) over ε. Performing the transformation in variables on the CME (2) and taking the limit of large Ω (see [8] for technical details of the general procedure), one finds that the CME in equation (2) can be approximated by the GFPE: Note that truncating terms on the right hand side of the above equation to order Ω 0 leads to the LNA, a Fokker-Planck equation with a drift term linear in ε and a diffusion coefficient which is independent of ε. This equation admits a Gaussian solution for the distribution of fluctuations about the solution of the rate equations. A closed formula for the non-Gaussian solution of the GFPE for any deterministically monostable single species system, truncated to any order, has been derived recently [14]. Applying this general formula, we find that the steady-state solution of the GFPE in equation (5) is given by: The first term π 0 ( ) is the Gaussian solution of the LNA and reads: where σ 2 = 3 4 φ s and φ s = k 1 /2k 2 . The latter is the steady-state solution of the rate equation in (4). The corrections beyond the LNA in (6) are given in terms of the functions ψ m ( ) which are proportional to Hermite polynomials H m as follows: The non-vanishing coefficients a (i) m in equation (6) read: (9) Figure 1 shows the solutions to the GFPE given in equation (6) truncated to orders Ω 0 (red), Ω −1 (blue) and Ω −2 (orange). The latter correspond to the GFPE with at most second-, fourth-and sixth-order derivatives, respectively. Note that figure 1 shows distributions in the variable n, which is obtained by reversing the transformation of variables from distributions in ε to distributions over the variable n .
Note that truncating the GFPE to higher-order leads to a better agreement of the solution of the GFPE with that of the exact solution of the CME for positive n. However this is achieved at the expense of negative probabilities for some negative n values (for other examples, such as the simple birth-death process, the higher-order GFPE solution leads to negative probabilities also for positive n values-see figure 1(a) in [14]). Note that all GFPE's (including the LNA) are not just defined for positive molecule numbers n but also for negative ones. This is in contrast to the CME which only allows transitions between positive molecule numbers. However for the LNA this is not a major cause of concern since the probability is positive and hence the LNA constitutes a stochastic process over the real domain; in contrast, due to negative probabilities, the higher-order GFPE's do not allow a probabilistic interpretation. Curiously, however, the means and variances of the continuous GFPE distributions computed over the open interval n = (−∞, ∞) are found to progressively approach the exact CME value as  (6) obtained by truncating van Kampen's SSE to orders Ω 0 (red, LNA), Ω −1 (blue) and Ω −2 (orange) for the dimerization reaction system in (1). The grey dots show the exact steady-state solution of the CME. The parameters used are k 1 = k 2 = 1 and Ω = 1. Note that while the LNA solution is positive for all n, the higher order solutions (blue, orange) are negative for some values of n. Hence while the LNA has a probabilistic interpretation, the higherorder GFPE solutions do not.
we include higher-order derivatives (see table 1). Hence we have a dilemma: despite being derived from a distribution that is not positive definite, these moments are closer to the true moments of the CME, when compared to the moments calculated from the positive-definite distribution of the LNA.
More generally one finds that for a fixed value of the system size Ω, consideration of higher-order derivatives of the GFPE leads to more accurate estimates of the moments until a certain order, after which the accuracy rapidly decreases; this threshold order increases with the system size and hence strictly speaking, the expansion is valid to all orders only in the limit Ω → ∞ [21]. This also implies that the SSE is divergent for finite Ω. This is typical of series derived from perturbation methods [22]. It is not possible to a priori know the threshold order however generally if the addition of a new term leads to either a large relative change in the value of a moment or if it leads to a negative value of the average of the molecule numbers raised to any power then clearly one should not add further terms to the series.
The accuracy of moments computed from higher order truncations of the GFPE (below the threshold order) leads one to surmise that it may be possible to derive them using an approach which bypasses the use of the GFPE and hence the knotty issue of negative probabilities. The rest of this paper is dedicated to derive such an approach.

The system-size expansion for a general chemical system
We start by reviewing the derivation of the SSE to higher orders introduced in [7]. Consider a system of N chemical species that interact through a set of R chemical reactions where the r th reaction has the form: Here the index r takes values between 1 and R, and X i denotes species i, i = 1, ..., N . The integers s ir and h ir are the stoichiometric coefficients and k r is the macroscopic reaction rate. We can write the corresponding CME in compact form using step operators as [1]: where n = (n 1 , ..., n N ), n i is the molecule number of species X i , E x i is a step operator returning f r ( n, Ω)P( n, t) but with n i replaced by n i + x, and we have defined the stoichiometric matrix S ir = h ir − s ir . The microscopic rate function f r ( n, Ω) depends on the type of the rth chemical reaction. The probability that a reaction occurs in a time interval [t, t + dt) somewhere in the volume Ω is given by Ωf r ( n, Ω)dt. Table 1. Comparison of the mean and variance in molecule numbers calculated from the exact CME distribution and the solutions of the GFPE truncated to orders Ω 0 (LNA), Ω −1 and Ω −2 (the orders of Ω −1/2 and Ω −3/2 do not contribute to the moments). Parameters are as in figure 1 where we also show the corresponding distributions. Note that the mean and variance become more accurate with increasing order. The exact values are computed using formulae (68) and (69) in [7]. We assume mass action kinetics, and specifically we consider reactions involving at most two reactant molecules, i.e. N i=1 s ir 2 ∀r, since reactions involving three or more molecules are relatively rare. This means that we can write the functions f r ( n, t) in (11) in the following form: where we explicitly split up the term into contributions from different types of possible reactions. The first and second term represent zeroth and first order (unimolecular) reactions, respectively. The third term corresponds to a second order (bimolecular) reaction, where the two reactant molecules are of different type, while the last term corresponds to a bimolecular reaction involving two identical reactant molecules. Note that for a given reaction r, only one of the coefficients ∆ r , β k r , γ k s r and ζ k k r in equation (12) is non-zero, depending on the type of reaction, and it is equal to k r .
To derive the SSE equations one makes the following ansatz: Note that this is the same ansatz that we made in equation (3) for the single species system, but this time it is applied to the concentrations of each species. The φ i are the solution of the deterministic rate equations. Next, we transform the CME in equation (11) from the variables n to the variables , leading to a GFPE for the distribution Π( , t). In order to perform this transformation, we transform the time derivative, expand the step operators E x i in powers of Ω −1/2 (for the explicit expression see [7]) and transform the functions f r accordingly. This leads to: where L (i) are differential operators-the reader is referred to [6] for their explicit definitions. Now by the law of mass action, one can deduce that the deterministic rate equations for the system (10) are: where φ is the vector of macroscopic concentrations and f r ( φ) is the macroscopic rate function of the r th reaction. In the macroscopic limit, terms of order Ω 1/2 in equation (14) dominate. However by the assumption that the CME must agree with the rate equations in the macroscopic limit, one finds that terms of order Ω 1/2 on both sides of the transformed CME disappear. Hence the transformed CME equation (14) reads: The form of the transformed CME in equation (16) Substitution of this ansatz in equation (16) and collecting terms of equal orders leads to timeevolution equations for Π j ( , t). From these equations we can subsequently obtain time-evolution equations for the pseudo-moments [ k m ..
Finally the equations of the moments of the concentrations are obtained from the van Kampen ansatz in equation (13) and are given by: where we truncated to order Ω −2 . Note that the short hand notation (i ↔ k) stands for all the expressions of the same form as the one proceeding the notation but with i and k interchanged.
Note that here we used the fact that: which follows from equation (17) and the definition of the pseudo-moments given below it. This procedure can also be generalized to chemical systems with non-mass action kinetics and those reactions involving three or more reactant molecules [13].
In what follows, we shall use the Einstein summation convention where all twice repeated indices in product terms are understood to be summed over 1 to N. This notation will be used throughout the paper when it is necessary to write equations in compact form. To simplify the presentation, we define the following coefficients: Because of the quadratic form of f r ( φ) (which stems from allowing at most bimolecular reactions), the only non-zero J coefficients are those involving at most second-order derivatives, i.e. those with at most two superscripts. These are: where δ k ,p is the Kronecker delta. The time-evolution equations for the pseudo-moments [.] j are needed to compute the moments in concentrations. These equations are obtained as detailed just before equation (18) and can be found in the appendix.
In summary, the estimates of the first, second and third moments of the concentrations up to order Ω −2 according to the SSE are generally given by equations (18)-(20) together with the solution of the closed set of differential equations given by equations (A.1)-(A.10). Equations for the fourth and higher-order moments in the concentrations can be derived similarly. For more details on the derivation we refer the reader to [7,13].

An alternative route
As we saw in the last section, the final output of the SSE as given by equations (18)-(20) is a series expansion of the moments in powers of Ω −1/2 . The open question is whether this expansion is accurate given that it is derived from the GFPE in equation (16) which has higher than second-order derivatives and hence admits no probabilistic interpretation according to Pawula's theorem. To solve this apparent issue, we next formulate an alternative derivation which relies on a direct expansion of the moment equations obtained from the CME in the inverse system size, thereby avoiding the use of a GFPE. The procedure is as follows.
The moment equations are obtained by multiplying the CME in equation (11) on both sides with n i . . . n l and subsequently summing over all allowed values of the molecule numbers. The moment equations of the first four moments read: S kr n i f r ( n, Ω) + S ir n k f r ( n, Ω) + S kr S ir f r ( n, Ω) , where f r ( n, Ω) is the microscopic rate function as defined in equation (12). We now assume that the moments can be written as a series in powers of Ω −1/2 . For the first three moments, for instance, we write where we introduced the coefficients α, ρ and λ. Note that we have here omitted terms of higher order than Ω −2 . The subscripts and superscripts indicate expansion orders and species indices, respectively. Substituting equations (28) into (24) and equating equal order terms on both sides of the equation, one obtains: where δ i,j is the Kronecker delta and α k j = 0 if j < 0. Similarly, one can obtain equations for the coefficients of higher order moments. For instance substituting equations (28) into (25) and equating equal order terms on both sides of the equation, we obtain time-evolution equations for the coefficients ρ ik j of the second order moments: Note that it is here understood that ρ ks j = 0 if j < 0. Equations for the coefficients of third and higher order moments can be similarly derived. As we will show in the following section, equations (29) and (30) agree exactly with the ones obtained in the previous section using the conventional SSE.

Equivalence of the two methods
Comparing equation (28) with (18)- (20) one obtains the following relationship between the coefficients of the system-size expansion and of the alternative method: (31) Equivalence of the two methods is then proved if, given these relations, we can obtain the time-dependent equations for the coefficients of the alternative method, i.e. equations (29) and (30), by starting from the time-dependent equations for the coefficients of the SSE, i.e. equations (A.1)-(A. 10). This is what we seek to establish next.

Equivalence at the deterministic level
It follows directly from the previous relations that the deterministic rate equations given by equations (15) correspond to (29) with j = 0.

Equivalence at the first moment level
We start by recalling the time evolution equation for [ i ] j derived using the SSE (see equations (A.1)-(A.4)): The next step is to substitute the explicit expression of the matrices J as given by equation (23). This leads to the equations: (37) These are exactly the same as equation (29) with j = 1, . . . , 4.

Equivalence at the second moment level
Using the SSE we previously derived the expression for [ i k ] 0 which is given by: Writing [ i k ] 0 and [ w k ] 0 in terms of the new coefficients using equation (31), we obtain: (39) This equation can be expanded as: (40) Next, we substitute the time derivative of α i k (which we derived earlier and is given by equation (37)) in the right hand side of the above equation, leading to: To write the final explicit equation for ρ ik 2 , one needs to insert the expressions for the J matrices given by equations (22) and (23). After some algebraic manipulation most of the terms inside the expression cancel leading to the final form: where we used the expressions for the λ coefficients given in equation (31). Note that equation (42) is precisely one and the same as equation (30) with j = 2. By a completely analogous derivation to the one above, starting from the SSE timeevo lution equations for [ i k ] 1 and [ i k ] 2 as given by equations (A.6) and (A.7), one can derive equation (30) with j = 3 and j = 4, respectively.
Hence in summary we have here proved that the SSE equations for the first two moments of the concentrations up to order Ω −2 are precisely the same as those obtained from the alternative route. The algebra becomes formidable to perform by hand for third and higher-order moments and for higher powers in Ω, but we have carried this out with Mathematica and have verified agreement of the SSE and the alternative method for these cases as well.
We finish this section by noting the difference between our method and moment-closure approximations [12,[23][24][25][26]. Essentially the latter close moment equations by postulating an underlying distribution such as Gaussian or log-normal but in contrast in our method no such artificial closure is used. Rather the equations for the moments close by themselves due to the nature of the perturbation expansion. Hence our procedure is unique and does not suffer from the ad-hoc nature of moment-closure approximations. A comparison of the SSE (to which our method is formally equivalent) and of common moment-closure approximations can be found in [27].

Summary and conclusion
In this paper we have introduced a novel way to approximately calculate the moments of concentrations in the chemical master equation. This approach is useful when traditional exact methods become computationally demanding, e.g. when one is using the Gillespie algorithm and the number of reactions fired per unit time is large or when one is directly solving the CME equations and the transition matrix has a large dimensionality. Our approach developed in section 4 is based on a Taylor series expansion of the moments of molecule numbers in powers of a small parameter-the inverse square root of the system size. Substituting this expansion in the exact equations for the moments derived from the CME, one obtains timeevolution equations for the coefficients of the Taylor series. We showed in section 5 that this approach leads to the same result as the SSE (stated in general in section 3), but in a much simpler and straightforward way. Our approach has one further and most important advantage: unlike the SSE, the moments are not calculated from a GFPE with third and higher-order coefficients and hence it circumvents the problem with the probabilistic interpretation of the GFPE. The fact that it nevertheless agrees with the SSE does prove that the latter is a trustworthy device to approximate the moments.
Note that according to our previously introduced shorthand notation, the terms J il [ k ] 0 + (k ↔ i) + (i ↔ l) in equation (A.8) stand for J il Finally we note that the last equation involves the zeroth order coefficient of the fourth moment, which satisfies the time-evolution equation: (A.10)