Approximate probability distributions of the master equation
AApproximate probability distributions of the master equation
Philipp Thomas ∗ School of Mathematics and School of Biological Sciences, University of Edinburgh
Ramon Grima † School of Biological Sciences, University of Edinburgh
Master equations are common descriptions of mesoscopic systems. Analytical solutions to theseequations can rarely be obtained. We here derive an analytical approximation of the time-dependentprobability distribution of the master equation using orthogonal polynomials. The solution is givenin two alternative formulations: a series with continuous and a series with discrete support both ofwhich can be systematically truncated. While both approximations satisfy the system size expansionof the master equation, the continuous distribution approximations become increasingly negativeand tend to oscillations with increasing truncation order. In contrast, the discrete approximationsrapidly converge to the underlying non-Gaussian distributions. The theory is shown to lead toparticularly simple analytical expressions for the probability distributions of molecule numbers inmetabolic reactions and gene expression systems.
I. INTRODUCTION
Master equations are commonly used to describe fluc-tuations of particulate systems. In most instances, how-ever, the number of reachable states is so large that theircombinatorial complexity prevents one from obtaininganalytical solutions to these equations. Explicit solutionsare known only for certain classes of linear birth-deathprocesses [1], under detailed balance conditions [2], orfor particularly simple examples in stationary conditions[3]. Considerable effort has been undertaken to approx-imate the solution of the master equation under moregeneral conditions including time-dependence and condi-tions lacking detailed balance [4, 5].A common technique addressing this issue was givenby van Kampen in terms of the system size expansion(SSE) [6]. The method assumes the existence of a spe-cific parameter, termed the system size, for which themaster equation approaches a deterministic limit as itsvalue is taken to infinity. The leading order term of thisexpansion describes small fluctuations about this limitin terms of a Gaussian probability density, called thelinear noise approximation (LNA). This approximationhas been widely applied in biochemical kinetics [7], butalso in the theory of polymer assembly [8], epidemics [9],economics [10], and machine learning [11]. The benefitof the LNA is that it yields generic expressions for theprobability density. Its deficiency lies in the fact that,strictly speaking, it is valid only in the limit of infinitesystem size. Hence one generally suspects that its pre-dictions become inaccurate when one studies fluctuationsthat are not too small compared to the mean and there-fore implying non-Gaussian statistics.Higher order terms in the SSE have been employed tocalculate non-Gaussian corrections to the LNA for the ∗ [email protected] † [email protected] first few moments [12–14]; alternative methods are basedon moment closure [15]. It is however the case that theknowledge of a limited number of moments does not al-low to uniquely determine the underlying distributionfunctions. Reconstruction of the probability distributiontherefore requires additional approximations such as themaximum entropy principle [16] or the truncated momentgenerating function [17] which generally yield differentresults. While the accuracy of these repeated approxi-mations remains unknown, analytical expressions for theprobability density can rarely be obtained, or might noteven exist [18]. A systematic investigation of the distri-butions implied by the higher order terms in the SSE,without resorting to moments, is therefore still missing.We here analytically derive, for the first time, a closed-form series expansion of the probability distribution un-derlying the master equation. We proceed by outliningthe expansion of the master equation in Section I andbriefly review the solution of the leading order termsgiven by the LNA in Section II. While commonly theSSE is truncated at this point, we show that the higherorder terms can be obtained using an asymptotic expan-sion of the continuous probability density. The resultingseries is given in terms orthogonal polynomials and canbe truncated systematically to any desired order in theinverse system size. Analytical expressions are given forthe expansion coefficients.Thereby we establish two alternative formulations ofthis expansion: a continuous and a discrete one both sat-isfying the expansion of the master equation. We showthat for linear birth-death processes, the continuous ap-proximation often fails to converge reasonably fast. Incontrast, the discrete approximation introduced in Sec-tion III accurately converges to the true distribution withincreasing truncation order. In Section IV, we show thatfor nonlinear birth-death processes, renormalization isrequired for achieving rapid convergence of the series.Our analysis is motivated by the use of simple examplesthroughout. Using a common model of gene expression,we conclude in Section VI that the new method allows to a r X i v : . [ c ond - m a t . s t a t - m ec h ] O c t predict the full time-dependence of the molecule numberdistribution. II. SYSTEM SIZE EXPANSION
As a starting point, we focus on the master equationformulation of biochemical kinetics. We therefore con-sider a set of R chemical reactions involving a singlespecies confined in a well-mixed volume Ω. Note thatfor chemical systems the system size coincides with thereaction volume. We denote by S r the net-change in themolecule numbers in the r th reaction and by γ r ( n, Ω) theprobability per unit time for this reaction to occur. Theprobability of finding n molecules in the volume Ω at time t , denoted by P ( n, t ), then obeys the master equationd P ( n, t )d t = R (cid:88) r =1 (cid:18) E − S r − (cid:19) γ r ( n, Ω) P ( n, t ) , (1)where E − S r is the step operator defined as E − S r g ( n ) = g ( n − S r ) for any function g ( n ) of the molecule numbers[6]. Note that throughout the article, deterministic initialconditions are assumed. The system size expansion nowproceeds by separating the instantaneous concentrationinto a deterministic part, given by the solution of the rateequations [ X ], and a fluctuating part (cid:15) , n Ω = [ X ] + Ω − / (cid:15), (2)which is van Kampen’s ansatz. The expansion of themaster equation can be summarized in three steps: (i) Using Eq. (2) one expands the step operator E − S r γ r ( n, Ω) P ( n, t ) = γ r ( n − S r , Ω) P ( n − S r , t )= e − Ω − / S r ∂ (cid:15) γ r (Ω[ X ] + Ω / (cid:15), Ω) P (Ω[ X ] + Ω / (cid:15), t ) , (3)where ∂ (cid:15) denotes ∂∂(cid:15) . (ii) Next, the probability for the molecule numbers iscast into a probability density Π( (cid:15), t ) for the fluctuationsusing van Kampen’s ansatz,Π( (cid:15), t ) = Ω / P (Ω[ X ] + Ω / (cid:15), t ) , (4)which is essentially a change of variables. Note thatthis step implicitly assumes a continuous approximationΠ( (cid:15), t ) of the probability distribution as thought in theoriginal derivation of van Kampen [6]. (iii) It remains to expand the propensity about thedeterministic limit γ r (Ω[ X ] + Ω / (cid:15), Ω) = ∞ (cid:88) k =0 Ω − k/ (cid:15) k k ! ∂ k γ r (Ω[ X ] , Ω) ∂ [ X ] k . (5)Note that γ r (Ω[ X ] , Ω) is just the propensity evaluated atthe macroscopic concentration and hence it must depend explicitly on Ω. We assume that the propensity possessesa power series in the inverse volume γ r (Ω[ X ] , Ω) = Ω ∞ (cid:88) s =0 Ω − s f ( s ) r ([ X ]) . (6)For mass-action kinetics, for instance, the propensity isgiven by γ r ( n, Ω) = Ω − (cid:96) r k r (cid:96) r ! (cid:0) n(cid:96) r (cid:1) , where (cid:96) r is the re-action order of the r th reaction. Using the Taylor ex-pansion of the binomial coefficient, we have f (0) r ([ X ]) = k r [ X ] (cid:96) r , f ( s ) r ([ X ]) = k r [ X ] (cid:96) r − s S (cid:96) r ,(cid:96) r − s , and f ( s ) r = 0 for s ≥ (cid:96) r , where S denotes the Stirling numbers of the firstkind. Note also that effective propensities being deducedfrom mass action kinetics have an expansion similar toEq. (6). The Michaelis-Menten propensity γ r ( n, Ω) =Ω k r nn + K Ω [19], for instance, has f (0) r ([ X ]) = k r [ X ][ X ]+ K and f ( s ) r ([ X ]) = 0 for s > − / , we find (cid:18) ∂∂t − Ω / d[ X ]d t ∂∂(cid:15) (cid:19) Π( (cid:15), t )= (cid:18) − Ω / R (cid:88) r =1 S r f (0) r ([ X ]) ∂∂(cid:15) + N (cid:88) k =0 Ω − k/ L k (cid:19) Π( (cid:15), t )+ O (Ω − ( N +1) / ) . (7)Equating terms to order Ω / yields the deterministicrate equation d[ X ]d t = R (cid:88) r =1 S r f (0) r ([ X ]) . (8)The higher order terms in the expansion of the masterequation can be written out explicitly L k = (cid:100) k/ (cid:101) (cid:88) s =0 k − s − (cid:88) p =1 D k − p − s − p,s p !( k − p − s − − ∂ (cid:15) ) p (cid:15) k − p − s − , (9)where (cid:100)·(cid:101) denotes the ceiling value and the coefficientsare given by D qp,s = R (cid:88) r =1 ( S r ) p ∂ q [ X ] f ( s ) r ([ X ]) , (10)which depend explicitly on the solution of the rate equa-tion (8). Note that in the following the abbreviation D qp = D qp, is used. III. EXPANSION OF THE CONTINUOUSPROBABILITY DENSITY
We here study the time-dependent solution of the par-tial differential equation approximation of the masterequation, Eq. (7). We therefore expand the probabil-ity density of Eq. (4),Π( (cid:15), t ) = N (cid:88) j =0 Ω − j/ π j ( (cid:15), t ) + O (Ω − ( N +1) / ) , (11)which also allows the expansion of the time-dependentmoments to be deduced in closed-form. Finally we re-cover the stationary solution as a particular case. A. Linear Noise Approximation
Substituting Eq. (11) into Eq. (7) and equating termsto order O (Ω ) we find (cid:18) ∂∂t − L (cid:19) π = 0 , (12)where L = − ∂ (cid:15) J (cid:15) + ∂ (cid:15) D is a Fokker-Planck operatorwith linear coefficients, and J = D is the Jacobian ofthe rate equation. The probability density of fluctuationsabout the macroscopic concentration, described by (cid:15) , isgiven by a centered Gaussian π ( (cid:15), t ) = 1 (cid:112) πσ ( t ) exp (cid:18) − (cid:15) σ ( t ) (cid:19) , (13)which acquires time-dependence via its variance σ ( t ).The latter satisfies ∂σ ∂t = 2 J ( t ) σ + D ( t ) , (14)which is the familiar LNA result [6]. In the followingwe will drop the time-dependence of the coefficients forconvenience of notation. B. Higher order terms
Substituting Eq. (11) into Eq. (7), rearranging theremaining terms, and equating terms to order Ω − j/ , wefind (cid:18) ∂∂t − L (cid:19) π j ( (cid:15), t ) = L π j − + . . . + L j π = j (cid:88) k =1 L k π j − k ( (cid:15), t ) . (15)This system of partial differential equations can be solvedusing the eigenfunction approach. We consider (cid:18) ∂∂t − L (cid:19) Ψ m = λ m Ψ m , (16)which is solved by λ m = − m J and Ψ m = ψ m ( (cid:15), t ) π ( (cid:15), t )with ψ m ( (cid:15), t ) = π − ( − ∂ (cid:15) ) m π = 1 σ m H m (cid:16) (cid:15)σ (cid:17) . (17)The functions H m denote the Hermite orthogonal poly-nomials which are given explicitly in Appendix A. Toverify the solution of the eigenvalue problem, we setΨ m +1 = ( − ∂ (cid:15) )Ψ m and observe that ( ∂ t − L )Ψ m +1 = −J Ψ m +1 − ∂ (cid:15) ( ∂ t − L )Ψ m . Using this in Eq. (16), weobtain λ m +1 = ( −J + λ m ) from which the result followsbecause λ = 0 and Ψ = π .Using the completeness of the eigenfunctions, we canwrite π j ( (cid:15), t ) = (cid:80) ∞ m =0 a ( j ) m ( t ) ψ m ( (cid:15), t ) π ( (cid:15), t ). We verifyin Appendix B that the j th order term in the expansioninvolves only the first N j = 3 j eigenfunctions. The con-tinuous SSE approximation is consequently given by theasymptotic expansionΠ( (cid:15), t ) = π ( (cid:15), t ) (cid:18) N (cid:88) j =1 Ω − j/ N j (cid:88) m =1 a ( j ) m ( t ) ψ m ( (cid:15), t ) (cid:19) + O (Ω − ( N +1) / ) , (18)for which the coefficients can be determined us-ing the orthogonality of the functions ψ m , i.e., σ n n ! (cid:82) d (cid:15) ψ n ( (cid:15), t ) ψ m ( (cid:15), t ) π ( (cid:15), t ) = δ m,n C. The equation for the expansion coefficients
The coefficients a ( j ) n are now determined by insertingthe expansion of π j into Eq. (15), multiplying the re-sult by σ n n ! (cid:82) d (cid:15) ψ n ( (cid:15), t ), and performing the integration.Using Eq. (16), the left hand side of Eq. (15) becomes σ n n ! (cid:88) m (cid:90) d (cid:15) ψ n ( (cid:15), t ) (cid:18) ∂∂t − L (cid:19) a ( j ) m ψ m ( (cid:15), t ) π ( (cid:15), t )= (cid:18) ∂∂t − n J (cid:19) a ( j ) n . (19)The calculation of terms in the summation on the righthand side of Eq. (15) is greatly simplified by defining theintegral I αβmn = σ n n ! α ! β ! (cid:90) d (cid:15) ψ n ( (cid:15), t )( − ∂ (cid:15) ) α (cid:15) β ψ m ( (cid:15), t ) π ( (cid:15), t ) , (20)which yields σ n n ! (cid:90) d(cid:15)ψ n ( (cid:15), t ) L k ψ m ( (cid:15), t ) π ( (cid:15), t )= (cid:100) k/ (cid:101) (cid:88) s =0 k − s − (cid:88) p =1 D k − p − s − p,s I p,k − p − s − mn . (21)Using Eqs. (19) and (21) in Eq. (15), we find that the co-efficients satisfy the following set of ordinary differentialequations (cid:18) ∂∂t − n J (cid:19) a ( j ) n = j (cid:88) k =1 N j − k (cid:88) m =0 a ( j − k ) m (cid:100) k/ (cid:101) (cid:88) s =0 k − s − (cid:88) p =1 D k − p − s − p,s I p,k − p − s − mn , (22)where we have assumed a ( j ) n = 0 for n > N j . Explicitly,the non-zero integrals are given by I αβmn = σ β − α + n − m α ! min( n − α,m ) (cid:88) s =0 (cid:18) ms (cid:19) × ( β + α + 2 s − ( m + n ) − β + α + 2 s − ( m + n ))!( n − α − s )! , (23)and zero for odd ( α + β ) − ( m + n ). Here (2 k − (2 k )!2 k k ! isthe double factorial. Along with Eq. (23), in AppendixB we verify the following two important properties ofthe asymptotic series solution given deterministic initialconditions: (i) We have N j = 3 j and hence Eq. (22)indeed yields a finite number of equations, and (ii) a ( j ) n vanishes for all times when ( n + j ) is odd.Finally, we note that D qp,s and I pqmn are generally time-dependent because they are functions of the solution ofthe rate equation and the LNA variance. Explicit expres-sions for the approximate probability density can now beevaluated to any desired order. D. Moments of the distribution
The solution for the probability density enables one toderive closed-form expressions for the moments. Theseare obtained by multiplying Eq. (18) by (cid:82) d (cid:15) (cid:15) β and per-forming the integration using Eq. (A4) of Appendix B.We find (cid:104) (cid:15) β (cid:105) = N (cid:88) j =0 Ω − j/ (cid:98) β/ (cid:99) (cid:88) k =0 β !2 k k ! σ k a ( j ) β − k + O (Ω − ( N +1) / ) , (24) where a ( j )0 = δ ,j and (cid:98)·(cid:99) denotes the floor value. Inparticular, it follows that mean and variance are given by (cid:104) (cid:15) (cid:105) = (cid:80) Nj =1 Ω − j/ a ( j )1 + O (Ω − ( N +1) / ) and (cid:104) (cid:15) (cid:105) = σ +2 (cid:80) Nj =1 Ω − j/ a ( j )2 + O (Ω − ( N +1) / ).It is now evident that the coefficients of the expansionare intricately related to the system size expansion of thedistribution moments. Naturally, one may seek to invertthis relation. Indeed, as we show in Appendix C, giventhe expansion for a finite set moments, the coefficientsin Eq. (18) can be uniquely determined. In particular,to construct the probability density to order Ω − j/ onerequires the expansion of the first 3 j moments to thesame order. Thus the problem of moments provides anequivalent route of systematically constructing solutionsto the master equation. E. Solution in stationary conditions
Of particular interest is the expansion of the proba-bility density under stationary conditions. Implicitly, weassume here that the rate equation, Eq. (8), has a sin-gle asymptotically stable fixed point, and hence the LNAvariance is given by σ = D / ( − J ). Setting the time-derivative on the left hand side of Eq. (22) to zero, wefind that the coefficients of Eq. (18) can be expressed interms of lower order ones a ( j ) n = − n J j (cid:88) k =1 (cid:100) k/ (cid:101) (cid:88) s =0 k − s − (cid:88) p =1 ×D k − p − s − p,s j − k ) (cid:88) m =0 a ( j − k ) m I p,k − p − s − mn . (25)For example, truncating after terms of order Ω − , weobtainΠ( (cid:15) ) = π ( (cid:15) ) + Ω − / (cid:16) a (1)1 ψ ( (cid:15) ) + a (1)3 ψ ( (cid:15) ) (cid:17) π ( (cid:15) )+ Ω − (cid:16) a (2)2 ψ ( (cid:15) ) + a (2)4 ψ ( (cid:15) ) + a (2)6 ψ ( (cid:15) ) (cid:17) π ( (cid:15) )+ O (Ω − / ) . (26)The non-zero coefficients to order Ω − / are given by a (1)1 = − σ D J − D , J ,a (1)3 = − σ D J − σ D J − D J , (27) FIG. 1. (Color online)
Linear birth-death process.
Weconsider the reaction system (29) in stationary conditions.(A) We compare the exact Poisson distribution (gray) tothe continuous SSE approximation [Eq. (18) together withEqs. (25) and (30)] truncated after Ω (LNA, blue line), Ω − (green), and Ω − -terms (red) for parameter values k = 0 . k = 1 and Ω = 1 giving half a molecule on average. We ob-serve that the continuous approximation becomes increasinglynegative and tends to oscillations with increasing truncationorder. (B) In contrast the discrete approximation shows nooscillations, and the overall agreement with the exact Poissondistribution (gray bars) improves with increasing truncationorder. while those to order Ω − are a (2)2 = − a (1)1 (cid:32) D , J + D J + 3 σ D J (cid:33) − a (1)3 D J− D , J − σ D , J − σ D J − σ D J ,a (2)4 = − a (1)1 (cid:18) D J + σ D J + σ D J (cid:19) − D J − σ D J− σ D J − σ D J − a (1)3 (cid:32) D , J + 3 D J + 7 σ D J (cid:33) ,a (2)6 = 12 ( a (1)3 ) . (28)The accuracy of this distribution approximation is stud-ied through an example in the following. F. The continuous approximation fails under lowmolecule number conditions
We now study the SSE solution for a linear birth-deathprocess, i.e., its propensities depend at most linearly onthe molecular populations. Specifically, we consider thesynthesis and decay of a molecular species X , ∅ k − (cid:42)(cid:41) − k X. (29)The master equation is constructed using S = +1, γ = Ω k , S = − γ = k n , and R = 2 in Eq. (1).The exact stationary solution of the master equation is aPoisson distribution with mean Ω[ X ] where [ X ] = k /k .The coefficients in Eq. (10) are then given by D mn = δ m, k + ( − n k ( δ m, [ X ] + δ m, ) , (30)and D mn,s = 0 for s >
0. The leading order correctionsto the LNA given by Eqs. (26-28) lead to very compactexpressions for the expansion coefficients and are givenby a (1)3 = [ X ]6 , a (2)4 = [ X ]24 , a (2)6 = [ X ]
72 (31)and a (1)1 = a (2)2 = 0.Though the continuous approximation is expected toperform well at large values of Ω, we are particularlyinterested in its performance when the value of Ω is de-creased. Since the expansion is carried out at constantaverage concentration, low values of Ω typically implylow numbers of molecules and non-Gaussian distribu-tions. In Fig. 1A we show that for parameters yieldinghalf a molecule on average, the continuous approxima-tion obtained in this section, given by Eq. (18) togetherwith Eqs. (25) and (30), is unsatisfactory since as higherorders are taken into account, one observes large oscil-lations in the tails of the distribution. In the followingsection we show that the disagreement arises due to theassumption that the support of the distribution is con-tinuous rather than discrete as implied by the masterequation. IV. DISCRETE APPROXIMATION OF THEPROBABILITY DISTRIBUTION
The aim of this paragraph is to establish a discrete for-mulation of the distribution approximations.. To clarifythis issue, we note that the exact characteristic function G ( k, t ) = (cid:80) ∞ n =0 e ikn P ( n, t ) is a 2 π -periodic function, andhence can be inverted as follows P ( n, t ) = (cid:90) π − π d k π e − ikn G ( k, t ) . (32)We now associate our continuous approximation, Eq.(18), with this characteristic function, i.e., G ( k, t ) = (cid:82) ∞−∞ d (cid:15) e ik Ω([ X ]+Ω − / (cid:15) ) Π( (cid:15), t ). Substituting this togetherwith Eq. (11) into Eq. (32) one establishes a connectionformula between these discrete and continuous approxi-mations via the convolution P ( n, t ) = N (cid:88) j =0 Ω − j/ (cid:90) ∞−∞ d (cid:15) K ( n − Ω[ X ] − Ω / (cid:15) ) π j ( (cid:15), t )+ O (Ω − ( N +1) / ) , (33)with kernel K ( s ) = (cid:90) π − π d k π e − iks = sin( πs ) πs . The convolution can be used to define the derivatives ofthe discrete probability via ∂ n P ( n, t ) = (cid:90) ∞−∞ d (cid:15) K ( n − Ω[ X ] − Ω / (cid:15) )(Ω − / ∂ (cid:15) )Π( (cid:15), t ) , (34)and hence it satisfies E − S j P ( n, t ) = (cid:82) ∞−∞ d (cid:15) K ( n − Ω[ X ] − Ω / (cid:15) ) e − Ω − / ∂ (cid:15) S j Π( (cid:15), t ), as well as γ j ( n, Ω) P ( n, t ) = (cid:82) ∞−∞ d (cid:15) K ( n − Ω[ X ] − Ω / (cid:15) ) γ j (Ω[ X ]+Ω / (cid:15), Ω)Π( (cid:15), t ) foranalytic γ j . It then follows from the fact that P ( n, t ) andΩ / Π(Ω − / ( n − Ω[ X ]) , t ) have the same characteristicfunction expansion, that (i) both approximations possessthe same asymptotic expansion of their moments, andthat (ii) they satisfy the same expansion of the masterequation.For example, to leading order Ω , Eq. (33) replacesthe conventional continuous LNA estimate, π given byEq. (13), with a discrete approximation P ( n, t ) = 12 e − y √ π Σ (cid:20) erf (cid:18) iy + π Σ √ (cid:19) − erf (cid:18) iy − π Σ √ (cid:19)(cid:21) , (35)where y = n − Ω[ X ], Σ = Ω σ is the LNA’s estimatefor the variance of molecule numbers, and erf is the errorfunction defined by erf( x ) = √ π (cid:82) x e − t d t .Associating the Ω − j/ -term of Eq. (18) with π j in Eq.(33), higher order approximations can now be obtainedfrom P ( n, t ) = P ( n, t )+ N (cid:88) j =1 Ω − j/ j (cid:88) m =1 a ( j ) m (cid:16) − Ω / ∂ n (cid:17) m P ( n, t )+ O (Ω − ( N +1) / ) . (36) The above follows from the definition of the eigenfunc-tions, Eq. (17), and using the derivative property of theconvolution given after Eq. (34). Note that the coeffi-cients in this equation are exactly the same as given inEq. (18) and hence are determined by Eq. (22). Onecan verify two limiting cases: (i) as Σ → X ] be-ing integer-valued, then P ( n ) = K ( n − Ω[ X ]) = δ n, Ω[ X ] is just the Kronecker delta as required for determinis-tic initial conditions; (ii) as Ω → ∞ with y/ Σ constant,the probability distribution P reduces to the density π given by Eq. (13) and hence it follows that in this limitthe continuous and discrete series give the same results. A. The discrete approximation performs well forlinear birth-death processes
For the linear birth-death process in the previous sec-tion, in Fig. 1B we show that the discrete approximationgiven by Eq. (36) with Eq. (31) is in good agreementwith the true distribution when truncated after terms oforder Ω − and shows no oscillations. This agreement isremarkable given the compact form of the solution givenby Eq. (31) and (36). The approximation is almost in-distinguishable from the exact result when the series istruncated after Ω − -terms using Eqs. (25) and (30) inEq. (36). We hence conclude that the discrete seriesapproximates better the underlying distribution of themaster equation than the continuous approximation. B. The discrete approximation fails for non-linearbirth processes
Next, we turn our attention to the analysis of nonlin-ear birth-death processes, i.e., a process whose propensi-ties depend nonlinearly on the number of molecules. Aparticular feature of such processes is that the LNA es-timates for mean and variances are generally no longerexact, but agree with those of the true distribution onlyin the limit of large system size [13].Exemplary, we here consider a simple metabolic re-action confined in a small subcellular compartment ofvolume Ω with substrate input, ∅ h −→ S, (37a) S + E h − (cid:42)(cid:41) − h C h −→ E. (37b)The reactions describe the input of substrate molecules S and their catalytic conversion by enzyme species E via the enzyme-substrate complex C . The SSE of theaverage concentrations correcting the macroscopic rateequations have been extensively studied [12]. Since ourtheory applies to a single species only, we here considera reduced model in which reaction (37b) is modelled viaan effective propensity: this gives S = +1, γ = Ω k , FIG. 2. (Color online)
Nonlinear birth-death process.
A metabolic reaction with Michaelis-Menten kinetics, scheme (37),is studied using the reduced model described in Sec. IV B. The exact stationary distribution is a negative binomial (shown ingray). (A) The discrete SSE approximation given by Eq. (36) with Eq. (25) and (30) is shown in the low molecule numberregime ( k /k = 0 .
25, 1 molecule on average) when truncated after Ω (blue), Ω − / (green) and Ω − -terms (red dots). Weobserve that the expansion tends to oscillations and negative values of probability as the truncation order is increased. (B)Similar oscillations are observed for moderate molecule numbers ( k /k = 0 .
9, 27 molecules on average) for the discrete seriestruncated after Ω (blue), Ω − / (green) and Ω − -terms (red lines). In (C) and (D) we show the approximations correspondingto the same parameters used in (A) and (B), respectively, but obtained using the renormalization procedure given by Eq.(41) with Eq. (42) as described in the main text. The renormalized approximations avoid oscillations and are in excellentagreement with the true probability distributions (gray bars). We note that for the cases (B) and (D) the continuous anddiscrete approximations give essentially the same result. The remaining parameters are given by Ω = 10 and K = 0 . and S = − γ = Ω k nn +Ω K . This simplification isvalid when the enzyme-substrate association is in rapidequilibrium, which holds when [ E T ] (cid:28) K and h (cid:28) h where [ E T ] is the total enzyme concentration [19]. Theparameters in the reduced model are related to those inthe developed model by k = h , k = h [ E T ], and K = h /h . This reduced master equation is solved exactlyby a negative binomial distribution [20].The system size coefficients are obtained from Eq.(10), and are given by D mn = δ m, k + ( − n k ∂ m ∂ [ X ] m [ X ] K + [ X ] , (38)and D mn,s = 0 for s >
0. In Fig. 2A and 2B, we con-sider two parameter sets corresponding to low and mod-erate numbers of substrate molecules, respectively. Weobserve that in contrast to the linear case, the discreteapproximation of the nonlinear birth-death process tendsto oscillate with increasing truncation order. This issueis addressed in the following section.
V. RENORMALIZATION OF NONLINEARBIRTH-DEATH PROCESSES
Van Kampen’s ansatz, Eq. (2), bears the particularlysimple interpretation that for linear birth-death processes (cid:15) denotes the fluctuations about the average given by thesolution of the rate equation [ X ]. As noted in the pre-vious example, for nonlinear birth-death processes theseestimates are only approximate. Their asymptotic seriesexpansions will therefore require additional terms thatcompensate for the deviations of the LNA from the trueconcentration mean and variance. It would therefore bedesirable to find an approximation for nonlinear processesthat yields more accurate mean and variance than theLNA. For instance by rewriting van Kampen’s ansatz as n Ω = [ X ] + Ω − / (cid:104) (cid:15) (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) mean + Ω − / ¯ (cid:15) (cid:124) (cid:123)(cid:122) (cid:125) fluctuations . (39)Here, ¯ (cid:15) = (cid:15) − (cid:104) (cid:15) (cid:105) denotes a centered variable that quan-tifies the fluctuations about the true average which is apriori unknown. These estimates can however be approx-imated using the SSE beforehand, and the asymptotic ex-pansion of the distributions can then be performed aboutthese new estimates. This idea is called renormalizationand makes use of the fact that the terms correcting meanand variances can be summed exactly. As we show in thefollowing the resummation allows to better control theconvergence by effectively reducing the number of termsin the summation while at the same time it retains theaccuracy of the expansion.The system size expansion of the moments, Eq. (24),yields the following estimates for mean and variance ofthe fluctuations (cid:104) (cid:15) (cid:105) = N (cid:88) j =0 Ω − j/ a ( j )1 + O (Ω − ( N +1) / ) , (40a)¯ σ = σ + N (cid:88) j =1 Ω − j/ σ j ) + O (Ω − ( N +1) / ) , (40b)respectively, where ¯ σ j ) = 2( a ( j )2 − B j, ( { χ ! a ( χ )1 } j − χ =1 ) /j !)and B j,n are the partial Bell polynomials [21].The renormalization procedure amounts to replacing y by ¯ y = ( n − Ω[ X ] − Ω / (cid:104) (cid:15) (cid:105) ), Σ by ¯Σ = Ω¯ σ in Eq.(35) and associating a new Gaussian ¯ P ( n ) with theseestimates. The renormalized expansion is then given by P ( n, t ) = ¯ P ( n, t )+ N (cid:88) j =1 Ω − j/ j (cid:88) m =1 ¯ a ( j ) m (cid:16) − Ω / ∂ n (cid:17) m ¯ P ( n, t )+ O (Ω − ( N +1) / ) , (41)where the renormalized coefficients can be calculatedfrom the bare ones using¯ a ( j ) m = j (cid:88) k =0 3 k (cid:88) n =0 a ( k ) n κ ( j − k ) m − n , (42)and κ ( n ) j = 1 n ! (cid:98) j/ (cid:99) (cid:88) m =0 ( − ( j + m ) n − m (cid:88) k = j − m (cid:18) nk (cid:19) × B k,j − m (cid:16)(cid:110) χ ! a ( χ )1 (cid:111)(cid:17) B n − k,m (cid:18)(cid:26) χ !2 ¯ σ χ ) (cid:27)(cid:19) , (43)where again B k,n ( { x χ } ) denote the partial Bell polyno-mials [21]. The result is verified at the end of this section.Note that the renormalized series has generally less non-zero coefficients since by construction ¯ a ( j )1 = ¯ a ( j )2 = 0.Note that for linear birth-processes, mean and varianceare exact to order Ω (LNA), and hence for this caseexpansion (36) coincides with Eq. (41).For example, truncating after Ω − -terms, from Eq.(40) it follows that (cid:104) (cid:15) (cid:105) = Ω − / a (1)1 + O (Ω − / ) and ¯ σ = σ +Ω − (2 a (2)2 − ( a (1)1 ) )+ O (Ω − / ). Using Eq. (42)the renormalized coefficients can be expressed in termsof the bare ones¯ a (1)1 = 0 , ¯ a (1)3 = a (1)3 , (44a)¯ a (2)2 = 0 , ¯ a (2)4 = a (2)4 − a (1)1 a (1)3 , ¯ a (2)6 = a (2)6 . (44b)This result can for instance be used to renormalize thestationary solution using the bare coefficients given inSec. III E, Eqs. (27-28). The non-zero renormalizedcoefficients evaluate to¯ a (1)3 = − σ D J + σ D J + D J , (45a)¯ a (2)4 = − D J − σ D J − σ D J − σ D J− ¯ a (1)3 (cid:18) D J + 3 σ D J (cid:19) , ¯ a (2)6 = 12 (¯ a (1)3 ) . (45b)Note that for linear birth-death processes D mn,s = 0 for s > m >
1, and hence the above equations reduceto Eqs. (27-28).
A. The renormalized approximation performs wellfor nonlinear birth-death processes
For the metabolic reaction (37), mean and variancecan be obtained to be (cid:104) (cid:15) (cid:105) = Ω − / ς + O (Ω − ), ¯ σ = σ + Ω − ς ( ς + 1) + O (Ω − ), where ς = [ X ] /K is thereduced substrate concentration and σ = Kς ( ς + 1).Substituting now Eq. (38) into Eqs. (45), we obtain theexpansion coefficients¯ a (1)3 = σ ς + 1) , (46a)¯ a (2)4 = σ
24 (6 ς ( ς + 1) + 1) , ¯ a (2)6 = 12 (¯ a (1)3 ) , (46b)which determine the renormalized series expansion to or-der Ω − . Using Eq. (25), (38) and (42) we can give thenext order terms to order Ω − / analytically¯ a (3)3 = ¯ a (1)3 K , ¯ a (3)5 = ¯ a (1)3
20 (12 ς ( ς + 1) + 1) , ¯ a (3)7 = ¯ a (1)3 ¯ a (2)4 , ¯ a (3)9 = 16 (¯ a (1)3 ) . (46c)In Fig. 2C and 2D we compare the renormalized approx-imation given by Eq. (41) with the respective bare ap-proximations in Fig. 2A and 2B. We observe that therenormalization technique avoids oscillations and eventhe simple analytical approximation given by Eqs. (46) isin reasonable agreement with the exact result. We notethat the asymptotic approximations shown in C and Dare almost indistinguishable for higher truncation orders. B. Proof of the renormalization formula
The renormalized coefficients can in principle be ob-tained by matching the expansions given by Eq. (36) and(41) via their characteristic functions. For conveniencewe consider the characteristic function of the series (18) G ( k ) = G ( k ) ∞ (cid:88) j =1 Ω − j/ j (cid:88) n =1 a ( j ) n ( ik ) n , (47)with G ( k ) = e − ( kσ ) / being the characteristic functionsolution of the LNA π ( (cid:15) ) and we have omitted the ex-plicit time-dependence to ease the notation. We are nowlooking for a different expansion with corrected estimatesfor the mean and variance.¯ G ( k ) = ¯ G ( k ) ∞ (cid:88) j =1 Ω − j/ j (cid:88) n =1 ¯ a ( j ) n ( ik ) n , (48)Note that ¯ G ( k ) = e ik (cid:104) (cid:15) (cid:105) e − ( k ¯ σ ) / is the characteristicfunction for a Gaussian random variable with mean (cid:104) (cid:15) (cid:105) and variance ¯ σ given by Eqs. (40).Equating now Eq. (47) and (48), we find1+ ∞ (cid:88) j =1 Ω − j/ j (cid:88) n =1 ¯ a ( j ) n ( ik ) n = G ( k )¯ G ( k ) ∞ (cid:88) j =1 Ω − j/ j (cid:88) n =1 a ( j ) n ( ik ) n . (49)Expanding the prefactor in the above equation in powersof k and then in Ω, we have G ( k )¯ G ( k ) = ∞ (cid:88) j =0 ( ik ) j κ j = ∞ (cid:88) n =0 Ω − n/ n (cid:88) j =0 ( ik ) j κ ( n ) j , (50)from which Eq. (42) follows, which expresses the new co-efficients ¯ a ( j ) n in terms of the bare ones a ( j ) n . It remains toderive an explicit expression for the κ ( n ) j . The expansionin powers of ( ik ) yields κ j = (cid:98) j/ (cid:99) (cid:88) m =0 ( − ( j + m ) (cid:104) (cid:15) (cid:105) j − m ( j − m )! (cid:16) ¯ σ − σ (cid:17) m m ! . (51)We now expand the first term in inverse powers of Ωusing the partial Bell polynomials1( j − m )! (cid:32) ∞ (cid:88) n =1 Ω − n/ a ( n )1 (cid:33) j − m = ∞ (cid:88) n =1 Ω − n/ n ! n (cid:88) k =0 δ j − m,k B n,k (cid:16)(cid:110) χ ! a ( χ )1 (cid:111)(cid:17) , (52) and similarly for the second term1 m ! (cid:32) ∞ (cid:88) n =1 Ω − n/ ¯ σ n ) (cid:33) m = ∞ (cid:88) n =0 Ω − n/ n ! n (cid:88) k =0 δ m,k B n,k (cid:18)(cid:26) χ !2 ¯ σ χ ) (cid:27)(cid:19) . (53)Using the above expansions in Eq. (51) and rearrangingin powers of Ω − / , Eq. (43) for the coefficients κ ( n ) j follows.Finally, one associates with the centered variable ¯ (cid:15) = (cid:15) − (cid:104) (cid:15) (cid:105) , a Gaussian ¯ π (¯ (cid:15) ) with variance ¯ σ . It thenfollows from inverting Eq. (48) that Π(¯ (cid:15) ) = ¯ π (¯ (cid:15) ) + (cid:80) Nj =1 Ω − j/ (cid:80) jn =1 ¯ a ( j ) n ψ n (¯ (cid:15) ) ¯ π (¯ (cid:15) ) + O (Ω − ( N +1) / ). As-sociating now the Ω − j/ -term of this equation with π j inEq. (33), the discrete series for P ( n, t ) given by Eq. (41)follows. VI. APPLICATION
The models studied so far have been useful to developthe method. It remains however to be demonstrated thatit remains accurate in cases where analytical solution isnot feasible, as for instance, for out-of-steady-state andnon-detailed balance systems. We here consider the syn-thesis of a protein P which is degraded through an en-zyme ∅ h −→ M h −→ ∅ , M h −→ M + P, (54a) P + E h − (cid:42)(cid:41) − h C h −→ E, (54b)where M denotes the transcript, E the enzyme and C complex species as has been studied in Ref. [22]. Sinceour theory applies only to a single species, we consider thelimiting case in which the protein dynamics representsthe slowest timescale of the system. It has be shown [23]that when species M is degraded much faster than theprotein P , the protein synthesis (54a) reduces to the tran-sition S = + z , γ = Ω k in which z is a random variablefollowing the geometric distribution ϕ ( z ) = b (cid:16) b b (cid:17) z with average b , which is called the burst approximation.Similarly to the metabolic reaction studied in Sec. IV,the enzymatic degradation process (54b) can be reducedto S = − γ = Ω k n Ω K + n with a nonlinear depen-dence on the protein number n . The master equationdescribing the protein number is then given bydd t P ( n ) =Ω ∞ (cid:88) z =0 ( E − z − k ϕ ( z ) P ( n )+ Ω( E +1 − k n Ω K + n P ( n ) . (55)0 FIG. 3. (Color online)
Predicting transient distributionsof gene expression.
The dynamics of protein synthesiswith enzymatic degradation, scheme (54), is studied usingthe burst approximation (55). (A) We compare the time-dependence of the renormalized discrete approximations toexact stochastic simulations at times 1, 2, and 14 min . Theoverall shape (mode, skewness, distribution tails) of the sim-ulated distributions (bars) is in excellent agreement with theseries approximation when truncated after Ω − -terms (solidlines) but not when only Ω are taken into account (dashedlines). This agreement is also observed for the first two mo-ments shown in the inset: while the Ω − -approximation (bluesolid line) agrees with the moment dynamics of the simu-lated distributions (dots) of the reduced model (55), the Ω -approximation underestimates the mean (gray solid line) andvariance by 25%. The area within one standard deviation ofthe mean obtained from simulations is shown in blue, theboundary obtained from the approximations are shown asdashed lines (Ω grey, Ω − blue). (B) Despite the good agree-ment shown in (A) we found that at very short times (12 s –blue solid line) the series truncated after Ω − -terms tendsto oscillations which quickly disappear for later times (24 s –green, 48 s – red solid line). See main text for discussion. Pa-rameters are k Ω = 15 min − , k Ω = 100 min − , K Ω = 20,Ω = 100, and b = 5. Histograms were obtained from 10 , The relation between the parameters in the reduced andthe developed model are given by k = h h /h , b = h /h , k = h [ E T ], K = h /h , where [ E T ] denotes thetotal enzyme concentration. This description involvescountably many reactions: one for the degradation ofthe protein, and one for each value of z . Therefore, thereactions cannot obey detailed balance in steady state.The system size coefficients now follow from Eq. (10), and are given by D mn = δ m, k (cid:104) z n (cid:105) ϕ + ( − n k ∂ m ∂ [ X ] m [ X ] K + [ X ] , (56)and D mn,s = 0 for s >
0, where (cid:104) z n (cid:105) ϕ = (cid:80) ∞ z =0 z n ϕ ( z ) = b Li − n ( b b ) denotes the average over the geometricdistribution in terms of the polylogarithm function [24].The deterministic equation is given byd[ X ]d t = k b − k [ X ] K + [ X ] , (57)which follows from the expression for D . Using the Ja-cobian J = D and D in Eq. (14), we find that theLNA variance obeys ∂σ ∂t = − k K ([ X ] + K ) σ + k b (1 + 2 b ) + k [ X ] K + [ X ] . (58)The ODEs given by Eq. (57) and (58) are integratednumerically and the solution is used in Eq. (35) fromwhich the leading order approximation follows. Higherorder approximations are now be obtained by using Eq.(56) in (22) which govern the time-evolution of the coef-ficients a ( j ) m ( t ) and using the result in Eq. (41) and (42).We assume deterministic initial conditions with zero pro-teins meaning a ( j ) m (0) = δ m, δ j, . In Fig. 3A we comparethe time-evolution obtained by the leading order approx-imation P and Eq. (41) truncated after the Ω − -term.The latter distributions are in excellent agreement withthe distributions sampled using the stochastic simulationalgorithm [25]. In particular, unlike the leading orderapproximation, these describe well mode, skewness, andtails of the distribution. We note that also the meanand variance of these distribution approximations are inexcellent agreement as verified in inset of Fig. 3A.Despite the overall good agreement, in Fig. 3B weshow that there are discrepancies at very short timeswhere and, again, the distribution approximations tendto oscillations. Motivated by this numerical observation,we speculate that this behavior of the expansion is due atemporal boundary layer as commonly observed in singu-lar perturbation expansions [26]. Theoretically, the layermust be located at times of the same order as the expan-sion parameter, i.e., t = (Ω K ) − / min ≈ s , coincidingwith the simulation in Fig. 3B. This suggests that ourapproach does only describe the outer solution. Furtheranalysis would be required to investigate also the innersolution which is beyond the scope of this article. VII. DISCUSSION
We have here presented an approximate solutionmethod for the probability distribution of the master1equation. The solution is given in terms of an asymp-totic series expansion that can be truncated systemati-cally to any desired order in the inverse system size. Forbiochemical systems with large numbers of molecules, wehave derived a continuous series approximation that ex-tends van Kampen’s LNA to higher orders in the SSE. Inlow molecule number conditions, we have found that thiscontinuous approximation becomes inaccurate. Instead,in most practical situations the prescribed discrete distri-bution approximations incorporating higher order termsin the SSE better capture the underlying solution of themaster equation. While the terms to order Ω − have beengiven explicitly, we found that for the examples studiedhere up to Ω − or Ω − -terms had to be taken into accountto accurately characterize these non-Gaussian distribu-tions. We note, however, that the asymptotic expansioncannot generally guarantee the positivity of the probabil-ity law. These undulations are particularly pronouncedin the short-time behavior of the expansion studied inSec. VI, which our theory does not describe.Previous means of solving the master equation haveeither been numerical in nature [27] or have focused onthe inverse problem, i.e., reconstruction of the proba-bility density from the moments. While a numericalsolution for the master equation of a single species israther straightforward, we expect our procedure to be-come computationally advantageous when generalized tothe multivariate case where numerical solution is usuallyprohibitive because of combinatorial explosion.Methods based on moments typically require approx-imations such as moment closure [16] and also requirethe prior assumption of the first few moments contain-ing all information on the probability distribution. Con-versely, using the system size expansion, we have hereobtained the probability distribution directly from themaster equation without the need to resort to moments .This method enjoys the particular advantage over previ-ous ones that the first few terms of this expansion canbe written down explicitly as a function of the rate con-stants and for any number of reactions. For small modelswe have demonstrated that the procedure leads to par-ticularly simple expressions for the non-Gaussian distri-butions. This development could prove particularly valu-able for parameter estimation of biochemical reactions inliving cells. ACKNOWLEDGMENTS
It is a pleasure to thank Claudia Cianci and DavidSchnoerr for careful reading of the manuscript.
Appendix A: Useful properties of the Hermitepolynomials
We here briefly review some properties of the Hermiteorthogonal polynomials. The polynomials can be defined in terms of the derivatives of a centered Gaussian π withvariance σ , H n (cid:16) (cid:15)σ (cid:17) = π − ( (cid:15) )( − σ∂ (cid:15) ) n π ( (cid:15) ) . (A1)An explicit formula is H n (cid:16) (cid:15)σ (cid:17) = (cid:98) n/ (cid:99) (cid:88) k =0 (cid:18) n k (cid:19) ( − k (cid:16) (cid:15)σ (cid:17) n − k (2 k − . (A2)These functions are orthogonal n ! (cid:82) ∞−∞ d (cid:15) H m (cid:0) (cid:15)σ (cid:1) H n (cid:0) (cid:15)σ (cid:1) π ( (cid:15) ) = δ nm , with respectto the Gaussian measure π . The derivative satisfies( σ∂ (cid:15) ) m H n (cid:16) (cid:15)σ (cid:17) = n !( n − m )! H n − m (cid:16) (cid:15)σ (cid:17) . (A3)Since these polynomials are complete, every function f ( (cid:15) )in L ( R , π ) (not necessarily positive) can be expandedas f ( (cid:15) ) = (cid:80) ∞ n =0 b n H n (cid:0) (cid:15)σ (cid:1) π ( (cid:15) ) , where the coefficientsare given by b n = n ! (cid:82) d (cid:15) H n (cid:0) (cid:15)σ (cid:1) f ( (cid:15) ). We note because H (cid:0) (cid:15)σ (cid:1) = 1 and π is normalized, we must have b = 1 if (cid:82) d (cid:15) f ( (cid:15) ) = 1. Appendix B: Explicit derivation of Eq. (23) and theproperties of the expansion coefficients
Changing variables (cid:15) = xσ and letting ˜ I αβmn = σ α − β + m − n I αβmn , the integral (20) can be written˜ I αβmn = 1 n ! α ! β ! (cid:90) d xH n ( x ) ( − ∂ x ) α x β H m ( x ) π ( x ) , (A1)where π ( x ) is a centered Gaussian with unit variance.Using partial integration, property (A3), and the relation H α ( x ) H β ( x ) = α ! β ! min( α,β ) (cid:88) s =0 H α + β − s ( x ) s !( α − s )!( β − s )! , (A2)given in Ref. [28], one obtains˜ I αβmn = 1 α ! β ! min( n − α,m ) (cid:88) s =0 (cid:18) ms (cid:19) (cid:82) d x x β H m + n − α − s ( x ) π ( x )( n − α − s )! . (A3)The remaining integral can now be evaluated as the mo-ments of the unit Gaussian which yields (cid:90) d x x b H a ( x ) π ( x ) = b !( b − a )! (cid:90) d x x b − a π ( x )= b !( b − a )! ( b − a − . (A4)2for even ( b − a ) ≥ I αβmn = 1 α ! min( n − α,m ) (cid:88) s =0 (cid:18) ms (cid:19) × ( β + α + 2 s − ( m + n ) − β + α + 2 s − ( m + n ))!( n − α − s )! , (A5)for even ( α + β ) − ( m + n ) and zero otherwise. Note thatthe above quantity is strictly positive. Note also that theargument of the double factorial is taken to be positiveand hence the summation is non-zero only if α + β +2 min( n − α, m ) ≥ m + n and hence for even β = 2 k wehave n = m + α ± l, while for odd β = (2 k + 1) we have n = m + α ± (2 l + 1), with l = 0 , . . . , k .The integral formula can be used to verify two impor-tant properties of the solution of Eq. (22) given determin-istic initial conditions: (i) We have N j = 3 j and henceEq. (22) indeed yields a finite number of equations. (ii)The coefficients a ( j ) n for which ( n + j ) is odd vanish at alltimes.To verify property (i), let N j be the index of the highesteigenfunction required to order Ω − j/ . Using Eq. (22)one can show that a ( j ) N j ∼ a ( j − N j − I p, − pN j − ,N j for p ∈ { , , } .By virtue of the properties given after Eq. (23), we find N j = N j − + 3. Since for deterministic initial conditionswe have N = 0, it follows that N j = 3 j .Finally, we verify property (ii). To the summationin Eq. (22) there contribute only terms for which I p,k − p − s − mn is non-zero. Hence, by the condition givenafter Eq. (23), k − ( m + n ) is an even number. Consider-ing the equation for a ( j ) n for which n + j is even, it followsthat in the summation on the right hand side of Eq. (22)there appear only coefficients for which m + ( j − k ) iseven. Conversely, for n + j being odd then same holdsfor m + ( j − k ). Hence the pairs of equations for a ( j ) n for which ( j + n ) is even or odd are mutually uncoupled.For deterministic initial conditions, only terms with j + n even differ from zero initially from which the result fol-lows. Appendix C: Solution to the problem of momentsusing the system size expansion
Having obtained the moment expansion in terms of thecoefficients a ( j ) n , it would be desirable to invert this rela-tion and the coefficients in terms of the expansion of themoments. This can be derived using the completenessof the Hermite polynomials, and writing the probabil-ity density as Π( (cid:15) ) = (cid:80) ∞ n =0 b n H n (cid:0) (cid:15)σ (cid:1) π ( (cid:15) ) , where the b n = n ! (cid:82) d(cid:15)H n (cid:0) (cid:15)σ (cid:1) Π( (cid:15) ) can be expressed in terms of themoments using Eq. (A2), as follows b n = 1 n ! (cid:98) n/ (cid:99) (cid:88) k =0 (cid:18) n k (cid:19) ( − k (cid:104) (cid:15) n − k (cid:105) σ n − k (2 k − . (A1)Assuming now that the moments can be expanded in aseries in powers of Ω, i.e., (cid:104) (cid:15) β (cid:105) = N (cid:88) j =0 Ω − j/ [ (cid:15) β ] j + O (Ω − ( N +1) / ) , (A2)the b n can be matched to the coefficients a n in Eq.(18) using σ n b n = (cid:80) Nj =0 Ω − j/ a ( j ) n + O (Ω − ( N +1) / ), fromwhich one obtains a ( j ) n = 1 n ! (cid:98) n/ (cid:99) (cid:88) k =0 (cid:18) n k (cid:19) ( − σ ) k (2 k − (cid:15) n − k ] j , (A3)with [ (cid:15) ] j = δ j, . The above formula relates the expan-sion of the moments to the expansion of distribution func-tions. It is now evident that the system size expansion ofthe distribution can be constructed from the system sizeexpansion for a finite set of moments.Specifically, to order Ω − / the non-zero coefficientsevaluate to a (1)1 = [ (cid:15) ] , a (1)3 = 13! (cid:0) [ (cid:15) ] − σ [ (cid:15) ] (cid:1) (A4)while the coefficients to order Ω − are given by a (2)2 = 12 [ (cid:15) ] , a (2)4 = 14! (cid:0) [ (cid:15) ] − σ [ (cid:15) ] (cid:1) ,a (2)6 = 16! (cid:0) σ [ (cid:15) ] − σ [ (cid:15) ] + [ (cid:15) ] (cid:1) . (A5)A different series is obtained using the Edgeworth expan-sion which instead of using the system size expansion ofthe moments, Eq. (A2), proceeds by scaling the cumu-lants by a size parameter.3 [1] T. Jahnke and W. Huisinga, J Math Biol , 1 (2007).[2] H. Haken, Phys Lett A , 443 (1974); N. G. Van Kam-pen, ibid . , 333 (1976).[3] R. M. Mazo, J Chem Phys , 4244 (1975); J. Peccoudand B. Ycart, Theor Popul Biol , 222 (1995); P. Bokes,J. R. King, A. T. Wood, and M. Loose, J Math Biol ,829 (2012).[4] R. G¨ortz and D. Walls, Z Phys B Con Mat , 423 (1976);G. Haag and P. H¨anggi, ibid . , 411 (1979).[5] D. Schnoerr, G. Sanguinetti, and R. Grima, J ChemPhys , 024103 (2014).[6] N. G. Van Kampen, Adv Chem Phys , 245 (1976); Stochastic Processes in Physics and Chemistry. , 3rd ed.(Elsevier, Amsterdam, 1997).[7] J. Elf and M. Ehrenberg, Genome Res , 2475 (2003).[8] A. Melbinger, L. Reese, and E. Frey, Phys Rev Lett ,258104 (2012); J. Szavits-Nossan, K. Eden, R. J. Morris,C. E. MacPhee, M. R. Evans, and R. J. Allen, ,098101 (2014).[9] G. Rozhnova and A. Nunes, Phys Rev E , 041922(2009).[10] M. Aoki, Modeling aggregate behavior and fluctuationsin economics (Cambridge University Press, Cambridge,2001).[11] T. Heskes, J Phys A: Math Gen , 5145 (1994).[12] R. Grima, Phys Rev Lett , 218103 (2009);P. Thomas, A. V. Straube, and R. Grima, J Chem Phys , 195101 (2010).[13] R. Grima, P. Thomas, and A. V. Straube, J Chem Phys , 084103 (2011).[14] C. Cianci, F. Di Patti, D. Fanelli, and L. Barletti,Eur Phys J Special Topics , 5 (2012); P. Thomas,C. Fleck, R. Grima, and N. Popovic, J Phys A: MathTheor , 455007 (2014).[15] S. Engblom, Appl Math Comput , 498 (2006);R. Grima, J Chem Phys , 154105 (2012); A. Ale,P. Kirk, and M. P. Stumpf, , 174101 (2013).[16] V. Sotiropoulos and Y. N. Kaznessis, Chem Eng Sci ,268 (2011); P. Smadbeck and Y. N. Kaznessis, Proc NatlAcad Sci , 14261 (2013); A. Andreychenko, L. Mi-keev, and V. Wolf, ACM Trans Model Comput Simul , 12 (2015).[17] C. Cianci, F. Di Patti, and D. Fanelli, Europhys Lett , 50011 (2011).[18] T. M. Cover and J. A. Thomas, Elements of informa-tion theory , 2nd ed. (John Wiley & Sons, Hoboken, NewYersey, 2012).[19] P. Thomas, A. V. Straube, and R. Grima, J Chem Phys , 181103 (2011); K. R. Sanft, D. T. Gillespie, andL. R. Petzold, IET Syst Biol , 58 (2011).[20] J. Paulsson, O. G. Berg, and M. Ehrenberg, Proc NatlAcad Sci , 7148 (2000).[21] L. Comtet, Advanced Combinatorics (Springer, Nether-lands, 1974); The partial Bell polynomials are defined as B n,k ( { x χ } ) = (cid:80) (cid:48) n ! j ! ..j n − k +1 ! (cid:0) x (cid:1) j .. (cid:16) x n − k +1 ( n − k +1)! (cid:17) j n − k +1 , where the summation (cid:80) (cid:48) is such that j + . . . + j n − k +1 = k and j + 2 j + . . . + ( n − k + 1) j n − k +1 = n and { x χ } denotes the sequence ( x , x , x , . . . ). These are availablein Mathematica via the function BellY.[22] P. Thomas, H. Matuschek, and R. Grima, Computationof biochemical pathway fluctuations beyond the linearnoise approximation using iNA. in Bioinformatics andBiomedicine (BIBM), 2012 IEEE International Confer-ence on (IEEE, 2012) pp. 1–5.[23] V. Shahrezaei and P. S. Swain, Proc Natl Acad Sci ,17256 (2008).[24] The polylogarithm is defined by Li s ( x ) = (cid:80) ∞ k =1 x k k s andimplemeted in Mathematica by PolyLog[ x, s ].[25] D. T. Gillespie, J Phys Chem , 2340 (1977); Ex-act realizations of the process described by Eq. (55) areobtained using the Gillespie’s algorithm with the incre-ments of the synthesis reaction being replaced by inde-pendently and identically geometrically distributed inte-gers of mean b .[26] F. Verhulst, Methods and applications of singular pertur-bations (Springer, New York, 2006).[27] B. Munsky and M. Khammash, J Chem Phys ,044104 (2006).[28] N. N. Lebedev,