1

Computing the projected reachable set of switched
affine systems: an application to systems biology

arXiv:1705.00400v1 [] 1 May 2017

Francesca Parise, Maria Elena Valcher and John Lygeros

Abstract—A fundamental question in systems biology is what
combinations of mean and variance of the species present in
a stochastic biochemical reaction network are attainable by
perturbing the system with an external signal. To address this
question, we show that the moments evolution in any generic
network can be either approximated or, under suitable assumptions, computed exactly as the solution of a switched affine
system. Motivated by this application, we propose a new method
to approximate the reachable set of switched affine systems. A
remarkable feature of our approach is that it allows one to easily
compute projections of the reachable set for pairs of moments of
interest, without requiring the computation of the full reachable
set, which can be prohibitive for large networks. As a second
contribution, we also show how to select the external signal in
order to maximize the probability of reaching a target set. To
illustrate the method we study a renown model of controlled gene
expression and we derive estimates of the reachable set, for the
protein mean and variance, that are more accurate than those
available in the literature and consistent with experimental data.

I. I NTRODUCTION
One of the most impressive results achieved by synthetic
biology in the last decade is the introduction of externally
controllable modules in biochemical reaction networks. These
are biochemical circuits that react to external signals, as for
example light pulses [1], [2], [3] or concentration signals [4],
[5], allowing researchers to influence and possibly control the
behavior of cells in vivo. To fully exploit these tools, it is
important to first understand what range of behaviors they
can exhibit under different choices of the external signal. For
deterministic systems, this amounts to computing the set of
states that can be reached by the controlled system trajectories
starting from a known initial configuration [6], [7]. Since
chemical species are often present in low copy numbers
inside the cell, biochemical reaction networks can however be
inherently stochastic [8]. In other words, if we apply the same
signal to a population of identical cells, then every cell will
have a different evolution (with different likelihood), requiring
a probabilistic analysis.
If we interpret each cell has an independent realization, we
can then study the effect of the external signal on a population
of cells by characterizing how such a signal influences the
F. Parise is with the Laboratory for Information and Decision
Systems, MIT, Cambridge, MA: parisef@mit.edu, M.E.
Valcher is with the Department of Information Engineering, University
of Padova, Italy: meme@dei.unipd.it and J. Lygeros is
with the Automatic Control Laboratory, ETH, Zurich, Switzerland:
lygeros@control.ee.ethz.ch. We thank M. Khammash for
allowing us to perform the experiments in the CTSB Laboratory, ETH, and
J. Ruess and A.M. Argeitis for their help in collecting the data of Fig. 5a).
This work was supported by the SNSF grant number P2EZP2 168812.

moments of the underlying stochastic process. Specifically, in
this paper we pose the following question:
“What combinations of moments of the stochastic process
can be achieved by applying the external signal?”
This approach is motivated for example by biotechnology
applications, where one would like to control the average
behavior of the cells in large populations, instead of each
cell individually. More on the theoretical side, this perspective
can be useful to investigate fundamental questions on noise
suppression in biochemical reaction networks, as in [9].
The cornerstone of our approach is the observation that
while the number of copies in each cell is stochastic, the
evolution of the moments is deterministic and can either
be described or approximated by a switched affine system.
Consequently, the above question can be reformulated as a
reachability problem in the moment space. Computing the
exact reachable set of a switched affine system is in general
far from trivial, see [10], [11]. We thus start our analysis by
proposing a new method to approximate the reachable set of a
switched affine system. This is an extension of the hyperplane
method for linear systems suggested in [12] and is of interest
on its own. We then show how to apply the proposed approach
to biochemical reaction networks by distinguishing two cases:
1) If all the reactions follow the laws of mass action kinetics
and are at most of order one, the system of moments
equations is switched affine. Consequently, for this class
of networks, the above question can be solved by directly
applying the newly suggested hyperplane method in the
moments space;
2) For all other reaction networks the moments equations
are in general non-closed (i.e., the evolution of mean and
variance depends on higher order moments). We show
however that the evolution of the probability of being in
a given state can be described by an infinite dimensional
switched system and that the desired moments can be
computed as the output of such system. We then show:
i) How to approximate such an infinite dimensional
system with a finite dimensional one, by extending the
finite state projection method [13] to controllable networks, ii) How to compute the reachable set of the finite
dimensional system by applying the newly suggested
hyperplane method in the probability space, and iii) How
to recover an approximation of the original reachable set
from the reachable set of the finite dimensional system.
In the last part of the paper, we change perspective and,
instead of focusing on population properties, we consider the
behaviour of a single cell (i.e., a single realization of the

2

process), given a fixed initial condition or an initial probability
distribution. Such perspective has been commonly employed
for the case without external signals, see e.g. [14], [15], [16],
[13]. Our objective is to show how the external signal can be
used to control single cell realizations by posing the following
question
“What external signal should be applied to maximize the
probability that the cell trajectory reaches a prespecified
subset of the state space at the end of the experiment?”
We show that such a problem can be addressed by using
similar tools as those derived for the population analysis.
Comparison with the literature: A vast literature has been
devoted to the analysis of the reachable set of piecewiseaffine systems in the context of hybrid systems, see e.g. [10],
[17], [18], [19], [20], [21], [22] among many. Our results
are different because we exploit the specific structure of the
problem at hand, that is, the fact that the switching signal is
a control variable and that the dynamics in each mode are
autonomous and affine. In other words, we consider switched
affine systems for which the switching signal is the only
control action. We also note that many different methods have
been proposed in the literature to compute the reachable set
of generic nonlinear systems. Among these there are level set
methods [23], ellipsoidal methods [24] and sensitivity based
methods [25]. For example, we became aware at the time
of submission that the authors of [26] extended our previous
works [27], [28] by suggesting the use of ellipsoidal methods.
It is important to stress that the choice of a method that
scales well with the system size is essential in our context,
since biochemical networks are typically very large. Moreover,
biologists are often interested in analyzing the behavior of
only a few chemical species of the possibly many involved
in the network. Consequently, one is usually interested in
computing the projection of the reachable set (which is a
high-dimensional object) on some low-dimensional space of
interest. The hyperplane method that we propose stands out
in this respect since, by using a method tailored for switched
systems, it allows one to compute directly the projections of
the reachable set, without requiring the computation of the
full high-dimensional reachable set first. We thus avoide the
curse of dimensionality that characterises all the previously
mentioned methods. We note that part of the results of this
paper appeared in our previous works [28], [29]. Specifically,
in [28] we first suggest the use of the hyperplane method to
compute the reachable set of biochemical networks with linear
moment equations, which we then adapted in [28] to the case
of switched affine moment equations. As better detailed in
Section IV-A, the assumptions made both in [28] and [29]
do not allow for bimolecular reactions, which are instead
present in the vast majority of biochemical networks. The key
contribution of this paper is the generalisation of our analysis
to any biochemical network by using the approach described
in point 2) above. The analysis of single cell realizations is
also entirely new.
Outline: In Section II we present the hyperplane method.
In Section III-A we review how to compute the hyperplane
constants for linear systems, while in Section III-B we propose

a new procedure for switched affine systems. In Section IV
we introduce stochastic biochemical reaction networks and the
controlled chemical master equation (CME). Additionally, we
recap how to derive the moments equations from the CME
(Section IV-A) and we derive an extension of the finite state
projection method to controlled biochemical networks (Section
IV-B). In Section V we show how to compute the reachable
set of biochemical networks and in Section VI we derive the
results on single cell realizations. Section VII illustrates our
theoretical results on a gene expression case study.
Notation: Given a < b ∈ N, we set N[a, b] := {a, a +
1, . . . , b}. Given a set S, the symbol ∂S denotes its boundary,
conv(S) its convex hull and |S| its cardinality. For a vector
x ∈ Rn , xp := [x]p denotes its pth component, |x| :=
[|x1 |> , . . . , |xn |> ]> and kxk∞ := maxp=1,2,...,n |xp | denotes
the infinity norm. 1 denotes a vector of all ones. Given two
random variables Z1 , Z2 , we denote by V[Z1 ] and V[Z1 , Z2 ]
their variance and covariance, respectively.
II. R EACHABILITY TOOLS
A. The reachable set and the hyperplane method
Consider the n-dimensional nonlinear control system
ẋ(t) = f (x(t), σ(t)),

t ≥ 0,

(1)

where x is the n-dimensional state and σ the m-dimensional
input function. Set a final time T > 0 and let S be the set of
admissible input functions that we assume to be a subset of
the set of all measurable functions that map [0, T ] into Rm .
We assume that the function f : Rn × Rm → Rn is such that,
for every initial condition x(0) ∈ Rn and every input function
σ ∈ S, the solution of (1), denoted by x(t; x(0), σ), t ≥ 0, is
well defined and unique at every time t ≥ 0. The reachable
set of system (1) at time T is defined as the set of all states
x ∈ Rn that can be reached at time T , starting from x(0), by
using an admissible input function σ ∈ S.
Definition 1 (Reachable set at time T ). The reachable set at
time T > 0 from x(0) = x0 , for system (1) with admissible
input set S, is
RT (x0 ) := {x ∈ Rn | ∃ σ ∈ S : x = x(T ; x0 , σ)}.

(2)

From now on we will assume that the set RT (x0 ) is
compact, since this will be the case for all the systems of
interest analysed in the following. Computing such a reachable
set for nonlinear systems is in general a very difficult task. For
the case of linear systems with bounded inputs a method to
construct an outer approximation of RT (x0 ) as the intersection
of a family of half-spaces that are tangent to its boundary (see
Fig. 1) was proposed in [12].
We present here a generalisation of this method to system (1). For a given direction c ∈ Rn , let us define
vT (c) :=

max

c> x,

x∈RT (x0 )

(3)

where, for simplicity, we omitted the dependence of vT (c) on
the initial condition x0 . Let
HT (c) := {x ∈ Rn | c> x = vT (c)}

(4)

3

required to obtain a good characterisation of the reachable
set. In the next subsection we show how to avoid this curse
of dimensionality, in cases when only the projection of the
reachable set on a plane of interest is needed.
B. The output reachable set
Let the output of system (1) be
y(t) = Lx(t),

(7)

p×n

Fig. 1. Illustration of the hyperplane method for a convex reachable set
RT (x0 ) (in blue). The external parallelogram is the outer approximation, the
region in between the dotted lines is the inner approximation.

be the corresponding hyperplane. By definition of the constant
vT (c), the associated half-space
HT (c) := {x ∈ Rn | c> x ≤ vT (c)}

(5)

is a superset of RT (x0 ). We note that if ∂RT (x0 ) is smooth,
then HT (c) is the tangent plane to ∂RT (x0 ). By evaluating
the above hyperplanes and half-spaces for various directions,
one can construct an outer approximation of the reachable set,
as illustrated in the next theorem. If the reachable set is convex
then an inner approximation can also be derived.
Theorem 1 (The hyperplane method [12]). Given system (1),
an initial condition x0 ∈ Rn , a fixed time T > 0, an integer
number D ≥ 2, and a set of D directions C := {c1 , . . . , cD },
define the half-spaces HT (cd ) as in (5), for d = 1, . . . , D.
1) The set
D
d
Rout
T (x0 ) := ∩d=1 HT (c )

is an outer approximation of the reachable set RT (x0 )
at time T starting from x0 .
2) If the set RT (x0 ) is convex and for each d = 1, 2, . . . , D,
we select a (tangent) point
x?T (cd ) ∈ RT (x0 ) ∩ HT (cd )

(6)

then the set
Rin
T (x0 )

:= conv

{x?T (cd ), d

= 1, 2, . . . , D}



is an inner approximation of the reachable set RT (x0 )
at time T starting from x0 .

Remark 1. We note that by construction the outer approximation Rout
T (x0 ) is a convex object. Specifically, when the
number of hyperplanes tends to infinity Rout
T (x0 ) coincides
with the convex hull of RT (x0 ). Similarly, for any set RT (x0 ),
the set Rin
T (x0 ) is an inner approximation of the convex hull
of RT (x0 ). However, the inner approximation of the convex
hull of a set is an inner approximation of the set itself only if
such set is convex, as assumed in the previous theorem. 
The main advantage of this method is that hyperplanes
are very easy objects to handle and visualise. The main
disadvantage is that the higher the dimension n of the state
space, the higher in general is the number of directions D

for L ∈ R
, and the output reachable set be the set of all
output values that can be generated at time T from x(0) = x0 ,
by using an admissible input function σ ∈ S.
Definition 2 (Output reachable set at time T ). The output
reachable set RyT (x0 ) from x0 at time T > 0, for system (1)
with admissible input set S and output as in (7), is
RyT (x0 ) := {y ∈ Rp | ∃ x ∈ RT (x0 ) : y = Lx}.
For simplicity, in the following we restrict our discussion
to the case of a two-dimentional output vector, that is
>

l1 x(t)
y(t) = Lx(t) = >
∈ R2 ,
(8)
l2 x(t)
for some l1 , l2 ∈ Rn , the generalization to higher dimentions
is however immediate. Note that, for any pair of indices
i, j ∈ {1, . . . , n}, i 6= j, one can recover the projection of
the reachable set RT (x0 ) onto an (xi , xj )-plane of interest by
imposing l1 = ei and l2 = ej . The two-dimentional output
vector case can therefore be applied to study the relation
between the mean behavior of two species or between mean
and variance of a single species in large biochemical networks.
In the following theorem we show that inner and outer
approximations of RyT (x0 ) can be efficiently computed by
selecting only hyperplanes that are perpendicular to the plane
of interest.
Theorem 2 (Projection on a two dimensional subspace).
Consider system (1), with output (8) and initial condition
x0 ∈ Rn . Let T > 0 be a fixed time, D ≥ 2 an integer number
and choose D values γ d ∈ R. Set cd := l2 − γ d l1 ∈ Rn and
HTy (γ d ) := {y ∈ R2 | y2 ≤ γ d y1 + vT (cd )},
where vT (cd ) is as in (3). Set yT? (γ d ) := Lx?T (cd ), where
x?T (cd ) is defined as in (6). Then the set
y
d
Ry,out
(x0 ) := ∩D
d=1 HT (γ )
T

is an outer approximation of
convex then the set

RyT (x0 ).

Moreover, if RT (x0 ) is

?
d
Ry,in
T (x0 ) := conv {yT (γ ), d = 1, 2, . . . , D}

is an inner approximation of

(9)



(10)

RyT (x0 ).

for any ȳ ∈ RyT (x0 )
ȳ > = [l1> x̄, l2> x̄]. By



Proof: By definition,
there exists
an x̄ ∈ RT (x0 ) such that
Theorem 1,
for any direction cd it holds that RT (x0 ) ⊂ HT (cd ). Consequently, x̄ ∈ RT (x0 ) implies x̄ ∈ HT (cd ). By substituting the
definition of cd given in the statement we get
x̄ ∈ HT (cd ) ⇔ (cd )> x̄ ≤ vT (cd ) ⇔
⇔ (l2 − γ d l1 )> x̄ ≤ vT (cd ) ⇔ l2> x̄ ≤ γ d l1> x̄ + vT (cd ).

4

The last inequality implies ȳ > = [l1> x̄, l2> x̄] ∈ HTy (γ d ).
Consequently, RyT (x0 ) ⊆ HTy (γ d ) for any γ d and therefore
RyT (x0 ) ⊆ Ry,out
(x0 ). If RT (x0 ) is convex, then RyT (x0 )
T
is convex as well. The points yT? (γ d ) belong to RyT (x0 ) by
construction. Consequently, by convexity, it must hold that
y
Ry,in
T (x0 ) ⊆ RT (x0 ).
III. C OMPUTING THE TANGENT HYPERPLANES
The success of the hyperplane method hinges on the possibility of efficiently evaluating, for any given direction c, the
constant vT (c) in (3). Note that this problem is equivalent to
the following finite time optimal control problem
vT (c) := max c> x(T )

(11)

σ∈S

s.t.

ẋ(t) = f (x(t), σ(t)),

∀t ∈ [0, T ],

x(0) = x0 .
In the rest of this section, we aim at solving (11). To this end,
we start by recalling the linear case, for which the hyperplane
method was originally derived in [12].
A. Linear systems with bounded input
The hyperplane method was originally proposed for linear
systems with bounded inputs
ẋ(t) = Ax(t) + Bσ(t),
n

n×n

n×m

(12)
m

where x(t) ∈ R , A ∈ R
, B ∈R
and σ(t) ∈ R .
Since biological signals are non-negative and bounded, we
here make the following assumption on the input set S.
Assumption 1. The input function σ belongs to the admissible set SΣ := {σ | σ(t) ∈ Σ, ∀t ∈ [0, T ]}, where
Σ = Σ1 × . . . × Σm . Moreover, there exist σ̄r > 0, r ∈
N[1, m], such that either (a) every set Σr is the interval
Σcr := [0, σ̄r ] (continuous and bounded input set), or (b)
for every
 set Σr there exists 2 ≤ qr < +∞ such that
Σdr := 0 = σr1 < σr2 < . . . < σrqr = σ̄r ⊂ R≥0 (finite input
set). We set Σc := Σc1 × . . . × Σcm , Σd := Σd1 × . . . × Σdm , and
denote by SΣc and SΣd the corresponding admissible sets.
In the case of a continuous and bounded input set, i.e. under
Assumption 1-(a), it was shown in [12] that it is possible to
solve the control problem in (11) in closed form by using the
Maximum Principle [30].
Proposition 1 (Tangent hyperplanes for linear systems with
bounded and continuous inputs). Consider system (12) and
suppose that Assumption 1-(a) holds. Define the following
admissible input function, expressed component-wise for every
rth entry, r = 1, . . . , m, as


if c> eA(T −t) br > 0;
σ̄r
σr? (t) : = 0
(13)
if c> eA(T −t) br < 0;


0 ≤ σr ≤ σ̄r if c> eA(T −t) br = 0;
where br denotes the rth column of B. Then

RT 
Pm
vT (c) = c> eAT x0 + r=1 σ̄r 0 c> eA(T −t) br + dt, (14)

where [g(t)]+ denotes the positive part of the function, namely
[g(t)]+ = g(t) when g(t) > 0 and zero otherwise. Suppose
additionally that the pair (A, br ) is reachable, for every r ∈
N[1, m]. Then there exists no interval [τ1 , τ2 ], with 0 ≤ τ1 <
τ2 ≤ T , such that c> eA(T −t) br = 0 for every t ∈ [τ1 , τ2 ].
Consequently, a tangent point can be obtained as
RT
x?T (c) := eAT x0 + 0 eA(T −t) Bσ ? (t)dt.
(15)

The proof follows the same lines as [12, Lemma 2.1 and
Theorem 2.1] and is omitted for the sake of brevity.
By using the explicit characterisation given in Proposition 1
together with Theorems 1 and 2, one can efficiently construct
both an inner and an outer approximation of the (output)
reachable set for linear systems with continuous and bounded
input set Σc , as summarised in the next corollary. Therein we
also show how the same result can be extended to finite input
sets Σd .
Corollary 1 (The hyperplane method for linear systems).
Consider system (12) and suppose that either Assumption 1(a) or Assumption 1-(b) holds. Let vT (cd ) and x?T (cd ) be
in
computed as in (14) and (15). Then Rout
T (x0 ) and RT (x0 )
y,out
y,in
(RT (x0 ) and RT (x0 ), resp.) as defined in Theorem 1
(Theorem 2, resp.) are outer and inner approximations of
RT (x0 ) (of RyT (x0 ), resp.).
Proof: In the case of continuous and bounded input,
that is, under Assumption 1-(a), the reachable set RT (x0 )
is convex and the statement is a trivial consequence on
Theorems 1 and 2 and Proposition 1. We here show that the
same result holds also under Assumption 1-(b). The proof of
this second part follows from the fact that the reachable set
RcT (x0 ), obtained by using the continuous input set Σc , and
the reachable set RdT (x0 ), obtained by using the discrete input
set Σd , coincide. To prove this, let Rbb
T (x0 ) be the reachable
:=
{0,
σ̄
}
for
any
r, that is, the set of
set obtained using Σbb
r
r
vertices of Σc . Consider now an arbitrary point x̄ ∈ RcT (x0 ),
which is a compact set. By definition there exists an admissible
input function in Σc that steers x0 to x̄ in time T . Since Σc
is a convex polyhedron, by [31, Theorem 8.1.2], system (12)
with input set Σc has the bang-bang with bound on the number
of switchings (BBNS) property. That is, for each x̄ ∈ RcT (x0 )
there exists a bang-bang input function in Σbb that reaches x̄
in the same time T with a finite number of discontinuities.
c
Thus x̄ ∈ Rbb
T (x0 ). Since this is true for any x̄ ∈ RT (x0 ),
c
bb
bb
d
c
we get RT (x0 ) ⊆ RT (x0 ). From Σ ⊆ Σ ⊆ Σ we get
d
c
Rbb
T (x0 ) ⊆ RT (x0 ) ⊆ RT (x0 ), concluding the proof.
B. Switched affine systems
In this section, we propose an extension of the hyperplane
method to the case of a switched affine system of the form
ẋ(t) = Aσ(t) x(t) + bσ(t) ,

(16)

where the switching signal σ(t) ∈ N[1, I] is the input
function, I ≥ 2 is the number of modes, x(t) ∈ Rn and
Ai ∈ Rn×n , bi ∈ Rn for all i ∈ N[1, I]. We make the
following assumption.

5

Assumption 2. The switching signal σ(t) switches K times
within the finite set N[1, I] at fixed switching instants 0 = t0 <
. . . < tK+1 = T , that is, σ ∈ SIK , where

following optimisation problem
vT (c) := max

xk ,zik ,γik

s.t.
SIK

:= {σ | σ(t) = ik ∈ N[1, I], ∀t ∈ [tk , tk+1 ), k ∈ N[0, K]}.

For every k ∈ N[0, K] and i ∈ N[1, I] we define Āki :=
R (t −t )
e
and b̄ki = [ 0 k+1 k eAi τ dτ ]bi . Moreover, we
set xk := x(tk ). Note that under Assumption 2 the reachable
set of system (16) consists of a finite number of points that
can be computed by solving the state equations for each
possible switching signal. Since the cardinality of the set
SIK grows exponentially with K, this approach is however
computationally infeasible even for small systems. We here
show that, on the other hand, the hyperplane constants defined
in (11) can be computed by solving a mixed integer linear
program (MILP), thus allowing us to exploit the sophisticated
software that has been developed to solve large MILPs in the
last years.
Ai (tk+1 −tk )

Proposition 2 (Tangent hyperplanes for switched affine systems). Consider system (16) and suppose that Assumption 2
holds. Take a vector M ∈ Rn such that M ≥ |xk | componentwise for all k ∈ N[0, K]. Then
vT (c) = max

xk ,zik ,γik

c> xK+1

(17)

c> xK+1

(19)

zik+1 = (Āki xk + b̄ki )γik , ∀i ∈ Σ,
PI
k
∀k ∈ N[0, K],
i=1 γi = 1,
PI
k
xk = i=1 zi , ∀k ∈ N[1, K + 1],
x0 ∈ R assigned.

Finally, by using the big-M method in [33, Eq. (5b)], the first
equality constraint in the optimization problem (19) can be
equivalently replaced by
zik+1 ≤ (Āki xk + b̄ki ) + M(1 − γik ),

zik+1 ≥ −Mγik ,

zik+1 ≥ (Āki xk + b̄ki ) − M(1 − γik ),

zik+1 ≤ Mγik ,

leading to the equivalent reformulation given in (17).
We summarize our results on the hyperplane method for
switched affine systems in the next corollary, which is an
immediate consequence of Proposition 2 and Theorems 1, 2.
Corollary 2 (The hyperplane method for switched affine
systems). Given system (16), let x0 ∈ Rn be the initial state
and suppose that Assumption 2 holds. Let vT (cd ) be computed
y,out
(x0 ) as defined in
as in (17). Then Rout
T (x0 ) and RT
Theorems 1 and 2 are outer approximations of RT (x0 ) and
RyT (x0 ), respectively.


Note that in the case of switched affine systems it is not
possible to recover an inner approximation, since there is
no guarantee in general that the reachable set is convex. By
zik+1 ≥ (Āki xk + b̄ki ) − M(1 − γik ),
computing the convex hull of the points xK+1 in (17) for each
k+1
k+1
k
k
zi
≥ −Mγi , zi
≤ Mγi ,
direction c one could however recover an inner approximation
zik ∈ Rn , ∀k ∈ N[1, K + 1], ∀i ∈ N[1, I], of the convex hull of R (x ).
T
0
γik ∈ {0, 1}, ∀k ∈ N[0, K], ∀i ∈ N[1, I],
PI
IV. C ONTROLLED STOCHASTIC BIOCHEMICAL REACTION
xk = i=1 zik ∈ Rn , ∀k ∈ N[1, K + 1],
PI
NETWORKS
k
∀k ∈ N[0, K],
i=1 γi = 1,
A biochemical reaction network is a system comprising S
x0 ∈ Rn assigned.
molecular species Z1 , ..., ZS that interact through R reactions.
Let Z(t) = [Z1 (t), ..., ZS (t)]> be the vector describing the
Proof: To prove the statement we follow a procedure number of molecules present in the network for each species
similar to the one in [32, Section IV.A]. Under Assumption 2 at time t, that is, the state of the network at time t. Since
the switching signal σ(t) is such that σ(t) = ik , ∀t ∈ each reaction r is a stochastic event [8], Z(t) is a stochastic
[tk , tk+1 ), ∀k ∈ N[0, K]. Therefore, the finite time optimal process. In the following, we always use the upper case to
denote a process and the lower case to denote its realizations.
control problem in (11) can be rewritten as
For example, z = [z1 , ..., zS ]> denotes a particular realization
of the state Z(t) of the stochastic process at time t.
vT (c) = max
c> xK+1
(18)
A typical reaction r ∈ N[1, R] can be expressed as
ik ∈{1,...,I}
s.t.

zik+1 ≤ (Āki xk + b̄ki ) + M(1 − γik ),

s.t. xk+1 = Ākik xk + b̄kik

∀k ∈ N[0, K]

x0 ∈ R assigned.
Let us introduce the binary variables γik ∈ {0, 1} defined
so that, for each i ∈ N[1, I] and k ∈ N[0, K], γik = 1
if and only if the system is in mode i in the time interval
[tk , tk+1 ). Moreover, let us introduce a copy of the state
vector for each possible update of the system in each possible
mode: zik+1 = (Āki xk + b̄ki )γik . Then (18) is equivalent to the

0
0
ν1r
Z1 + . . . + νSr
ZS

−→

00
00
ν1r
Z1 + . . . + νSr
ZS , (20)

0
0
00
00
where ν1r
, . . . , νSr
∈ N and ν1r
, . . . , νSr
∈ N are the coefficients that determine how many molecules for each species
are respectively consumed and produced by the reaction. The
net effect of each reaction can thus be summarized with
the stoichiometric vector νr ∈ NS , whose components are
00
0
νsr
− νsr
for s = 1, . . . , S. We say that
is of order
PaS reaction
0
k if it involves k reactant units (i.e., s=1 νsr
= k) and we
distinguish two classes of reactions:

6

-uncontrolled reactions that happen, in the infinitesimal interval [t, t + dt], with probability
αr (θr , z)dt := θr · hr (z) · dt,

(21)

where hr (z) is a given function of the available molecules z
and θr ∈ R≥0 is the so-called rate parameter;
- controlled reactions for which there exists an external signal
ur (t) such that the reaction fires at time t with probability
ur (t) · αr (θr , z)dt.

(22)

counts can theoretically grow unbounded. Consequently, the
controlled CME in (23) is a system of infinitely many coupled
ordinary differential equations that cannot be solved, even
for very simple systems. Several analytical and computational
methods have been proposed in the literature to circumvent this
difficulty, see [34], [35], [36] for a comprehensive review. In
the following we limit our discussion to two methods: moment
equations [37] and finite state projection (FSP) [13].
A. The moment equations

In the following we refer to αr (θr , z) as the propensity of the
reaction and without loss of generality we assume that the con-
trolled reactions are the first Q ones. If hr (z) := ΠSs=1 νz0s
sr
we say that reaction r follows the laws of mass action kinetics
as derived in [8]. Our analysis can however be applied to
generic functions hr (z), allowing us to model different types
of kinetics, as the Michaelis-Menten [34, Section 7.3].
To illustrate the following results, we consider a model of
gene expression as running example.

We start by considering the case when all the reactions
follow the laws of mass action kinetics and are at most of
order one. In this case for each reaction r the propensity hr (z)
is affine in the molecule counts vector z and one can show
that the moments equations are closed (i.e., the dynamics of
moments up to any order k do not depend on higher order
moments), see for example [38]. Let x≤2 (t) be a vector whose
components are the moments of Z(t) up to second order. From
[38, Equations (6) and (7)] one gets

Example 1 (Gene expression reaction network). Consider a
biochemical network consisting of two species, the mRNA (M )
and the corresponding protein (P ), and the following reactions

ẋ≤2 (t) = A(u(t))x≤2 (t) + b(u(t)).

∅

α1 (kr ,z)

−−−−−−−−→
α2 (γr ,z)

−−−−−−−−→

M

M
∅

M
P

α3 (kp ,z)

−−−−−−−−→
α4 (γp ,z)

−−−−−−−−→

M +P
∅

where the parameters kr and kp are the mRNA and protein
production rates, while γr and γp are the mRNA and protein
degradation rates, respectively. The empty set notation is used
whenever a certain species is produced or degrades without
involving the other species. In this context, Z = [M, P ]> ,
z = [m, p]> , θ = [θ1 , θ2 , θ3 , θ4 ]> := [kr , γr , kp , γp ]> and the
stoichiometric matrix is


1 −1 0 0
ν := [ν1 , ν2 , ν3 , ν4 ] =
.
0 0 1 −1
In the case of mass action kinetics the propensities αr (θr , z)
can be further specified as α1 (kr , z) = kr , α2 (γr , z) = γr ·
m, α3 (kp , z) = kp · m, α4 (γp , z) = γp · p.

Note that since the propensity of each reaction depends
only on the current state of the system, the process Z(t) is
Markovian. Let p(t, z) := P[Z(t) = z] be the probability that
the realization of the process Z at time t is z. Following the
same procedure as in [8] one can derive a set of equations,
known as chemical master equation (CME), describing the
evolution of p(z, t) as a function of the external signal u(t)

ṗ(z, t) =

Q
X

[p(z − νr , t)αr (θr , z − νr ) − p(z, t)αr (θr , z)] ur (t)

r=1

+

R
X

[p(z − νr , t)αr (θr , z − νr ) − p(z, t)αr (θr , z)] , ∀z ∈ NS .

r=Q+1

(23)

Since the previous set of equations depends on the external
signal u we refer to it as the controlled CME. Typical biochemical reaction networks involve many different species, whose

(24)

Example 2. Consider the gene expression model of Example 1. Assume that the reactions follow the mass action kinetics
and that an external input signal influencing the first reaction,
that is the mRNA production, is available (as in [1], [2], [3],
[4], [5]), so that α1 (kr , z) := kr · u(t). Set
x≤2 := [E[M ], E[P ], V[M ], V [M, P ], V[P ]]> .
Then the moments evolution over time is expressed as
ẋ≤2 (t) = Ax≤2 (t) + Bu(t),

(25)

where


−γr 0
0
0
0
 kp −γp 0
0
0 


0
0 ,
A =  γr 0 −2γr
 0
0
kp −(γr + γp ) 0 
kp γ p
0
2kp
−2γp

 
kr
0
 
B =  kr  .
0
0


Since the input u(t) may appear in the entries of the A
matrix, the moment equations (24) are in general nonlinear.
To overcome this issue we introduce the following assumption
on the external signal u(t).
Assumption 3. The external signal u(t) can switch at most
K times within the set Σd , as defined in Assumption 1, at
preassigned switching instants 0 = t0 < . . . < tK+1 = T .
Assumption 3 imposes that the number of switchings and
their timing during a given experiment is fixed a priori.
This assumption can be motivated by the fact that changes
in the external stimulus are costly and/or stressful for the
cells. Moreover, it is trivially satisfied if the stimulus can
only be changed simultaneously with some fixed events, such
as culture dilution or measurements. The great advantage of
Assumption 3 is that, as illustrated in the following remark, it
allows us to rewrite the nonlinear moment equations (24) as a
switched affine system so that the theoretical tools described
in Section III-B can be applied.

7

Remark 2. The set Σd has finite cardinality I := Πm
r=1 qr
and we can enumerate its elements as ui , i ∈ N[1, I].
Consequently, for any fixed external signal u(t) satisfying
Assumption 3 we can construct a sequence of indices in N[1, I]
such that, at any time t, σ(t) = i if and only if u(t) = ui .
Such switching sequence σ satisfies Assumption 2.

B. The finite state projection
Let us introduce a total ordering {z j }∞
j=1 in the set of
all possible state realizations z ∈ NS . For the system in
Example 1, we could for instance use the mapping
z 1 = (0, 0), z 2 = (1, 0), z 3 = (0, 1), z 4 = (2, 0),
z 5 = (1, 1), z 6 = (0, 2), z 7 = (3, 0), z 8 = (2, 1), . . .
where (m, p) denotes the state with m mRNA copies and p
proteins (see Fig. 2).

Note that
 while
 the full matrix Fσ(t) is stochastic, the reduced
matrix Fσ(t) J is substochastic. Consequently, the probability
mass is in general not preserved in (28) (i.e. 1> P̄J (t) may
decrease with time). From now on, we denote by P (T ; σ)
and P̄J (T ; σ) the solutions at time T of system (27) and system (28), respectively, when the switching signal σ is applied.
The dependence on the initial conditions P (0) and PJ (0) is
omitted to keep the notation compact. As in the uncontrolled
case, the truncated system (28) is a good approximation of
the original system (27) if most of the probability mass lies in
J. However in the controlled case we need to guarantee that
this happens for all possible switching signals. This intuition
is formalized in the following assumption.
Assumption 4. For a given finite set of state indices J, an
initial condition PJ (0), a given tolerance ε > 0 and a finite
instant T > 0,

1> P̄J (T ; σ) ≥ 1 − ε, ∀σ ∈ SIK .

(29)

Note that Assumption 4 holds if and only if
1 − ε ≤ min

σ∈SIK

1> P̄J (T ; σ)



s.t. P̄˙J (t; σ) = Fσ(t) J P̄J (t; σ), P̄J (0) = PJ (0).

Fig. 2.

State space for the gene expression system of Example 1.

Following the same steps as in [13] and setting1 Pj (t) :=
p(z j , t), the controlled CME in (23) can be rewritten as the
nonlinear infinite dimensional system
Ṗ (t) = F (u(t))P (t),

(27)

with switching signal σ(t) constructed from u(t) as detailed
i
in Remark 2, I = Πm
r=1 qr modes and matrices Fi := F (u ).
Note that system (27) can also be thought of as a Markov
chain with countably many states z j ∈ NS and time-varying
transition matrix Fσ(t) .
As in the FSP method for the uncontrolled CME [13], one
can try to approximate the behavior of the infinite Markov
chain in (27) by constructing a reduced Markov chain that
keeps track of the probability of visiting only the states
indexed in a suitable set J. To this end, let us define the
reduced order system


P̄˙J (t) = Fσ(t) J P̄J (t), P̄J (0) = PJ (0),
(28)
where PJ (0) is the subvector of P (0) corresponding to the
indices in J, and [F ]J denotes the submatrix of F obtained
by selecting only the rows and columns with indices in J.
1 Not

Proposition 3 (FSP for controlled CME). If Assumptions 2
and 4 hold, then for every switching signal σ ∈ SIK , it holds

(26)

where P (t) is an infinite dimensional vector with entries in
[0, 1]. If the signal
u(t) satisfies Assumption 3, then (26) can be rewritten as
an infinite dimensional linear switched system
Ṗ (t) = Fσ(t) P (t),

This problem has the same structure as (11). Therefore, as
illustrated in Section III-B, Assumption 4 can be checked by
solving the MILP (17) for the switched affine system (28) by
setting c = 1 and M = 1. Under Assumption 4, the following
relation between the solutions of (27) and (28) holds.

to be confused with the symbol used to denote the amount of protein.

Pj (T ; σ) ≥ P̄j (T ; σ),

∀j ∈ J

kPJ (T ; σ) − P̄J (T ; σ)k1 ≤ ε.
Proof: This result has been proven in [13] for linear
systems. We extend it here to the case of switched systems
with K switchings. Note that for any i ∈ N[1, I], Fi := F (ui )
has non-negative off diagonal elements [13]. Hence, using the
same argument as in [13, Theorem 2.1] it can be shown that
for any index set J, and any τ ≥ 0
[exp(Fi τ )]J ≥ exp([Fi ]J τ ) ≥ 0,

∀i ∈ 1, . . . , I.

Consider an arbitrary switching signal σ ∈ SIK . We have
PJ (T ; σ) = [ΠK
k=0 exp(Fik (tk+1 − tk )) · P (0)]J
≥
≥

ΠK
k=0 [exp(Fik (tk+1 − tk ))]J
ΠK
k=0 exp([Fik ]J (tk+1 − tk ))

(30)

· PJ (0)

· PJ (0) = P̄J (T ; σ).
P
Moreover, from 1 =
j=1 Pj (T ; σ) ≥
j∈J Pj (T ; σ) =
1> PJ (T ; σ) and Assumption 4, we get
P∞

1> P̄J (T ; σ) ≥ 1 − ε ≥ 1> PJ (T ; σ) − ε.

(31)

Combining (30) and (31) yields 0 ≤ 1> PJ (T ; σ) −
1> P̄J (T ; σ) ≤ ε, thus kPJ (T ; σ) − P̄J (T ; σ)k1 ≤ ε.

8

V. A NALYSIS OF THE REACHABLE SET
We here show how the reachability tools of Sections II
and III can be applied to the moment equation and FSP
reformulations derived in Sections IV-A and IV-B, under
different assumptions. Fig. 3 presents a conceptual scheme
of this section.

Fig. 3.

Conceptual scheme for the reachable set analysis of biochemical networks.

A. Reachable set of networks with affine propensities via
moment equations
The methods developed in Sections II and III can be applied
to the moments equations in (24) to approximate the desired
projected reachable set. To illustrate the proposed procedure,
we distinguish two cases depending on whether the external
signal u(t) influences reactions of order zero or one.
1) Linear moments equations: We start by considering the
case when all and only the reactions of order zero are con>
trolled, so that hr (z) = 1 for r ∈ N[1, Q] and hr (z) = νr0 z
for r ∈ N[Q + 1, R]. This is the simplest scenario since the
system of moment equations given in (24) becomes linear
ẋ≤2 (t) = Ax≤2 (t) + Bu(t),

(32)

see [38, Equations (6) and (7)]. Consequently, the theoretical
results of Section III-A can be applied to (32) by setting
σ(t) ≡ u(t). If the external signal u ≡ σ satisfies Assumption
1, both inner and outer approximations of the reachable set
can be computed by using Corollary 1.
2) Switched affine moments equations: If reactions of order
one are controlled then the external input u(t) appears also in
the entries of the A matrix and system (24) is nonlinear. To
overcome this issue we exploit Assumption 3. Specifically, let
σ(t) be the switching signal associated with u(t) as described
in Remark 2. Then (24) can be equivalently rewritten as the
switched affine system
ẋ≤2 (t) = Aσ(t) x≤2 (t) + bσ(t) ,

(33)

with matrices Ai := A(ui ), bi := b(ui ), for all i ∈ N[1, I].
Consequently, the theoretical results of Section III-B can be
applied to (33) and an outer approximation of the reachable
set can be computed by using Corollary 2.

B. Reachable set of networks with generic propensities via
finite state projection
If the network contains reactions of order higher than one or
if the reactions do not follow the laws of mass action kinetics,
then hr (z) might be non-affine. In such cases, the arguments
illustrated in the previous subsection cannot be applied. We
here show how the FSP approximation of the CME derived in
Section IV-B can be used to overcome this problem.
Firstly note that, from system (27), one can compute the
evolution of the uncentered moments of Z(t), as a linear
function of P (t). 2 For example, if we let zsj be the amount of
species Zs in the state z j , then the mean E[Zs ] of any species

>
s can be obtained as l> P (t), by setting l := zs1 , zs2 , . . . ,
and the second uncentered moment E[Zs2 ] can be obtained as

>
l> P (t), by setting l := (zs1 )2 , (zs2 )2 , . . . . Consequently
the desired projected reachable set coincides with the output
reachable set of the infinite dimensional linear switched system (27) with linear output

>
y(t) = l1 , l2 P (t),
(34)
where l1 and l2 are the infinite vectors associated with any
desired pair of moments. Note that l1 and l2 are non-negative.
Example 1 (cont.) With the ordering introduced at the
beginning of the section, the uncentered protein moments up
to order two can be computed as the output of (27) by setting

>
l1 = 0 0 1 0 1 2 0 1 . . .
,
(35)


>
l2 = 0 0 1 0 1 4 0 1 . . .
.
Let lj1 and lj2 be the j-th components of the vectors l1 and l2 ,
respectively, as defined in (34). For a given species of interest
s and set J, we denote by
P
y1 (t; σ) :=

j∈J
P

lj1 · Pj (t; σ)

j∈J

Pj (t; σ)

P
, y2 (t; σ) :=

j∈J
P

lj2 · Pj (t; σ)

j∈J

Pj (t; σ)
(36)

the moments associated with l1 and l2 conditioned on
the fact that Z(t) is in J and the switching signal σ
is applied. For example if one is interested in the mean
and second order moment of a specific species Zs (t) we
get y1 (t; σ) = E [Zs (t)
 | Z(t) ∈ J, σ(·)] and y2 (t; σ) =
E Zs2 (t) | Z(t) ∈ J, σ(·) . The aim of this section is to obtain an outer approximation of the output reachable set of
the infinite system (27) with the nonlinear output (36), by
using computations involving only the finite dimensional system (28). To this end, we define the two entries of the output
of the finite dimensional system as
P
ȳ1 (t; σ) := j∈J lj1 · P̄j (t; σ) =: (¯l1 )> P̄J (t; σ)
P
(37)
ȳ2 (t; σ) :=
l2 · P̄j (t; σ) =: (¯l2 )> P̄J (t; σ).
j∈J j

Theorem 3. Suppose Assumptions 3 and 4 hold. Let RyT (x0 )
be the output reachable set at time T > 0 of system (27)
with output (36). Choose D values γ d ∈ R and set cd :=
(¯l2 ) − γ d (¯l1 ) ∈ Rn , with ¯l1 , ¯l2 as in (37). Set
HTy (γ d ) := {w ∈ R2 | w2 ≤ γ d w1 + v̄T (cd ) + δ(γ d )},
2 The reachable set for the centered moments can be immediately computed
from the reachable set of the uncentered ones, since there is a bijective relation
between the set of centered and uncentered moments up to any desired order.

9

where v̄T (cd ) is the constant that makes the hyperplane
HT (cd ) in (5) tangent to the reachable set of the finite
system (28) (i.e. v̄T (cd ) can be computed as in (17)) and
δ(γ d ) :=

2ε
1−ε

Therefore for every signal σ and every γ d it holds
y2 (T ; σ) ≤ γ d y1 (T ; σ) + v̄T (cd ) + δ(γ d ) and consequently
[y1 (T ; σ), y2 (T ; σ)]> ∈ HTy (γ d ).

· (max{0, −γ d } · k¯l1 k∞ + k¯l2 k∞ ),

with ε as in Assumption 4. Then the set Ry,out
(x0 ) :=
T
y
y
d
∩D
{H
(γ
)}
is
an
outer
approximation
of
R
(x
).

0
d=1
T
T
Proof: Firstly note that if the external signal u satisfies
Assumption 3 then the corresponding switching signal σ(t)
(constructed as in Remark 2) satisfies Assumption 2. Let
R̄yT (x0 ) be the output reachable set of the finite dimensional
system (28) with output (37). Proposition 2 guarantees that for
any direction cd the constant v̄T (cd ) that makes
H̄Ty (γ d ) := {w ∈ R2 | w2 ≤ γ d w1 + v̄T (cd )}
tangent to R̄yT (x0 ) can be computed by solving the MILP (17)
for system (28). The main idea of the proof is to show that if
we shift the halfspace H̄Ty (γ d ) by a suitably defined constant
δ(γ d ) we can guarantee that the original reachable set RyT (x0 )
is a subset of the shifted halfspace HTy (γ d ) defined in the
statement. The result then follows since Ry,out
(x0 ) is defined
T
as the intersection of hyperspaces containing RyT (x0 ).
To derive the constant δ(γ d ) we start by focusing on the
first component of the output and for simplicity we will
omit the dependence on (T ; σ) in Pj , P̄j , y and ȳ. Take any
switching signal σ ∈ SIK . By taking into account the following
conditions: (1) lj1 ≥ 0 for all j P
∈ J; (2) Pj ≥ P̄j for all j ∈ J,
due to Proposition 3, and (3) j∈J Pj ≤ 1, we get y1 ≥ ȳ1 .
Consequently, at time t = T we have
P
P
l1 ·Pj
|y1 − ȳ1 | = y1 − ȳ1 = Pj∈J jPj − j∈J lj1 · P̄j
j∈J
P
P
lj1 ·Pj
−
≤ j∈J
lj1 · P̄j
1−ε
 P j∈J

P
ε
1
1
= 1 + 1−ε
j∈J lj · Pj −
j∈J lj · P̄j
P
P
ε
1
1
= 1−ε
j∈J lj · Pj +
j∈J lj · (Pj − P̄j )


P
P
ε
≤ k¯l1 k∞ 1−ε
j∈J Pj +
j∈J (Pj − P̄j )


ε
2ε
≤ k¯l1 k∞ 1−ε
+ kPJ − P̄J k1 ≤ k¯l1 k∞ 1−ε
,
P
P
where we used
j∈J Pj ≥
j∈J P̄j ≥ 1 − ε (due to
Assumption 4), and Pj ≥ P̄j , kPJ − P̄J k1 ≤ ε (following from
2ε
.
Proposition 3). To summarize, ȳ1 ≤ y1 ≤ ȳ1 + k¯l1 k∞ 1−ε
2ε
2
¯
Similarly, it can be proven that ȳ2 ≤ y2 ≤ ȳ2 + kl k∞ 1−ε
.
Consider any pair (y1 , y2 ) ∈ RyT (x0 ) and the associated
pair (ȳ1 , ȳ2 ) ∈ R̄yT (x0 ) (i.e. the two output pairs obtained
from (27) and (28) when the same σ is applied). Note that
(ȳ1 , ȳ2 ) ∈ R̄yT (x0 ) implies (ȳ1 , ȳ2 ) ∈ H̄Ty (γ d ) for any γ d .
The previous relations then imply that if γ d ≥ 0,
2ε
2ε
y2 ≤ ȳ2 + k¯l2 k∞ 1−ε
≤ γ d ȳ1 + v̄T (cd ) + k¯l2 k∞ 1−ε
≤ γ d y1 + v̄T (cd ) + k¯l2 k∞ 2ε = γ d y1 + v̄T (cd ) + δ(γ d ).
1−ε

d

VI. A NALYSIS OF SINGLE CELL REALIZATIONS
The previous analysis focused on characterising what combinations of moments of the stochastic biochemical reaction
network are achievable by using the available external input.
In this section, we change perspective and instead of looking
at population properties we focus on single cell trajectories.
Specifically, we are interested in characterising the probability
that a single realization of the stochastic process will satisfy a
specific property at the final time T (e.g. the number of copies
of a certain species is higher/lower than a certain threshold)
when starting from an initial condition P (0). Note that we can
start either deterministically from a given state z i (by setting
P (0) = ei ) or stochastically from any state according to a
generic vector of probabilities P (0). To define the problem
let us call T the target set, that is, the set of all indices i
associated with a state z i in the Markov chain (26) that satisfies
the desired property. Note that this set might be of infinite
size. We restrict our analysis to external signals satisfying
Assumption 3, so that we can map the external signal u to the
switching signal σ, as detailed in Remark 2. For a fixed signal
σ the solution of (27) immediately allows one to compute the
probability that the state at time T belongs to T , and thus
has the desired property, as PT (σ) := 1>
T P (T ; σ) where 1T
is an infinite vector that has the ith component equal to 1 if
i ∈ T and 0 otherwise. Our objective is to select the switching
signal σ(t) (and thus the external signal u(t)) that maximizes
the probability PT (σ).3 That is, we aim at solving
PT? := max PT (σ),
σ∈SIK

σ ? := arg max PT (σ),
σ∈SIK

(38)

where I is the cardinality of Σd as by Remark 2. Note that
PT (σ) in (38) is computed according to P (T ; σ) which is an
infinite dimentional vector. In the next theorem we show how
to overcome this issue and approximately solve (38) by using
the FSP approach of Proposition 3 and the reformulation as
MILP given in Proposition 2. To this end, let
σ̄ ? := arg max P̄T (σ).
σ∈SIK

(39)

where P̄T (σ) := 1̄>
T P̄J (T ; σ) is the probability that the final
state of the reduced Markov chain (28) belongs to T ∩ J at
time T given the switching signal σ, and 1̄T is a vector of
size |J| that has 1 in the positions corresponding to states of
J that belong also to T , and 0 otherwise.
Theorem 4. Suppose that Assumptions 3 and 4 hold. Then
PT (σ̄ ? ) ≥ PT? − 2ε.

On the other hand, when γ < 0
2ε
2ε
y2 ≤ ȳ2 + k¯l2 k∞ 1−ε
≤ γ d ȳ1 + v̄T (cd ) + k¯l2 k∞ 1−ε
≤ γ d y1 + v̄T (cd ) + (k¯l2 k∞ − γ d k¯l1 k∞ ) 2ε
1−ε

= γ d y1 + v̄T (cd ) + δ(γ d ).

Moreover (39) can be solved by solving the MILP in (17) for
system (28) with c = 1̄T and M = 1.

3 Note that one can use the same tools to maximize the probability of
avoiding a given set D by maximizing the probability of being in T = Dc .

10

Proof: Under Assumption 3 and 4, for any set T and any
signal σ, we get
P
P
P
P
1>
Pi ≤ i∈T ∩J Pi + i∈J
TP =
/ Pi ≤
i∈T ∩J Pi +ε
Pi∈T
P
≤ i∈T ∩J P̄i + i∈T ∩J |Pi − P̄i | + ε
P
≤ i∈T ∩J P̄i + kPJ − P̄J k1 + ε = 1̄>
T P̄ + 2ε,
P
P
P
>
and 1>
TP =
i∈T Pi ≥
i∈T ∩J Pi ≥
i∈T ∩J P̄i = 1̄T P̄ ,
where we used Assumption 4 and Proposition 3 and we
omitted (T ; σ) for simplicity. To sum up, for each σ,
P̄T (σ) ≤ PT (σ) ≤ P̄T (σ) + 2ε.
By imposing σ = σ ? we get PT? = PT (σ ? ) ≤ P̄T (σ ? )+2ε ≤
P̄T (σ̄ ? )+2ε. By imposing σ = σ̄ ? we get P̄T (σ̄ ? ) ≤ PT (σ̄ ? ).
Combining the last two inequalities we get the desired bound.
The last result can be proven as in Proposition 2. Note that P̄
is a vector of probabilities, hence we can set M = 1.
VII. T HE GENE EXPRESSION NETWORK CASE STUDY
To illustrate our method we consider again the gene expression model of Example 1 and determine what combinations of
the protein mean and variance are achievable starting from the
zero state, under different assumptions on the external signal.
A. Single input
Consider the gene expression model with one external
signal and reactions following the mass action kinetics, as
described in Example 2. In this case, the moments equations
are linear and the protein mean and variance can be obtained
by assuming as output matrix for the linear system (25)


0 1 0 0 0
L :=
.
0 0 0 0 1
Depending on the experimental setup, the external signal u(t)
may take values in the set Σd := {0, 1}, if the input is of
the ON-OFF type [1], [2], [3], [5], or in the interval Σc :=
[0, 1], if the input is continuous [4]. Corollary 1 guarantees
the validity of the following results both for Σd and Σc . The
problem of computing an outer approximation of the reachable
set of this system was studied in [27] using ad hoc methods. In
Fig. 4 we compare the outer approximation obtained therein
(magenta dashed/dotted line) with the inner (solid red) and
outer (dashed blue) approximations that we obtained using
the methods for linear moment equations of Section V-A1.
We used the parameters kr = 0.0236, γr = 0.0503, kp =
0.18, γp = 0.0121 (all in units of min−1 ) and set T = 360
min. Figure 4 shows that the outer approximation computed
using the hyperplane method is more accurate than the one
previously obtained in the literature. Moreover, since inner and
outer approximations practically coincide, this method allows
one to effectively recover the reachable set.
B. Single input and saturation
As second case study we consider again Example 2, but we
now assume that not all the reactions follow the laws of mass
action kinetics. Specifically, we are interested in investigating
how the reachable set changes if we assume that the number

Fig. 4.
Comparison of the inner (red solid) and outer (blue dashed)
approximations of the reachable set for the protein mean and variance,
according to model (25), computed using the hyperplane method and the
outer approximation computed according to [27] (magenta dashed/dotted).

of ribosomes in the cell is limited and consequently we impose
a saturation to the translation propensity. Following [39], we
assume that the translation rate follows the Michaelis-Menten
kinetics so that
α3 (kp , z) = k̃p ·

a·m
b+a·m

instead of α3 (kp , z) = kp · m.

For the simulations we impose k̃p = 0.7885, b = 0.06, a =
0.02, so that the maximum reachable protein mean is the
same as in the case without saturation analysed in the previous subsection. The corresponding propensity function is
illustrated in Fig. 5a). All the other propensities are assumed
as in Section VII-A. Note that in this case the propensities
are not affine. Consequently, we estimate the reachable set by
using the FSP approach derived in Theorem 3. Specifically
we consider as set J the indices corresponding to states with
less than 6 mRNA copies and 40 protein copies. By assuming
T = 360 min and that u can switch any 30 minutes in the
set Σd = {0, 1}, we obtain an error ε = 2.84 · 10−4 . Fig. 5b)
shows the comparison of the reachable sets obtained for the
cases with and without saturation. From this plot it emerges
that, for the chosen values of parameters, saturation leads to
a decrease of variability in the population.
a)

b)

Fig. 5. Comparison of the reachable set for the protein mean and variance, according
to the model in Example 2, when all the reactions follows the mass action kinetics (as
in Fig 4) and when the translation is saturated. Figure a): h3 (z) when the translation
reaction follows the mass action kinetics (dashed blue) or the Michaelis Menten kinetics
(dashed dotted green). Figure b): comparison of the outer approximations of the reachable
sets in the two cases. The blue dashed line is as in Fig 4. The grey line is the outer
approximation of the reachable set of the FSP system (28), the green dashed dotted line
is the outer approximation of the original system (27) according to Theorem 3.

11

a)

C. Fluorescent protein and the two inputs case
Consider again Example 1, but now assume that:
1) mRNA production and degradation can both be controlled, so that the vector of propensities is α(z)h = [kri·
(t)
u1 (t), γr ·m·u2 (t), kp ·m, γp ·p]> and u(t) := uu21 (t)
;
2) the protein P can mature into a fluorescent protein F
according to the additional maturation and degradation
reactions
α5 (kf ,z)

−−−−−→

P

F,

F

α6 (γp ,z)

−−−−−→

∅,

where α5 (kf , z) := kf ·p, α6 (γp , z) := γp ·f and kf >
0 is the maturation rate. For simplicity, the degradation
rate of F is assumed to be the same as that of P ;
3) the fluorescence intensity I(t) of each cell can be measured and is proportional to the amount of fluorescence
protein, that is, I(t) = rF (t) for a fixed scaling parameter
r > 0.
Since all the propensities are affine, the system describing the
evolution of means and variances of the augmented network
is
ẋ≤2 (t) = Af (u(t))x≤2 (t) + bf (u(t)),
(40)

b)

c)

where the state vector x≤2 (t) and Af (u(t)), bf (u(t)) are
x≤2= [E[M ], E[P ], E[F ], V[M, P ], V[M, F ], V [P ], V[P, F ], V [F ]]>
d (u (t))
0
0
0
0
0 0 0
1

2

kp
0
kp
0
kp
0
0





f 
A =




d2
γp
0
0
(γp + kf )
−γp
γp


bf = kr u1 (t)

0

0
0
0
d3
0
0
0 d4 (u2 (t))
0
0
γp
d5 (u2 (t))
0
2kp
0
0
0
kp
kf
0
0

0 0
0 0
0 0
0 0
d6 0
γp d7
0 2γp

0

0

0

0

0
0

0
,
0

0

0
d8
>
0 0
Fig. 6.

with d1 (u2 (t)) = −γr u2 (t), d2 = −(γp + kf ), d3 = −kf ,
d4 (u2 (t)) = −(γr u2 (t) + γp + kf ), d5 (u2 (t)) = −(γr u2 (t) + kf ),
d6 = −2(γp + kf ), d7 = −(2kf + γp ), d8 = −2kf .

System (40) depends on the parameter vector θ =
[kr , γr , kp , γp , kf , r] (for more details see [40, Supplementary
Information pg. 16]). For the parameters we use the MAP
estimates identified in [28] (all in min−1 )
kr
kf

= 0.0236
= 0.0212

γr = 0.0503 kp = 178.398
γp = 0.0121 r−1 = 646.86

(41)

and we set
f

L :=



0
0

0
0

r
0

0 0 0 0
0 0 0 0

0
r2


,

(42)

to compute the mean and variance reachable set for the
fluorescence intensity.
Our first aim is to compare the reachable set of such extended model with experimental data, when only one external
signal (“1in”) is available. In the case of one input, (40) is a
linear system and the methods of Section V-A1 can be applied.
Fig. 6a) shows the estimated reachable set compared with the
real data collected in [2].
Our second goal is to investigate how the reachable set
changes when both mRNA production and degradation are

Output reachable set of system (40) with output as in (42) and parameters as
in (41). Figure a) [1 external signal]: Comparison between the inner (red contour) and
outer (blue lines) approximation of the output reachable set, when the set of possible
modes is Σ1in , and the measured data. Different colors refers to data collected in different
experiments. Figure b) [2 external signals]: Outer approximation of the output reachable
set, when the set of possible modes is Σ2in . The two green dots represent the outputs
when u(t) = [1, 0.5]> ∀t and u(t) = [1, 1]> ∀t, respectively. The black crosses
represent the output for random signals in Σ2in . Figure c) [Comparison]: The red solid
line is the outer approximation obtained for u(t) ∈ Σ1in , as in Fig. a), the blue dashed
line the one for u(t) ∈ Σ2in , as in Fig. b).

controlled (“2in”), as studied in [41]. Note that in this case,
system (40) is nonlinear. We therefore set T = 300 min
and assume that switchings can occur every 20 min, so that
Assumption 3 is satisfied with K = 15 and use the hyperplane
method as described in Section V-A2 with input sets
  
   

0
0
1
1
Σ2in :=
,
,
,
, so that I = 4,
1
0.5
1
0.5
   
0
1
Σ1in :=
,
, so that I = 2,
1
1
respectively. Note that we set the minimum input for the
mRNA degradation to 0.5 > 0 to avoid unboundedness. With
these input choices it is intuitive that the largest possible state
is reached when the mRNA production is at its maximum and
the mRNA degradation is at its minimum. Therefore, in the
1 ] ∀t)
MILPs we can use the bounds M = x (T ; 0, u(t) = [ 0.5

12

for the case of two inputs and M = x(T ; 0, u(t) = [ 11 ] ∀t), for
the case of one input. Fig. 6b) shows the output reachable set
for the case of two inputs. The simulation time for computing
the outer approximation with the hyperplane method was 5.6
hrs. Computing the exact reachable set by simulating all the
possible switching signals, assuming that one simulation takes
10−4 sec and neglecting the time needed to enumerate all
possible signals, would take 29.8 hrs. The black crosses in
Fig. 6b) are obtained by simulating the output of the system for
5000 randomly constructed input signals. This simulation
illustrates that random approaches might lead to significantly
under estimate the reachable set. Fig. 6c) shows a comparison
of the reachable sets obtained in Fig. 6a) and b) when the
input set is Σ1in and Σ2in , respectively.
VIII. C ONCLUSION
In the paper we have: i) proposed a method to approximate
the projected reachable set of switched affine systems with
fixed switching times, ii) extended the FSP approach to
controllable networks, iii) illustrated how these new theoretical
tools can be used to analyse generic networks both from a
population and single cell perspective and iv) provided an
extensive gene expression case study using both in silico
and in vivo data. Even though our analysis is motivated by
biochemical reaction networks, our results can actually be
applied to study the moments of any Markov chain with
transitions rates that switch among I possible configurations
at K fixed instants of times. Our results hold both in case
of finite and infinite state space. Moreover, while we have
assumed here that cells are identical, we showed in [29] that
also in the case of heterogeneous population one can derive
equations describing the moments evolution. The reachable set
of such populations can be obtained, as described in this paper,
by applying Corollary 1 or 2 to such system.
R EFERENCES
[1] 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,
vol. 29, pp. 1114–1116, 2011.
[2] J. Ruess, F. Parise, A. Milias-Argeitis, M. Khammash, and J. Lygeros,
“Iterative experiment design guides the characterization of a lightinducible gene expression circuit,” National Academy of Sciences of the
USA, vol. 112, no. 26, pp. 8148–8153, 2015.
[3] E. J. Olson, L. A. Hartsough, B. P. Landry, R. Shroff, and J. J.
Tabor, “Characterizing bacterial gene circuit dynamics with optically
programmed gene expression signals.” Nature Methods, vol. 11, pp.
449–455, 2014.
[4] J. Uhlendorf, A. Miermont, T. Delaveau, G. Charvin, F. Fages, S. Bottani, G. Batt, and P. Hersen, “Long-term model predictive control of gene
expression at the population and single-cell levels.” National Academy
of Sciences of the USA, vol. 109, no. 35, pp. 14 271–14 276, 2012.
[5] F. Menolascina, M. Di Bernardo, and D. Di Bernardo, “Analysis, design
and implementation of a novel scheme for in-vivo control of synthetic
gene regulatory networks,” Automatica, pp. 1265–1270, 2011.
[6] G. Batt, H. De Jong, M. Page, and J. Geiselmann, “Symbolic reachability
analysis of genetic regulatory networks using discrete abstractions,”
Automatica, vol. 44, no. 4, pp. 982–989, 2008.
[7] N. Chabrier and F. Fages, “Symbolic model checking of biochemical
networks,” in Computational Methods in Systems Biology, C. Priami,
Ed. Springer, 2003, pp. 149–162.
[8] D. T. Gillespie, “A rigorous derivation of the chemical master equation,”
Physica A: Statistical Mechanics and its Applications, pp. 404–425,
1992.

[9] I. Lestas, G. Vinnicombe, and J. Paulsson, “Fundamental limits on the
suppression of molecular fluctuations,” Nature, pp. 174–178, 2010.
[10] Z. Sun and S. S. Ge, Switched linear systems: Control and design.
Springer Science & Business Media, 2006.
[11] C. Altafini, “The reachable set of a linear endogenous switching system,”
Systems & control letters, vol. 47, no. 4, pp. 343–353, 2002.
[12] T. J. Graettinger and B. H. Krogh, “Hyperplane method for reachable
state estimation for linear time-invariant systems,” Journal of optimization theory and applications, vol. 69, no. 3, pp. 555–588, 1991.
[13] B. Munsky and M. Khammash, “The finite state projection algorithm for
the solution of the chemical master equation.” The Journal of chemical
physics, vol. 124, no. 4, p. 044104, 2006.
[14] C. Baier, B. Haverkort, H. Hermanns, and J.-P. Katoen, “Modelchecking algorithms for continuous-time Markov chains,” IEEE Trans.
on Software Engineering, vol. 29, no. 6, pp. 524–541, 2003.
[15] A. Abate, J.-P. Katoen, J. Lygeros, and M. Prandini, “Approximate model
checking of stochastic hybrid systems,” European Journal of Control,
vol. 16, no. 6, pp. 624–641, 2010.
[16] H. El-Samad, S. Prajna, A. Papachristodoulou, J. Doyle, and M. Khammash, “Advanced methods and algorithms for biological networks analysis,” IEEE, vol. 94, no. 4, pp. 832–853, 2006.
[17] R. Alur, C. Courcoubetis, T. A. Henzinger, and P.-H. Ho, “Hybrid
automata: An algorithmic approach to the specification and verification
of hybrid systems,” in Hybrid systems. Springer, 1993, pp. 209–229.
[18] R. Alur, T. A. Henzinger, G. Lafferriere, and G. J. Pappas, “Discrete
abstractions of hybrid systems,” Proceedings of the IEEE, vol. 88, no. 7,
pp. 971–984, 2000.
[19] X. D. Koutsoukos and P. J. Antsaklis, “Safety and reachability of piecewise linear hybrid dynamical systems based on discrete abstractions,”
Discrete Event Dynamic Systems, vol. 13, no. 3, pp. 203–243, 2003.
[20] R. Ghosh and C. Tomlin, “Symbolic reachable set computation of piecewise affine hybrid automata and its application to biological modelling:
Delta-notch protein signalling,” Systems Biology, pp. 170–183, 2004.
[21] L. Habets, P. J. Collins, and J. H. van Schuppen, “Reachability and
control synthesis for piecewise-affine hybrid systems on simplices,”
IEEE Trans. on Automatic Control, vol. 51, no. 6, pp. 938–948, 2006.
[22] A. Hamadeh and J. Goncalves, “Reachability analysis of continuoustime piecewise affine systems,” Automatica, pp. 3189–3194, 2008.
[23] I. M. Mitchell and J. A. Templeton, “A toolbox of Hamilton-Jacobi
solvers for analysis of nondeterministic continuous and hybrid systems,”
in Hybrid systems: Computation and control, M. Morari and L. Thiele,
Eds. Springer, 2005, pp. 480–494.
[24] A. Kurzhanski and I. Valyi, Ellipsoidal calculus for estimation and
control. Springer, 1997.
[25] A. Donzé and O. Maler, “Systematic simulation using sensitivity analysis,” in Hybrid Systems: Computation and Control, A. Bemporad,
A. Bicchi, and G. Buttazzo, Eds. Springer, 2007, pp. 174–189.
[26] E. Lakatos and M. Stumpf, “Control mechanisms for stochastic biochemical systems via computation of reachable sets,” bioRxiv, 2016.
[27] F. Parise, M. E. Valcher, and J. Lygeros, “On the reachable set of the
controlled gene expression system,” in IEEE CDC, 2014, pp. 4597–4604.
[28] ——, “ On the use of hyperplane methods to compute the reachable set
of controlled stochastic biochemical reaction networks,” in IEEE CDC,
2015, pp. 1259–1264.
[29] ——, “Reachability analysis for switched affine autonomous systems
and its application to controlled stochastic biochemical reaction networks,” in IEEE CDC, 2016.
[30] D. Liberzon, Calculus of variations and optimal control theory: A
concise introduction. Princeton University Press, 2011.
[31] H. J. Sussmann, “Lie Brackets, Real Analyticity and Geometric Control,”
in Differential Geometric Control Theory, R. W. Brockett, R. S. Millman,
and H. J. Sussmann, Eds. Birkhauser Boston, 1983, vol. 27, pp. 1–116.
[32] C. Seatzu, D. Corona, A. Giua, and A. Bemporad, “Optimal control
of continuous-time switched affine systems,” IEEE Trans. on Automatic
Control, vol. 51, no. 5, pp. 726–741, 2006.
[33] A. Bemporad and M. Morari, “Control of systems integrating logic,
dynamics, and constraints,” Automatica, pp. 407–427, 1999.
[34] D. J. Wilkinson, Stochastic modeling for systems biology. CRC, 2011.
[35] J. Goutsias and G. Jenkinson, “Markovian dynamics on complex reaction
networks,” Physics Reports, vol. 529, no. 2, pp. 199–264, 2013.
[36] J. Ruess, “Moment-based methods for the analysis and identification of
stochastic models of biochemical reaction networks.” Doctoral dissertation, ETH Zurich, 2014.
[37] J. Hespanha, “Moment closure for biochemical networks,” in IEEE
International Symposium on Communications, Control and Signal Processing, 2008, pp. 142–147.

13

[38] C. H. Lee, K.-H. Kim, and P. Kim, “A moment closure method for
stochastic reaction networks,” The Journal of Chemical Physics, vol.
130, no. 13, p. 134107, 2009.
[39] R. Brockmann, A. Beyer, J. J. Heinisch, and T. Wilhelm, “Posttranscriptional expression regulation: what determines translation rates?” PLoS
Comput Biol, vol. 3, no. 3, p. e57, 2007.
[40] J. Ruess, A. Milias-Argeitis, and J. Lygeros, “Designing experiments to
understand the variability in biochemical reaction networks.” Journal of
the Royal Society Interface, vol. 10, no. 88, p. 20130588, 2013.
[41] C. Briat and M. Khammash, “Computer control of gene expression:
Robust setpoint tracking of protein mean and variance using integral
feedback,” in IEEE CDC, 2012, pp. 3582–3588.

Francesca Parise Francesca Parise was born in
Verona, Italy, in 1988. She received the B.Sc. and
M.Sc. degrees (cum Laude) in Information and Automation Engineering from the University of Padova,
Italy, in 2010 and 2012, respectively. She conducted
her master thesis research at Imperial College London, UK, in 2012. She graduated from the Galilean
School of Excellence, University of Padova, Italy,
in 2013. She defended her PhD at the Automatic
Control Laboratory, ETH Zurich, Switzerland in
2016 and she is currently a Postdoctoral researcher
at the Laboratory for Information and Decision Systems, M.I.T., USA.
Her research focuses on identification, analysis and control of complex
systems, with application to distributed multi-agent networks and systems
biology.

Maria Elena Valcher Maria Elena Valcher received
the Laurea degree and the PhD from the University
of Padova, Italy. Since January 2005 she is full
professor of Automatica at the University of Padova.
She is author/co-author of 76 papers appeared
in international journals, 92 conference papers, 2
text-books and several book chapters. Her research
interests include multidimensional systems theory,
polynomial matrix theory, behavior theory, convolutional coding, fault detection, delay-differential
systems, positive systems, positive switched systems
and Boolean control networks.
She has been involved in the Organizing Committees and in the Program
Committees of several conferences. In particular, she was the Program Chair
of the CDC 2012. She was in the Editorial Board of the IEEE Transactions
on Automatic Control (1999-2002), Systems and Control Letters (2004-2010)
and she is currently in the Editorial Boards of Automatica (2006-today),
Multidimensional Systems and Signal Processing (2004-today), SIAM J. on
Control and Optimization (2012-today), European Journal of Control (2103today) and IEEE Access (2014-today).
She was Appointed Member of the CSS BoG (2003); Elected Member of
the CSS BoG (2004-2006; 2010-2012); Vice President Member Activities of
the CSS (2006-2007); Vice President Conference Activities of the CSS (20082010); CSS President (2015). She is presently a Distinguished Lecturer of the
IEEE CSS and IEEE CSS Past President. She received the 2011 IEEE CSS
Distinguished Member Award and she is an IEEE Fellow since 2012.

John Lygeros John Lygeros completed a B.Eng. degree in electrical engineering in 1990 and an M.Sc.
degree in Systems Control in 1991, both at Imperial
College of Science Technology and Medicine, London, U.K.. In 1996 he obtained a Ph.D. degree from
the Electrical Engineering and Computer Sciences
Department, University of California, Berkeley. During the period 1996-2000 he held a series of research
appointments at the National Automated Highway
Systems Consortium, Berkeley, the Laboratory for
Computer Science, M.I.T., and the Electrical Engineering and Computer Sciences Department at U.C. Berkeley. Between 2000
and 2003 he was a University Lecturer at the Department of Engineering,
University of Cambridge, U.K., and a Fellow of Churchill College. Between
2003 and 2006 he was an Assistant Professor at the Department of Electrical
and Computer Engineering, University of Patras, Greece. In July 2006 he
joined the Automatic Control Laboratory at ETH Zurich, first as an Associate
Professor, and since January 2010 as a Full Professor. Since 2009 he is
serving as the Head of the Automatic Control Laboratory and since 2015
as the Head of the Department of Information Technology and Electrical
Engineering. His research interests include modelling, analysis, and control of
hierarchical, hybrid, and stochastic systems, with applications to biochemical
networks, automated highway systems, air traffic management, power grids
and camera networks. He teaches classes in the area of systems and control
both at the undergraduate and at the graduate level at ETH Zurich, notably
the 4th semester class Signals and Systems II, which he delivers in a flipped
classroom format. John Lygeros is a Fellow of the IEEE, and a member of the
IET and the Technical Chamber of Greece. He served as an Associate Editor
of the IEEE Transactions on Automatic Control and in the IEEE Control
Systems Society Board of Governors; he is currently serving as the Treasurer
of the International Federation of Automatic Control.

