A Convex Approach to Steady State Moment Analysis for Stochastic Chemical Reactions
aa r X i v : . [ q - b i o . M N ] A ug A Convex Approach to Steady State MomentAnalysis for Stochastic Chemical Reactions
Yuta Sakurai and Yutaka Hori ∗† Abstract
Model-based prediction of stochastic noise in biomolecular reactions often resorts to ap-proximation with unknown precision. As a result, unexpected stochastic fluctuation causes aheadache for the designers of biomolecular circuits. This paper proposes a convex optimizationapproach to quantifying the steady state moments of molecular copy counts with theoreticalrigor. We show that the stochastic moments lie in a convex semi-algebraic set specified bylinear matrix inequalities. Thus, the upper and the lower bounds of some moments can becomputed by a semidefinite program. Using a protein dimerization process as an example, wedemonstrate that the proposed method can precisely predict the mean and the variance of thecopy number of the monomer protein.
A grand challenge of synthetic biology is to build and control layers of artificial biomolecularreaction networks that enable complex tasks by utilizing naturally existing reaction resources inmicro-scale organisms such as E. coli. In control engineering community, many theoretical toolswere developed over the last decade to enable model-guided design of biomolecular circuits includinga fold-change detector [1], oscillators [2–5], an event detector [6], and a light-controlled feedbackcircuit for setpoint control [7, 8]. Other successfully implemented examples include logic gates [9],a genetic memory [10], and a communication system between cells [11], to name a few. Given theelementary modules of artificial bimolecular networks, the next step stone is to assemble the circuitparts to build multiple layers of biomolecular networks that can provide practical functions.As is the case with mechanical and electrical engineering, guaranteeing the performance of in-dividual circuit modules at high precision is a key to successfully building large-scale and robustbiomolecular circuits. Biosystems engineering, however, has a unique challenge, in contrast withother fields of engineering, that the signal-to-noise ratio is so low that the heterogeneity of circuitstates between biological cells is almost impossible to avoid. A major source of the heterogeneityis the stochastic chemical reactions that are caused by the low copy nature of reactive molecules ina small-volume reactor, that is, a biological cell. In other words, the events of molecular collisionthat fire reactions are more accurately captured by stochastic process rather than continuous and ∗ This work was supported in part by JSPS KAKENHI Grant Number JP16H07175, Okawa Foundation ResearchGrant under grant number 16-10, Keio Gijuku Academic Development Funds and Research Grant of Keio Leading-edge Laboratory of Science and Technology. c (cid:13) † Y. Sakurai and Y. Hori are with Department of Applied Physics and Physico-Informatics, Keio University,Japan. [email protected] , [email protected] eterministic process governed by ordinary differential equations (ODEs). This situation raisesa strong need for the theoretical tools that can rigorously certify the performance of stochasticchemical reactions using statistical norms such as mean and variance of molecular copy numbers.The random fluctuation of the copy number of molecules can be regarded as a Markov processthat follows the chemical master equation (CME) [12]. Despite a linear ODE, the CME is oftenhard to solve analytically because of the infinite dimensionality of the equation. A typical solutionto this problem is to use Monte Carlo simulations [13]. This approach, however, requires highcomputational time. Moreover, a strict error bound is hard to obtain. To resolve these issues, thefinite state project [14] provides a systematic way to truncate the equation with an error bound.Although useful for analyzing transient dynamics, the FSP still suffers from high dimensionalityif one wants to evaluate the probability distribution near the steady state, which is often theoperation point of interest.A more direct approach to evaluating the statistics of biomolecular circuits is to use the momentequation, an infinite dimensional ODE derived from the CME. Recently developed moment closureapproaches [15,16] allow for approximately solving the ODE by order reduction technique. The cen-tral idea is to approximate the high order moments in the equation using the lower order momentsto derive a finite order closed ODE. Using this method, we can efficiently and directly computethe statistics of the molecular copy numbers of biomolecular circuits. Although appealing, thisapproach requires assuming the underlying probability distribution, which is not apriori in prac-tice. As a result, biocircuit designers have to engineer biomolecular circuits based on approximatedstatistics with unknown precision.The objective of this paper is to propose an algebraic approach to rigorously quantifying thesteady state statistics of stochastic biomolecular reactions without assuming underlying probabilitydistributions. Specifically, we formulate a semidefinite program that computes the upper and thelower bounds of statistical values of molecular copy numbers using the moment approach [ ? ].During the review process, the authors were informed that the same approach has been exploredindependently by two other research groups [ ? , ? ]. In comparison with these works, this paperpresents an optimization algorithm that directly computes the upper bound of the second order central moment, or variance, without running multiple optimizations for computing raw moments,allowing for direct and tighter quantification of the upper bound. We illustrate the proposedmethod using a stochastic protein dimerization process.The organization of this paper is as follows. In Section II, we introduce the CME and derivethe moment equation. Then, the optimization problem for moment computation is formulatedin Section III. In section IV, we analyze the mean and the variance of a stochastic dimerizationprocess. Finally, Section V summarizes our findings and concludes this paper.The following notations are used in this paper. N := { , , , · · · } . N + := { , , , · · · } . R + := { x ∈ R | x ≥ } . deg( p ( x )) is the degree of the polynomial p ( x ). For multivariate polynomials p ( x ) = Q nj =1 x p j j , deg( p ( x )) := P nj =1 p j . We consider a chemical reaction network consisting of n ∈ N + species of molecules and denotethe copy number of each molecular species by x = [ x , x , · · · , x n ] T ∈ N n . There are r types ofreactions in the reaction network, and let s i = [ s i , s i , · · · , s in ] T ∈ Z n be the increment of themolecular copy x by the i -th chemical reaction ( i = 1 , , · · · , r ).hemical reactions are caused by collisions of molecules. When the reaction volume is sufficientlylarge, the copy number of molecules is so large that the collision event occurs almost continuouslyin time. On the other hand, in small reactor systems such as microbes, the collision events tend tobe stochastic as the copy number of molecules is small. This results in the stochastic fluctuationof x . It is known that stochastic chemical reactions can be modeled by a Markov process. Morespecifically, we define P x ( t ) as the probability that there are x molecules at time t . The stochasticchemical reactions can be modeled by the following Chemical Master Equation (CME) [12], dP x ( t ) dt = r X i =1 { w i ( x − s i ) P x − s i ( t ) − w i ( x ) P x ( t ) } , (1)where w i ( x ) is the transition rate associated with the i -th chemical reaction and is defined by w i ( x ) = lim ∆ t → P x + s i ( t + ∆ t ) − P x ( t )∆ t ( i = 1 , , · · · , r ) . (2)The transition rate w i ( x ) is characterized by the rate of the reaction i ( i = 1 , , · · · , r ) We assumethat all of the r reactions are elementary as non-elementary reactions can always be decomposedinto elementary reactions. For elementary reactions, the transition rate w i ( x ) is a polynomialof x i ’s, and we utilize this fact in the following subsection to derive the equations of stochasticmoments.As x is the vector of non-negative positive integers N n , the CME (1) is a linear but an infinitedimensional ODE with respect to P x ( t ). An analytic solution to the CME is, thus, hardly obtainedexcept for some simple examples. In the next subsection, we derive the equations of momentdynamics to directly compute the statistical values without computing the distribution P x ( t ). In this section, we derive the model of moment dynamics based on the CME (1). First, we definea stochastic moment by m α := E n Y j =1 x α j j = ∞ X x =0 ∞ X x =0 · · · ∞ X x n =0 n Y j =1 x α j j P x ( t ) , (3)where α := [ α , α , · · · , α n ] ∈ N n . We refer to the sum of the entries of α as the order of themoment and denote by deg( m α ) := P j α j with a little abuse of notation.The stochastic moments carry all the information of the probability distribution, and some ofthe low order moments are useful to define the design specification of stochastic biomolecularcircuits. Figure 1 illustrates an example of the distribution of a molecular copy number x of somebiomolecular reaction. In this case, the analytic form of the distribution is not known, and thus itis more appropriate to specify the design goal by statistical norms such mean E [ x ] = m and thestandard deviation p E [( x − E [ x ]) ] = p m − m of the distribution using the moments.To derive the dynamics of stochastic moments, we multiply the model (1) by Q j x α j j and calculatethe expected value by taking sum over N n to obtain the following differential equation. ddt X x n Y j =1 x α j j P x ( t ) = X x n Y j =1 x α j j r X i =1 { w i ( x − s i ) P x − s i ( t ) − w i ( x ) P x ( t ) } , (4)igure 1: An example of probability distribution P x ( t ) and its relationship with stochastic momentswhere we use P x to abbreviate the sum over the positive orthant. X x := ∞ X x =0 ∞ X x =0 · · · ∞ X x n =0 . (5)As w i ( x ) is a polynomial of x i ( i = 1 , , · · · , n ), we can rewrite the equation (4) as ddt m α = X β X β · · · X β n A α , β m β ( α ∈ N n ) , (6)where A α , β are the coefficients of the monomials of the polynomial r X i =1 X x { ( x + s i ) α ( x + s i ) α · · · ( x n + s in ) α n − x α x α · · · x α n n } w i ( x ) . (7)Note that the range of the summation in (6) depends on the degree of the polynomial (7). In manycases, the higher order moments are necessary to characterize the low order moments, i.e., thereare moments m β satisfying deg( m β ) > deg( m α ) in the right-hand side of (6). More specifically,when the reactions are elementary, it is reasonable to assume that the polynomial order of thetransition rates w i ( x ) are at most two since most practical reactions are unimolecular or bimolecularreactions [17]. When deg( w i ( x )) = 0 or 1, m β in the right-hand side of equation (6) satisfies0 ≤ deg( m β ) ≤ deg( m α ) , (8)implying that the moment equation is closed, i.e , the dynamics of the i -th order moment can bewritten by the i -th or lower order moments. On the other hand, when one or more w i ( x )’s arequadratic, i.e., deg( w i ( x )) = 2, 0 ≤ deg( m β ) ≤ deg( m α ) + 1 . (9)Thus, the moment equation (6) becomes an infinite dimensional coupled linear ODE, which is hardto solve analytically. Moment Analysis of Stochastic Chemical Reactions
When we design stochastic chemical reactions using a genetic circuit, the specification of thecircuit is often given by a set of statistical constraints such as the mean and the variance of thecopy numbers of the molecules. In this paper, we consider the following problem.
Problem 1.
For a given CME (1), find the lower and the upper bounds of the mean and thevariance of molecular copy number x at steady state.To solve this problem, let dm α /dt = 0 in the equation (6), and consider the subset of linearequations by limiting the equations up to the µ -th order moments. Specifically,0 = X β X β · · · X β n A α , β m β ( α ∈ A µ ⊂ N n ) , (10)where the set A µ is defined by A µ := ( α ∈ N n (cid:12)(cid:12)(cid:12)(cid:12) ≤ n X i α i ≤ µ ) . (11)In what follows, we refer to µ as the truncation order. When deg( w i ( x )) = 0 or 1, the number ofequations is the same as the number of variables (see (8)), and thus we can compute the uniquesolution (if the equations are not degenerated). In the more general case where the transitionrates w i ( x ) are quadratic, the right-hand side of (10) depends on the higher order moments, i.e., deg( m β ) = µ + 1, as shown in (9) Consequently, the equations become underdetermined, and wecan only conclude that the solution lies on a certain hyperplane unless we consider an infinitenumber of equations.To further specify the solution space of the moments m β , we utilize the fact that m β must bethe moments of some (probability) measure defined on the positive orthant. In particular, we use aso-called moment condition to constrain m β and formulate a semidefinite program solving Problem1. m β to be moments Let x p be the vector that consists of all of the monomial bases satisfying deg( Q j x p j j ) = p , and( H ) p,q be a matrix of the form ( H ) p,q := E h x p ( x q ) T i , (12)where E [ x p ( x q ) T ] represents the entry-wise expected value of the matrix x p ( x q ) T . In a similarmanner, we define ( H k ) p,q by( H k ) p,q := E h x k x p ( x q ) T i ( k = 1 , , · · · , n ) . (13)It should be noted that the entries of the matrices ( H ) p,q and ( H k ) p,q can be written with m β satisfying deg( m β ) = ( p + q for ( H ) p,q p + q + 1 for ( H k ) p,q . (14)sing the matrices ( H ) p,q and ( H k ) p,q , we define the following block Hankel matrices H ( γ ) ( { m β } )and H ( γ ) k ( { m β } ) ( k = 1 , , · · · , n ). H ( γ ) ( { m β } ) := ( H ) , ( H ) , · · · ( H ) ,γ ( H ) , ( H ) , · · · ( H ) ,γ ... ... . . . ...( H ) γ , ( H ) γ , · · · ( H ) γ ,γ , (15) H ( γ ) k ( { m β } ) := ( H k ) , ( H k ) , · · · ( H k ) ,γ ( H k ) , ( H k ) , · · · ( H k ) ,γ ... ... . . . ...( H k ) γ , ( H k ) γ , · · · ( H k ) γ ,γ , (16)where γ and γ are determined as follows. γ = (cid:26) ( µ + 1) / µ is odd) µ/ µ is even) , (17) γ = (cid:26) ( µ − / µ is odd) µ/ µ is even) . (18)Using H ( γ ) ( { m β } ) and H ( γ ) k ( { m β } ), the following proposition provides linear matrix inequality(LMI) conditions that the moments of some (probability) measure on R n + must satisfy. Proposition 1. [18] Consider a sequence ( m β ) β ∈A µ +1 with the set A µ +1 defined by (11). Thesequence ( m β ) β ∈A µ +1 constitutes moments of some measure defined on the positive orthant R n + := { x ∈ R n | x k ≥ k = 1 , , · · · , n ) } only if H ( γ ) ( { m β } ) (cid:23) , (19) H ( γ ) k ( { m β } ) (cid:23) k = 1 , , · · · , n ) . (20)The conditions (19) and (20) are called a moment condition and localized moment conditions,respectively. Combining Proposition 1 with (10), the following theorem specifies semi-algebraicconditions that the moments of the stochastic chemical reactions must satisfy. Theorem 1.
Consider the stochastic chemical reaction modeled by the equation (1). Let m ∗ β denote the steady state moments of the random variable x . For a given truncation order µ , themoments m ∗ β satisfy the following conditions.0 = X β X β · · · X β n A α , β m ∗ β ( α ∈ A µ ⊂ N n ) ,H ( γ ) ( { m ∗ β } ) (cid:23) , (21) H ( γ ) k ( { m ∗ β } ) (cid:23) k = 1 , , · · · , n ) , where A α , β are the coefficients of the polynomial (7).Theorem 1 implies that the moments of the stochastic chemical reaction governed by (1) lie in thesemi-algebraic set (21). Thus, Problem 1 can be recast as a relaxation problem over the convexsemi-algebraic set. In particular, we show, in the following section, that the bounds of the meanand the variance can be computed by semidefinite programming. emark 1. In the special case of n = 1 random variable, the conditions (19) and (20) with γ , γ → ∞ become both necessary and sufficient for the sequence ( m β ) β ∈ N to be the moments ofsome measure on R + [18]. The associated problem is called Stieltjes moment problem named afterthe renowned mathematician Thomas Stieltjes. Similar conditions are available for the momentsof measures defined on various domains such as [0 , R and semi-algebraic sets (see [18, 19] forexample). Theorem 1 asserts that the lower and the upper bounds of the steady state statistics, say f ( m α ),can be calculated by solving the following optimization problems, respectively.min f ( m α )subject to (21) , (22)max f ( m α )subject to (21) . (23)Note that the constrains are convex. In what follows, we show that the bounds of mean andvariance, the two most important statistics of stochastic chemical reactions, can be formulated assemidefinite programs. Mean
The objective function f ( m α ) is obviously linear, and thus the problem falls into the classof semidefinite programming. For example, if one wants to compute the mean of the copy numberof the i -th molecule, the objective function should be f ( m α ) = E [ x i ] = m e i (24)with e i representing the row vector whose i -th entry is 1 and 0 elsewhere.In general, the gap between the upper and the lower bounds decreases monotonically with theincrease of the truncation order µ . Thus, there is a tradeoff between computational cost and theconservativeness of the bounds. Upper bound of variance
By definition, the variance of the copy number of the i -th moleculeis given by f ( m α ) = e i E (cid:2) ( x − E [ x ])( x − E [ x ]) T (cid:3) e T i = m e i − m e i . (25)The second term of the objective function f ( m α ) = m e i − m e i is quadratic, and the optimizationproblem (23) is seemingly not semidefinite programming. However, using the Schur complement,we can convert the optimization (23) into an equivalent semidefinite programming as follows. Theorem 2.
The solution of the optimization problem (23) with the objective function (25) isthe same as that of the following optimization problem.max V subject to (21) (26) (cid:20) m e i − V m e i m e i (cid:21) (cid:23) O. The proof of this theorem can be found in Appendix A. = = ( − 1) = Figure 2: Genetic circuit producing monomer and dimer proteins
Lower bound of variance
The lower bound of variance can be also computed by semidefiniteprogramming with additional relaxation. Specifically, we first compute the upper bound of themean value m e i , and then substitute the result into m e i of the equation (25) to obtain a lin-ear objective function. The estimated lower bound of variance monotonically increases with thetruncation order µ since the upper bound of the mean value m e i monotonically decreases.As the optimization problems shown in this subsection use convex relaxation, there may be thecase where the truncation order µ needs to be unrealistically large to reduce the conservativenessof the bounds to a practically useful level. To the authors’ experience, however, one will onlyneed to use µ = 6 to 12-th order moments to obtain a practically useful bounds for small reactionnetworks with n = 3 to 5 molecules. The general quantitative analysis, however, is still an openquestion. Remark 2.
The moment closure approaches [15, 16] allow for approximating the high order mo-ments, m β with deg( m β ) = µ + 1, using the low order moments, m β with deg( m β ) ≤ µ , to solvethe equation (10). This approach, however, is based on assumptions on the underlying distribution P x ( t ). Thus, the quantification of the approximation error is not easy. On the other hand, the pro-posed semidefinite programs are advantageous in that they can compute the upper and the lowerbounds of the statistics of stochastic chemical reactions with theoretical rigor without assumingan underlying distribution. To illustrate the proposed semidefinite programming approach, we analyze a protein dimerizationprocess, which is one of the simplest examples of stochastic biomolecular reactions. Consider thechemical reactions shown in Fig. 2, where there are three species of molecules, DNA, monomerprotein and dimer protein, and r = 3 reactions. In Fig. 2, the protein monomer is first producedfrom DNA, and then it is either degraded or dimerized. As a result, the copy number of proteinmonomer reaches a steady state at t → ∞ .In what follows, our goal is to analyze the steady state mean and variance of the copy numberof the protein monomer. Let the copy number of monomer protein and its α -th order moment bedefined by x and m α := E [ x α ] = ∞ X x =0 x α P x ( t ) , respectively. It should be noticed that m = 1 by definition. Let D T denote the total copy numberof DNA molecules and assume that D T is a constant. The reaction rates of monomer production,degradation and dimerization w i ( x ) ( i = 1 , ,
3) are then defined by the polynomials in Table 1.able 1: Transition rates w i ( x ) in Fig. 2Reaction i Transition rate w i ( x )1 k D T k x k x ( x − k . − Transcription and translation rate k ln(2) /
20 copy − · min − Protein (monomer) degradation k . − · min − Dimerization rate D T
50 copy Total copy of DNAThe stoichiometry of the reactions, or the increment of the copy number of monomer protein byeach reaction, is given by s = 1 , s = − , s = − . Suppose the values of the parameters are given by Table 2. First, we derive the moment equation(10) by expanding the polynomial (7). Specifically, let the truncation order be µ = 3. Then, wehave = . . . . . . . − . . . . . . − . . . . . . − . m m m m m . (27)It should be noted that the low order moments m , m and m depend on the higher order moment m since the transition rate w ( x ) is quadratic in x (see Table 1). The equation (27) implies thatthe moments of the stochastic reactions are in the nullspace of the matrix. To further constrainthe feasibility set, we use the moment condition and the localized moment condition, which aregiven by m m m m m m m m m (cid:23) O, (cid:20) m m m m (cid:21) (cid:23) O, (28)respectively. In practice, the truncation order µ needs to be determined so that the gap between thelower and the upper bounds is sufficiently small. This might require computing the optimizationproblems (22) and (23) multiple times.To analyze the mean copy number of protein monomer, let the objective function be f ( m α ) = E [ x ] = m and solve the optimization problems (22) and (23). Figure 3 (A) illustrates the computedlower and the upper bounds of the mean monomer protein copy number for each truncation order µ . The gap between the bounds monotonically decreases with increasing µ , and with µ = 5, wecan conclude 13 . ≤ E [ x ] ≤ . , (29)which is a sufficient resolution in practice. Order M e a n o f p r o t e i n c o p y m o l e c u l e Maximum ValueMinimum Value
Order V a r i a n c e o f p r o t e i n c o p y m o l e c u l e Maximum ValueMinimum Value (A) (B)
Figure 3: Upper and lower bounds of mean copy number and variance.Next, to calculate the upper bound of the variance, let f ( m α ) = E [( x − E [ x ]) ] = m − m ,and solve the optimization problem (26). The lower bound is also computed using the relaxationdescribed in Section III. The optimization result in Fig. 3 (B) shows that the variance is boundedby 9 . ≤ E [( x − E [ x ]) ] ≤ . . (30)It should be noticed that these results are mathematically strict and thus enable rigorous quantifi-cation of biocircuit performance using the statistical norms. Remark 3.
The optimization problem was solved with SeDuMi 1.32 on MATLAB 2016b. Toavoid numerical issues, we normalized the variables m i by constants 20 i in the above numericalexample. It should be noted that optimization problems associated with the moment matricestend to be numerically unstable as reported in [ ? ]. In this paper, we have formulated the optimization problem to compute the bounds of the steadystate statistics of stochastic chemical reactions in genetic circuits. First, we have introduced thesteady state moment equation, which is a set of underdetermined linear equations. To restrictthe solution space of the moment equation, our key idea is to use the LMI conditions, whichthe moments of some measure must satisfy. Consequently, we have obtained a convex relaxationproblem that can be solved by semidefinite programming. A distinctive feature of the proposedapproach is that it can provide mathematically rigorous upper and lower bounds of the statisticsfor any stochastic chemical reactions modeled by the CME. To demonstrate the method, we haveanalyzed the protein dimerization process and have obtained tight bounds of the mean and thevariance.Although we have only considered the steady state moments, we can extend the proposed ap-proach to the analysis of transient moments. The authors are currently working to build a stochasticbiocircuit design tool based on the specifications of transient statistics.
References [1] S. Guo, Y. Hori, and R. M. Murray, “Systematic design and implementation of a novel syn-thetic fold-change detector biocircuit in vivo,” tech. rep., California Institute of Technology,2014.2] M. B. Elowitz and S. Leibler, “A synthetic oscillatory network of transcriptional regulators,”
Nature , vol. 403, no. 6767, pp. 335–338, 2000.[3] Y. Hori, T.-H. Kim, and S. Hara, “Existence criteria of periodic oscillations in cyclic generegulatory networks,”
Automatica , vol. 47, no. 5, pp. 1203–1209, 2011.[4] Y. Hori, M. Takada, and S. Hara, “Biochemical oscillations in delayed negative cyclic feedback:existence and profiles,”
Automatica , vol. 49, no. 9, pp. 2581–2590, 2013.[5] H. Niederholtmeyer, Z. Sun, Y. Hori, E. Yeung, A. Verpoorte, R. M. Murray, and S. J. Maerkl,“Rapid cell-free forward engineering of novel genetic ring oscillators,” eLife , vol. 4, p. e09771,2015.[6] V. Hsiao, Y. Hori, P. W. K. Rothemund, and R. M. Murray, “A population-based temporallogic gate for timing and recording chemical events,”
Molecular Systems Biology , vol. 12,no. 869, 2016.[7] A. Milias-Argeitis, S. Summers, J. Stewart-Ornstein, I. Zuleta, D. Pincus, H. El-Samad,M. Khammash, and J. Lygeros, “In silico feedback for in vivo regulation of a gene expressioncircuit,”
Nature Biotechnology , vol. 29, no. 12, pp. 1114–1116, 2011.[8] A. Milias-Argeitis, M. Rullan, S. K. Aoki, P. Buchmann, and M. Khammash, “Automatedoptogenetic feedback control for precise and robust regulation of gene expression and cellgrowth,”
Nature Communications , vol. 7, no. 12, pp. 1114–1116, 2011.[9] T.-S. Moon, C. Lou, A. Tamsir, B. C. Stanton, and C. A. Voigt, “Genetic programs constructedfrom layered logic gates in single cells,”
Nature , vol. 491, no. 7423, pp. 249–253, 2012.[10] L. Yang, A. A. K. Nielsen, J. Fernandez-Rodriguez, C. J. McClune, M. T. Laub, T. K. Lu, andC. A. Voigt, “Permanent genetic memory with > Nature Method , vol. 11,no. 12, pp. 1261–1266, 2014.[11] T. Danino, O. Mondrag´on-Palomino, L. Tsimring, and J. Hasty, “A synchronized quorum ofgenetic clocks,”
Nature , vol. 463, no. 7279, pp. 326–330, 2010.[12] D. T. Gillespie, “A rigorous derivation of the chemical master equation,”
Physica A , vol. 188,no. 1–3, pp. 404–425, 1992.[13] D. T. Gillespie, “A general method for numerically simulating the stochastic time evolutionof coupled chemical reactions,”
Journal of Computational Physics , vol. 22, no. 4, pp. 403–434,1976.[14] B. Munsky and M. Khammash, “The finite state projection algorithm for the solution of thechemical master equation,”
Journal of Chemical Physics , vol. 124, no. 4, p. 044104, 2006.[15] Y.-B. Zhao, J. Kim, and J. P. Hespanha, “Hybrid moment computation algorithm for bio-chemical reaction networks,” in
Proceedings of IEEE Conference on Decision and Control ,pp. 1693–1698, 2010.[16] A. Singh and J. P. Hespanha, “Approximate moment dynamics for chemically reacting sys-tems,”
IEEE Transactions on Automatic Control , vol. 56, no. 2, pp. 414–418, 2011.[17] E. T. Denisov, O. M. Sarkisov, and G. I. Likhteshtein,
Chemical kinetics: fundamentals andnew developments . Elsevier, 2003.18] H. J. Landau,
Moments in Mathematics . American Mathematical Society, 1987.[19] J. A. Shohat and J. D. Tamarkin,
The problem of moments . American Mathematical Society,1943.[20] S. P. Boyd and L. Vandenberghe,
Convex Optimization . Cambridge University Press, 2004.
A Proof of Theorem 2.
Proof
The optimization problem (23) is equivalent to the following problem.max V subject to f ( m α ) ≥ V (21) . As the objective function is f ( m α ) = m e i − m e i , the inequality f ( m α ) ≥ V in the constraint canbe equivalently written by (cid:20) m e i − V m e i m e i (cid:21) (cid:23) O (31)using Schur compliment [20]. The optimal value of (23) is thus equal to that of (26).(31)using Schur compliment [20]. The optimal value of (23) is thus equal to that of (26).