Analytical derivation of moment equations in stochastic chemical kinetics

Chemical Engineering Science 66 (2011) 268–277

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Analytical derivation of moment equations in stochastic chemical kinetics Vassilios Sotiropoulos, Yiannis N. Kaznessis  Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Ave. SE, Minneapolis, MN 55455, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 27 May 2010 Received in revised form 5 October 2010 Accepted 9 October 2010 Available online 26 October 2010

The master probability equation captures the dynamic behavior of a variety of stochastic phenomena that can be modeled as Markov processes. Analytical solutions to the master equation are hard to come by though because they require the enumeration of all possible states and the determination of the transition probabilities between any two states. These two tasks quickly become intractable for all but the simplest of systems. Instead of determining how the probability distribution changes in time, we can express the master probability distribution as a function of its moments, and, we can then write transient equations for the probability distribution moments. In 1949, Moyal defined the derivative, or jump, moments of the master probability distribution. These are measures of the rate of change in the probability distribution moment values, i.e. what the impact is of any given transition between states on the moment values. In this paper we present a general scheme for deriving analytical moment equations for any N-dimensional Markov process as a function of the jump moments. Importantly, we propose a scheme to derive analytical expressions for the jump moments for any N-dimensional Markov process. To better illustrate the concepts, we focus on stochastic chemical kinetics models for which we derive analytical relations for jump moments of arbitrary order. Chemical kinetics models are widely used to capture the dynamic behavior of biological systems. The elements in the jump moment expressions are a function of the stoichiometric matrix and the reaction propensities, i.e. the probabilistic reaction rates. We use two toy examples, a linear and a non-linear set of reactions, to demonstrate the applicability and limitations of the scheme. Finally, we provide an estimate on the minimum number of moments necessary to obtain statistical significant data that would uniquely determine the dynamics of the underlying stochastic chemical kinetic system. The first two moments only provide limited information, especially when complex, non-linear dynamics are involved. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Multiscale models Master equation Model reduction Statistical thermodynamics Mathematical modelling Computational chemistry

1. Introduction Many stochastic processes in physics and chemistry can be considered to follow the Markov property; the movement of the system in the available state space depends only on the previous state and not on the path that the stochastic process has followed before. Any stochastic process that obeys the Markov property is a Markov process. The underlying probability distribution, that is the probability of finding the system at any possible state at a certain time, is governed by a single partial differential equation (PDE) called the master equation (ME). If an initial condition is known and given the transitional probabilities, i.e. the probabilities of the system transitioning from any state to any other state, ME uniquely determines the probability distribution at any later time (van Kampen, 1992; Gardiner, 2004).

 Corresponding author.

E-mail addresses: [email protected] (V. Sotiropoulos), [email protected] (Y.N. Kaznessis). 0009-2509/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2010.10.024

The solution of ME requires enumerating the states and finding the transition matrix, a task that is tractable only for the simplest of systems. For more complex systems, in lieu of a solution of ME, kinetic Monte Carlo techniques are often used to sample the underlying probability distribution. There is a growing community of researchers that develops computationally efficient and accurate stochastic simulation algorithms (Salis and Kaznessis, 2005a, 2005b; Gibson and Bruck, 2000; Cao et al., 2006; Rathinam et al., 2003; Cao et al., 2005; Samant and Vlachos, 2005; Munsky and Khammash, 2006; Sotiropoulos and Kaznessis, 2008; Sotiropoulos et al., 2009; Salis et al., 2006; Hill et al., 2008). Alternatively, instead of sampling the probability distribution all important information for the system’s behavior can be had through the moments of the probability distribution (Gillespie, 1992). The first moment relates to the mean of the probability distribution, the second to the variance, the third to skewness and the fourth to kurtosis. Higher order moments may also contain important information on systems dynamics especially when these systems are complex and highly non-linear. In this work we start from the master equation and derive a system of ordinary differential equations (ODE) that describe the

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

dynamics of the probability distribution moments. We express the transient moment dynamics in terms of the derivative or jump moments (Moyal, 1949). Jump moments are measures of the rate of change in the probability distribution moment values, i.e. they quantify the impact of any given transition between states on the moment values. What is new in the current work is a scheme to derive analytical expressions for the jump moments for any N-dimensional Markov process. These expressions result in an infinite linear system of ODEs that describes the moment dynamics of any such process. We then apply this scheme to stochastic chemical kinetics. Stochastic chemical kinetics have emerged in recent years as an appropriate modeling formalism for biological systems that are away from the thermodynamic limit (Kaznessis, 2006; Sotiropoulos and Kaznessis, 2007; Kaznessis, 2007; Tomshine and Kaznessis, 2006; Tuttle et al., 2005; Salis and Kaznessis, 2006; Ramalingam et al., 2009). The idea of using moment equations to predict dynamics of systems in the stochastic chemical kinetics regime has also been considered early on by McQuarrie (1967). In recent years, publications have proposed alternative ways to derive the moment equations (Hespanha and Singh, 2005; Goutsias, 2007; Gillespie, 2009; Lee et al., 2009). Our work complements this rich literature, providing useful analytical relations for the jump moments in stochastic chemical kinetics. These relations are then used to provide analytical equations for the probability distribution moments. The paper is organized as follows. First we briefly present background information on Markov processes and their governing equations. We then define the jump moments of Markov processes and we use this definition to derive the general form of the moment equations. Then we concentrate on the stochastic chemical kinetics regime and we derive analytical relations for the jump moments. We use examples to discuss the applicability of the derived equations. Finally we discuss the drawbacks of any moment based method that mainly stem from the fact that the resulting ODE system is infinite dimensional. We discuss literature momentclosure schemes and speculate on how the analytical expressions we develop may assist in the development of accurate closing schemes. We also argue that higher order moments, at least up to order six, are necessary for biologically-relevant chemical kinetics systems. This is contrary to the widely held belief that only the mean and variance are important in stochastic chemical kinetics models.

The probability distribution function of any stochastic process with state vector X ðtÞ ¼ ðX1 ðtÞ, . . . ,XN ðtÞÞ and initial condition X ðt ¼ 0Þ ¼ X 0 obeying the Markov property is governed by a differential difference equation in probability known as the master equation (ME) (van Kampen, 1992) Z @PðX ,tÞ ¼ ½TðX =X uÞPðX u,tÞTðX u=X ÞPðX ,tÞ dX u ð2Þ @t PðX ,tÞ is the probability of the system being at state X at time t. TðX =X uÞ is the transition probability per unit time for the system to jump from state X u to state X . In principle, the ME uniquely determines the probability PðX ,tÞ of the system being at a state X ¼ X ðtÞ at time t 40. Note that using integrals in the ME we have made the assumption that variable X ðtÞ is continuous. An alternative description of the Markov process probability distribution is obtained through the Kramers–Moyal expansion (van Kampen, 1992). 1 N X @PðX ,tÞ ð1Þm X @m ¼ ½a ðX ÞPðX ,tÞ @t m! j ,...,j ¼ 1 @Xj1    @Xjm m m¼1 1

2.1. Markov processes Consider a Markov process, X(t), a stochastic process that for any n successive set of times, t1 o t2 o    otn , obeys the following property, known as the Markov property P1jn1 ðXn ,tn jX1 ,t1 , . . . ,Xn1 ,tn1 Þ ¼ P1j1 ðXn ,tn jXn1 ,tn1 Þ

ð1Þ

where Xk is the state of the system at time tk. P1jn1 denotes the conditional probability that the outcome at time tn is Xn provided that the outcomes were X1, X2,y,Xn  1 at times t1, t2,y,tn  1, respectively. P1j1 denotes the conditional probability that the outcome at time tn is Xn provided that the outcome at time tn  1 was Xn  1, irrespective of the previous outcomes. The Markov property states that the transition of the system to state Xn at time tn depends only on the previous state Xn  1 at time tn  1 and therefore is independent of the history of the process, i.e. X1, t1, y, Xn  2, tn  2.

ð3Þ

m

where a m are m-order tensors and were originally called derivative or jump moments by Moyal (1949) and later referred to as propagator moment functions (Gillespie, 1992). Jump moments are measures of the rate of change in the probability distribution moment values, i.e. they quantify the impact of any given transition between states on the moment values.

2.2. Moments and jump moments Let us consider an N-dimensional Markov process X ðtÞ ¼ ðX1 ðtÞ, . . . ,XN ðtÞÞ with probability distribution PðX ðtÞ,tÞ. The m-th moment of the stochastic variable X is defined as Z /Xim S ¼ Xim PðX ,tÞ dX , i ¼ 1, . . . ,N ð4Þ The first order moments of any Xi is usually referred to as the mean while second order moments relate to the variance of the probability distribution through the relation, varðXi Þ ¼ /Xi2 S/Xi S2 The joint moments are defined through the following relations Z /X1m1 X2m2    Xnmn S ¼ ½X1m1 X2m2    Xnmn PðX ,tÞ dX m1 ,m2 , . . . mn ¼ 1,2,3, . . .

2. Theory

269

ð5Þ

where the sum m ¼ m1 þm2 þ    mn is the order of the joint moment. The set of the m-th order moments uniquely determines the underlying probability distribution. The existence of higher order moments depends on how fast PðX ,tÞ approaches zero as X goes to infinity. On the other hand, jump moments relate to the transition probability per unit time rather than the probability distribution. Based on Moyal’s work jump moments are defined through the following relations (Moyal, 1949) Z aim ðX Þ ¼ ðXi uXi Þm TðX u=X Þ dX u, i ¼ 1, . . . ,N ð6Þ Index m refers to the order of the moment and index i to the element within the tensor. Analogously the joint jump moments are Z ðX Þ ¼ ½ðXi uXi Þmi ðXj uXj Þmj   TðX u=X Þ dX u ð7Þ ai,j mi þ mj þ  Jump moments are tensors and the indexes i,j, . . . denote the element of the m ¼ mi þ mj þ    order tensor. Jump moments of

270

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

2,1 order 2 and greater are symmetric tensors, i.e. a1,2 2 ¼ a2 . The same is true for regular moments. In the simple case of an onedimensional Markov process jump moments are simply scalars.

where by identity a0 ðX Þ ¼ 0, ðmjÞ denotes the binomial coefficient and index notation is used to symbolize the following operation index ¼ f1, . . . ,1 ,2, . . . ,2 , . . . ,N, . . . ,N g |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflffl{zfflfflfflffl} |fflfflfflffl{zfflfflfflffl} m1

2.3. Derivation of moment equations

m2

j1

Given the definitions for both regular and jump moments we can now start deriving the moment equations, i.e. the system of ODEs that describes moment dynamics. Starting with the first moment we have Z Xi PðX ,tÞ dX ð8Þ /Xi S ¼ From the defining equation we have the exact identity for the time derivative of the first moment Z @PðX ,tÞ d/Xi S dX ¼ Xi ð9Þ @t dt substituting Eq. (2) in the last equation we have Z Z d/Xi S ¼ Xi ½TðX =X uÞPðX u,tÞTðX u=X ÞPðX ,tÞ dX u dX dt ZZ ¼ ½Xi TðX =X uÞPðX u,tÞXi TðX u=X ÞPðX ,tÞ dX u dX

thus we have ZZ d/Xi S ¼ ðXi uXi ÞTðX u=X ÞPðX ,tÞ dX u dX dt

ð10Þ

ð11Þ

ð12Þ

the right hand side (RHS) of Eq. (12) can be recast in the following form Z d/Xi S ¼ ai1 ðX ÞPðX ,tÞ dX ð14Þ dt where the RHS term denotes the average of ai1 ðX Þ. The ODE equation describing the dynamics of the first moment of variable Xi is then simply

m1 j1

m2 j2

m 1 X d/X ðmÞ S ¼ vecðam Þ þ P m /X ðjÞ  vecðamj ÞS j dt j¼1

/X1m1 j1 X2m2 j2    XNmN jN aindex j1 þ j2 þ jN ðX ÞS ð16Þ

mN ðt ¼ 0ÞS /X1m1 X2m2    XNmN St ¼ 0 ¼ /X1m1 ðt ¼ 0ÞX2m2 ðt ¼ 0Þ    XN0

ð19Þ

with initial conditions /X ðmÞ St ¼ 0 ¼ /X ðmÞ ðt ¼ 0ÞS where X ðmÞ ¼ X  X     X |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} m

and where vecðÞ symbolizes the vectorization of any matrix or tensor in general into a column vector, which is obtained by stacking the columns of the matrix or tensor on top of one another. Finally P m is the sum of ðmjÞ ¼ m!=ðj!ðmjÞ!Þ permutation matrices j

(Nm  Nm), defined as follows: m

Pm ¼

ðjÞ X k¼1

p

m,k

where each permutation matrix, p

j1 ,j2 ,...,jN ¼ 0

ð18Þ

Alternatively, the multivariate case can be written in efficient matrix form using Kronecker product notation (Brewer, 1978). Kronecker products, ðÞ, have interesting implementation properties on computers and thus may improve computational efficiency. We can write

ð15Þ

Following a similar way of thinking and carrying out slightly more complicated algebraic calculations, which we omit in the present section for brevity but we include in Appendix A for completeness, we can write the moment equations for second, third and fourth order moment dynamics. Using inductive reasoning, from Eqs. (15), (A.5), (A.14) and (A.23) we derive the general moment equation for the mth order moment of an N-dimensional Markov process as a function of joint jump moments ! ! ! m1 ,m 2 ,...,mN X m2 m1 mN d/X1m1 X2m2    XNmN S ¼  j1 j2 jN dt

with initial conditions

ð17Þ

mN jN

In other words, the index notation refers to the appropriate element of the jth order tensor, j ¼ j1 þ j2 þ    jN . To our knowledge there have been no other attempts to deduce the general moment equation from jump moments in the case of a N-dimensional Markov process. In the case of a 1-dimensional Markov process the above equation reduces to the following simpler relation (Gillespie, 1992) ! m X m d/X m S ¼ /X mj aj ðX ÞS j dt

j

i ¼ 1, . . . ,N

jN

/X m St ¼ 0 ¼ /X m ðt ¼ 0ÞS

Recalling the definition of the jump moment in Section 2.2 and in particular that of the first moment Z ai1 ðX Þ ¼ ðXi uXi ÞTðX u=X Þ dX u ð13Þ

d/Xi S ¼ /ai1 ðX ÞS, dt

j2

¼ f1, . . . ,1 ,2, . . . ,2 , . . . ,N, . . . ,N g |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflffl{zfflfflfflffl} |fflfflfflffl{zfflfflfflffl}

j¼1

Noticing that the integration over X and X u runs over the same domain we can interchange indexes in the last equation, i.e. Xi TðX =X uÞPðX u,tÞ ¼ Xi uTðX u=X ÞPðX ,tÞ

mN

f1, . . . ,1 ,2, . . . ,2 , . . . ,N, . . . ,N g |fflfflfflffl{zfflfflfflffl} |fflfflfflffl{zfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl}

p

m,k

¼

m,k

, is defined as (Brewer, 1978)

N X

NN NN ENN i1 i1  Ei2 i2     Eim im i1 ,...,im ¼ 1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} m

where Epq is a p  q matrix (elementary matrix) that has a ‘‘1’’ in ij the (i,j) position and zeros everywhere else and k denotes the different permutation sets for the subscripts of all Epq ij . The methodology to choose the appropriate permutation sets is m m detailed in Appendix B. Note that Pm 0 ¼ IN m N m , i.e. the N  N identity matrix. The system of Eqs. (16) or its equivalent (19), completely characterizes the moment dynamics for any given N-dimensional Markov process. The only prerequisite is the knowledge of analytical relations for the jump moments. These relations depend on the underlying physics of the problem as it will become evident in the following section. In particular, we derive analytical relations for the components of the jump moments tensors in the case where the

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

underlying Markov process stems from a stochastic chemical kinetics model. The final equations may appear cumbersome but they are very simple to generate with the help of a computer program. Even though we have made the assumption that variables X ðtÞ are continuous, this needs not be the case. The derivation of Eq. (16) could be carried out by starting from the discrete ME (van Kampen, 1992). The main difference in the derivation is that instead of integrals, summation signs over all possible states would be used. 2.4. Jump moments and stochastic chemical kinetics Consider a system of N distinct chemical species, Xi (i¼1,y,N), participating in M chemical reactions in a well-mixed volume V N X i¼1

kj

rij Xi !

N X

pji Xi ,

j ¼ 1, . . . ,M

ð20Þ

i¼1

Consider also n j the stoichiometric vector associated with the jth reaction 2 j 3 p r j 6 1 17 7 ð21Þ nj ¼ 6 4 ^ 5, j ¼ 1, . . . ,M pjN rNj Such systems of reactions are widely used to model biological interactions, such as transcription, translation, degradation, regulation and protein–protein interactions (Wolf and Arkin, 2002; Kaznessis, 2007; St-Pierre and Endy, 2008; Tuttle et al., 2005; Sotiropoulos and Kaznessis, 2007). Dilute and sparse species populations render the traditional continuous-deterministic modeling approach false. Instead, stochastic chemical kinetics models are more appropriate to describe systems of reactions that are far from the thermodynamic limit (Elowitz et al., 2002). Therefore system (20) is frequently considered as a Markov process where a chemical master equation (CME) governs the evolution of the probability distribution. Under the stochastic chemical kinetics regime reaction rates become reaction propensities, aj ðX Þ. These are the probabilistic equivalents and are defined as follows (Gillespie, 1977)

aj ðX Þ ¼ kj cj ðX Þ, cj ðX Þ ¼

N Y

Xi ! j j i ¼ 1 ri !ðXi ri Þ!

ð22Þ

Each constant kj represents the mesoscopic reaction rate of the jth reaction. In what follows we derive analytical relations for the jump moments for any stochastic chemical kinetics model. Starting with the CME, we can determine the Kramers–Moyal expansion of the probability distribution (cf. Eq. (3)). If we truncate the expansion and retain only the first two terms then the resulting equation is the well-known Fokker–Plank equation (FPE) (van Kampen, 1992) i @PðX ,tÞ @ 1 @ @ h ¼ ½a ðX ÞPðX ,tÞ þ : a ðX ÞPðX ,tÞ ð23Þ 2 @X 1 2 @X @X @t where a and a 1 are the first two jump moments tensors or using 2 the terminology most often referred to them the drift vector and diffusion tensor, respectively. Instead of solving the Fokker–Planck equations it is usually more convenient to sample the underlying probability distribution generating ensembles of trajectories obtained as solutions of the corresponding chemical Langevin equations (CLE) or systems of CLEs (Gillespie, 2000). From the work of Gillespie we know that for systems of chemical reactions under the Markov process regime the corresponding chemical Langevin equation is (Gillespie, 2000) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dX ¼ n a ðX Þ dt þ n D ð a ðX ÞÞ dW ð24Þ

271

where n corresponds to the N  M stoichiometric matrix, a ðX Þ corresponds to the M  1 propensities vector and the notation pffiffiffi D ð F Þ denotes the diagonal matrix whose (i,i)th element are the only nonzero elements and their value equals the square root of the ith component of vector F . Given any FPE we can deduce the corresponding system of CLEs and vice versa. This direct relation between the two allows us to infer through inductive reasoning analytical relations for the jump moment tensors elements. From Eq. (24) we obtain the following representation for the FPE @PðX ,tÞ @ 1 @ @ ¼ ½n a ðX ÞPðX ,tÞ þ : ½n D ða ðX ÞÞn T PðX ,tÞ @t @X 2 @X @X

ð25Þ

Comparing Eqs. (23) and (25) we retain the following relations for the first two jump moment tensors a 1 ðX Þ ¼ n a ðX Þ a ðX Þ ¼ n D ða ðX ÞÞn T

ð26Þ

2

or using summation notation each of the tensor elements is defined as follows: ai1 ðX Þ ¼

M X

nik ak ðX Þ

k¼1

aij2 ðX Þ ¼

M X

nik njk ak ðX Þ

ð27Þ

k¼1

We notice that the difference between the two jump moments relations relies on an additional component of the stoichiometric vector. Following the discussion in the previous section we infer the relations for the third and fourth jump moments, where a3 ðX Þ and a4 ðX Þ are N  N  N and N  N  N  N tensors, respectively. aijl 3 ðX Þ ¼

M X

nik njk nlk ak ðX Þ

k¼1

aijlm 4 ðX Þ ¼

M X

nik nlk njk nmk ak ðX Þ

ð28Þ

k¼1

In general, for stochastic chemical kinetics models we infer that given the stoichiometric matrix and the reaction propensities vector all jump moments can be defined analytically through the following recursive formula m indeces

zfflffl}|fflffl{ M X ij    l am ðX Þ ¼ nik njk    nlk ak ðX Þ |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} k¼1

ð29Þ

m nus

The formula is relatively simple and intuitive. Substituting the last equation into Eq. (16) returns analytical relations for the moment equations. A first obvious comment is that the linear or non-linear character of the underlying reaction networks manifests itself through the reaction propensities which in turn render jump moments linear or non-linear functions of the chemical species concentrations. The impact of the non-linear jump moment relations will become evident in the examples following. In general, a similar approach may lead to analytical relations for jump moments for any given Markov process. In the supplementary material, we present two toy examples, a linear and a non-linear reaction network, to highlight the analytical form of the moment equations for stochastic chemical kinetic systems.

272

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

3. Discussion 3.1. Infinite dimensional moment equations For any given non-linear model, whether it is a chemical kinetics model or an electrical circuit model, jump moments will be nonlinear functions of the state vector, hence the system of ODEs describing the moment evolution will be infinite dimensional. The general form of the ODEs describing the moment dynamics of a stochastic chemical kinetics model up to a desired order can be summarized through the following matrix equation dY ¼ AY þA Y h þB h dt

ð30Þ

where Y is a vector containing the moment elements of interest and A is a square matrix with appropriate dimensions. Its elements depend on the kinetic constants and the stoichiometric matrix of the reactions set. Matrix A may be sparse depending on the connectivity of the system. Vector B contains the terms that depend on any zeroth order reactions and are always constant as the corresponding reaction rates are independent from the state vector. Finally, the term A Y h denotes the dependence of lower h

order moments on higher order moments. Vector Y h contains higher order moments and matrix A depends on the reaction h

constants and stoichiometric matrix. A Y h vanishes only when the h

underlying system is linear. 3.2. Moment closure schemes As mentioned earlier Eq. (30) is not solvable when the system of reactions is non-linear, as the dependence on higher moments renders the system infinite. In this section we briefly discuss existing approaches, called moment closure techniques, that attempt to bypass the problem and obtain approximate solutions. The main idea behind any moment closure scheme is to approximate the infinite dimensional system of ODEs depicted in Eq. (30) with a finite system of the form dY t ¼ A Y t þ f h ðY t Þ þ B dt

ð31Þ

where Y t is a truncated version of the infinite dimensional vector Y (cf. Eq. (30)). The last equation explicitly states that the term A Y h h

in Eq. (30) should be substituted with a vector containing linear or non-linear functions of Y , f h ðY t Þ, that would closely approximate the effect of the higher order moment terms. In other words, the effect of higher order moment terms is approximated through functions of lower order moment terms. This approximation renders the system finite but also non-linear as all the developed approaches consider f to be non-linear functions. At the same time the approximation introduces an error in the moment values with its significance depending on how well f follows A Y h . h

Numerous attempts have been made to develop such closure schemes (Thakur et al., 1978; Gillespie, 2009; Krishnarajah et al., 2005; Hespanha (2006) and Singh, 2005). A summary of the existing moment closure techniques in stochastic chemical kinetics system can be found in Singh and Hespanha (2006). To our knowledge there is no satisfactory solution to the moment closure scheme with broad applicability and the topic is still being actively researched.

higher than the first two order moments are necessary to capture the dynamic behavior. ¨ model that under specific Consider the well-studied Schlogl values for the kinetic rates exhibits bistability, similar to the bistable switch (Gardner et al., 2000) and or the lphage infection (St-Pierre and Endy, 2008). Consider the network of reactions shown in Table 1. Kinetic data and initial conditions are taken from Gunawan et al. (2005). The kinetic values depicted in Table 1 render the system ¨ model enables us to determine bistable. The use of the Schlogl the minimum number of moments needed to reconstruct the underlying probability distribution within a reasonable error. Using Hy3S (Salis et al., 2006) and SynBioSS (Hill et al., 2008), a suite of multiscale algorithms, we generate 105 independent trajectory trials in the interval [0,20] s. We should note that Hy3S is a publicly available software suite for simulating chemical reaction networks. Hy3S is available with Open Licenses at sourceforge.net. We use Matlab to compute the underlying probability distribution of species X at three different time points. The shapes and characteristics of the corresponding probability distributions are depicted in Fig. 1a. Next we compute the first eight moments of X at the three different time points using the data obtained from the stochastic simulation. Their values are shown in Table 2. Using the numerical values of the moments we can reconstruct the corresponding probabilities using an algorithm based on the maximum entropy principle (Mohammad-Djafari, 2001). In order to determine the minimum number of moments we first reconstruct the distributions using only the first two moments and then we keep increasing the number of moments by two and up to order eight. In Figs. 1b–d the comparison between the actual and reconstructed probabilities are depicted. Reconstructed probabilities using only the first two moments expectedly fail to capture the bimodality of the distribution especially in the cases of t ¼4 and 20 s. Using four moments seems to capture the essential futures of the distributions, i.e. the two modes, but it fails to weight in the relative peaks of the two modes. Six moments produce adequate results, especially when the separation between the two modes is distinct (cf. Fig. 1d), but results are inferior when there is no distinct separation (cf. Figs. 1b and c). Still results are relatively accurate. Finally, using the first eight moments produces distributions that are almost identical to the actual distribution with minor, if any, deviations. Overall, it is evident that using the first two moments is not adequate to reconstruct the probability distribution. Our analysis shows that use of at least the first six, or in some cases eight, moments is needed for adequate reconstruction of the probability distribution.

Table 1 ¨ model. Reactions and parameters for the Schlogl Set of reactionsa k1

A þ 2X!3X k2

3X!2X þ A k3

B!X k4

X!B

3.3. Importance of higher moments in non-linear chemical kinetics

a

Initial valuesc

k1 ¼3  10  7

[X]0 ¼247

k2 ¼10

4

2

 10

k3 ¼ 10  3

[A]¼ 105 2

k4 ¼3.5  10

[B] ¼2  105

A and B are buffer species. For 1st order reactions the units are s  1, for 2nd order reactions the units are molecules  1 s  1 and for 3rd order reactions the units are molecules  2 s  1. c Number of molecules. b

In this section we determine the probability distribution of a simple non-linear system of chemical reaction to illustrate how

Mesoscopic reaction ratesb

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

5

273

x 10−3

Probability Distribution

Probability Distribution

Actual

0.02 0.015 0.01 0.005 0 0 200 400 Numb er of M ole

2 600 800 cules of X

4 ) e (s Tim

20

4 Moments 6 Moments 8 Moments

2

0 0

200

400

600

800

Number of Molecules of X

0.012

0.015 Actual 2 Moments 4 Moments 6 Moments 8 Moments

0.008

Probability Distribution

Probability Distribution

2 Moments

4

0.004

0

Actual 2 Moments 4 Moments 6 Moments 8 Moments

0.01

0.005

0 0

200

400

600

800

Number of Molecules of X

0

200

400

600

800

Number of Molecules of X

¨ model, depicted in Table 1, using the moment data in Table 2. (a) Probability distribution of species X Fig. 1. Actual and reconstructed probability distributions for the Schlogl at times t¼ 2, 4 and 20 s. (b) Comparison between the actual and reconstructed probability distributions using different moment sets at time t ¼2 s. (c) Comparison between the actual and reconstructed probability distributions using different moment sets at time t ¼ 4 s. (d) Comparison between the actual and reconstructed probability distributions using different moment sets at time t ¼ 20 s.

Table 2 ¨ model. First eight moment values for X in the Schlogl Moments

Value at t ¼2 s

Value at t ¼4 s

Value at t ¼20 s

/XS /X 2 S /X 3 S /X 4 S /X 5 S /X 6 S /X 7 S /X 8 S

268.9871 9.2903  104 3.7966  107 1.7179  1010 8.2743  1012 4.1528  1015 2.1465  1018 1.1347  1021

295.6759 1.3364  105 7.0001  107 3.8403  1010 2.1488  1013 1.2170  1016 6.9578  1018 4.0100  1021

301.6574 1.4733  105 8.1413  107 4.6203  1010 2.6470  1013 1.5262  1016 8.8511  1018 5.1617  1021

We should stress that there is no general approach to choose the number of moments necessary to accurately capture the probability distribution. This choice remains system dependent and is one driven by empiricism.

Acknowledgements This work was supported by a Grant from the National Institutes of Health (American Recovery and Reinvestment Act Grant GM086865), and Grants from the National Science Foundation (CBET-0425882 and CBET-0644792). Computational support from the Minnesota Supercomputing Institute (MSI) is gratefully acknowledged. This work was also supported by the National Computational Science Alliance under TG-MCA04N033.

4. Conclusions A new derivation of the moment equations for any N-dimensional Markov process using the definition of jump moments was presented. The applicability of the scheme is general and leads to analytical relations given that the functional form of the jump moments is known. Focusing on stochastic chemical kinetics models, we derived analytical relations for the elements of any jump moment tensor, demonstrating the ease of equation setup. The elements of the analytical relations are a function of the stoichiometric matrix and the reaction propensities, i.e. the probabilistic reaction rates. Using a simple example of non-linear kinetics we illustrated how higher than second order moments are important to accurately reconstruct the probability distribution for biologically-relevant chemical kinetics systems. Certainly, then, more effort is warranted on the development of moment closure schemes.

Appendix A. Derivation of second through fourth order moment equations In this section of the Appendix we derive the second through fourth order moment equations. For the second order moment equations, the starting point is the defining Eq. (5) from which the following identity is derived for the time derivative Z @PðX ,tÞ d/Xi Xj S dX , i,j ¼ 1, . . . ,N ¼ Xi Xj ðA:1Þ @t dt Analogously to the derivation of the first moment equations, we substitute Eq. (2) in the last equation and by interchanging notation using the same arguments as previously (cf. Eq. (11)). Then ZZ d/Xi Xj S ¼ ðXi uXj uXi Xj ÞTðX u=X ÞPðX ,tÞ dX u dX ðA:2Þ dt

274

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

ZZ

Through algebraic manipulations we can rearrange the terms in the last equation. In particular we use the following identity

þ

ðXi uXj uXi Xj Þ ¼ ðXi uXi ÞðXj uXj Þ þ Xj ðXi uXi Þ þ Xi ðXj uXj Þ

þ

Substituting the term in the RHS of Eq. (A.2) yields ZZ d/Xi Xj S ¼ ðXi uXi ÞðXj uXj ÞTðX u=X ÞPðX ,tÞ dX u dX dt ZZ þ Xj ðXi uXi ÞTðX u=X ÞPðX ,tÞ dX u dX ZZ þ Xi ðXj uXj ÞTðX u=X ÞPðX ,tÞ dX u dX

ZZ

ðA:3Þ

ZZ þ

Xi Xj ðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX Xi Xl ðXj uXj ÞTðX u=X ÞPðX ,tÞ dX u dX Xj Xl ðXi uXi ÞTðX u=X ÞPðX ,tÞ dX u dX

Using the definitions of the joint jump moments (cf. Eq. (7)) and that of the averages the RHS eventually becomes ðA:4Þ

d/Xi Xj Xl S jl ij li ¼ /aijl 3 ðX ÞS þ /Xi a2 ðX ÞS þ /Xj a2 ðX ÞS þ /Xl a2 ðX ÞS dt þ /Xi Xj al1 ðX ÞS þ/Xj Xl ai1 ðX ÞS þ /Xl Xi aj1 ðX ÞS, i,j,l ¼ 1, . . . ,N

Invoking the definition of the joint jump moments (cf. Eq. (7)) and also that of averages the RHS becomes d/Xi Xj S ¼ /aij2 ðX ÞS þ /Xi aj1 ðX ÞS þ /Xj ai1 ðX ÞS, dt

i,j ¼ 1, . . . ,N

ðA:5Þ

where a ðX Þ is a second order tensor, whereas a 1 ðX Þ is a first 2

order tensor. In the trivial case where i¼j the last equation simplifies to d/Xi Xi S ¼ /aii2 ðX ÞSþ 2/Xi ai1 ðX ÞS, dt

i ¼ 1, . . . ,N

ðA:6Þ

In the derivation of the third moment equations we start again with the defining moment equation (cf. Section 2.2) Z Xi Xj Xl PðX ,tÞ dX ðA:7Þ /Xi Xj Xl S ¼ and the time derivative is defined as follows: Z @PðX ,tÞ d/Xj Xi Xl S dX ¼ Xi Xj Xl @t dt

ðA:13Þ

ðA:8Þ

ðA:14Þ

where a ðX Þ is a third order tensor. In the trivial case where i ¼j¼l 3 the last equation simplifies to d/Xi3 S ii 2 i ¼ /aiii 3 ðX ÞSþ 3/Xi a2 ðX ÞS þ 3/Xi a1 ðX ÞS dt

i ¼ 1, . . . ,N ðA:15Þ

For the fourth order moments the starting point is again the moment defining equation (cf. Section 2.2) Z Xi Xj Xl Xm PðX ,tÞ dX ðA:16Þ /Xi Xj Xl Xm S ¼ and the time derivative is defined as follows: d/Xj Xi Xl Xm S ¼ dt

Z

Xi Xj Xl Xm

@PðX ,tÞ dX @t

ðA:17Þ

substituting Eq. (2) in the last equation we have

substituting Eq. (2) in the last equation we have Z Z d/Xj Xi Xl S ¼ Xj Xi Xl ½TðX =X uÞPðX u,tÞTðX u=X ÞPðX ,tÞ dX u dX dt ZZ ¼ ½Xj Xi Xl TðX =X uÞPðX u,tÞXj Xi Xl TðX u=X ÞPðX ,tÞ dX u dX

d/Xj Xi Xl Xm S ¼ dt

Z

Xj Xi Xl Xm

ZZ

¼

Z

½TðX =X uÞPðX u,tÞTðX u=X ÞPðX ,tÞ dX u dX

½Xj Xi Xl Xm TðX =X uÞPðX u,tÞXj Xi Xl Xm TðX u=X ÞPðX ,tÞ dX u dX

ðA:18Þ ðA:9Þ Noticing that the integration over X and X u runs over the same domain we can interchange indexes in the last equation, i.e. Xj Xi Xl TðX =X uÞPðX u,tÞ ¼ Xj uXi uXl uTðX u=X ÞPðX ,tÞ thus we have ZZ d/Xj Xi Xl S ¼ ðXj uXi uXl uXj Xi Xl ÞTðX u=X ÞPðX ,tÞ dX u dX dt

ðA:10Þ Xj Xi Xl Xm TðX =X uÞPðX u,tÞ ¼ Xj uXi uXl uXm uTðX u=X ÞPðX ,tÞ ðA:11Þ

Skipping through some tedious calculations we use the following relation to transform the RHS of Eq. (A.11) ðXj uXi uXl uXi Xj Xl Þ ¼ ðXi uXi ÞðXj uXj ÞðXl uXl Þ þXi ðXj uXj ÞðXl uXl Þ þ Xj ðXi uXi ÞðXl uXl Þ þXl ðXi uXi ÞðXj uXj Þ þ Xi Xj ðXl uXl Þ þXi Xl ðXj uXj Þ þ Xj Xl ðXi uXi Þ substituting yields ZZ d/Xi Xj Xl S ¼ ðXi uXi ÞðXj uXj ÞðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX dt ZZ þ Xi ðXj uXj ÞðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX ZZ þ Xj ðXi uXi ÞðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX ZZ þ Xl ðXi uXi ÞðXj uXj ÞTðX u=X ÞPðX ,tÞ dX u dX

Similarly to the third order moments case the integration in Eq. (A.18) over X and X u runs over the same domain thus we can interchange indexes, i.e. ðA:19Þ

thus we have d/Xj Xi Xl S ¼ dt

ZZ

ðXj uXi uXl uXm uXj Xi Xl Xm ÞTðX u=X ÞPðX ,tÞ dX u dX

ðA:20Þ

Skipping through some tedious calculations we use the following relation to transform the RHS of Eq. (A.11) ðA:12Þ

ðXj uXi uXl uXm uXi Xj Xl Xm Þ ¼ ðXi uXi ÞðXj uXj ÞðXl uXl ÞðXm uXm Þ þ Xi ðXj uXj ÞðXl uXl ÞðXm uXm Þ þ Xj ðXi uXi ÞðXl uXl ÞðXm uXm Þ þ Xl ðXi uXi ÞðXj uXj ÞðXm uXm Þ þ Xm ðXi uXi ÞðXj uXj ÞðXl uXl Þ þ Xi Xj ðXl uXl ÞðXm uXm Þ þ Xi Xl ðXj uXj ÞðXm uXm Þ þ Xi Xm ðXj uXj ÞðXl uXl Þ þ Xj Xl ðXi uXi ÞðXm uXm Þ

þXj Xm ðXi uXi ÞðXl uXl Þ þXl Xm ðXi uXi ÞðXj uXj Þ þXi Xj Xl ðXm uXm Þ þ Xi Xj Xm ðXl uXl Þ þXi Xl Xm ðXj uXj Þ þXj Xl Xm ðXi uXi Þ

ðA:21Þ

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

The number of the necessary permutation matrices to construct

substituting yields

P m depends on both m and j, as the binomial coefficient denotes in

ZZ

d/Xi Xj Xl S ¼ dt

275

j

ðXi uXi ÞðXj uXj ÞðXl uXl ÞðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX

Eq. (B.1). The choice of the appropriate ðmjÞ permutation matrices

ZZ

in turn depends on the terms that P m multiplies, i.e. /X ðjÞ 

þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ ZZ þ

Xi ðXj uXj ÞðXl uXl ÞðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX Xj ðXi uXi ÞðXl uXl ÞðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX

j

vecðamj ÞS. These terms are column vectors ðNm  1Þ and each of i

,...,im

jþ1 the Nm elements can be depicted as /Xi1 Xi2    Xij amj

2

Xl ðXi uXi ÞðXj uXj ÞðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX

S.

3

^

6 7 ij þ 1 ,...,im S7 /X ðjÞ  vecðamj ÞS ¼ 6 4 /Xi1 Xi2    Xij amj 5

Xm ðXi uXi ÞðXj uXj ÞðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX

ðB:3Þ

^

Xi Xj ðXl uXl ÞðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX Xi Xl ðXj uXj ÞðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX

The idea behind choosing the appropriate permutation matrices, p , is simple and straight forward. Through Eq. (B.2), m,k

Xi Xm ðXj uXj ÞðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX

the k permutation matrices are identified by keeping the first subscript of the elementary matrices constant and changing the second subscript according to the order in which subscript indices

Xj Xl ðXi uXi ÞðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX

i

Xj Xm ðXi uXi ÞðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX

,...,im

jþ1 appear in the /Xi1 ,i2 ,...,ij amj

S terms. Note that we are only

interested in unique subscript permutations sets. The methodology presented above will become evident in the following examples. We consider a system with N ¼2 and we construct the P m matrices for the cases m¼2 through m ¼4.

Xl Xm ðXi uXi ÞðXj uXj ÞTðX u=X ÞPðX ,tÞ dX u dX Xi Xj Xl ðXm uXm ÞTðX u=X ÞPðX ,tÞ dX u dX

j

Xi Xj Xm ðXl uXl ÞTðX u=X ÞPðX ,tÞ dX u dX

 m¼2 According to Eq. (19) we have

Xi Xl Xm ðXj uXj ÞTðX u=X ÞPðX ,tÞ dX u dX Xj Xl Xm ðXi uXi ÞTðX u=X ÞPðX ,tÞ dX u dX

ðA:22Þ

Using the definitions of the joint jump moments (cf. Eq. (7)) and that of the averages the RHS eventually becomes d/Xi Xj Xl Xm S jlm lmi ¼ /aijlm 4 ðX ÞS þ/Xi a3 ðX ÞSþ /Xj a3 ðX ÞS dt ijl lm þ /Xl amij 3 ðX ÞS þ /Xm a3 ðX ÞS þ/Xi Xj a2 ðX ÞS

d/X ð2Þ S ¼ vecða2 Þ þ P 2 /X  vecða1 ÞS 1 dt The column vector /X  vecða1 ÞS has entries of the form /Xi1 ai12 S or /Xi2 ai11 S. Note that these two permutation of indices i1 and i2 are the only ones producing unique terms. Therefore the corresponding permutation matrices, p , are m,k

ij jl þ /Xj Xl ami 2 ðX ÞS þ /Xl Xm a2 ðX ÞS þ/Xm Xi a2 ðX ÞS

/Xi1 ai12 S

:p

2,1

2 X

¼

/Xi2 ai11 S : p

þ /Xm Xi Xj al1 ðX ÞS,

i,j,l,m ¼ 1, . . . ,N

2,2



E22 i2 i2

i1 ,i2 ¼ 1

m þ /Xj Xm ali2 ðX ÞS þ /Xl Xi ajm 2 ðX ÞS þ/Xi Xj Xl a1 ðX ÞS

þ /Xj Xl Xm ai1 ðX ÞS þ/Xl Xm Xi aj1 ðX ÞS

E22 i1 i1

2 X

¼

22 E22 i1 i2  Ei2 i1

i1 ,i2 ¼ 1

ðA:23Þ

The sum of the above two terms yields the desired P 2 matrix. In 1

where a ðX Þ is a fourth order tensor. In the trivial case where 4

i¼j ¼l ¼m the last equation simplifies to d/Xi4 S dt

P2 ¼ 1

iii 2 ii ¼ /aiiii 4 ðX ÞS þ4/Xi a3 ðX ÞS þ 6/Xi a2 ðX ÞS

þ4/Xi3 ai2 ðX ÞS,

the general case where we have N state variables then

i ¼ 1, . . . ,N

ðA:24Þ

Appendix B. Appropriate permutation sets for P m

N X

N X

EiNN  ENN i2 i2 þ 1 i1

i1 ,i2 ¼ 1

EiNN  ENN i2 i1 1 i2

ðB:4Þ

i1 ,i2 ¼ 1

 m¼3 According to Eq. (19) we have d/X ð3Þ S ¼ vecða3 Þ þ P 3 /X  vecða2 ÞS þ P 3 /X ð2Þ  vecða1 ÞS 1 2 dt

j

Recall that P m is defined as follows: j

m

Pm ¼ j

ðjÞ X

p

k¼1

ðB:1Þ

m,k

where each permutation matrix, p

p

m,k

¼

N X

m,k

unique terms. For instance, the term /Xi1 ai23 i2 S is equivalent to ðB:2Þ

m

Epq ij

1

of the form /Xi1 ai22 i3 S or /Xi2 ai23 i1 S or /Xi3 ai21 i2 S. Note that these three permutation of indices i1–i3 are the only ones producing

, is defined as

NN NN ENN i1 i1  Ei2 i2     Eim im i1 ,...,im ¼ 1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

Starting with P 3 , the column vector /X  vecða2 ÞS has entries

is a p  q matrix (elementary matrix) that has a ‘‘1’’ in where the (i,j) position and zeros everywhere else and k denotes the different permutation sets for the subscripts of all Epq . ij

/Xi1 ai22 i3 S and therefore is not unique. Then the corresponding permutation matrices, p , are m,k

/Xi1 ai22 i3 S : p

3,1

¼

2 X i1 ,i2 ¼ 1

22 22 E22 i1 i1  Ei2 i2  Ei3 i3

276

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

/Xi2 ai23 i1 S : p

/Xi3 ai21 i2 S : p

3,2

3,3

2 X

¼

The sum of the above ð41Þ ¼ 4 terms yields the desired P 4

22 22 E22 i1 i2  Ei2 i3  Ei3 i1

1

i1 ,i2 ¼ 1 2 X

¼

matrix. In the general case where we have N state variables then 22 22 E22 i1 i3  Ei2 i1  Ei3 i2

P4 ¼

i1 ,i2 ¼ 1

1

The sum of the above ð31Þ ¼ 3 terms yields the desired P 3 matrix. 1

1

N X

NN NN ENN i1 i1  Ei2 i2  Ei3 i3 þ

i1 ,i2 ¼ 1

NN NN ENN i1 i2  Ei2 i3  Ei3 i1

NN NN ENN i1 i3  Ei2 i1  Ei3 i2

ðB:5Þ

N X

þ

NN NN NN ENN i1 i3  Ei2 i4  Ei3 i1  Ei4 i2

NN NN NN ENN i1 i4  Ei2 i1  Ei3 i2  Ei4 i3

/Xi1 Xi2 ai13 S

/Xi2 Xi3 ai11 S

Using similar arguments we can also conclude that the above permutation matrices are those that apply in the case of P 4 .

/Xi3 Xi1 ai12 S.

or or entries of the form Note that these three permutation of indices i1–i3 are the only

3

Therefore

ones producing unique terms. For instance, the term /Xi1 Xi3 ai12 S

P4 ¼ P4 3

is equivalent to /Xi3 Xi1 ai12 S and therefore is not unique. Then the corresponding permutation matrices, p

3,1

3,2

3,3

¼

2 X

m,k

2

/X ð2Þ  vecða2 ÞS and the unique entries of column are of the form /Xi1 Xi2 ai23 i4 S or /Xi2 Xi3 ai24 i1 S or /Xi3 Xi4 ai21 i2 S or /Xi4 Xi1 ai22 i3 S or /Xi2 Xi4 ai23 i1 S or /Xi3 Xi1 ai22 i4 S. Note that there are no other unique sequences i1, i2, i3 and i4 that would create unique column entries. On the other hand, for example we could have substituted

i1 ,i2 ¼ 1

¼

2 X

22 22 E22 i1 i2  Ei2 i3  Ei3 i1

/Xi2 Xi4 ai23 i1 S with /Xi2 Xi4 ai23 i1 S, but we could not use both as they are dependent. This dependence is even more evident if we recall that the tensors am are symmetric. Then the corresponding permutation matrices, p , are

i1 ,i2 ¼ 1

¼

2 X

ðB:8Þ

1

In the case of P 4 the corresponding column vector is

, are

22 22 E22 i1 i1  Ei2 i2  Ei3 i3

22 22 E22 i1 i3  Ei2 i1  Ei3 i2

i1 ,i2 ¼ 1

m,k

The sum of the above ð32Þ ¼

3

3 terms yields the desired P matrix,

/Xi1 Xi2 ai23 i4 S : p

2

which if you notice is equal to P 3 . This also implies that in the

4,1

¼

1

3

2

1

P ¼P

/Xi2 Xi3 ai24 i1 S : p

ðB:6Þ

 m¼4 According to Eq. (19) we have

/Xi3 Xi4 ai21 i2 S : p

d/X ð4Þ S ¼ vecða4 Þ þ P 4 /X  vecða3 ÞS 1 dt þ P 4 /X ð2Þ  vecða2 ÞS þ P 4 /X ð3Þ  vecða1 ÞS 2

/Xi4 Xi1 ai22 i3 S : p

4,2

4,3

4,4

/Xi2 Xi4 ai23 i1 S : p

Starting with P , the column vector /X  vecða3 ÞS has unique

¼

¼

/Xi4 ai31 i2 i3 S. Then the corresponding permutation matrices, p

4,5

m,k

/Xi3 Xi1 ai22 i4 S : p

,

4,6

are

/Xi2 ai33 i4 i1 S : p

/Xi3 ai34 i1 i2 S

:p

/Xi4 ai31 i2 i3 S : p

4,2

4,3

4,4

22 22 22 E22 i1 i2  Ei2 i3  Ei3 i4  Ei4 i1

2 X

22 22 22 E22 i1 i3  Ei2 i4  Ei3 i1  Ei4 i2

i1 ,i2 ¼ 1

¼

2 X

22 22 22 E22 i1 i4  Ei2 i1  Ei3 i2  Ei4 i3

i1 ,i2 ¼ 1

¼

¼ ¼

2 X i1 ,i2 ¼ 1 2 X

¼

22 22 22 E22 i1 i2  Ei2 i4  Ei3 i3  Ei4 i1

¼

2 X

22 22 22 E22 i1 i3  Ei2 i1  Ei3 i2  Ei4 i4

i1 ,i2 ¼ 1

The sum of the above ð42Þ ¼ 6 terms yields the desired P 4 matrix. In E22 i1 i1



E22 i2 i2



E22 i3 i3



2

E22 i4 i4

the general case where we have N state variables then N X

P4 ¼ P4 þ

22 22 22 E22 i1 i2  Ei2 i3  Ei3 i4  Ei4 i1

2

1

i1 ,i2 ¼ 1 2 X

2 X i1 ,i2 ¼ 1

1

4,1

2 X i1 ,i2 ¼ 1

entries of the form /Xi1 ai32 i3 i4 S or /Xi2 ai33 i4 i1 S or /Xi3 ai34 i1 i2 S or

:p

22 22 22 E22 i1 i1  Ei2 i2  Ei3 i3  Ei4 i4

3

4

/Xi1 ai32 i3 i4 S

2 X i1 ,i2 ¼ 1

general case for N state variables 3

ðB:7Þ

i1 ,i2 ¼ 1

2

/Xi3 Xi1 ai12 S : p

NN NN NN ENN i1 i2  Ei2 i3  Ei3 i4  Ei4 i1

i1 ,i2 ¼ 1

Continuing with P 3 , the column vector /X ð2Þ  vecða1 ÞS has

/Xi2 Xi3 ai11 S : p

N X

þ

i1 ,i2 ¼ 1

/Xi1 Xi2 ai13 S : p

N X i1 ,i2 ¼ 1

i1 ,i2 ¼ 1

N X

þ

N X

NN NN EiNN  ENN i2 i2  Ei3 i3  Ei4 i4 1 i1

i1 ,i2 ¼ 1

þ

In the general case where we have N state variables then

P3 ¼

N X

þ 22 22 22 E22 i1 i3  Ei2 i4  Ei3 i1  Ei4 i2

NN NN NN ENN i1 i2  Ei2 i4  Ei3 i3  Ei4 i1

i1 ,i2 ¼ 1 2 X

NN NN NN ENN i1 i3  Ei2 i1  Ei3 i2  Ei4 i4

ðB:9Þ

i1 ,i2 ¼ 1

i1 ,i2 ¼ 1

¼

2 X i1 ,i2 ¼ 1

22 22 22 E22 i1 i4  Ei2 i1  Ei3 i2  Ei4 i3

The same methodology can be applied to construct any given

P m matrix. Though as m increases the terms in the sum of Eq. (B.1) j

V. Sotiropoulos, Y.N. Kaznessis / Chemical Engineering Science 66 (2011) 268–277

increase which implies that p

m,k

increase too. The idea is to find all

the permutations sets of i1, i2,y,im that uniquely describe the elements of the column vector, /X ðjÞ  vecðamj ÞS and then substitute them in Eq. (B.2).

Appendix C. Supplementary data Supplementary data associated with this article can be found in the online version at doi:10.1016/j.ces.2010.10.024.

References Brewer, J., 1978. Kronecker products and matrix calculus in system theory. IEEE Transactions on Circuits and Systems 25, 772–781. Cao, Y., Gillespie, D.T., Petzold, L.R., 2005. The slow-scale stochastic simulation algorithm. The Journal of Chemical Physics 122, 014116. Cao, Y., Gillespie, D.T., Petzold, L.R., 2006. Efficient step size selection for the tauleaping simulation method. Journal of Chemical Physics 124, 44109. Elowitz, M.B., Levine, A.J., Siggia, E.D., Swain, P.S., 2002. Stochastic gene expression in a single cell. Science 297, 1183–1186. Gardiner, C., 2004. Handbook of Stochastic Methods for Physics, Chemistry and Natural Sciences, Synergetics third ed Springer, Berlin. Gardner, T.S., Cantor, C.R., Collins, J.J., 2000. Construction of a genetic toggle switch in escherichia coli. Nature 403, 339–342. Gibson, M.A., Bruck, J., 2000. Efficient exact stochastic simulation of chemical systems with many species and many channels. Journal of Physical Chemistry A 104, 1876. Gillespie, C.S., 2009. Moment-closure approximations for mass-action models. Systems Biology, IET 3, 52. Gillespie, D., 1992. Markov Processes: An Introduction for Physical Scientists. Academic Press. Gillespie, D.T., 1977. Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry 81, 2340–2361. Gillespie, D.T., 2000. The chemical langevin equation. Journal of Chemical Physics 113, 297. Goutsias, J., 2007. Classical versus stochastic kinetics modeling of biochemical reaction systems. Biophysical Journal 92, 2350–2365. Gunawan, R., Cao, Y., Petzold, L., Doyle, R.F.J., 2005. Sensitivity analysis of discrete stochastic systems. Biophysical Journal 88, 2530–2540. Hespanha, J.P., Singh, A., 2005. Stochastic models for chemically reacting systems using polynomial stochastic hybrid systems. International Journal of Robust and Nonlinear Control 15, 669–689. Hill, A.D., Tomshine, J.R., Weeding, E.M.B., Sotiropoulos, V., Kaznessis, Y.N., 2008. Synbioss: the synthetic biology modeling suite. Bioinformatics 24, 2551–2553. Kaznessis, Y.N., 2006. Multi-scale models for gene network engineering. Chemical Engineering Science 61, 940. Kaznessis, Y.N., 2007. Models for synthetic biology. BMC Systems Biology 1, 47. Krishnarajah, I., Cook, A., Marion, G., Gibson, G., 2005. Novel moment closure approximations in stochastic epidemics. Bulletin of Mathematical Biology 67, 855–873.

277

Lee, C.H., Kim, K., Kim, P., 2009. A moment closure method for stochastic reaction networks. The Journal of Chemical Physics 130, 134107. McQuarrie, D.A., 1967. Stochastic approach to chemical kinetics. Journal of Applied Probability 4, 413. Mohammad-Djafari, A., 2001. A matlab program to calculate the maximum entropy distributions. Moyal, J., 1949. Stochastic processes and statistical physics. Journal of the Royal Statistical Society. Series B (Methodological) 11, 150–210. Munsky, B., Khammash, M., 2006. The finite state projection algorithm for the solution of the chemical master equation. Journal of Chemical Physics 124, 44104. Ramalingam, K., Tomshine, J., Maynard, J., Kaznessis, Y., 2009. Forward engineering of synthetic bio-logical and gates. Biochemical Engineering Journal 47, 38–47. Rathinam, M., Petzold, L.R., Yang, C., Gillespie, D.T., 2003. Stiffness in stochastic chemically reacting systems: the implicit tau-leaping method. Journal of Chemical Physics 119, 12784. Salis, H., Kaznessis, Y., 2005a. Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. Journal of Chemical Physics 122, 54103. Salis, H., Kaznessis, Y.N., 2005b. An equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks. Journal of Chemical Physics 123, 214106. Salis, H., Kaznessis, Y.N., 2006. Computer-aided design of modular protein devices: Boolean and gene activation. Physical Biology 3, 295–310. Salis, H., Sotiropoulos, V., Kaznessis, Y.N., 2006. Multiscale hy3s: hybrid stochastic simulation for supercomputers. BMC Bioinformatics 7, 93. Samant, A., Vlachos, D.G., 2005. Overcoming stiffness in stochastic simulation stemming from partial equilibrium: a multiscale monte carlo algorithm. Journal of Chemical Physics 123, 144114. Singh, A., Hespanha, J. Moment closure techniques for stochastic models in population biology. In: American Control Conference, 2006, p. 6. Sotiropoulos, V., Kaznessis, Y.N., 2007. Synthetic tetracycline-inducible regulatory networks: computer-aided design of dynamic phenotypes. BMC Systems Biology 1, 7. Sotiropoulos, V., Kaznessis, Y.N., 2008. An adaptive time step scheme for a system of stochastic differential equations with multiple multiplicative noise: chemical langevin equation, a proof of concept. The Journal of Chemical Physics 128, 014103. Sotiropoulos, V., Contou-Carrere, M., Daoutidis, P., Kaznessis, Y.N., 2009. Model reduction of multiscale chemical langevin equations: a numerical case study. Computational Biology and Bioinformatics, IEEE/ACM Transactions on 6, 470–482. St-Pierre, F., Endy, D., 2008. Determination of cell fate selection during phage lambda infection. Proceedings of the National Academy of Sciences 105, 20705. Thakur, A.K., Rescigno, A., DeLisi, C., 1978. Stochastic theory of second-order chemical reactions. The Journal of Physical Chemistry 82, 552–558. Tomshine, J., Kaznessis, Y.N., 2006. Optimization of a stochastically simulated gene network model via simulated annealing. Biophysical Journal 91, 3196–3205. Tuttle, L., Salis, H., Tomshine, J., Kaznessis, Y., 2005. Model-driven designs of an oscillating gene network. Biophysical Journal 89, 3873–3883. van Kampen, N., 1992. Stochastic Processes in Physics and Chemistry. North Holland. Wolf, D.M., Arkin, A.P., 2002. Fifteen minutes of fim: control of type 1 pili expression in e. coli. Omics 6, 91–114.