Variance reduction for antithetic integral control of stochastic
reaction networks

arXiv:1711.08291v1 [math.OC] 22 Nov 2017

Corentin Briat∗, Ankit Gupta∗, Mustafa Khammash†
Department of Biosystems Science and Engineering,
ETH-Zürich, Switzerland

Abstract
The antithetic integral feedback motif recently introduced in [6] is known to ensure robust perfect adaptation
for the mean dynamics of a given molecular species involved in a complex stochastic biomolecular reaction
network. However, it was observed that it also leads to a higher variance in the controlled network than that
obtained when using a constitutive (i.e. open-loop) control strategy. This was interpreted as the cost of the
adaptation property and may be viewed as a performance deterioration for the overall controlled network. To
decrease this variance and improve the performance, we propose to combine the antithetic integral feedback
motif with a negative feedback strategy. Both theoretical and numerical results are obtained. The theoretical
ones are based on a tailored moment closure method allowing one to obtain approximate expressions for the
stationary variance for the controlled network and predict that the variance can indeed be decreased by
increasing the strength of the negative feedback. Numerical results verify the accuracy of this approximation
and show that the controlled species variance can indeed be decreased, sometimes below its constitutive level.
Three molecular networks are considered in order to verify the wide applicability of two types of negative
feedback strategies. The main conclusion is that there is a trade-off between the speed of the settling-time of
the mean trajectories and the stationary variance of the controlled species; i.e. smaller variance is associated
with larger settling-time.

Author summary
Homeostasis, the ability of living organisms to regulate their internal state, is of fundamental importance
for their adaptation to environmental changes and their survival. This is the reason why complex regulatory
genetic networks evolved and allowed for the emergence of more and more complex organisms. Recently, the
theoretical study of those regulatory networks using ideas and concepts from control theory and the design of
novel ones have gained a lot of attention. Synthetic regulatory circuits are seen as elementary building blocks
for the design of complex networks that need to incorporate some regulating elements to be fully functional.
This is for instance the case in metabolic engineering where the production of biomolecules, such as drugs
or biofuels, need to be optimized and tightly regulated. A particular circuit, the so-called antithetic integral
controller, is now known to ensure homeostasis even when regulatory circuits are subject to randomness.
However, it is also known that this circuit increases variability in the network. The effects of a correcting
negative feedback loop on the variance are discussed here and it is shown that variability can be reduced
this way. Notably, we show that there is a tradeoff between speed of the network and variability.
∗ These

two coauthors contributed equally
author: mustafa.khammash@bsse.ethz.ch

† Corresponding

1

Introduction
The design and implementation of artificial in-vivo biomolecular controllers have become very popular [4, 6,
9,11,20,25] because of their potential applications for the tight and robust control of gene expression [6], the
optimization of metabolic networks for the efficient production of biomolecules [10, 30] , or the development
of new treatments for certain genetic diseases [31]. Indeed, many of the instances of those problems can
be interpreted from an homeostatic point of view in the sense that they may all be solved by achieving or
restoring homeostasis in the corresponding genetic network using synthetic regulatory circuits [6,10,27,30,31].
In this regard, those problems essentially reduce to the design and the implementation of robust and reliable
regulatory circuits that can optimize an inefficient network or correct a malfunctioning one – an observation
which strongly suggests that ideas from control theory and control engineering [2] could be adapted to
biochemical control problems [6, 12, 16]. A cornerstone in control theory and engineering is the so-called
integral controller that can ensure precise constant set-point regulation for a regulated variable in a given
system. Such mechanism, where the action onto the controlled system is depending on the integral of
the deviation of the regulated variable from the desired set-point, is to be contrasted with the so-called
proportional controller where the system is simply actuated proportionally to the deviation of the regulated
variable from the desired set-point. Unlike integral control, the latter one is unable to achieve robust constant
set-point regulation for the controlled variable and to reject constant disturbances. In other words, integral
control has the capacity of ensuring perfect adaptation for the regulated variable. The downside, however,
is that it has a destabilizing effect on the dynamics (emergence of oscillations or even diverging trajectories)
of the overall controlled system which can be then compensated by adjoining a proportional action, thus
giving rise to the so-called Proportional-Integral (PI) controller [1].
Based on the strength of these facts, an integral controller referred to as the antithetic integral controller
was proposed in [6] for the control of the mean level of some molecular species of interest in a given biochemical reaction network. This controller is implementable in terms of elementary biochemical reactions with
mass-action kinetics, making it practically implementable in-vivo using, for instance, sigma- and anti-sigmafactors [20]. This controller theoretically works in both the deterministic and the stochastic settings. In the
latter setting, it was notably shown that, under some reasonable conditions, the ergodicity properties of the
controlled network are independent of the parameters of the antithetic integral controller – a surprising key
property that has no counterpart in the deterministic setting and that dramatically simplifies its implementation. A drawback, however, is the increase of the stationary variance of the regulated species compared
to the constitutive variance that would be obtained by using a static open-loop strategy, even though the
latter one would be unable to ensure regulation and perfect adaptation for the mean level of the regulated
species. This phenomenon is seemingly analogous to the destabilizing behavior of the deterministic integral
controller mentioned in the previous paragraph. This variance increase can hence be interpreted as the price
to pay for perfect adaptation at the mean species level.
The goal of this paper is to investigate the effect of adding a negative feedback to the antithetic integral
motif in a way akin, yet different, to deterministic PI controllers. As discussed above, adding a proportional
action in the deterministic setting compensates for the destabilizing effect of the integrator. Comparatively,
it may seem reasonable to think that, in the stochastic setting, a proportional action could have an analogous
effect and would result in a decreased variance for the controlled variable (this is, for instance, what happens
when considering certain linear systems driven by white noise). In fact, it has been shown that negative
feedback at a transcriptional level in a gene expression network leads to a variance reduction in the protein
levels; see e.g. [5, 18, 24, 29] and the references therein.
Two types of negative feedback are considered in the present paper: the first one consists of an ON/OFF
proportional action whereas the second one is governed by a repressing Hill function. First we theoretically
prove using a tailored moment closure method that, in a gene expression network controlled with an antithetic
integral controller, the stationary variance in the protein copy number can be decreased by the use of a
negative feedback. In this specific case, the steady-state variance is decreasing monotonically as a function
of the strength of the negative feedback. An immediate consequence is that it is theoretically possible to
set the steady-state variance to a level that lies below the constitutive steady-state variance, which is the
value of the steady-state variance that would have been obtained using a constitutive (i.e. open-loop) control
2

strategy. The theoretical prediction will also be observed by exact numerical predictions using Gillespie’s
algorithm (Stochastic Simulation Algorithm - SSA [13]). A caveat, however, is that setting the gain of the
negative feedback very high will likely result in a very low steady-state variance but may also result in a
regulation error for the mean of the controlled species and in a loss of ergodicity for the overall controlled
network. In this regard, reducing the steady-state variance below its constitutive level may not always be
physically possible. Finally, it is also emphasized that a low stationary variance for the controlled species
is often associated with higher settling-time for the controlled species. Hence, there is a tradeoff between
variability and fast dynamics/small settling-time. The two negative feedback actions also exhibit quite
different behaviors. Indeed, while the ON/OFF proportional feedback seems to be efficient at reducing the
stationary variance through an increase of its gain, the dynamics of the mean gets first improved by reducing
the settling-time but then gets dramatically deteriorated by the appearance of a fast initial transient phase
followed by a very slow final one resulting then in a high settling-time. On the other hand, the Hill controller
leads to very homogeneous mean dynamics for different feedback strength but the steady-state variance is
also much less sensitive and does not vary dramatically. It is then argued that those differences may find
an explanation by reasoning in a deterministic point of view. The ON/OFF controller (an error-feedback)
introduces a stable zero in the dynamics of the closed-loop network which is small in magnitude when the
gain of the negative feedback is high. When this is the case, the zero is close to the origin and the closed-loop
dynamics will almost contain a derivative action, whence the fast initial transient phase. On the other hand,
the Hill negative feedback (an output-feedback) does not introduce such a zero in the closed-loop dynamics,
which may explain the homogeneity of the mean trajectories. Another possible reason is that the effective
proportional gain (which will be denoted by β) is much less sensitive to changes in the feedback strength
than the ON/OFF controller.
Approximate equations for the stationary variance are then obtained in the general unimolecular network
case. The obtained expressions shed some light on an interesting connection between the covariances of the
molecular species involved in the stochastic reaction network and the stability of a deterministic linear system
controlled with a standard PI controller, thereby unveiling an unexpected, yet coherent, bridge between
the stochastic and deterministic settings. Applying this more general framework to the a gene expression
network with protein maturation allows one to reveal that the steady-state variance may not be necessarily
a monotonically decreasing function of the negative feedback strength. In spite of this, the same conclusions
as in the gene expression network hold: the variance can sometimes be decreased below its constitutive
level but this may also be accompanied with a loss of ergodicity. The same qualitative conclusions for the
transient of the mean dynamics and the properties of the controller also hold in this case.
Even though the proposed theory only applies to unimolecular networks, stochastic simulations are
performed for a gene expression network with protein dimerization; a bimolecular network. Once again,
the same conclusions as in for previous networks hold with the difference that the constitutive variance level
is unknown in this case due to openness of the moment equations. These results tend to suggest that negative
feedback seems to operate in the same way in bimolecular networks as in unimolecular networks.

Reaction networks
Let us consider a stochastic reaction network (X, R) involving d molecular species X1 , . . . , Xd interacting
through K reaction channels R1 , . . . , RK defined as
Rk :

d
X
i=1

ρk

l
ζk,i
Xi −−−→

d
X

r
ζk,i
Xi , k = 1, . . . , K

(1)

i=1

r
r
l
l
where ρk ∈ R>0 is the reaction rate parameter and ζkr = col(ζk,1
, . . . , ζk,d
), ζkl = col(ζk,1
, . . . , ζk,d
) are
the left/right stoichiometric vectors of the reaction Rk . The corresponding stoichiometric vector is hence
given by ζk := ζkr − ζkl ∈ Zd indicating that when this reaction fires, the state jumps
from x to x + ζk .

The stoichiometric matrix S ∈ Zd×K of this reaction network is defined as S := ζ1 · · · ζK . When the
Qd
xi !
kinetics is mass-action, the propensity function λk of the reaction Rk is given by λk (x) = ρk i=1 (xi −ζ
l )! .
n,i

3

Under the well-mixed assumption, the above network can be described by a continuous-time Markov process
(X1 (t), . . . , Xd (t))t≥0 with the d-dimensional nonnegative lattice Zd≥0 as state-space; see e.g. [3].

The regulation/perfect adaptation problems and antithetic integral control
Let us consider here a stochastic reaction network (X, R). The regulation problem consists of finding another
reaction network (i.e. a set of additional species and additional reactions) interacting with (X, R) in a way
that makes the interconnection well-behaved (i.e. ergodic) and such that the mean of some molecular species
X` for some given ` ∈ {1, . . . , d} converges to a desired set-point (given here by µ/θ for some µ, θ > 0) in a
robust way; i.e. irrespective of the values of the parameters of the network (X, R).
It was shown in [6] that, under some assumptions on the network (X, R), the antithetic integral controller
defined as
η
µ
θX`
kZ1
X1 ,
Z2 , Z1 + Z2 −−−→ ∅, ∅ −−−→
(2)
∅ −−−→ Z1 , ∅ −−−→
{z
} |
{z
} |
{z
} |
{z
}
|
comparison
actuation
reference measurement
solves the above regulation problem. This regulatory network consists of two additional species Z1 and Z2 ,
and four additional reactions. The species Z1 is referred to as the actuating species as it is the species that
governs the rate of the actuation reaction which produces the actuated species X1 at a rate proportional to
Z1 . The species Z2 is the sensing species as it is produced at a rate proportional to the controlled species
X` through the measurement reaction. The first reaction is the reference reaction as it encodes part of the
set-point µ/θ whereas the third reaction is the comparison reaction that compares the population of the
controller species and annihilates them accordingly, thereby closing negatively the loop while, and at the
same time correlating the populations of the controller species. The comparison (or titration) reaction is the
crucial element of the above controller network and, to realize such a reaction, one needs to rely on intrinsic
strongly binding properties of certain molecules such as sigma- and anti-sigma-factors [6] or small RNAs and
RNAs [19, 25, 32].

Variance amplification in antithetic integral control
We discussed above about the convergence properties of the mean level of the controlled species X` when
network (X, R) is controlled with the antithetic integral controller (2). However, it was remarked in [6]
that while the mean of X` converges to the desired steady-state, the stationary variance of the controlled
species could be much larger than its constitutive value that would be obtained by simply considering a naive
constitutive production of the species X1 that would lead to the same mean steady-state value µ/θ. This
was interpreted as the price to pay for having the perfect adaptation property for the controlled species. To
illustrate this phenomenon, let us consider the following gene expression network:
k

kp

γr

γp

r
∅ −−−→
X1 , X1 −−−→ X1 + X2 , X1 −−−→ ∅, X2 −−−→ ∅

(3)

where X1 denotes mRNA and X2 denotes protein. The objective here is to control the mean level of the
protein by acting at a transcriptional level using the antithetic controller (2); hence, we set kr = kZ1 . Using
a tailored moment closure method, it is proved in the SI that the stationary variance VarIπ (X2 ) for the
protein copy number is approximately given by the following expression


kkp
kp
1+
+
µ
γr + γp
γr γp 
 , k > 0, k/η  1.
VarIπ (X2 ) ≈ 
(4)


kθkp
θ
1−
γr γp (γr + γp )
The rationale for the assumption k/η  1 is that it allows for closing the moments equation (which is open
because of the presence of the comparison reaction) and obtain a closed-form solution for the stationary

4

variance. On the other hand, the constitutive (i.e. open-loop) stationary variance VarOL
π (X2 ) for the protein
copy number obtained with the constitutive strategy
kr =

µ γr γp
θ kp

(5)

is given by
VarOL
π (X2 ) =

µ
θ


1+

kp
γr + γp


.

(6)

It is immediate to see that the ratio
VarIπ (X2 )
≈
VarOL
π (X2 )

kkp
γr γp (γr + γp )
, k, θ > 0, k/η  1
kθkp (kp + γr + γp )
1−
γr γp (γr + γp )
1+

(7)

is greater than 1 for all k, θ > 0 such that the denominator is positive. Note that the above formula is not
valid when k = 0 or θ = 0 since this would result in an open-loop network for which set-point regulation
could not be achieved. This expression is also a monotonically increasing function of the gain k, a fact that
was numerically observed in [6]. This means that choosing k very small will only result in a small increase
of the stationary variance of the controlled species when using an antithetic integral feedback. However, this
will very likely result in very slow dynamics for the mean of the controlled species.
Finally, it is important to stress that while this formula is obviously not valid when the denominator is
nonpositive, we know from [6] that in the case of the gene expression network, the closed-loop network will
be ergodic with converging first and second-order moments for all k > 0 and all θ > 0 (assuming that the
ratio µ/θ is kept constant). This inconsistency stems from the fact that the proposed theoretical approach
relies on a tailored moment closure approximation that will turn out to be connected to the Hurwitz stability
of a certain matrix that may become unstable when the gain k of the integrator is too large. This will be
elaborated more in the following sections.

Negative feedback action
We will consider in this paper two types of negative feedback action. The first one, referred to as the
ON/OFF proportional feedback, is essentially theoretical and cannot be exactly implemented, but it may be
seen as a local approximation of some more complex (e.g. nonlinear) repressing function. It is given by the
reaction
F (X` )
∅ −−−→ X1
(8)
together with the propensity function F (X` ) = Kp max{0, µ − θX` } where Kp is the so-called feedback
gain/strength. It is similar to the standard proportional feedback action used in control theory with the
difference that a regularizing function, in the form of a max function, is involved in order the restrict the
propensity function to nonnegative values. Note that this controller can still be employed for the in-silico
control of single-cells using a stochastic controller as, in this case, we would not be restricted anymore to
mass-action, Hill or Michaelis-Menten kinetics. This was notably considered in the case of in-silico population
control in [7, 8, 14].
The second type of negative feedback action, referred to as the Hill feedback, consists of the reaction (8)
but involves the non-cooperative repressing Hill function F (X` ) = Kp /(1 + X` ) as propensity function. This
type of negative feedback is more realistic as such functions have empirically been shown to arise in many
biochemical, physiological and epidemiological models; see e.g. [22].
In both cases, the total rate of production of the molecular species X1 can be expressed as the sum
kZ1 + F (X` ) which means that, at stationarity, we need to have that Eπ [kZ1 + F (X` )] = u∗ where u∗ is
equal to the value of the constitutive (i.e. deterministic) production rate for X1 for which we would have
that Eπ [X` ] = µ/θ. Noting now that for both negative feedback functions, we will necessarily have that
5

Eπ [F (X` )] > 0, then this means that if the gain Kp is too large, it may be possible that the mean of the
controlled species do not converge to the desired set-point implying, in turn, that the overall controlled
network will fail to be ergodic. This will be notably the case when Eπ [F (X` )] > u∗ . In particular, when
F (X` ) = Kp max{0, µ − θX` }, a very conservative sufficient condition for the closed-loop network to be
ergodic is that Kp < u∗ /µ whereas when F (X` ) = Kp /(1 + X` ), this condition becomes Kp < u∗ . These
conditions can be determined by considering the worst-case mean value of the negative feedback strategies;
i.e. Kp µ and Kp , respectively.

Results
Invariants for the antithetic integral controller
We describe some important invariant properties of the antithetic integral controller (2) which are independent of the parameters of the controlled network under the assumption that these invariants exist; i.e. they
are finite. Those invariants are given by
Covπ (X` , Z1 − Z2 ) =

µ
,
θ

(9)

µ
,
η

(10)

Eπ (Z12 Z2 ) =

µ
(1 + Eπ (Z1 ))
η

(11)

Eπ (Z1 Z22 ) =

µ + θEπ (X` Z2 )
η

(12)

Eπ (Z1 Z2 ) =

and

play an instrumental role in proving all the theoretical results of the paper. Interestingly, we can notice that
Covπ (X` , Z1 − Z2 ) = Eπ [X` ], which seems rather coincidental. From the second invariant we can observe
that, if η  µ, then Eπ (Z1 Z2 ) ≈ 0, which indicates that the values taken by the random variable Z2 (t)
will be most of the time equal to 0. Note that it cannot be Z1 (t) to be mostly taking zero values since Z1
is the actuating species whose mean must be nonzero (assuming here that the natural production rates of
the molecular species in the controlled network are small). Similarly, setting η large enough in the third
expression will lead to a similar conclusion. Note that Eπ (Z1 ) is independent of η here and only depends
on the set-point µ/θ, the integrator gain k and the parameters of the network which is controlled. The
last expression again leads to similar conclusions. Indeed, if η is sufficiently large, then Eπ (X` Z2 ) ≈ 0 and,
hence, Eπ (Z1 Z22 ) ≈ 0 which implies that Z2 (t) needs to be most of the time equal to 0. These properties
will be at the core of the moment closure method used to obtain an approximate closed-form formula for
the covariance matrix for the closed-loop network.

An approximate formula for the stationary variance of the controlled species
Let us assume here that the open-loop network (X, R) is mass-action and involves, at most, unimolecular
reactions. Hence, the vector of propensity functions can be written as
λ(x) = W x + w0

(13)

for some nonnegative matrix W ∈ RK×d and nonnegative vector w0 ∈ RK . It is proved in the SI that, under
the assumption k/η  1, we can overcome the moment closure problem arising from the presence of the

6

annihilation reaction in the antithetic controller and show that the exact stationary covariance matrix of the
network given by


I
PI
CovP
π (X, X) Covπ (X, Z)
, Z := Z1 − Z2
I
I
CovP
VarP
π (Z, X)
π (Z)
is approximatively given by the matrix Σ solving the Lyapunov equation
RΣ + ΣRT + Q = 0
where

(14)



SW − βe1 eT` ke1
,
−θeT`
0 

SDS T + ce1 eT1
0
Q =
,
0
2µ
D = diag(W Eπ [X] + w0 ),

µ
1
T
−1
+
e
(SW
)
Sw
,
c = − T
0
`
e` (SW )−1 e1 θ
I
CovP
π (F (X` ), X` )
β = −
.
I
CovP
π (X` , X` )
R

=

Note that since the function F is decreasing then the effective proportional gain, β, is always a positive
constant and seems to be mostly depending on Kp but does not seem to change much when k varies (see
e.g. Figure 19 and Figure 20 in the appendix). It can also be seen that for the Lyapunov equation to have a
positive definite solution, we need that the matrix R be Hurwitz stable; i.e. all its eigenvalues have negative
real part. In parallel of that, it is known from the results in [6] that the closed-loop network will remain
ergodic when β = 0 even when the matrix R is not Hurwitz stable. In this regard, the formula (14) can only
be valid when the parameters β and k are such that the matrix R is Hurwitz stable. When this is not the
case, the formula is out its domain of validity and is meaningless. The stability of the matrix R is discussed
in more details in the SI.

Connection to deterministic proportional-integral control
Interestingly, the matrix R coincides with the closed-loop system matrix of a deterministic linear system
controlled with a particular proportional-integral controller. To demonstrate this fact, let us consider the
following linear system
ẋ(t) = SW x(t) + e1 u(t)
(15)
y(t) = eT` x(t)
where x is the state of the system, u is the control input and y is the measured/controlled output. We
propose to use the following PI controller in order to robustly steers the output to a desired set-point µ/θ
β
u(t) = (µ − θy(t)) + k
θ

Z

t

(µ − θy(s))ds

(16)

0

where θ is the sensor gain, β/θ is the proportional gain and k is the integral gain. The closed-loop system
is given in this case by

 

 β 
ẋ(t)
SW − βe1 eT` ke1 x(t)
=
(17)
+ θ µ
˙
I(t)
−θeT`
0
I(t)
1
where we can immediately recognize the R matrix involved in the Lyapunov equation (14).

7

Example - Gene expression network
We present here the results obtained for the gene expression network (3) using the two negative feedback
actions. In particular, we will numerically verify the validity of the formula (4) and study the influence of
the controller parameters on various properties of the closed-loop network. The matrix R is given in this
case by


−γr −β k
R =  kp −γp 0 .
(18)
0
−θ 0
It can be shown that the above matrix is Hurwitz stable (i.e. all its eigenvalues are located in the open left
half-plane) if and only if the parameters k, β > 0 satisfy the inequality
1−

βkp
kθkp
+
> 0.
γr γp (γr + γp ) γr γp

(19)

Hence, given k > 0, the matrix R will be Hurwitz stable for any sufficiently large β > 0 illustrating the
stabilizing effect of the proportional action. When the above condition is met, then the closed-loop stationary
I
variance VarP
π (X2 ) of the protein copy number is approximately given by the expression


kp
kkp
βkp
1+
+
+
µ
γr + γp
γr γp
γr (γr + γp ) 
I

.
(20)
VarP
π (X2 ) ≈ Σ22 =


βkp
kθkp
θ
+
1−
γr γp (γr + γp ) γr γp
For any fixed k > 0 such that (19) is satisfied, the closed-loop steady-state variance is a monotonically
decreasing function of β. As a consequence, there will exist a βc > 0 such that


µ
kp
Σ22 <
1+
(21)
θ
γr + γp
for all β > βc . In particular, when β → ∞, then we have that
Σ22 →

µ
µ γp
< .
θ γr + γp
θ

(22)

We now analyze the results obtained with the antithetic integral controller combined with an OF/OFF
proportional feedback. The first step is the numerical comparison of the approximate formula (20) with
the stationary variance computed using 106 SSA simulations with the parameters kp = 2, γr = 2, γp = 7,
µ = 10, θ = 2 and η = 100. The absolute value of the relative error between the exact and the approximate
stationary variance of the protein copy number for several values for the gains k and Kp is depicted in Figure
1. We can observe there that the relative error is less than 15% except when k is very small where the relative
error is much larger. However, in this latter case, the mean trajectories do not have time to converge to
their steady state value and, therefore, what is depicted in the figure for this value is not very meaningful.
In spite of that, we can observe that the approximation is reasonably accurate.
We now look at the performance of the antithetic integral controller combined with an OF/OFF proportional feedback. Figure 2 depicts the trajectories of the mean protein copy number while Figure 3 depicts
the trajectories of the variance of the protein copy number, both in the case where k = 3. Regarding the
mean copy number, we can observe that while at the beginning increasing Kp seems to improve the transient
phase, then the dynamics gets more and more abrupt at the start of the transient phase as the gain Kp
continues to increase and gets slower and slower at the end of the transient phase, making the means very
slow to converge to their set-point. On the other hand, we can see that the stationary variance seems to
be a decreasing function of the gain Kp . More interestingly, when the gain Kp exceeds 20, the stationary
variance becomes smaller than its constitutive value. Figure 4 helps at establishing the influence of the gains
k and Kp onto the stationary variance of the protein copy number. We can see that, for any k, increasing
8

Kp reduces the stationary variance while for any Kp , reducing k reduces the variance, as predicted by the
approximate formula (20). Hence, a suitable choice would be to pick k small and Kp large. We now compare
this choice for the parameters with the one that would lead to a small settling-time for the mean dynamics;
see Figure 5. We immediately see that a small k is not an option if one wants to have fast mean dynamics.
A sweet spot in this case would be around the right-bottom corner where the settling-time is the smallest.
Interestingly, the variance is still at a quite low level even if sometimes higher than the constitutive value.
We now perform the same analysis for the antithetic integral controller combined with the Hill feedback
and first verify the accuracy of the approximate formula (20). We can observe in Figure 6 that the formula
is very accurate in this case. To explain this, it is important to note that the gains Kp in both controllers
are not directly comparable, only the values for the parameter β are. For identical Kp ’s, the value of β for
the ON/OFF proportional feedback is much larger than for the Hill feedback (see Figure 19 and Figure 20
in the appendix). The Figure 1 and Figure 6 all together simply say that the formula is very accurate when
β is small.
We now look at the performance of the antithetic integral controller combined with a Hill feedback.
Similarly to as previously, Figure 7 depicts the trajectories of the mean protein copy number while Figure 8
depicts the trajectories of the variance of the protein copy number, both in the case where k = 3. Regarding
the mean copy number, we can observe than the dynamics are much more homogeneous than in the previous
case and that increasing Kp reduces the overshoot and, hence, the settling-time. This can again be explained
by the fact that β is much smaller in this case. Similarly, the spread of the variances is much tighter than when
using the other negative feedback again because of the fact that β is small in this case. This homogeneity
is well illustrated in Figure 9 and Figure 10 where we conclude on the existence of a clear tradeoff between
settling-time and stationary variance.
As can been seen in Figure 2 and Figure 7, the mean dynamics are quite different and it would be
interesting to explain this difference in terms of control theoretic ideas. A first explanation lies in the
sensitivity of the parameter β in terms of the feedback strength Kp . In the case of the ON/OFF proportional
feedback, this sensitivity is quite high whereas it is very low in the case of the Hill feedback (see Figure 19
and Figure 20 in the appendix). This gives an explanation on why the mean trajectories are very different in
the case of the ON/OFF proportional feedback for different values of Kp while the mean trajectories are very
close to each other in the case of the Hill feedback. A second explanation lies in the type of feedback in use.
Indeed, the ON/OFF proportional feedback is an error-feedback and, when combined with the antithetic
integral controller, may introduce a stable zero in the mean dynamics. On the other hand, the Hill feedback
is an output-feedback that does not seem to introduce such a zero. When increasing the negative feedback
gain Kp , this zero moves towards the origin. Once very close to the origin, this zero will have an action in the
closed-loop mean dynamics that is very close to a derivative action, leading then to abrupt initial transient
dynamics. A theoretical basis for this discussion is developed in more details in the SI.

9

Figure 1: Absolute value of the relative error between the exact stationary variance of the protein copy
number and the approximate formula (20) when the gene expression network is controlled with the antithetic
integral controller (2) and an ON/OFF proportional controller.

10

Figure 2: Mean trajectories for the protein copy number when the gene expression network is controlled
with the antithetic integral controller (2) with k = 3 and an ON/OFF proportional controller. The set-point
value is indicated as a black dotted line.

11

Figure 3: Variance trajectories for the protein copy number when the gene expression network is controlled
with the antithetic integral controller (2) with k = 3 and an ON/OFF proportional controller. The stationary
constitutive variance is depicted in black dotted line.

12

Figure 4: Stationary variance for the protein copy number when the gene expression network is controlled
with the antithetic integral controller (2) and an ON/OFF proportional controller.

13

Figure 5: Settling-time for the mean trajectories for the protein copy number when the gene expression
network is controlled with the antithetic integral controller (2) and an ON/OFF proportional controller.

14

Figure 6: Absolute value of the relative error between the exact stationary variance of the protein copy
number and the approximate formula (20) when the gene expression network is controlled with the antithetic
integral controller (2) and a Hill controller.

15

Figure 7: Mean trajectories for the protein copy number when the gene expression network is controlled
with the antithetic integral controller (2) with k = 3 and a Hill controller. The set-point value is indicated
as a black dotted line.

16

Figure 8: Variance trajectories for the protein copy number when the gene expression network is controlled
with the antithetic integral controller (2) with k = 3 and a Hill controller. The stationary constitutive
variance is depicted in black dotted line.

17

Figure 9: Stationary variance for the protein copy number when the gene expression network is controlled
with the antithetic integral controller (2) and a Hill controller.

18

Figure 10: Settling-time for the mean trajectories for the protein copy number when the gene expression
network is controlled with the antithetic integral controller (2) and a Hill controller.

Example - Gene expression network with protein maturation
The results obtained in the previous section clearly only hold for the gene expression network and it would be
quite hasty to directly generalize those results to more complex unimolecular networks. This hence motivates
the consideration of a slightly more complicated example, namely, the gene expression network involving a
protein maturation reaction given by
k

kp

γr

γp

r
∅ −−−→
X1 , X1 −−−→ X1 + X2 , X1 −−−→ ∅, X2 −−−→ ∅
0
kp

γp0

(23)

X2 −−−→ X3 , X3 −−−→ ∅
where, as before, X1 denotes mRNA, X2 denotes protein and, now, X3 denotes the mature protein. In this
case, the goal is to control the average mature protein copy number by, again, acting at a transcriptional
level. As this network is still unimolecular, the proposed framework remains valid. In particular, the matrix
R is given by


−γr
0
−β k
 kp −(γp + kp0 )
0
0

R=
(24)
0
0
 0
kp
−γp 0
0
0
−θ 0
and is Hurwitz stable provided that the two following conditions are satisfied
β<


1
(γr + γp + γp0 + kp0 )(γr γp + γr γp . + γp γp0 + γr kp0 + γp0 kp0 ) − γr γp0 (γp + kp0 )
kp kp0

(25)

and
kp2 kp02 β 2 + σ1 β + σ0 < 0
19

(26)

where

= −kp kp0 (γr + γp0 + γp + kp0 )(γr γp + γr γp0 + γp γp0 + γr kp0 + γp0 kp0 )
+2γr γp0 kp kp0 (γp + kp0 ),
= −γr γp0 (γp + kp0 )(γr + γp + γp0 + kp0 )(γr γp + γr γp0 + γp γp0 + γr kp0 + γp0 kp0 )
+γr2 γp 22 (γp + kp0 )2 + kkp kp0 θ(γr + γp + γp0 + kp0 )2 .

σ1
σ0

(27)

Considering, for instance, the following parameters kp = 1, γr = 2, γp = 1, kp0 = 3, γp0 = 1, µ = 10, θ = 2
and η = 100, the above conditions reduce to
β < 30
(28)
and
9β 2 − 246β + 294k − 720 < 0.
The intersection of these conditions yield the stability conditions
√
√


41 − 7 49 − 6k 41 + 7 49 − 6k
,
∩ (0, ∞).
k ∈ (0, 49/6) and β ∈
3
3

(29)

(30)

It can be verified that for values on the boundary of at least one of those intervals, the matrix R has
eigenvalues on the imaginary axis. Standard calculations on the moments equation show that the open-loop
variance is given by


kp0 + γr + γp + γp0
µ
0
VarOL
(X
)
=
1
+
k
k
.
(31)
3
p p
π
θ
(γr + γp0 )(γr + γp + kp0 )(γp + γp0 + kp0 )
With the numerical values for the parameters previously given, the open-loop variance is approximately
equal to 37/6 ≈ 6.1667. The closed-loop variance, however, is approximately given by


ζk
θ
ζβ
ζkβ
OL
(X
)
+
Var
k
+
β
+
kβ
3
π

µ
ζd
ζd
ζd
I
µ

VarP
(32)
π (X3 ) ≈ Σ33 =


ξβ 2 2
ξk
ξβ
θ
1+ k+ β+
β
ξd
ξd
ξd
where

= γr γp0 (γr + γp0 )(γp + kp0 )(γr + γp + kp0 )(γp + γp0 + kp0 )
= −kp kp0 θ(γr + γp + γp0 + kp0 )2
= kp kp0 (γr2 γp + γr2 γp0 + γr2 kp0 + γr γp2 + γr γp γp0 + 2γr γp kp0 + γr γp02
+γr γp0 kp0 + γr kp02 + γp2 γp0 + γp γp02 + 2γp γp0 kp0 + γp02 kp0 + γp0 kp02 )
= −kp2 kp02

(33)

ξd
kp kp0 (γr2 γp + γr2 γp0 + γr2 kp0 + γr γp2 + 2γr γp γp0 + 2γr γp kp0 + γr γp02
+2γr γp0 kp0 − θγr γp0 + γr kp02 + γp2 γp0 + γp γp02 + 2γp γp0 kp0 − θγp γp0
+γp02 kp0 − θγp02 + γp0 kp02 − θγp0 kp0 )
= γp0 kp kp0 (γr2 + γr γp + γr kp0 + γp0 γr + γp2 + 2γp kp0 + γp0 γp + kp02 + γp0 kp0 )
= −kp2 kp02 .

(34)

ξd
ξk
ξβ
ξβ 2
and
ζd
ζk

ζβ
ζkβ

=
=

An expression that is more complex than, yet very similar to, the formula (20) obtained for the simple gene
expression network. For the considered set of parameter values, the approximated variance is a nonmonotonic
function of the parameter β as it can be theoretically observed in Figure 21 in the appendix. It turns out
that this behavior can also be observed in the numerical simulations depicted in Figure 23 in the appendix
where we can see that the variance exhibits this monotonic behavior. However, it should also be pointed
out that this increase is accompanied with the emergence of an tracking error for the mean dynamics (see
Figure 23 in the appendix) and a loss of ergodicity for the overall controlled network as emphasized by
diverging mean dynamics for the sensing species (see Figure 25 in the appendix). This contrasts with the
20

gene expression case where the variance was a monotonically decreasing function of β. Regarding the mean
dynamics, we can see that increasing Kp and, hence, β, to reasonable levels improves the settling-time as
depicted in Figure 11 for the special case of k = 3. However, this is far from being the general case since
the settling-time can exhibit a quite complex behavior for this network (see 14). The stationary variance
depicted in Figure 13 exhibits here a rather standard and predictive behavior where a small k and a large
Kp both lead to its reduction. Similar conclusions can be drawn when the network is controlled with a Hill
negative feedback controller; see Figure 26, Figure 26, Figure 26 and Figure 26 in the appendix.

Figure 11: Mean trajectories for the mature protein copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF
proportional controller. The set-point value is indicated as a black dotted line.

21

Figure 12: Variance trajectories for the mature protein copy number when the gene expression network
with protein maturation is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF
proportional controller. The stationary constitutive variance is depicted in black dotted line.

22

Figure 13: Stationary variance for the mature protein copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) and an ON/OFF proportional
controller.

23

Figure 14: Settling-time for the mean trajectories for the mature protein copy number when the gene
expression network with protein maturation is controlled with the antithetic integral controller (2) and an
ON/OFF proportional controller.

Example - Gene expression network with protein dimerization
The proposed theory is only valid for unimolecular networks but, in spite of that, it is still interesting to see
whether similar conclusions could be obtained for a network that is not unimolecular. This motivates the
consideration of the following gene expression network with protein dimerization:
kp

k

γp

γr

r
∅ −−−→
X1 , X1 −−−→ X1 + X2 , X1 −−−→ ∅, X2 −−−→ ∅

k

γd

γ0

(35)

d
d
X2 + X2 −−−→
X3 , X3 −−−→ X2 + X2 , X3 −−−→
∅

where, as before, X1 denotes mRNA, X2 denotes protein but, now, X3 denotes a protein homodimer. In this
case, the Lyapunov equation (14) is not valid anymore because of the presence of the dimerization reaction
but we can still perform stochastic simulations. The considered parameter values are given by kp = 1, γr = 2,
γp = 1, kd = 3, γd = γd0 = 1, µ = 10, θ = 2 and η = 100. We can see in Figure 15, Figure 16, Figure 17,
Figure 18.

24

Figure 15: Mean trajectories for the homodimer copy number when the gene expression network with protein
dimerization is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF proportional
controller. The set-point value is indicated as a black dotted line.

25

Figure 16: Variance trajectories for the homodimer copy number when the gene expression network with
protein dimerization is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF
proportional controller. The stationary constitutive variance is depicted in black dotted line.

26

Figure 17: Stationary variance for the homodimer copy number when the gene expression network with
protein dimerization is controlled with the antithetic integral controller (2) and an ON/OFF proportional
controller.

27

Figure 18: Settling-time for the mean trajectories for the homodimer copy number when the gene expression
network with protein dimerization is controlled with the antithetic integral controller (2) and an ON/OFF
proportional controller.

Discussion
Adjoining a negative feedback strategy to the antithetic integral controller was shown to reduce the stationary
variance for the controlled species, an effect that was expected from previous studies and predicted by the
obtained theoretical results. The structure of the negative feedback strategy was notably emphasized to have
important consequences on the magnitude of the variance reduction. Indeed, the ON/OFF controller can
be used to dramatically reduce the variance while still preserving the ergodicity of the closed-loop network.
This can be explained mainly because the proportional effective gain β is very sensitive to changes in the
feedback strength Kp and can reach reasonably large values (still smaller than Kp ); see Fig. 19 in the
appendix. The preservation of the ergodicity property for the closed-loop network comes from the fact that
Eπ [Kp max{0, µ−θX` }] remains smaller than the value of the nominal stationary control input (the constant
input for which the stationary mean of the controlled species equals the desired set-point) for a wide range
of values for Kp . Regarding the mean dynamics, this feedback leads to a decrease of the settling-time but
also also leads to abrupt transient dynamics for large values of Kp because of the presence of a stable zero
in the mean closed-loop dynamics that is inversely proportional to β (which is very sensitive to changes
in Kp in this case and which can reach high values). Unfortunately, this controller cannot be implemented
in-vivo because it does not admit any reaction network implementation. However, it can still be implemented
in-silico for the stochastic single-cell control for the control of cell populations using, for instance, targeted
optogenetics; see e.g. [26]. On the other hand, the Hill feedback, while being practically implementable, has
a much less dramatic impact on the stationary variance and on the mean dynamics. The first reason is that
the effective proportional gain β is less sensitive with respect to changes in Kp and remains very small even
when Kp is large; see Fig. 20. The absence of zero does not lead to any abrupt transient dynamics even for

28

large values for Kp but this may also be due to the fact that β always remains small as opposed to as in
the ON/OFF proportional feedback case. A serious issue with this feedback is that ergodicity can be easily
lost since Eπ [Kp /(1 + X` )] becomes very quickly larger than the value of the nominal control input as we
increase Kp . The properties of both feedback strategies are summarized in Table 1.
To prove the main theoretical results, a tailored closure method had to be developed to deal with the
bimolecular comparison reaction. A similar one has also been suggested in [23] for exactly the same purpose.
These methods rely on the assumption that the molecular count of the controller species Z2 is, most of the
time, equal to 0, a property that is ensured by assuming that k/η  1. This allowed for the simplification
and the closure of the moment equations. The theory was only developed for unimolecular networks because
of the problem solvability. However, the extension of those theoretical results to more general reaction
networks, such as bimolecular networks, is a difficult task mainly because of the moment closure problem
that is now also present at the level of the species of the controlled network. In this regard, this extension is,
at the moment, only possible using existing moment closure methods (see e.g. [17, 21, 28]) which are known
to be potentially very inaccurate and would then compromise the validity of the obtained approximation.
We believe that obtaining accurate and general theoretical approximations for the stationary variance for
bimolecular networks is currently out of reach. It is also unclear whether the obtained qualitative and
quantitative results still hold when the assumption k/η  1 on the controller parameters is not met.
Interestingly, the results obtained in the current paper provide some interesting insights on an unexpected
connection between deterministic PI control and its stochastic analogues. In particular, it is possible to
observe that the destabilizing effect of deterministic integral control is analogous to the variance increase
due to the use of the stochastic antithetic integral controller. In a similar way, the stabilizing property of
deterministic proportional controllers is the deterministic analogue of the property of variance decrease of
the stochastic proportional controller; see Table 2).
The controller considered in this paper is clearly analogous to PI controllers. A usual complemental
element is the so-called derivative action (or a filtered version of it) in order to add an anticipatory effect
to the controller and prevent high overshoot; see [1]. So far, filtered versions of the derivative action have
been proposed in a deterministic setting. Notably, the incoherent feedforward loop locally behaves like a
filtered derivative action. More recently, a reaction network approximating a filtered derivative action was
proposed in [15] in the deterministic setting. It is unclear at the moment whether a stochastic version for
the derivative action can be found but it is highly possible that such a stochastic derivative action can be
implemented in terms of elementary reactions.
The negative feedback strategy considered here is an ideal/simplified one. Indeed, it was assumed in this
paper that the controlled species was directly involved in the negative feedback. However, it is very likely
that, the controlled species may not be directly usable in the feedback, that intermediary species may be
involved (e.g. a gene expression network is involved in the feedback) or that the feedback is in terms of a
species upstream the controlled species (for instance feedback uses a protein while the controlled species is
the corresponding homodimer). The theory may be adapted to deal with such cases as long as the controlled
network is unimolecular. It is however expected that the same qualitative behavior will be observed. The
reason for that is that in unimolecular networks, species cooperate in the sense that they act positively on
each other. Hence, decreasing the variance of one species will also decrease the variance of all the species
that are created from it. For instance, in a gene expression network, if the mRNA variance is decreased, the
protein variance will decrease as well, and vice-versa.
Finally, the implementation of such negative feedback loops is an important, yet elusive, task. It is
unclear at the moment how in-vivo experiments could be conducted. Preliminary experimental results to
validate the theoretical/computational ones could be obtained using optogenetics and single-cell control for
population control. In-vivo experiments will certainly require a lot more effort.

29

Table 1: Effects of the different feedback strategies on the mean dynamics and the stationary variance.
ON/OFF Proportional Feedback
Hill Feedback
Ergodicity
robust (+)
fragile (-)
β
very sensitive (+)
poorly sensitive (-)
wide range (+)
small range
Mean Dynamics
reduce settling-time (+)
reduce settling-time (+)
zero dynamics (-)
no zero dynamics (+)
Stationary variance
dramatic reduction (++)
slight reduction (+)
Table 2: The effects of the proportional and integral actions on the dynamics of a system in both the
deterministic and stochastic setting.
Integral action
Proportional action
Deterministic
regulation (+)
no regulation (-)
Setting
destabilizing (-)
stabilizing (+)
Stochastic
regulation (+)
no regulation (-)
Setting
increases variance (-) decreases variance (+)

Supplementary figures for the gene expression network

Figure 19: Evolution of β as a function of the gains k and Kp in a gene network controlled with an antithetic
integral controller combined with an ON/OFF proportional feedback.

30

Figure 20: Evolution of β as a function of the gains k and Kp in a gene network controlled with an antithetic
integral controller combined with a Hill feedback.

31

Supplementary figures for the gene expression network with protein
maturation

4

log10 ('33 )

3.5
3
2.5
2
1.5
1
0.5
0

5

10

15

20

25

30

Figure 21: Evolution of the logarithm of the predicted stationary variance for the mature protein copy
number as a function of the gain β. We can observe a nonmonotonic behavior due to the fact that the
stability region for R in terms of the parameters k and β is bounded. The black dashed line is the logarithm
of the open-loop variance.

32

Figure 22: Mean trajectories for the mature protein copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF
proportional controller. The set-point value is indicated as a black dotted line.

Figure 23: Variance trajectories for the mature protein copy number when the gene expression network
with protein maturation is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF
proportional controller. The stationary constitutive variance is depicted in black dotted line.

33

Figure 24: Mean trajectories for the actuating species copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF
proportional controller.

Figure 25: Mean trajectories for the sensing species copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) with k = 3 and an ON/OFF
proportional controller.

34

Figure 26: Mean trajectories for the mature protein copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) with k = 3 and a Hill negative
feedback controller. The set-point value is indicated as a black dotted line.

Figure 27: Variance trajectories for the mature protein copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) with k = 3 and a Hill negative
feedback controller. The stationary constitutive variance is depicted in black dotted line.

35

Figure 28: Stationary variance for the mature protein copy number when the gene expression network with
protein maturation is controlled with the antithetic integral controller (2) and a Hill feedback controller.

Figure 29: Settling-time for the mean trajectories for the mature protein copy number when the gene
expression network with protein maturation is controlled with the antithetic integral controller (2) and a
Hill negative feedback controller.

36

Supplementary figures for the gene expression network with protein
dimerization

Figure 30: Mean trajectories for the homodimer copy number when the gene expression network with protein
dimerization is controlled with the antithetic integral controller (2) with k = 3 and a Hill negative feedback.
The set-point value is indicated as a black dotted line.

37

Figure 31: Variance trajectories for the homodimer copy number when the gene expression network with
protein dimerization is controlled with the antithetic integral controller (2) with k = 3 and a Hill negative
feedback. The stationary constitutive variance is depicted in black dotted line.

Figure 32: Stationary variance for the homodimer copy number when the gene expression network with
protein dimerization is controlled with the antithetic integral controller (2) and aa Hill negative feedback.

38

Figure 33: Settling-time for the mean trajectories for the homodimer copy number when the gene expression
network with protein dimerization is controlled with the antithetic integral controller (2) and a Hill negative
feedback.

References
[1] K. J. Åström and T. Hägglund. PID Controllers: Theory, Design, and Tuning. Instrument Society of
America, Research Triangle Park, North Carolina, USA, 1995.
[2] P. Albertos and I. Mareels. Feedback and Control for Everyone. Springer, Berlin Heidelberg, Germany,
2010.
[3] D. Anderson and T. G. Kurtz. Continuous time Markov chain models for chemical reaction networks. In
H. Koeppl, D. Densmore, G. Setti, and M. di Bernardo, editors, Design and analysis of biomolecular circuits - Engineering Approaches to Systems and Synthetic Biology, pages 3–42. Springer Science+Business
Media, 2011.
[4] F. Annunziata, A. Matyjaszkiewicz, G. Fiore, C. S. Grierson, L. Marucci, M. di Bernardo, and N. J.
Savery. An orthogonal multi-input integration system to control gene expression in escherichia coli (in
press). ACS Synthetic Biology, 2017.
[5] A. Becskei and L. Serrano. Engineering stability in gene networks by autoregulation. Nature, 405:590–
593, 2000.
[6] C. Briat, A. Gupta, and M. Khammash. Antithetic integral feedback ensures robust perfect adaptation
in noisy biomolecular networks. Cell Systems, 2:17–28, 2016.
[7] C. Briat and M. Khammash. Computer control of gene expression: Robust setpoint tracking of protein
mean and variance using integral feedback. In 51st IEEE Conference on Decision and Control, pages
3582–3588, Maui, Hawaii, USA, 2012.
[8] C. Briat and M. Khammash. Integral population control of a quadratic dimerization process. In 52nd
IEEE Conference on Decision and Control, pages 3367–3372, Florence, Italy, 2013.
39

[9] C. Briat, C. Zechner, and M. Khammash. Design of a synthetic integral feedback circuit: dynamic
analysis and DNA implementation. ACS Synthetic Biology, 5(10):1108–1116, 2016.
[10] B. F. Cress, E. A. Trantas, F. Ververidis, R. J. Linhardt, and M. A. G. Koffas. Sensitive cells: enabling
tools for static and dynamic control of microbial metabolic pathways. Current Opinion in Biotechnology,
36:205–214, 2015.
[11] C. Cuba Samaniego and E. Franco. An ultrasensitive biomolecular network for robust feedback control.
In 20th IFAC World Congress, pages 11437–11443, 2017.
[12] D. Del Vecchio and R. M. Murray, editors. Biomolecular Feedback Systems. Princeton University Press,
2015.
[13] 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.
[14] C. Guiver, H. Logemann, R. Rebarber, A. Bill, B. Tenhumberg, D. Hodgson, and S. Townley. Integral
control for population management. Journal of Mathematical Biology, 70:1015–1063, 2015.
[15] W. Halter, Z. A. Tuza, and F. Allgöwer. Signal differentiation with genetic networks. In 20th IFAC
World Congress, pages 10938–10943, 2017.
[16] A. W. K. Harris, J. A. Dolan, C. L. Kelly, J. Anderson, and A. Papachristodoulou. Designing genetic
feedback controllers. IEEE Transactions on Biomedical Circuits and Systems, 9(4):475–484, 2015.
[17] J. P. Hespanha. Moment closure for biochemical networks. In 3rd International Symposium on Communications, Control and Signal Processing, pages 142–147, St. Julian’s, Malta, 2008.
[18] M. Kaern, T. C. Elston, W. J. Blake, and J. J. Collins. Stochasticity in gene expression: from theories
to phenotypes. Nature Reviews Genetics, 6(6):451–464, 2005.
[19] E. Levine, Z. Zhang, T. Kuhlman, and T. Hwa. Quantitative characteristics of gene regulation by small
rna. PLoS Biol, 5(9):e229, 2013.
[20] G. Lillacci, S. K. Aoki, D. Schweingruber, and M. Khammash. A synthetic integral feedback controller
for robust tunable regulation in bacteria. bioRxiv, 2017.
[21] P. Milner, C. S. Gillepsie, and D. J. Wilkinson. Moment closure approximations for stochastic kinetic
models with rational rate laws. Mathematical Biosciences, 231:99–104, 2011.
[22] J. D. Murray. Mathematical Biology Part I. An Introduction. 3rd Edition. Springer-Verlag Berlin
Heidelberg, 2002.
[23] N. Olsman, A.-A. Ania-Ariadna, F. Xiao, Y. P. Leong, J. Doyle, and R. Murray. Hard limits and
performance tradeoffs in a class of sequestration feedback systems. bioRxiv, 2017.
[24] J. Paulsson. Summing up the noise in gene networks. Nature, 427:415–418, 2004.
[25] Y. Qian and D. Del Vecchio. Realizing “integral control” in living cells: How to overcome leaky
integration due to dilution? bioRxiv, 2017.
[26] M. Rullan, D. Benzinger, G. W. Schmidt, A. Gupta, A. Milias-Argeitis, and M. Khammash. Optogenetic
single-cell control of transcription achieves mrna tunability and reduced variability. BioRxiv, 2017.
[27] L. Schukur and M. Fussenegger. Engineering of synthetic gene circuits for (re-)balancing physiological
processes inchronic diseases. WIREs Systems Biology and Medicine, 8:402–422, 2016.
[28] P. Smadbeck and Y. N. Kaznessis. A closure scheme for checmical master equations. Proc. Natl. Acad.
Sci. USA., 110(35):14261–14265, 2013.
40

[29] M. Thattai and A. van Oudenaarden. Intrinsic noise in gene regulatory networks. Proceedings of the
National Academy of Sciences, 98(15):8614–8619, 2001.
[30] N. Venayak, N. Anesiadis, W. R. Cluett, and R. Mahadevan. Engineering metabolism through dynamic
control. Current Opinion in Biotechnology, 34:142–152, 2015.
[31] H. Ye and M. Fussenegger. Synthetic therapeutic gene circuits in mammalian cells. FEBS Letters,
588(15):2537–2544, 2014.
[32] S. M. Yoo, D. Na, and S. Y. Lee. Design and use of synthetic regulatory small rnas to control gene
expression in escherichia coli. Nat. Protocols, 8(9):1694–1707, 2013.

41

