Shaping Pulses to Control Bistable
Monotone Systems Using Koopman
Operator
Aivar Sootla ∗ Alexandre Mauroy ∗∗ Jorge Goncalves ∗∗∗

arXiv:1605.06252v1 [cs.SY] 20 May 2016

∗

Montefiore Institute, University of Liège, B-4000 Liège, Belgium
(e-mail: asootla@ulg.ac.be).
∗∗
Luxembourg Centre for Systems Biomedicine, University of
Luxembourg, L-4367 Belvaux, Luxembourg (e-mail:
alexandre.mauroy@uni.lu)
∗∗∗
Luxembourg Centre for Systems Biomedicine, University of
Luxembourg, L-4367 Belvaux, Luxembourg (e-mail:
jorge.goncalves@uni.lu)

Abstract: In this paper, we further develop a recently proposed control method to switch
a bistable system between its steady states using temporal pulses. The motivation for using
pulses comes from biomedical and biological applications (e.g. synthetic biology), where it is
generally difficult to build feedback control systems due to technical limitations in sensing and
actuation. The original framework was derived for monotone systems and all the extensions
relied on monotone systems theory. In contrast, we introduce the concept of switching function
which is related to eigenfunctions of the so-called Koopman operator subject to a fixed control
pulse. Using the level sets of the switching function we can (i) compute the set of all pulses that
drive the system toward the steady state in a synchronous way and (ii) estimate the time needed
by the flow to reach an epsilon neighborhood of the target steady state. Additionally, we show
that for monotone systems the switching function is also monotone in some sense, a property
that can yield efficient algorithms to compute it. This observation recovers and further extends
the results of the original framework, which we illustrate on numerical examples inspired by
biological applications.
1. INTRODUCTION

optimally) but are also simple enough to be implemented
in an experimental setting.

In many applications, the use of a time-varying feedback
control signal is impeded by the limitations in sensing and
actuation. One of such applications is synthetic biology,
which aims to engineer and control biological functions in
living cells (Brophy and Voigt (2014)), and which is an
emerging field of science with applications in metabolic
engineering, bioremediation and energy sector (Purnick
and Weiss (2009)). Recently, control theoretic regulation
of protein levels in microbes was shown to be possible
by Milias-Argeitis et al. (2011); Menolascina et al. (2011);
Uhlendorf et al. (2012). However, the proposed methods are hard to automate due to physical constraints in
sensing and actuation (for example, using the techniques
from Levskaya et al. (2009); Mettetal et al. (2008)). In the
context of actuation, adding a chemical solution to the
culture is fairly straightforward, but in contrast, removing
a chemical from the culture is much more complicated
(this could be done through diluting, but would be time
consuming and hard to perform repeatedly). In regard to
these constraints, it is therefore desirable to derive control
policies which can not only solve a problem (perhaps not

The pioneering development in synthetic biology was the
design of the so-called genetic toggle switch by Gardner
et al. (2000), which is a synthetic genetic system (or a
circuit) of two mutually repressive genes LacI and TetR.
Mutual repression means that only one of the genes can
be activated or switched “on” at a time. The activated
gene expresses proteins within a cell, hence the number of
proteins expressed by the “on” gene is much higher than
the number of proteins of the “off” gene. This entails the
possibility of modeling this genetic circuit by a bistable
dynamical system. Since toggle switches serve as one of the
major building blocks in synthetic biology, we set up the
control problem of switching from one stable fixed point to
another (or toggling a gene). Recently, Sootla et al. (2015,
2016) proposed to solve the problem using temporal pulses
u(t) with fixed length τ and magnitude µ:

1 0≤t≤τ,
u(t) = µh(t, τ )
h(t, τ ) =
(1)
0 t>τ.

1 A. Sootla holds an F.R.S–FNRS fellowship. This work was performed when A. Mauroy was with the University of Liège and held
a return grant from the Belgian Science Policy (BELSPO).

In the case of monotone systems (cf. Angeli and Sontag
(2003)), the set of all pairs (µ, τ ) allowing a switch (i.e. the
switching set) was completely characterized. In particular,
the boundary of this set, called the switching separatrix,
was shown to be monotone, a result which significantly
simplifies the computation of the switching set. However,

the contributions of Sootla et al. (2016) provide only a
binary answer (on whether a given control pulse switches
the system or not), but do not characterize the time needed
to converge to the steady state.
In this paper, we conduct a theoretical study which extends the results by Sootla et al. (2016) and provides
a temporal characterization of the effects of switching
pulses. To do so, we exploit the framework of the so-called
Koopman operator (cf. Mezić (2005)), which is a linear
infinite dimensional representation of a nonlinear dynamical system. In particular, we use the spectral properties
of the operator, focusing on the Koopman eigenfunctions
(i.e. infinite dimensional eigenvectors of the operator). We
first introduce the switching function, which we define as
a function of µ and τ related to the dominant Koopman
eigenfunction. Each level set of the switching function
characterizes a set of pairs (µ, τ ) describing control pulses
that drive the system synchronously to the target fixed
point. Hence, the switching function provides a temporal
characterization of the controlled trajectories. The switching separatrix introduced by Sootla et al. (2016) is interpreted in this framework as a particular level set of the
switching function. Furthermore, there is a direct relationship between the level sets of the switching function and
the so-called isostables introduced in Mauroy et al. (2013).
Since the switching function is defined through a Koopman
eigenfunction, it can be computed in the Koopman operator framework with numerical methods based on Laplace
averages. These methods can be applied to a very general
class of systems, but usually require extensive simulations.
However, we show that the switching function of monotone
systems is also monotone in some sense, so that its level
sets can be computed in a very efficient manner by using
the algorithm proposed in Sootla and Mauroy (2016b).
The key to reducing the computational complexity is to
exploit the properties of the Koopman eigenfunctions of a
monotone system.
The main contribution of this paper is to provide a theoretical framework that relates the Koopman operator to control problems. We note, however, that the eigenfunctions
of the Koopman operator can be estimated directly from
the observed data using dynamic mode decomposition
methods (cf. Schmid (2010); Tu et al. (2014)). Therefore
our results could potentially be extended to a data-based
setting, which would increase their applicability.
The rest of the paper is organized as follows. In Section 2,
we cover some basics of monotone systems theory and
Koopman operator theory. In Section 3, we review the
shaping pulses framework from Sootla et al. (2016) and
present the main results of this paper. We illustrate the
theoretical results on examples in Section 4.
2. PRELIMINARIES
Consider control systems in the following form
ẋ = f (x, u), x(0) = x0 ,
(2)
with f : D × U → Rn , u : R≥0 → U, and where
D ⊂ Rn , U ⊂ R and u belongs to the space U∞ of
Lebesgue measurable functions with values from U. We
define the flow map φ : R×D×U∞ → Rn , where φ(t, x0 , u)
is a solution to the system (2) with an initial condition

x0 and a control signal u. If u = 0, then we call the
system (2) unforced. We denote the Jacobian matrix of
f (x, 0) as J(x). If x∗ is a fixed point of the unforced
system, we assume that the eigenvectors of J(x∗ ) are
linearly independent, for the sake of simplicity. We denote
the eigenvalues of J(x∗ ) by λi .
Koopman Operator. Spectral properties of nonlinear dynamical systems can be described through an operatortheoretic framework that relies on the so-called Koopman
operator L = f T ∇, which is an operator acting on the
functions g : Rn → C (also called observables). We
limit our study of the Koopman operator to unforced
systems (2) (that is, with u = 0) on a basin of attraction
B ⊂ Rn of a stable hyperbolic fixed point x∗ (that is,
the eigenvalues λj of the Jacobian matrix J(x∗ ) are such
that <(λj ) < 0 for all j). In this case, the Koopman
operator admits a point spectrum and the eigenvalues λj
of the Jacobian matrix J(x∗ ) are also eigenvalues of the
Koopman operator. In the non-hyperbolic case, the analysis is more involved since the spectrum of the Koopman
operator may be continuous. The operator L generates a
semigroup acting on observables g
U t g(x) = g ◦ φ(t, x, 0),
(3)
where ◦ is the composition of functions and φ(t, x, 0) is
a solution to the unforced system for x ∈ B. Since the
operator is linear (cf. Mezic (2013)), it is natural to study
its spectral properties. In particular, the eigenfunctions
sj : B 7→ C of the Koopman operator are defined as the
functions satisfying Lsj = f T ∇sj = λj sj , or equivalently
U t sj (x) = sj (φ(t, x, 0)) = sj (x) eλj t ,
where λj ∈ C is the associated eigenvalue.

x ∈ B,

(4)

If the vector field f is a C 2 function, then the eigenfunctions sj are C 1 functions (Mauroy and Mezic (2016)). If
the vector field f is analytic and if the eigenvalues λj are
simple, the flow of the system can be expressed through
the following expansion (see e.g. Mauroy et al. (2013)):
n
X
φ(t, x, 0) = x∗ +
sj (x)vj eλj t +
(5)
j=1

X

vk1 ,...,kn sk11 (x) · · · sknn (x)e(k1 λ1 +...kn λn )t ,

k1 ,...,kn ∈N0
k1 +···+kn >1

where N0 is the set of nonnegative integers, λj , vj are
the eigenvalues and right eigenvectors of the Jacobian
matrix J(x∗ ), respectively, and the vectors vk1 ,...,kn are
the Koopman modes (see Mezić (2005); Mauroy and Mezic
(2016) for more details). Note that a similar expansion can
also be obtained if the eigenvalues are not simple (cf. Mezic
(2015)).
Throughout the paper we assume that λj are such that
0 > <(λ1 ) > <(λj ) for all j ≥ 2. In this case, the
eigenfunction s1 , which we call a dominant eigenfunction,
can be computed through the so-called Laplace average
ZT
1
∗
gλ (x) = lim
(g ◦ φ(t, x, 0))e−λt dt.
(6)
t→∞ T
0

For all g that satisfy g(x∗ ) = 0 and ∇g(x∗ ) · v1 6= 0, the
Laplace average gλ∗1 is equal to s1 (x) up to a multiplication
with a scalar. If we let g(x) = w1T (x − x∗ ), where w1 is the

left eigenvector of J(x∗ ) corresponding to λ1 , the limit
in (6) does not converge if x 6∈ B. Therefore, we do not
require the knowledge of B in order to compute s1 . The
other eigenfunctions sj (x) are generally harder to compute
using Laplace averages. The eigenfunction s1 can also be
estimated with linear algebraic methods (cf. Mauroy and
Mezic (2016)), or obtained directly from data by using
dynamic mode decomposition methods (cf. Schmid (2010);
Tu et al. (2014)).
The eigenfunction s1 (x) captures the dominant (i.e.
asymptotic) behavior of the unforced system. Hence the
boundaries ∂B α of the sets B α = {x | |s1 (x)| ≤ α}, which
are called isostables, are important for understanding the
dynamics of the system. It can be shown that trajectories
with initial conditions on the same isostable ∂B α1 converge
synchronously toward the fixed point, and reach other
isostables ∂B α2 , with α2 < α1 , after a time
 
1
α1
T =
ln
.
(7)
|<(λ1 )|
α2
In particular, for λ1 ∈ R, it follows directly from (5)
that the trajectories starting from ∂B α share the same
asymptotic evolution
φ(t, x, 0) → x∗ + v1 αeλ1 t , t → ∞ .
Note that isostables could also be defined when the system
is driven by an input u 6= 0, but they are here considered
only to describe the dynamics of the unforced system. A
more rigorous definition of isostables and more details can
be found in (Mauroy et al. (2013)).
In the case of bistable systems characterized by two
equilibria x∗ and x• with basins of attraction B(x∗ ) and
B(x• ), respectively, the Koopman operator admits two sets
of eigenfunctions s∗k and s•k . The eigenfunctions s∗k (resp.
s•k ) are related to the asymptotic convergence toward x∗
(resp. x• ). The dominant eigenfunctions s∗1 and s•1 define
two families of isostables, each of which is associated with
one equilibrium and lies in the corresponding basin of
attraction.
Monotone Systems and Their Spectral Properties. We will
study the properties of the system (2) with respect to a
partial order induced by positive cones in Rn . A set K is a
positive cone if R≥0 K ⊆ K, K + K ⊆ K, K ∩ −K ⊆ {0}. A
relation ∼ is called a partial order if it is reflexive (x ∼ x),
transitive (x ∼ y, y ∼ z implies x ∼ z), and antisymmetric
(x ∼ y, y ∼ x implies x = y). We define a partial order
K through a cone K ∈ Rn as follows: x K y if and only
if x − y ∈ K. We write x 6K y, if the relation x K y
does not hold. We also write x K y if x K y and x 6= y,
and x K y if x − y ∈ int(K). Similarly we can define a
partial order on the space of signals u ∈ U∞ : u K v if
u(t) − v(t) ∈ K for all t ≥ 0.
Systems whose flows preserve a partial order relation
K are called monotone systems. We have the following
definition.
Definition 1. The system ẋ = f (x, u) is called monotone
with respect to the cones Kx , Ku if φ(t, x, u) Kx φ(t, y, v)
for all t ≥ 0, and for all x Kx y, u Ku v.
The properties of monotone systems require additional
definitions. A function g : Rn → R is called increasing
with respect to the cone K if g(x) ≤ g(y) for all x K y.

Let [x, y]K denote the order-interval defined as [x, y]K =
{z|x K z K y}. A set A is called order-convex if, for all
x, y in A, the interval [x, y]K is a subset of A. A set M is
called p-convex if, for every x, y in M such that x K y and
every λ ∈ (0, 1), we have that λx + (1 − λ)y ∈ M . Clearly,
order-convexity implies p-convexity. If K = Rn≥0 we say
that the corresponding partial order is standard. Without
loss of generality, we will only consider the standard partial
order throughout the paper.
Proposition 2. (Angeli and Sontag (2003)). Consider the
control system (2), where the sets D, U are p-convex and
f ∈ C 1 (D×U). Then the system (2) is monotone on D×U∞
with respect to Rn≥0 , Rm
≥0 if and only if
∂fi
≥ 0, ∀ i 6= j, (x, u) ∈ cl(D) × U
∂xj
∂fi
≥ 0, ∀ i, j, (x, u) ∈ D × U.
∂uj
A generalization of this result to other cones can be
found in Angeli and Sontag (2003). We finally consider
the spectral properties of monotone systems that are
summarized in the following result. The proof can be found
in Sootla and Mauroy (2016b).
Proposition 3. Consider the system ẋ = f (x) with f ∈ C 2 ,
which admits a stable hyperbolic fixed point x∗ with a
basin of attraction B. Assume that <(λ1 ) > <(λj ) for all
j ≥ 2. Let v1 be a right eigenvector of the Jacobian matrix
J(x∗ ) and let s1 be an eigenfunction corresponding to λ1
(with v1T ∇s1 (x∗ ) = 1). If the system is monotone with
respect to Rn≥0 on int(B), then λ1 is real and negative.
Moreover, there exist s1 (·) and v1 such that s1 (x) ≥ s1 (y)
for all x, y ∈ B satisfying x  y, and v1  0.
This result shows that the sets Bα = {x||s1 (x)| ≤ α}
are order-convex for any α > 0 (cf. Sootla and Mauroy
(2016b)).
3. SHAPING PULSES TO SWITCH BETWEEN
FIXED POINTS
In this paper we consider the problem of switching between
two stable fixed points by using temporal pulses (1). We
formalize this problem by making the following assumptions:
A1. Let f (x, u) in (2) be continuous in (x, u) and C 2 in x
for every fixed u on Df × U.
A2. Let the unforced system (2) have two stable hyperbolic fixed points in Df , denoted by x∗ and x• , and
let Df = cl(B(x∗ ) ∪ B(x• )).
A3. For any u = µh(·, τ ) with finite µ and τ let φ(t, x∗ , u)
belong to int(Df ). Moreover, let there exist µ > 0,
τ > 0 such that lim φ(t, x∗ , µh(·, τ )) = x• .
t→∞

Assumption A1 guarantees existence and uniqueness of
solutions, while Assumption A2 defines a bistable system.
Note that in Sootla et al. (2016), the assumptions A1–A2
are less restrictive. That is, f (x, u) is Lipschitz continuous
in x for every fixed u, and the fixed points are asymptotically stable. Our assumptions are guided by the use of the
Koopman operator. Assumptions A1 and A2 guarantee
the existence of eigenfunctions s∗1 (x) and s•1 (x) that are
continuously differentiable on each basin of attraction.

Assumption A3 is technical and ensures that the switching
problem is feasible.
The goal of our control problem is to characterize the socalled switching set S defined as
n
o
S = (µ, τ ) ∈ R2>0 ∀t > τ : φ(t, x∗ , µh(·, τ )) ∈ B(x• ) .
(8)
It is shown in Sootla et al. (2016) that the set S is simply
connected and order-convex under some assumptions, a
property which is useful to obtain an efficient computational procedure. In particular, the boundary ∂S, called
the switching separatrix, is such that for all (µ, τ ), (ν, ξ)
in ∂S we cannot have that µ > ν and τ > ξ. The following
result sums up one of the theoretical contribution in Sootla
et al. (2016).
Proposition 4. Let the system ẋ = f (x, u) satisfy Assumptions A1–A3. The following conditions are equivalent:
(i) the set S is order-convex and simply connected;
(ii) let φ(τ1 , x∗ , µ1 h(·, τ1 )) belong to B(x• ), then the flow
φ(τ2 , x∗ , µ2 h(·, τ2 )) ∈ B(x• ) for all µ2 ≥ µ1 , τ2 ≥ τ1 .
Moreover, if the system ẋ = f (x, u) is monotone with
respect to K ×R on D ×U∞ and satisfies Assumptions A1–
A3, then (i)-(ii) hold.
We also note that the results in Sootla et al. (2016) were
extended to account for parametric uncertainty in the
vector field under additional constraints. In particular, it
is possible to estimate bounds on the switching set S.
Now, we proceed by providing an operator-theoretic point
of view on shaping pulses, which allows to study rates of
convergence to the fixed point.
Definition 5. Let S ⊆ R2>0 be a set of (µ, τ ) such that
φ(τ, x∗ , µ) ∈ B(x• ). We define the switching function
r : S 7→ C by
r(µ, τ ) = s•1 (φ(τ, x∗ , µ))
for all (µ, τ ) such that φ(τ, x∗ , µ) ∈ B(x• ).
The level sets ∂S α of |r| defined as
n
o
∂S α = (µ, τ ) ∈ S |r(µ, τ )| = α ,

α≥0

α

are reminiscent of the isostables ∂B , which are the level
sets of |s1 |. We also consider the sublevel sets of |r|
n
o
S α = (µ, τ ) ∈ R2>0 |r(µ, τ )| ≤ α
and it is straightforward toSshow that the switching set
S in (8) is equal to S ∞ = α≥0 S α . The level sets ∂S α
can therefore be seen as a generalization of the switching
separatrix ∂S = ∂S ∞ .
The level sets ∂S α capture the pairs (µ, τ ) such that the
trajectories φ(t, x∗ , µh(·, τ )) reach the isostable ∂B α at
time t = τ and cross the same isostables for all t ≥ τ .
This implies that the trajectories with the pairs (µ, τ ) on
the same level sets ∂S α will take the same time to converge
towards the fixed point x• when the control is switched off.
We can for instance estimate the time T needed to reach
the set B ε for a positive ε. It follows from (7) that, for
(µ, τ ) ∈ ∂S α , we have
α
1
T (µ, τ, ε) =
ln
,
(9)
|λ1 |
ε

where a negative T means that the trajectory is inside the
set B ε at time t = τ . Hence, the quantity Ttot = T + τ
is the time it takes to reach B ε if T is nonnegative. For
small enough ε, the function Ttot (µ, τ, ε) approximates the
amount of time required to reach a small neighborhood of
the fixed point x• .
In order to compute the switching function r, we can again
employ Laplace averages
ZT
1
r(µ, τ ) = lim
g ◦ φ(t, φ(τ, x∗ , µ), 0)e−λ1 t dt
T →∞ T
0

1
T →∞ T

ZT

= lim

g ◦ φ(t, x∗ , µh(·, τ ))e−λ1 (t−τ ) dt, (10)

τ

where λ1 is the dominant Koopman eigenvalue and g
satisfies g(x• ) = 0 and v1T ∇g(x• ) 6= 0. Note again that
the limit does not converge unless φ(τ, x∗ , µ) belongs to
B(x• ).
The computation of S α is not an easy task in general,
but certainly possible. However, additional assumptions on
the system simplify the computation of these sets. From
this point on we will assume that r(µ, τ ) has only real
values (i.e. s•1 ∈ R), which holds if the dominant Koopman
eigenvalue on B(x• ) is real (see Mauroy et al. (2013)). In
this case, the set ∂S α can be split into two sets
n
o
∂− S α = (µ, τ ) ∈ R2>0 r(µ, τ ) = −α ,
n
o
∂+ S α = (µ, τ ) ∈ R2>0 r(µ, τ ) = α .
If S α is order-convex (as it is shown below for the case of
monotone systems), then ∂− S α and ∂+ S α are the sets of
minimal and maximal elements of S α , respectively. That
is, if x  y for some x ∈ ∂+ S α (respectively, if x  y for
α
S) then y 6∈ S. This implies that ∂− S α and
some x ∈ ∂−
α
∂+ S are monotone maps, which significantly facilitates
computations of S α by applying the algorithm from Sootla
and Mauroy (2016b) with a minor modification.
Monotonicity also plays a role in the properties of the sublevel sets S α , as it does in the properties of the switching
separatrix. The main result of the section establishes that,
for monotone systems, the sets S α are order-convex and
∂S α are monotone maps.
Theorem 6. Let the system (2) satisfy Assumptions A1–
A3 and be monotone on D × U ∞ . Then
(i) the set S α is order-convex (with respect to the positive
orthant) for any non-zero α;
(ii) the set ∂+ S α is a monotone map, that is for all (µ1 , τ1 ),
(µ2 , τ2 ) ∈ ∂+ S α , if τ1 < τ2 then µ1 ≥ µ2 , and if µ1 < µ2
then τ1 ≥ τ2 . Moreover, the set ∂− S α is a graph of a
monotonically decreasing function for any finite non-zero
α;
(iii) if additionally φ(t, x, µ)  φ(t, x, ν) for all x, all
µ > ν ≥ 0 and all t > 0, then ∂− S α and ∂+ S α are graphs
of monotonically decreasing functions for any finite nonzero α.
The proof of Theorem 6 is in Appendix A. An interesting
detail is that the level sets ∂− S α are graphs of decreasing
functions. This implies that the switching separatrix ∂S ∞

Eight Species Generalized Repressilator. This system is an
academic example (cf. Strelkowa and Barahona (2010)),
where each of the species represses another species in a
ring topology. The corresponding dynamic equations for a
symmetric generalized repressilator are as follows:
p1
ẋ1 =
+ p4 − p5 x1 + u,
1 + (x8 /p2 )p3
p1
+ p4 − p5 x2 ,
ẋ2 =
(11)
1 + (x1 /p2 )p3
p1
ẋi =
+ p4 − p5 xi , ∀i = 3, . . . 8,
1 + (xi−1 /p2 )p3
where p1 = 100, p2 = 1, p3 = 2, p4 = 1, and p5 = 1. This
system has two stable equilibria x∗ and x• and is monotone
with respect to the cones Kx = Px R8 and Ku = R, where
Px = diag([1, − 1, 1, − 1, 1, − 1, 1, − 1]). We have
also x∗ Kx x• . It can be shown that the unforced system
is strongly monotone in the interior of R8≥0 for all positive
parameter values. One can also verify that there exist pulse
control signals u that switch the system from the state x∗
to the state x• .
The level sets ∂S α are depicted in Figure 1, where instead
of the values of the level sets we provide the time needed
to converge to ∂B 0.01 . As the reader may notice, the two
curves related to T = 5 (i.e. blue solid curves) lie close to
each other. They approximate the pairs (µ, τ ) that drive
the flow to the zero level set of s•1 (x). It is also noticeable
that the level sets ∂S α are less dense on the right of these
lines. This is explained by the fact that the flow is driven
by the pulse beyond the zero level set of s•1 (x) and has to
counteract the dynamics of the system.
The generalized repressilator is a monotone system, and
hence the premise of Theorem 6 is fulfilled. The level sets
∂S α in Figure 1 appear to be graphs of monotonically
decreasing functions, an observation which is consistent
with the claim of Theorem 6.
Toxin-antitoxin system. Consider the toxin-antitoxin system studied in Cataudella et al. (2013).
σT
1

Ṫ = 
−
T
[Af ][Tf ]
1+
(1 + β [T ]) (1 + βC [Tf ])
K0

M

f

σA

− ΓA A + u
Ȧ = 
[Af ][Tf ]
1 + K0
(1 + βM [Tf ])


[Af ][Tf ] [Af ][Tf ]2
ε[Ȧf ] = A − [Af ] +
+
KT
KT KT T


[Af ][Tf ]
[Af ][Tf ]2
ε[Ṫf ] = T − [Tf ] +
+2
,
KT
KT KT T
where A and T is the total number of toxin and antitoxin
proteins, respectively, while [Af ], [Tf ] is the number of free
toxin and antitoxin proteins. In Cataudella et al. (2013),
the authors considered the model with ε = 0. In order to
simplify our analysis we set ε = 10−6 . For the parameters

Pulse Magnitude

4. EXAMPLES

Time to reach
300

200

100

5

15

10

20

Pulse length

Fig. 1. The level sets of T =

1
|λ1 |

ln



|r(µ,τ )|
ε



, where

ε = 10−2 . The pairs (µ, τ ) on the same curve are
related to trajectories which converge synchronously
to the stable equilibrium.
Time to reach
40

Pulse Magnitude

can be approximated by a graph of a function by setting
α  0. We can also partially recover the results in Sootla
et al. (2016) by letting α → +∞. Note, however, that
∂− S ∞ is not necessarily a graph of a function, since strict
inequalities may no longer hold in the limit.

30

20

10

5

25

15

Pulse length

Fig. 2. The level sets of T =

1
|λ1 |

ln



|r(µ,τ )|
ε
−2

antitoxin system, where ε = 10



for the toxin-

.

σT = 166.28, K0 = 1, βM = βc = 0.16, σA = 102
ΓA = 0.2, KT = KT T = 0.3,
the system is bistable with two stable steady states:
x• = (27.1517 80.5151 58.4429 0.0877)
∗
x = (162.8103 26.2221 0.0002 110.4375) .
It can be verified that the system is not monotone with respect to any orthant, however, it was established in Sootla
and Mauroy (2016a) that it is eventually monotone. This
means that the flow satisfies the monotonicity property
after some initial transient.
We depict the level sets ∂S α in Figure 2, where it appears
that these sets are monotone curves although the system
does not satisfy the assumptions of Theorem 6. This could
be explained by the property of eventual monotonicity, but
we have not further investigated this case.
Lorenz System. Now we illustrate the level sets ∂S α in
the case where r is complex-valued. Consider the Lorenz
system
ẋ1 = σ(x2 − x1 ) + u
ẋ2 = x1 (ρ − x3 ) − x2 + u
ẋ3 = x1 x2 − βx3
with parameters σ = 10, ρ = 2, β = 8/3, which is
bistable but not monotone. Note that the Jacobian matrix
at the steady states has two complex conjugate dominant

Time to reach
Pulse Magnitude

20

15

10

5

0.5

1

1.5

2.5

2

3

Pulse length
1
|λ1 |
−2

Fig. 3. The level sets of T =
system, where ε = 10

ln



|r(µ,τ )|
ε



for the Lorenz

.

eigenvalues. In this case, Theorem 6 cannot be applied and
Figure 3 shows that the level sets ∂S α are not monotone. It
is however noticeable that the lower part of the switching
separatrix (black curve) seems to be monotone (but the
upper part is not monotone).
5. CONCLUSION
In this paper, we have further developed a recent study
on the problem of switching a bistable system between its
steady states with temporal pulses. We have introduced a
family of curves in the control parameter space, denoted
as ∂S α , which provide an information on the time needed
by the system to converge to the steady state. The sets
∂S α can be viewed as an extension of the switching
separatrix defined in the previous study. They are related
to the dominant eigenfunction of the Koopman operator,
a property that provides a method to compute them. In
the case of monotone systems, we have also shown that the
level sets ∂S α are characterized by strong (monotonicity)
properties.
Future research will investigate the topological properties
of the level sets ∂S α such as connectedness. Moreover,
characterizing the properties of the level sets ∂S α (and
the switching separatrix) in the case of non-monotone (e.g.
eventually monotone) systems is still an open question.
REFERENCES
Angeli, D. and Sontag, E. (2003). Monotone control
systems. IEEE Trans Autom Control, 48(10), 1684–
1698.
Brophy, J. and Voigt, C. (2014). Principles of genetic
circuit design. Nat methods, 11(5), 508–520.
Cataudella, I., Sneppen, K., Gerdes, K., and Mitarai, N.
(2013). Conditional cooperativity of toxin-antitoxin
regulation can mediate bistability between growth and
dormancy. PLoS Comput Biol, 9(8), e1003174.
Gardner, T., Cantor, C.R., and Collins, J.J. (2000). Construction of a genetic toggle switch in escherichia coli.
Nature, 403, 339–342.
Levskaya, A., Weiner, O.D., Lim, W.A., and Voigt, C.A.
(2009). Spatiotemporal control of cell signalling using a
light-switchable protein interaction. Nature, 461, 997–
1001.
Mauroy, A. and Mezic, I. (2016).
Global stability analysis using the eigenfunctions of the Koop-

man operator.
IEEE Tran Autom Control.
doi:
10.1109/TAC.2016.2518918. In press.
Mauroy, A., Mezić, I., and Moehlis, J. (2013). Isostables,
isochrons, and Koopman spectrum for the action–angle
representation of stable fixed point dynamics. Physica
D, 261, 19–30.
Menolascina, F., Di Bernardo, M., and Di Bernardo, D.
(2011). Analysis, design and implementation of a novel
scheme for in-vivo control of synthetic gene regulatory
networks. Automatica, Special Issue on Systems Biology,
47(6), 1265–1270.
Mettetal, J.T., Muzzey, D., Gomez-Uribe, C., and van
Oudenaarden, A. (2008). The Frequency Dependence of
Osmo-Adaptation in Saccharomyces cerevisiae. Science,
319(5862), 482–484.
Mezić, I. (2005). Spectral properties of dynamical systems,
model reduction and decompositions. Nonlinear Dynam,
41(1-3), 309–325.
Mezic, I. (2013). Analysis of fluid flows via spectral
properties of the Koopman operator. Annual Review
of Fluid Mechanics, 45, 357–378.
Mezic, I. (2015). On applications of the spectral theory of
the koopman operator in dynamical systems and control
theory. In IEEE Conf Decision Control, 7034–7041.
Milias-Argeitis, A., Summers, S., Stewart-Ornstein, J.,
Zuleta, I., Pincus, D., El-Samad, H., Khammash, M.,
and Lygeros, J. (2011). In silico feedback for in vivo
regulation of a gene expression circuit. Nat biotechnol,
29(12), 1114–1116.
Purnick, P. and Weiss, R. (2009). The second wave of
synthetic biology: from modules to systems. Nat. Rev.
Mol. Cell Biol., 10(6), 410–422.
Schmid, P.J. (2010). Dynamic mode decomposition of
numerical and experimental data. Journal of Fluid
Mechanics, 656, 5–28.
Sootla, A. and Mauroy, A. (2016a). Operator-theoretic
characterization of eventually monotone systems. Provisionally accepted for publication in Automatica.
http://arxiv.org/abs/1510.01149.
Sootla, A. and Mauroy, A. (2016b).
Properties of
isostables and basins of attraction of monotone systems.
In Proc Amer Control Conf (to appear).
http://arxiv.org/abs/1510.01153v2.
Sootla, A., Oyarzún, D., Angeli, D., and Stan, G.B. (2015).
Shaping pulses to control bistable biological systems. In
Proc Amer Control Conf, 3138 – 3143.
Sootla, A., Oyarzún, D., Angeli, D., and Stan, G.B. (2016).
Shaping pulses to control bistable systems: Analysis,
computation and counterexamples. Automatica, 63,
254–264.
Strelkowa, N. and Barahona, M. (2010). Switchable genetic oscillator operating in quasi-stable mode. J R Soc
Interface, 7(48), 1071–1082. doi:10.1098/rsif.2009.0487.
Tu, J.H., Rowley, C.W., Luchtenburg, D.M., Brunton,
S.L., and Kutz, J.N. (2014). On dynamic mode decomposition: Theory and applications. J Comput Dynamics,
1(2), 391 – 421.
Uhlendorf, J., Miermont, A., Delaveau, T., Charvin, G.,
Fages, F., Bottani, S., Batt, G., and Hersen, P. (2012).
Long-term model predictive control of gene expression
at the population and single-cell levels. Proc. Nat.
Academy Sciences, 109(35), 14271–14276.

Appendix A. PROOF OF THEOREM 6
Before we proceed with the proof of Theorem 6, we show
a technical result, which establishes that for monotone
systems the transient during switching between operating
points is always an increasing function.
Proposition 7. Let the system (2) be monotone on D×U∞ ,
then
φ(ξ + h, x∗ , µ)  φ(ξ, x∗ , µ)
(A.1)
for any nonnegative scalars h, ξ, µ.
Proof: The proof stems from a well-know result in monotone systems theory, which states that the flow cannot
increase (or decrease) on two disjoint time intervals. We
show this result for completeness. Due to monotonicity we
have
φ(h, x∗ , µ)  φ(h, x∗ , 0) = x∗
(A.2)
for all nonnegative h, µ. Therefore, the semigroup property
of the dynamical systems implies
φ(ξ+h, x∗ , µ) = φ(ξ, φ(h, x∗ , µ), µ)  φ(ξ, x∗ , µ), (A.3)
for any nonnegative scalars ξ, h, µ.
Proof of Theorem 6: (i) Let (µ1 , τ1 ), (µ2 , τ2 ) belong to
S α for some finite α > 0 and µ2 ≥ µ1 , τ2 ≥ τ1 . Then
r(µ1 , τ1 ), r(µ2 , τ2 ) are finite. Due to monotonicity and
Proposition 7, we have that
φ(τ1 , x∗ , µ1 h(·, τ1 ))  φ(τ1 , x∗ , µ2 h(·, τ2 )) 
φ(τ2 , x∗ , µ2 h(·, τ2 )),
which according to Proposition 3 implies that
s•1 (φ(τ1 , x∗ , µ1 h(·, τ1 ))) ≤ s•1 (φ(τ2 , x∗ , µ2 h(·, τ2 ))).
Using this property it is rather straightforward to show
that S α is order-convex.
(ii) The sets ∂+ S α and ∂− S α contain the maximal and
minimal elements, respectively, of the order-convex set S α .
Assume that (µ1 , τ1 ), (µ2 , τ2 ) ∈ ∂+ S α (or (µ1 , τ1 ), (µ2 , τ2 ) ∈
∂− S α ). Then we cannot have (µ1 , τ1 )  (µ2 , τ2 ). Hence,
if τ1 < τ2 , we must have µ1 ≥ µ2 and if µ1 <
µ2 , we must have τ1 ≥ τ2 . We prove the second
part of the statement by contradiction. Let τ1 < τ2 ,
µ1 ≤ µ2 and let s•1 (φ(τ1 , x∗ , µ1 h(·, τ1 ))) = −α and
s•1 (φ(τ2 , x∗ , µ2 h(·, τ2 ))) = −α, where α > 0. Due to monotonicity we have that
φ(t, x∗ , µ1 h(·, τ1 ))  φ(t, x∗ , µ2 h(·, τ2 )) ∀t ≥ 0,
which according to Proposition 3 entails
s•1 (φ(t, x∗ , µ1 h(·, τ1 ))) ≤ s•1 (φ(t, x∗ , µ2 h(·, τ2 ))). (A.4)
The flow φ(t, x∗ , µ1 h(·, τ1 )) converges to x• freely for all
t > τ1 , since h(t, τ1 ) = 0 for all t > τ1 . Negativity
of −α implies that s•1 (·) is growing along the trajectory
φ(τ1 , x∗ , µ1 h(·, τ1 )) and hence s•1 (φ(τ2 , x∗ , µ1 h(·, τ1 ))) >
−α. This, however, contradicts (A.4), since
s•1 (φ(τ2 , x∗ , µ1 h(·, τ1 ))) ≤ s•1 (φ(τ2 , x∗ , µ2 h(·, τ2 ))) = −α.
(iii) Let r(µ1 , τ1 ) = r(µ2 , τ2 ) = α and pick τ1 < τ2 . Assume
that µ1 ≤ µ2 . We have that
φ(τ1 , x∗ , µ1 h(·, τ1 ))  φ(τ1 , x∗ , µ2 h(·, τ1 )) 
φ(τ2 , x∗ , µ2 h(·, τ1 ))
where the first inequality follows from monotonicity and
the second follows from Proposition 7. Due to the condi-

tion φ(t, x, µ)  φ(t, x, ν) for all x, all µ > ν and all t > 0
in the premise, we have
φ(τ1 , x∗ , µ1 h(·, τ1 ))  φ(τ2 , x∗ , µ2 h(·, τ1 ))
= φ(τ2 − τ1 , φ(τ1 , x∗ , µ2 h(·, τ1 )), 0)
 φ(τ2 − τ1 , φ(τ1 , x∗ , µ2 h(·, τ1 )), µ2 )
= φ(τ2 , x∗ , µ2 h(·, τ2 ))
or equivalently r(µ1 , τ1 ) < r(µ2 , τ2 ). We arrive at a
contradiction, hence µ1 > µ2 , which implies that ∂+ S α
is a graph of a decreasing function.

