arXiv:1704.07722v1 [q-bio.MN] 25 Apr 2017

A Convex Approach to Steady State Moment
Analysis for Stochastic Chemical Reactions
Yuta Sakurai and Yutaka Hori∗

Abstract
Model-based prediction of stochastic noise in biomolecular reactions often resorts to approximation with unknown precision. As a result, unexpected stochastic fluctuation causes a
headache for the designers of biomolecular circuits. This paper proposes a convex optimization
approach to quantifying the steady state moments of molecular copy counts with theoretical
rigor. We show that the stochastic moments lie in a convex semi-algebraic set specified by
linear matrix inequalities. Thus, the upper and the lower bounds of some moments can be
computed by a semidefinite program. Using a protein dimerization process as an example, we
demonstrate that the proposed method can precisely predict the mean and the variance of the
copy number of the monomer protein.

1

Introduction

A grand challenge of synthetic biology is to build and control layers of artificial biomolecular
reaction networks that enable complex tasks by utilizing naturally existing reaction resources
in micro-scale organisms such as E. coli. In control engineering community, many theoretical
tools were developed over the last decade to enable model-guided design of biomolecular circuits
including a fold-change detector [7], oscillators [4, 8, 9, 16], an event detector [10], and a lightcontrolled feedback circuit for setpoint control [12, 13]. Other successfully implemented examples
include logic gates [14], a genetic memory [19], and a communication system between cells [2], to
name a few. Given the elementary modules of artificial bimolecular networks, the next step stone
is to assemble the circuit parts 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 individual circuit modules at high precision is a key to successfully building large-scale and robust
biomolecular circuits. Biosystems engineering, however, has a unique challenge, in contrast with
other fields of engineering, that the signal-to-noise ratio is so low that the heterogeneity of circuit
states between biological cells is almost impossible to avoid. A major source of the heterogeneity
is the stochastic chemical reactions that are caused by the low copy nature of reactive molecules in
a small-volume reactor, that is, a biological cell. In other words, the events of molecular collision
that fire reactions are more accurately captured by stochastic process rather than continuous and
deterministic process governed by ordinary differential equations (ODEs). This situation raises
a strong need for the theoretical tools that can rigorously certify the performance of stochastic
chemical reactions using statistical norms such as mean and variance of molecular copy numbers.
∗ Y.

Sakurai and Y. Hori are with Department of Applied Physics and Physico-Informatics, Keio University, Japan.

y.sakurai-5861@keio.jp, yhori@appi.keio.ac.jp. This work was supported in part by JSPS KAKENHI
Grant Number JP16H07175 and Keio Gijuku Academic Development Funds.

The random fluctuation of the copy number of molecules can be regarded as a Markov process
that follows the chemical master equation (CME) [6]. Despite a linear ODE, the CME is often
hard to solve analytically because of the infinite dimensionality of the equation. A typical solution
to this problem is to use Monte Carlo simulations [5]. This approach, however, requires high
computational time. Moreover, a strict error bound is hard to obtain. To resolve these issues, the
finite state project [15] provides a systematic way to truncate the equation with an error bound.
Although useful for analyzing transient dynamics, the FSP still suffers from high dimensionality
if one wants to evaluate the probability distribution near the steady state, which is often the
operation point of interest.
A more direct approach to evaluating the statistics of biomolecular circuits is to use the moment
equation, an infinite dimensional ODE derived from the CME. Recently developed moment closure
approaches [18,20] allow for approximately solving the ODE by order reduction technique. The central idea is to approximate the high order moments in the equation using the lower order moments
to derive a finite order closed ODE. Using this method, we can efficiently and directly compute
the statistics of the molecular copy numbers of biomolecular circuits. Although appealing, this
approach requires assuming the underlying probability distribution, which is not apriori in practice. As a result, biocircuit designers have to engineer biomolecular circuits based on approximated
statistics with unknown precision.
The objective of this paper is to propose a novel algebraic approach to rigorously quantifying the
steady state statistics of stochastic biomolecular reactions. Specifically, the main contribution of
this paper is the formulation of an optimization problem that calculates the upper and the lower
bounds of statistical values of molecular copy numbers in biomolecular circuits. We show that the
bounds of the mean and the variance can be obtained by semidefinite programming. In particular,
the proposed method allows for rigorous quantification of the statistics for any stochastic chemical
reactions modeled by the CME as it does not assume underlying probability distribution unlike the
moment closure approach. We illustrate the proposed method by analyzing a stochastic protein
dimerization process.
The organization of this paper is as follows. In Section II, we introduce the CME and derive
the moment equation. Then, the optimization problem for moment computation is formulated
in Section III. In section IV, we analyze the mean and the variance of a stochastic dimerization
process. Finally, Section V summarizes our findings and concludes this paper.
The following notations are used in this paper. N0 := {0, 1, 2, · · · }. N+ := {1, 2, 3, · · · }. R+ :=
{x ∈ R Q
| x ≥ 0}. deg(p(x)) is P
the degree of the polynomial p(x). For multivariate polynomials
p
n
n
p(x) = j=1 xj j , deg(p(x)) := j=1 pj .

2
2.1

Moment Dynamics of Stochastic Chemical Reactions
Model of stochastic chemical reactions

We consider a chemical reaction network consisting of n ∈ N+ species of molecules and denote
the copy number of each molecular species by x = [x1 , x2 , · · · , xn ]T ∈ Nn0 . There are r types of
reactions in the reaction network, and let si = [si1 , si2 , · · · , sin ]T ∈ Zn be the increment of the
molecular copy x by the i-th chemical reaction (i = 1, 2, · · · , r).
Chemical reactions are caused by collisions of molecules. When the reaction volume is sufficiently
large, the copy number of molecules is so large that the collision event occurs almost continuously
in time. On the other hand, in small reactor systems such as microbes, the collision events tend to
be stochastic as the copy number of molecules is small. This results in the stochastic fluctuation
of x. It is known that stochastic chemical reactions can be modeled by a Markov process. More

specifically, we define Px (t) as the probability that there are x molecules at time t. The stochastic
chemical reactions can be modeled by the following Chemical Master Equation (CME) [6],
r
dPx (t) X
{wi (x − si )Px−si (t) − wi (x)Px (t)} ,
=
dt
i=1

(1)

where wi (x) is the transition rate associated with the i-th chemical reaction and is defined by
wi (x) = lim

Px+si(t + ∆t) − Px(t)
(i = 1, 2, · · · , r).
∆t

∆t→0

(2)

The transition rate wi (x) is characterized by the rate of the reaction i (i = 1, 2, · · · , r) We assume
that all of the r reactions are elementary as non-elementary reactions can always be decomposed
into elementary reactions. For elementary reactions, the transition rate wi (x) is a polynomial
of xi ’s, and we utilize this fact in the following subsection to derive the equations of stochastic
moments.
As x is the vector of non-negative positive integers Nn0 , the CME (1) is a linear but an infinite
dimensional ODE with respect to Px (t). An analytic solution to the CME is, thus, hardly obtained
except for some simple examples. In the next subsection, we derive the equations of moment
dynamics to directly compute the statistical values without computing the distribution Px (t).

2.2

The dynamics of moments

In this section, we derive the model of moment dynamics based on the CME (1). First, we define
a stochastic moment by


∞ X
∞
∞ Y
n
n
X
X
Y
α
α
(3)
xj j Px (t),
···
xj j  =
mα := E 
x1 =0x2 =0

j=1

xn =0 j=1

where α := [α1 , α2 , · · · , αn ] ∈ Nn0 . P
We refer to the sum of the entries of α as the order of the
moment and denote by deg(mα ) := j αj with a little abuse of notation.
The stochastic moments carry all the information of the probability distribution, and some of
the low order moments are useful to define the design specification of stochastic biomolecular
circuits. Figure 1 illustrates an example of the distribution of a molecular copy number x of some
biomolecular reaction. In this case, the analytic form of the distribution is not known, and thus it
is more appropriatep
to specify the design
p goal by statistical norms such mean E[x] = m1 and the
standard deviation E[(x − E[x])2 ] = m2 − m21 of the distribution using the moments.
Q α
To derive the dynamics of stochastic moments, we multiply the model (1) by j xj j and calculate
the expected value by taking sum over Nn0 to obtain the following differential equation.




n


X
Y
d
α

xj j  Px (t)

dt  x
j=1

=

n
XY

x j=1

where we use

P

x

α
xj j

r
X

{wi (x−si )Px−si(t)−wi (x)Px (t)} ,

(4)

i=1

to abbreviate the sum over the positive orthant.
X
x

:=

∞ X
∞
X

x1 =0 x2 =0

···

∞
X

xn =0

.

(5)

Figure 1: An example of probability distribution Px (t) and its relationship with stochastic moments
As wi (x) is a polynomial of xi (i = 1, 2, · · · , n), we can rewrite the equation (4) as
XX
X
d
mα =
Aα,β mβ
···
dt
β1

β2

(α ∈ Nn0 ),

(6)

βn

where Aα,β are the coefficients of the monomials of the polynomial
r X
X
i=1

{(x1 + si1 )α1 (x2 + si2 )α2

x

αn
1 α2
· · · (xn + sin )αn − xα
1 x2 · · · xn } wi (x).

(7)

Note that the range of the summation in (6) depends on the degree of the polynomial (7). In many
cases, the higher order moments are necessary to characterize the low order moments, i.e., there
are 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 the
transition rates wi (x) are at most two since most practical reactions are unimolecular or bimolecular
reactions [3]. When deg(wi (x)) = 0 or 1, mβ in the right-hand side of equation (6) satisfies
0 ≤ deg(mβ ) ≤ deg(mα ),

(8)

implying that the moment equation is closed, i.e, the dynamics of the i-th order moment can be
written by the i-th or lower order moments. On the other hand, when one or more wi (x)’s are
quadratic, i.e., deg(wi (x)) = 2,
0 ≤ deg(mβ ) ≤ deg(mα ) + 1.

(9)

Thus, the moment equation (6) becomes an infinite dimensional coupled linear ODE, which is hard
to solve analytically.

3

Moment Analysis of Stochastic Chemical Reactions

3.1

Truncation of the moment equation

When we design stochastic chemical reactions using a genetic circuit, the specification of the
circuit is often given by a set of statistical constraints such as the mean and the variance of the
copy 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 the
variance 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 linear
equations by limiting the equations up to the µ-th order moments. Specifically,
XX
X
Aα,β mβ (α ∈ Aµ ⊂ Nn0 ) ,
(10)
0=
···
β1

βn

β2

where the set Aµ is defined by
Aµ :=

(

α∈

Nn0

0≤

n
X
i

)

αi ≤ µ .

(11)

In what follows, we refer to µ as the truncation order. When deg(wi (x)) = 0 or 1, the number of
equations is the same as the number of variables (see (8)), and thus we can compute the unique
solution (if the equations are not degenerated). In the more general case where the transition
rates wi (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 we
can only conclude that the solution lies on a certain hyperplane unless we consider an infinite
number of equations.
To further specify the solution space of the moments mβ , we utilize the fact that mβ must be
the moments of some (probability) measure defined on the positive orthant. In particular, we use a
so-called moment condition to constrain mβ and formulate a semidefinite program solving Problem
1.

3.2

Necessary condition for mβ to be moments

Q p
Let x be the vector that consists of all of the monomial bases satisfying deg( j xj j ) = p, and
(H)p,q be a matrix of the form
h
i
T
(H)p,q := E xp (xq ) ,
(12)
p

where E[xp (xq )T ] represents the entry-wise expected value of the matrix xp (xq )T . In a similar
manner, we define (Hk )p,q by
i
h
T
(k = 1, 2, · · · , n).
(13)
(Hk )p,q := E xk xp (xq )

It should be noted that the entries of the matrices (H)p,q and (Hk )p,q can be written with mβ
satisfying
(
p+q
for (H)p,q
deg(mβ ) =
.
(14)
p + q + 1 for (Hk )p,q

Using the matrices (H)p,q and (Hk )p,q , we define the following block Hankel matrices H (γ1 ) ({mβ })
(γ )
and Hk 2 ({mβ }) (k = 1, 2, · · · , n).


(H)0,0 (H)0,1 · · · (H)0,γ1
 (H)1,0 (H)1,1 · · · (H)1,γ1 


(15)
H (γ1 )({mβ }) :=  .
,
..
..
..

 ..
.
.
.
(H)γ1 ,0 (H)γ1 ,1 · · · (H)γ1 ,γ1


(Hk )0,0 (Hk )0,1
 (Hk )1,0 (Hk )1,1

(γ )
Hk 2 ({mβ }) := 
..
..

.
.
(Hk )γ2 ,0 (Hk )γ2 ,1


· · · (Hk )0,γ2
· · · (Hk )1,γ2 

,
..
..

.
.
· · · (Hk )γ2 ,γ2

where γ1 and γ2 are determined as follows.

(µ + 1)/2 (if µ is odd)
γ1 =
,
µ/2
(if µ is even)
γ2 =



(µ − 1)/2 (if µ is odd)
.
µ/2
(if µ is even)

(16)

(17)

(18)

(γ )

Using H (γ1 ) ({mβ }) and Hk 2 ({mβ }), the following proposition provides linear matrix inequality
(LMI) conditions that the moments of some (probability) measure on Rn+ must satisfy.
Proposition 1. [11] Consider a sequence (mβ )β∈Aµ+1 with the set Aµ+1 defined by (11). The
sequence (mβ )β∈Aµ+1 constitutes moments of some measure defined on the positive orthant Rn+ :=
{x ∈ Rn | xk ≥ 0 (k = 1, 2, · · · , n)} only if
H (γ1 ) ({mβ })  0,
(γ )
Hk 2 ({mβ })

0

(19)
(k = 1, 2, · · · , 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-algebraic
conditions 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 µ, the
moments m∗β satisfy the following conditions.
0=

XX
β1

H

(γ1 )

···

β2

({m∗β })

(γ )
Hk 2 ({m∗β })

X

Aα,β m∗β

(α ∈ Aµ ⊂ Nn0 ) ,

βn

 0,
0

(21)
(k = 1, 2, · · · , 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 the
semi-algebraic set (21). Thus, Problem 1 can be recast as a relaxation problem over the convex
semi-algebraic set. In particular, we show, in the following section, that the bounds of the mean
and the variance can be computed by semidefinite programming.

Remark 1. In the special case of n = 1 random variable, the conditions (19) and (20) with
γ1 , γ2 → ∞ become both necessary and sufficient for the sequence (mβ )β∈N0 to be the moments of
some measure on R+ [11]. The associated problem is called Stieltjes moment problem named after
the renowned mathematician Thomas Stieltjes. Similar conditions are available for the moments
of measures defined on various domains such as [0, 1], R and semi-algebraic sets (see [11, 17] for
example).

3.3

Semidefinite programming for mean and variance computation

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
subject to

f (mα )
(21),

(22)

max
subject to

f (mα )
(21).

(23)

Note that the constrains are convex. In what follows, we show that the bounds of mean and
variance, the two most important statistics of stochastic chemical reactions, can be formulated as
semidefinite programs.
Mean The objective function f (mα ) is obviously linear, and thus the problem falls into the class
of semidefinite programming. For example, if one wants to compute the mean of the copy number
of the i-th molecule, the objective function should be
(24)

f (mα ) = E[xi ] = mei

with ei 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 the
increase of the truncation order µ. Thus, there is a tradeoff between computational cost and the
conservativeness of the bounds.
Upper bound of variance By definition, the variance of the copy number of the i-th molecule
is given by


f (mα ) = ei E (x − E[x])(x − E[x])T eT
i
= m2ei − m2ei .

(25)

The second term of the objective function f (mα ) = m2ei − m2ei is quadratic, and the optimization
problem (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) is
the same as that of the following optimization problem.
max
subject to

V
(21)


m2ei − V
mei

The proof of this theorem can be found in Appendix A.

(26)
mei
1



 O.

1

=
2

1

3

=

=

3

( − 1)

2

Figure 2: Genetic circuit producing monomer and dimer proteins
Lower bound of variance The lower bound of variance can be also computed by semidefinite
programming with additional relaxation. Specifically, we first compute the upper bound of the
mean value mei , and then substitute the result into m2ei of the equation (25) to obtain a linear objective function. The estimated lower bound of variance monotonically increases with the
truncation order µ since the upper bound of the mean value mei monotonically decreases.
As the optimization problems shown in this subsection use convex relaxation, there may be the
case where the truncation order µ needs to be unrealistically large to reduce the conservativeness
of the bounds to a practically useful level. To the authors’ experience, however, one will only
need to use µ = 6 to 12-th order moments to obtain a practically useful bounds for small reaction
networks with n = 3 to 5 molecules. The general quantitative analysis, however, is still an open
question.
Remark 2. The moment closure approaches [18, 20] allow for approximating the high order moments, mβ with deg(mβ ) = µ + 1, using the low order moments, mβ with deg(mβ ) ≤ µ, to solve
the equation (10). This approach, however, is based on assumptions on the underlying distribution
Px (t). Thus, the quantification of the approximation error is not easy. On the other hand, the proposed semidefinite programs are advantageous in that they can compute the upper and the lower
bounds of the statistics of stochastic chemical reactions with theoretical rigor without assuming
an underlying distribution.

4

Application to a Stochastic Genetic Circuit

To illustrate the proposed semidefinite programming approach, we analyze a protein dimerization
process, which is one of the simplest examples of stochastic biomolecular reactions. Consider the
chemical reactions shown in Fig. 2, where there are three species of molecules, DNA, monomer
protein and dimer protein, and r = 3 reactions. In Fig. 2, the protein monomer is first produced
from DNA, and then it is either degraded or dimerized. As a result, the copy number of protein
monomer reaches a steady state at t → ∞.
In what follows, our goal is to analyze the steady state mean and variance of the copy number
of the protein monomer. Let the copy number of monomer protein and its α-th order moment be
defined by x and
∞
X
xα Px (t),
mα := E [xα ] =
x=0

respectively. It should be noticed that m0 = 1 by definition. Let DT denote the total copy number
of DNA molecules and assume that DT is a constant. The reaction rates of monomer production,
degradation and dimerization wi (x) (i = 1, 2, 3) are then defined by the polynomials in Table 1.

Table 1: Transition rates wi (x) in Fig. 2
Reaction i Transition rate wi (x)
1
k1 DT
k2 x
2
3
k3 x(x − 1)

Table 2: Parameters of chemical reactions in Fig. 2
Parameter
Value
Meaning
k1
0.7 min−1
Transcription and translation rate
k2
ln(2)/20 copy−1·min−1 Protein (monomer) degradation
0.1 copy−1·min−1
Dimerization rate
k3
DT
50 copy
Total copy of DNA

The stoichiometry of the reactions, or the increment of the copy number of monomer protein by
each reaction, is given by
s1 = 1, s2 = −1, s3 = −2.
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, we
have



 m
0.0 0.0
0.0
0.0 0.0  0 
 35.0 0.2 −0.2 0.0 0.0  m1 


0=
(27)
 35.0 69.6 0.7 −0.4 0.0  m2  .
 m3 
35.0 105.8 105.5 1.7 −0.6
m4

It should be noted that the low order moments m1 , m2 and m3 depend on the higher order moment
m4 since the transition rate w3 (x) is quadratic in x (see Table 1). The equation (27) implies that
the moments of the stochastic reactions are in the nullspace of the matrix. To further constrain
the feasibility set, we use the moment condition and the localized moment condition, which are
given by




m0 m1 m2
m1 m2
 m1 m2 m3   O,
 O,
(28)
m2 m3
m2 m3 m4

respectively. In practice, the truncation order µ needs to be determined so that the gap between the
lower and the upper bounds is sufficiently small. This might require computing the optimization
problems (22) and (23) multiple times.
To analyze the mean copy number of protein monomer, let the objective function be f (mα ) =
E[x] = m1 and solve the optimization problems (22) and (23). Figure 3 (A) illustrates the computed
lower 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, we
can conclude
13.26 ≤ E[x] ≤ 13.27,
which is a sufficient resolution in practice.

(29)

(A)

(B)
30

Variance of protein
copy molecule

Mean of protein
copy molecule

15
10
5
0

Maximum Value
Minimum Value

1

3

5

Order

7

9

Maximum Value
Minimum Value

25
20
15
10
5
0

1

3

5

7

9

Order

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])2 ] = m2 − m21 ,
and solve the optimization problem (26). The lower bound is also computed using the relaxation
described in Section III. The optimization result in Fig. 3 (B) shows that the variance is bounded
by
9.94 ≤ E[(x − E[x])2 ] ≤ 9.95.

(30)

It should be noticed that these results are mathematically strict and thus enable rigorous quantification of biocircuit performance using the statistical norms.
Remark 3. We have only shown the example of one variable case for the illustration purpose
in this section, but generalization to n random variables is straightforward. However, a critical
challenge for the multivariate case, say n ≥ 10, is the reduction of computation time. The computation burden mainly comes from the number of the variables, which is given by µ+n Cµ . To
overcome the limitation, it would be worth developing a method to reduce the variables that have
little effect on the calculation results since most variables do not appear in the linear constraint as
n gets larger. Development of the reduction technique is left for our future work.

5

Conclusion

In this paper, we have formulated the optimization problem to compute the bounds of the steady
state statistics of stochastic chemical reactions in genetic circuits. First, we have introduced the
steady state moment equation, which is a set of underdetermined linear equations. To restrict
the solution space of the moment equation, our key idea is to use the LMI conditions, which
the moments of some measure must satisfy. Consequently, we have obtained a convex relaxation
problem that can be solved by semidefinite programming. A distinctive feature of the proposed
approach is that it can provide mathematically rigorous upper and lower bounds of the statistics
for any stochastic chemical reactions modeled by the CME. To demonstrate the method, we have
analyzed the protein dimerization process and have obtained tight bounds of the mean and the
variance.
Although we have only considered the steady state moments, we can extend the proposed approach to the analysis of transient moments. The authors are currently working to build a stochastic
biocircuit design tool based on the specifications of transient statistics.

References
[1] S. P. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
[2] T. Danino, O. Mondragón-Palomino, L. Tsimring, and J. Hasty. A synchronized quorum of
genetic clocks. Nature, 463(7279):326–330, 2010.
[3] E. T. Denisov, O. M. Sarkisov, and G. I. Likhteshtein. Chemical kinetics: fundamentals and
new developments. Elsevier, 2003.
[4] M. B. Elowitz and S. Leibler. A synthetic oscillatory network of transcriptional regulators.
Nature, 403(6767):335–338, 2000.
[5] D. T. Gillespie. A general method for numerically simulating the stochastic time evolution of
coupled chemical reactions. Journal of Computational Physics, 22(4):403–434, 1976.
[6] D. T. Gillespie. A rigorous derivation of the chemical master equation. Physica A, 188(1–
3):404–425, 1992.
[7] S. Guo, Y. Hori, and R. M. Murray. Systematic design and implementation of a novel synthetic
fold-change detector biocircuit in vivo. Technical report, California Institute of Technology,
2014.
[8] Y. Hori, T.-H. Kim, and S. Hara. Existence criteria of periodic oscillations in cyclic gene
regulatory networks. Automatica, 47(5):1203–1209, 2011.
[9] Y. Hori, M. Takada, and S. Hara. Biochemical oscillations in delayed negative cyclic feedback:
existence and profiles. Automatica, 49(9):2581–2590, 2013.
[10] V. Hsiao, Y. Hori, P. W. K. Rothemund, and R. M. Murray. A population-based temporal
logic gate for timing and recording chemical events. Molecular Systems Biology, 12(869), 2016.
[11] H. J. Landau. Moments in Mathematics. American Mathematical Society, 1987.
[12] A. Milias-Argeitis, M. Rullan, S. K. Aoki, P. Buchmann, and M. Khammash. Automated
optogenetic feedback control for precise and robust regulation of gene expression and cell
growth. Nature Communications, 7(12):1114–1116, 2011.
[13] 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 expression circuit. Nature Biotechnology, 29(12):1114–1116, 2011.
[14] T.-S. Moon, C. Lou, A. Tamsir, B. C. Stanton, and C. A. Voigt. Genetic programs constructed
from layered logic gates in single cells. Nature, 491(7423):249–253, 2012.
[15] B. Munsky and M. Khammash. The finite state projection algorithm for the solution of the
chemical master equation. Journal of Chemical Physics, 124(4):044104, 2006.
[16] 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, 4:e09771, 2015.
[17] J. A. Shohat and J. D. Tamarkin. The problem of moments. American Mathematical Society,
1943.
[18] A. Singh and J. P. Hespanha. Approximate moment dynamics for chemically reacting systems.
IEEE Transactions on Automatic Control, 56(2):414–418, 2011.

[19] L. Yang, A. A. K. Nielsen, J. Fernandez-Rodriguez, C. J. McClune, M. T. Laub, T. K. Lu, and
C. A. Voigt. Permanent genetic memory with > 1 byte capacity. Nature Method, 11(12):1261–
1266, 2014.
[20] Y.-B. Zhao, J. Kim, and J. P. Hespanha. Hybrid moment computation algorithm for biochemical reaction networks. In Proceedings of IEEE Conference on Decision and Control, pages
1693–1698, 2010.

A

Proof of Theorem 2.

Proof The optimization problem (23) is equivalent to the following problem.
max
subject to

V
f (mα ) ≥ V
(21).

As the objective function is f (mα ) = m2ei − m2ei , the inequality f (mα ) ≥ V in the constraint can
be equivalently written by


m2ei − V mei
O
(31)
1
mei
using Schur compliment [1]. The optimal value of (23) is thus equal to that of (26).



