Properties of Isostables and Basins of Attraction of
Monotone Systems

arXiv:1510.01153v2 [math.OC] 22 Mar 2016

Aivar Sootla and Alexandre Mauroy

Abstract—In this paper, we investigate geometric properties
of monotone systems by studying their isostables and basins of
attraction. Isostables are boundaries of specific forward-invariant
sets defined by the so-called Koopman operator, which provides
a linear infinite-dimensional description of a nonlinear system.
First, we study the spectral properties of the Koopman operator
and the associated semigroup in the context of monotone systems.
Our results generalize the celebrated Perron-Frobenius theorem
to the nonlinear case and allow us to derive geometric properties
of isostables and basins of attraction. Additionally, we show that
under certain conditions we can characterize the bounds on the
basins of attraction under parametric uncertainty in the vector
field. We discuss computational approaches to estimate isostables
and basins of attraction and illustrate the results on two and four
state monotone systems.

I. I NTRODUCTION
In many applications, such as economics [1] or biology [2],
linear dynamical systems have states, which take only nonnegative values. For example, protein concentrations in biology
are always nonnegative. These systems are called positive
and have received a considerable attention in the context of
systems theory [3], [4], model reduction [5], [6], distributed
control [7], [8], etc. Analysis of such systems is facilitated by
employing the celebrated Perron-Frobenius theorem (cf. [9]),
which establishes strong spectral properties of the drift matrix
in a linear positive system.
In the nonlinear setting positive systems have a couple
of generalizations, with the most known being cooperative
monotone systems (cf. [2]). More specifically positive systems
generalize to cooperative with respect to the positive orthant
systems. Similarly to the linear case, these nonlinear systems
generate trajectories (or flows) which for every time t are
increasing functions in every argument with respect to the
initial state. With a slight abuse of notation, we will refer to
cooperative monotone systems simply as monotone. In [2], it
was briefly mentioned that the flow of a monotone system can
be seen as a positive operator. Hence the authors argued that
an operator extension of the Perron-Frobenius theorem, which
is called the Krein-Rutman theorem [10], can be applied.
However, the investigation into spectral properties of these
operators lacked due to absence of a well-developed theory
of spectral elements of such operators. This gap was filled by
Aivar Sootla is with Montefiore Institute, University of Liège, B-4000, Belgium asootla@ulg.ac.be. Alexandre Mauroy is with the Luxembourg
Centre for Systems Biomedicine, University of Luxembourg, L-4367 Belvaux,
Luxembourg alexandre.mauroy@uni.lu
A. Sootla holds an F.R.S–FNRS fellowship. Most of 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).

the development of the so-called Koopman operator and the
associated semigroup (cf. [11]–[13]).
The Koopman semigroup is a linear infinite-dimensional
representation of a nonlinear dynamical system. In many practical situations, the spectral elements of this semigroup such
as eigenfunctions (which can be seen as infinite dimensional
eigenvectors) and eigenvalues can be computed [14], [15]. As
in the linear case, the dominant eigenfunctions (i.e., eigenfunctions corresponding to the eigenvalue with the largest real
part) offer an insight into the dynamics of the system. In [14],
it is shown that the level sets of the dominant eigenfunction
(called isostables) contain the initial states of trajectories that
converge synchronously toward the equilibrium. The level
set at infinity is the boundary of the basin of attraction of
the equilibrium. To summarize, the dominant eigenfunctions
and eigenvalues contain the information about the asymptotic
behavior of the system.
In this paper, we first study the properties of the Koopman
semigroup associated with a monotone system. We show that
for stable systems monotone with respect to the nonnegative
orthant the dominant eigenfunction is an increasing function
in every argument for all states in the basin of attraction.
Hence we offer yet another version of the Perron-Frobenius
theorem for nonlinear systems (see [16] for other results on
the subject). We use our version of the Perron-Frobenius
theorem to study geometric concepts of isostables and basins
of attraction for monotone systems. Using this statement we
provide a straightforward proof of a known result stating that
a basin attraction of a monotone system is a union of the socalled order-intervals (cf. [17]–[19]). We refine this result by
showing that the sublevel sets of the dominant eigenfunctions
are unions of order-intervals, as well.
We proceed our theoretical development by considering
the basins of attraction of bistable monotone systems under
parametric uncertainty in the vector field. We show that if the
system is monotone with respect to parameter variations, then
it is possible to estimate inner and outer bounds on the basin
of attraction of the uncertain system. In this instance, instead
of the spectral theory we use monotone control systems theory.
Moreover, it appears that the behavior of the eigenfunctions
and isostables under parametric uncertainty is more complicated in comparison with the behavior of basins of attraction.
The rest of the paper is organized as follows. In Section II,
we cover the main properties of the Koopman operator and
monotone systems. In Subsection III-A we obtain spectral
properties of monotone systems, while deriving some properties of the isostables in Subsection III-B. We investigate the
behavior of basins of attraction of monotone systems under

parameter variations in Subsection III-C. We also discuss
different algorithms to compute isostables and the boundary of
the basins of attraction in Section IV. We conclude the paper
with numerical examples in Section V.
II. P RELIMINARIES

B. Koopman Operator
Spectral properties of nonlinear dynamical systems can be
described through an operator-theoretic framework that relies
on the so-called Koopman operator. The Koopman operator
associated with ẋ = f (x) is an operator acting on the functions
g : Rn → C (also called observables). The Koopman operator
generates the semigroup

Throughout the paper we consider parameter-dependent
systems in the following form
ẋ = f (x, p),

x(0) = x0 ,

(1)

where f : D×P → Rn , where D ⊂ Rn and P ⊂ Rm for some
integers n and m. We define the flow map φf : R × D × P →
Rn , where φf (t, x0 , p) is a solution to the system (1) with an
initial condition x0 and a parameter p. We assume that f (x, p)
is continuous in (x, p) on D ×P and Lipschitz in x on a every
compact subset of D for every fixed p. When it is clear from
the context, we will drop the parameter-dependence.

A. Monotonicity
We will study the properties of the system (1) with respect to
a partial order induced by cones. A set K is called a positive
cone if it is closed under addition and multiplication by a
nonnegative scalar and if we have −x 6∈ K for any x ∈ K. 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 will also write x K y if x K y and
x 6= y, and x K y if x − y ∈ int(K). We say that the
order is standard if K = Rn≥0 , which, with a slight abuse of
notation, we denote as  without a subscript. We call a set
{z|x K z K y} an order-interval induced by a cone K and
denote it as [x, y]K . A set A is called p-convex if for every x,
y in A such that x K y, and every λ ∈ (0, 1) we have that
λx + (1 − λ)y ∈ A. Systems whose flows preserve a partial
order relation K are called monotone systems.
Definition 1: The system is monotone with respect to the
cones Kx , Kp if φf (t, x, p) Kx φf (t, y, q) for all t ≥ 0,
and for all x Kx y and p Kp q. The system is strongly
monotone with respect to the cones Kx , Kp if it is monotone
and φf (t, x, p) Kx φf (t, y, q) holds for all t > 0 provided
x Kx y, p Kp q, and either x ≺Kx y or p ≺Kp q holds.
A certificate for monotonicity with respect to an orthant is
called Kamke-Müller conditions [20], where some generalizations of this result may also be found.
Proposition 1 ( [20]): Consider the system (1), where f is
differentiable in x and p and let the sets D, P be p-convex.
Then the system (1) is monotone on D × P with respect to
the standard partial orders if and only if
∂fi
≥ 0, ∀ i 6= j,
∂xj
∂fi
≥ 0, ∀ i, j,
∂pj

(x, p) ∈ cl(D) × P
(x, p) ∈ D × P.

U t g(x) = g ◦ φf (t, x),

(2)

where ◦ is the composition of functions and φf (t, x) is a
solution to the considered system. The Koopman semigroup is
linear [11], hence it is natural to study its spectral properties. In
particular, the eigenfunctions sj (x) of the Koopman semigroup
are defined as the functions Rn → C satisfying
U t sj (x) = sj (φf (t, x)) = sj (x) eλj t ,

(3)

and λj ∈ C is the associated eigenvalue. We can also obtain
a very useful expression
f (x)T ∇sj (x) = λj sj (x).

(4)

If f is analytic, and if the eigenvalues λj of the Jacobian of the
vector field at the equilibrium x∗ are distinct and <(λj ) < 0
for all j, then the flow of the system can be expressed (at least
locally) through the following expansion:
∗

φf (t, x) = x +

n
X

sj (x)vj eλj t +

(5)

j=1

X

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

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

where N0 is the space of nonnegative integers, vj are the right
eigenvectors corresponding to λj , the vectors vk1 ,...,kn are the
so-called Koopman modes (cf. [11], [15] for more details). In
the case of a linear system ẋ = Ax with matrix A having the
left eigenvectors wi , the eigenfunctions si (x) are equal to wiT x
and the expansion (5) has only the finite sum. A similar (but
lengthy) expansion can be obtained even if the eigenvalues
λj are not distinct and have linearly dependent eigenvectors
(cf. [21]).
Let λj be such that 0 > <(λ1 ) ≥ <(λj ), j 6= 1, then the
eigenfunction s1 , which we call dominant, can be computed
using the so-called Laplace average (cf. [14])
gλ∗1 (x)

1
= lim
t→∞ t

Zt

(g ◦ φf (s, x))e−λ1 s ds.

(6)

0

For all g that satisfy g(x∗ ) = 0 and v1T ∇g(x∗ ) 6= 0, the
Laplace average gλ∗1 is equal to s1 (x) up to a multiplication
with a scalar. Therefore, without loss of generality, we will
only consider so-called increasing functions as observables.
A function g(x) : Rn → R is called increasing with
respect to K if for all x, z such that x K z we have
g(x) ≤ g(z). The eigenfunctions sj (x) with j ≥ 2 (nondominant eigenfunctions) are generally harder to compute and
are not considered in the present study.
The eigenfunction s1 (x) captures the dominant (i.e. asymptotic) behavior of the system. Moreover for λ1 ∈ R, it follows

from (5) that the trajectories starting from the boundary ∂Bα
of the set Bα = {x||s1 (x)| ≤ α} share the same asymptotic
evolution
φf (t, x) → x∗ + v1 αeλ1 t ,

t → ∞.

Naturally, these sets are important for understanding the
dynamics of the system. The sets ∂Bα = {x||s1 (x)| = α} are
called isostables, and contain the initial conditions of trajectories that converge synchronously toward the equilibrium. In
the neighborhood of x∗ , the isostables are parallel hyperplanes
(if λ ∈ R). We will also employ the notations
o
o
n
n
∂+ Bα = x s1 (x) = α , ∂− Bα = x s1 (x) = −α
for real and nonnegative α. A more rigorous definition of
isostables and more information can be found in [14].
III. G EOMETRIC P ROPERTIES OF M ONOTONE S YSTEMS
A. Spectral Properties of Monotone Systems
We start by establishing the properties of the dominant
eigenfunctions and eigenvalues of monotone systems.
Proposition 2: Assume that the system ẋ = f (x) admits a
stable equilibrium x∗ with a basin of attraction B. Let λi be
the eigenvalues of the Jacobian at x∗ such that 0 > <(λi ) ≥
<(λj ) for all i ≤ j, and let vj be corresponding right eigenvectors. Let s1 (·) ∈ C 1 (B) be the dominant eigenfunction
with v1T ∇s1 (x∗ ) = 1. Then:
(i) if the system is monotone with respect to K on int(B),
then λ1 is real and negative, s1 (x) ≥ s1 (y) for all x, y
satisfying x K y, and v1 K 0. Additionally, we have that
s1 (x) > s1 (y) for all x K y.
(ii) if the system is strongly monotone with respect to K on
int(B), then λ1 is simple, real and negative, s1 (x) > s1 (y)
for all x, y satisfying x K y, and v1 K 0.
Both conditions (i) and (ii) are only necessary and not
sufficient, which is consistent with the linear case. The proof
of this result is similar to the proof of the main result in [22]
concerning spectral properties of the so-called eventually
monotone systems. If the vector field is analytic, then the
second condition is necessary and sufficient for so-called
strong eventual monotonicity.
We can view this result through the prism of positive
operator theory. An operator A : V → V acting on a normed
space V is called positive (or invariant) with respect to a cone
K ⊆ V , if AK ⊆ K. It can be shown that the Koopman
semigroup U t associated with a monotone dynamical system
is a positive operator with respect to the cone of increasing
functions at every time t. Hence the following result may be
seen as a nonlinear version of the Perron-Frobenius theorem
or a version of the Krein-Rutman theorem for the Koopman
operator. We note that a semigroup associated with a nonlinear
system was studied from the positive operator viewpoint in the
context of diffusion processes in [23].
There were other successful attempts to extend the PerronFrobenius theorem to nonlinear systems (cf. [16], [24], [25]),
which however are based on applications of Perron-Frobenius
arguments to nonlinear maps, while our approach is rather
based on their infinite-dimensional extension. Moreover, our

main contribution is that our version uses eigenfunctions,
which give additional insight into the dynamics of the system
and can be estimated using data-based, algebraic or simulation
methods (cf. [15]). Nevertheless, there are similarities, and an
investigation into a connection between these results is one of
the future work directions.
B. Isostables of Monotone Systems
First, we recall a known result on the geometry of basins
of attraction B and their boundary ∂B. This result is briefly
mentioned in [2], [17], [18] in the case of bistable monotone
systems. The following version requires monotonicity on the
basin of attraction only (for a complete proof see e.g. [26]).
Proposition 3: Let the system ẋ = f (x) have an asymptotically stable equilibrium x∗ with a domain of attraction B and
be monotone on B with respect to K. Then for any x, y ∈ B
the order-interval [x, y]K is a subset of B.
We proceed by showing similar properties of the sets Bα =
{x||s1 (x)| ≤ α} and isostables ∂Bα = {x||s1 (x)| = α}.
Proposition 4: Let the system ẋ = f (x) have a stable
hyperbolic equilibrium x∗ with a domain of attraction B and
be monotone with respect to K on B with s1 (·) ∈ C 1 (B).
Then the following statements
hold:
T
(i) for any T
x, y ∈ Bα B, the order-interval [x, y]K is a
subset of Bα B,
(ii) the manifolds ∂+ Bα , ∂− Bα do not contain x, y such
that x K y;
(iii) if the system is strongly monotone then the manifold
∂Bα does not contain x, y such that x K y.
The proof of this result is identical to the case of (strongly)
eventually monotone systems [22]. It is important to note that
the boundary ∂B can contain two points x, y such that x K
y, but in this case these points belong to different manifolds
∂− Bα , ∂+ Bα . The point (ii) can also be shown for α = ∞,
hence the boundaries ∂− B∞ , ∂+ B∞ do not contain x, y such
that x K y as well. The point (iii) does not generalize in a
straightforward manner to the case α = ∞. This is however
shown in [17], [18] under additional assumptions.
C. Basins of Attraction of Monotone Systems Subject to Parameter Variations
First we need to discuss the following assumptions.
A1. The system ẋ = f (x, p) is a bistable system on D with
two stable equilibria x∗ (p), x• (p) for all p. Assume also
that the system is monotone in x and p.
A2. Parameters p take values from a compact set P with
vectors pmax , pmin ∈ P such that
pmin  p  pmax

∀p ∈ P.

A3. The following holds for all p ∈ P:
\
\
x∗ (p) ∈
B(x∗ (q)), x• (p) ∈
B(x• (q))
q∈P

q∈P

A4. x• (pmin )  x∗ (pmax ).
Assumption A1 defines a monotone bistable system, which
is also monotone with respect to parameter variations. In order

ii) Let y ∈ B(x∗ (p1 )), due to monotonicity we have that

State

φ(t, y, p1 )  φ(t, y, p2 ),

State
Fig. 1.

A schematic depiction of Assumptions A3 and A4.

to simplify the presentation we consider only the standard
orders for monotonicity in the initial conditions and parameters
(which can be verified by means of Proposition 1), but our
results hold for other orderings as well. In order to obtain
any meaningful bounds on the basins of attraction, we need to
assume that the set P is compact and ordered, which is done
in Assumption A2. We also assume that parameter variations
are small enough so that no bifurcation occurs under these
parameter variations. Moreover, we assume that the parameter
variations do not affect the locations of the equilibria in
a significant way, which is done in Assumption A3. This
assumption and Assumption A4 are technical and perhaps
can be avoided1 , but they are straightforward to satisfy and
will simplify the presentation of results. The assumptions on
the location of stable equilibria and the basins of attraction
are schematically depicted in Figure 1. Assumptions A1, A2,
and A4 are straightforward to check, while Assumption A3
can only be verified after basins of attraction for pmin and
pmax are computed. Using these assumptions we can derive
the following result.
Theorem 1: Consider the system ẋ = f (x, p) under the
assumptions A1 – A4. Then B(x∗ (p1 )) ⊆ B(x∗ (p2 )) and
B(x• (p1 )) ⊇ B(x• (p2 )) for all p1  p2 . In particular, we
have
B(x∗ (pmin )) ⊆ B(x∗ (p)) ⊆ B(x∗ (pmax )) ∀p ∈ P
B(x• (pmin )) ⊇ B(x• (p)) ⊇ B(x• (pmax )) ∀p ∈ P

(7)

Proof: i) We need to show that for all p1 , p2 we have
that x• (p1 )  x∗ (p2 ). But first we will verify that x∗ (p1 ) 
x∗ (p2 ) for all p1  p2 . Due to monotonicity for all p1  p2
we have that
φ(t, x∗ (p1 ), p1 )  φ(t, x∗ (p1 ), p2 ),
where φ(t, x∗ (p1 ), p1 ) is equal to x∗ (p1 ) for all t,
while φ(t, x∗ (p1 ), p2 ) converges to x∗ (p2 ) since x∗ (p1 ) ∈
B(x∗ (p2 )) due to Assumption A3. Hence, x∗ (p1 )  x∗ (p2 )
for all p1 , p2 ∈ P such that p1  p2 . Similarly we can
show that x• (p1 )  x• (p2 ) for all p1  p2 . Together with
Assumption A4, this implies that x• (p1 )  x∗ (p2 ) for all p1 ,
p2 ∈ P.
1
modify Assumption A3 to require the sets
T For example,
Twe can
B(x∗ (q)),
B(x• (q)) to be non-empty, however, this modified
q∈P

q∈P

assumption is harder to check

for p1  p2 , where φ(t, y, p1 ) converges to x∗ (p1 ) with
t → ∞. Now, by point i), we conclude that φ(t, y, p2 ) cannot
converge to x• (p2 ), hence it converges to x∗ (p2 ). This implies
that B(x∗ (p1 )) ⊆ B(x∗ (p2 )) for all p1  p2 .
iii) The inclusion B(x• (p1 )) ⊇ B(x• (p2 )) follows from
the fact that cl(B(x• (p))) = cl(D\B(x∗ (p))) for all p ∈ P.
While (7) is now straightforward to show.
Theorem 1 implies that we can predict a direction of change
in the basins of attraction subject to parameter variations if
we check a couple of simple conditions. In many applications (e.g. synthetic biology), so-called toggle switches are
used, which are bistable systems. The larger the domain of
attraction, the more robust is the toggle switch with respect
to the corresponding equilibrium. Therefore this result can be
valuable for design purposes in applications, where we need
to ensure operational robustness of a stable equilibrium in a
toggle switch subject to intrinsic and/or exogenous noise.
IV. T OOLS FOR C OMPUTING I SOSTABLES AND BASINS OF
ATTRACTION
In control theory, a go-to approach to compute forwardinvariant sets (not only basins of attraction) of dynamical
systems is sum-of-squares (SOS) programming. This approach
can be applied to systems with polynomial vector fields and
also can provide estimates on the basins of attraction for such
systems. There exists a number different SOS-based methods
for estimation of basins of attraction, see e.g. [27] and [28]
and the references within.
Another option is to compute the function s1 (·), which
provides the isostables and the basin of attraction. In the case
of a polynomial vector field, we can formulate the computation
of s1 (x) as an infinite dimensional linear algebraic problem
using (4). Hence, we can provide an approximation of s1 (x)
using linear algebra by parametrizing s1 (x) with a finite
number of basis functions (cf. [15]). On another hand, we
can estimate s1 (x) directly from data using dynamic mode
decomposition methods (cf. [29], [30]). These two options
provide extremely cheap estimates of s1 . In fact, the algebraic
methods (as demonstrated in [15]) also provide good estimates
on basins of attraction. However, we cannot typically compute estimates with an excellent approximation quality, which
comes as trade-off for fast computations.
If the quality of approximation is very important we can
simulate a number of trajectories with initial points on a
mesh grid of the state-space and compute s1 (·) in these points
using Laplace averages (6). After that different interpolation
or machine learning methods can be applied to estimate the
dominant eigenfunction.
We stress that computing the eigenfunction using Laplace
averages is not necessarily more computationally difficult than
computing the basin of attraction using sum-of-squares. Moreover, the Laplace average method can be applied to systems
with non-polynomial vector fields and is well-suited to high
dimensions. This is due to the fact that the simulation time

of a nonlinear system can significantly be reduced by model
reduction (time-scale separation) methods. In contrast to SOSbased methods, however, the Laplace average method does not
generally provide a certificate, that is we cannot guarantee that
the estimated set is an inner or an outer approximation of a
basin of attraction.
In this paper, we use Laplace averages in our computational algorithm. We use an adaptive grid and thus lower
computational complexity by exploiting the facts that the
eigenfunction is increasing and only one level set is needed.
Finally, monotonicity allows us to guarantee that we compute
an inner and an outer approximation of the basin, hence
providing a certificate.
A. Algorithm for Computation of a Level Set of an Increasing
Function
First we derive an algorithm for computation of the level
sets of an arbitrary increasing with respect to K function
y = g(z). Our algorithm can be extended to an arbitrary
positive cone K in Rn , however, in order to simplify the
presentation we focus on the ordering induced by a positive orthant. We compute the level sets within a hypercube
B = {z|b0  z  b1 }. In particular, we compute a set of
points Mmin , which lies in {z|g(z) < α}, and a set of points
Mmax , which does not lie in {z|g(z) < α}. The set Mmin
(respectively, Mmax ) will not contain points x, y such that
x  y. Since g(z) is an increasing function this will allow to
build a piecewise constant inner approximation (respectively,
an outer approximation) of the sublevel set. In order to do so
we will use the following oracle
(
0 if g(z) < α,
O(z) =
(8)
1 otherwise.
We will generate the points on the face of B with zn = b0n
randomly, while choosing the point zn greedily. In order to
do so efficiently, we need to modify a hypercube B in such
a way that the curve y = g(z) does not intersect the faces
zn = b0n and zn = b1n . Therefore, first we need to increase
b1n sufficiently such that all vertices of the face with zn = b1n
return the value of the oracle O equal to one. At the same
time, we need to lower the face with zn = b0n significantly so
that all the vertices of this face return O equal to zero. Due
the constraints on the function we can have b0 = 0 and z  0,
which is typically the case in the state-space of biological
applications. In this case, we need to move the points b01 , . . .
b0n−1 and shrink B in order to achieve the same goal. We can
do so, for example, by computing the level set g(z) = α with
zn = b0n and a pick a point with the smallest norm.
Our algorithm exploits the monotonicity property of g in
the following manner. If a sample z j is such that O(z j ) = 0,
then for all w  z j we have O(w) = 0. Similarly, if a sample
z j is such that O(z j ) = 1, then for all w  z j we have
O(w) = 1. Therefore, we need to keep track of the largest (in
the order) samples z j with O(z j ) = 0, the set of which we
denote Mmin , and the smallest (in the order) samples z j with
O(z j ) = 1, the set of which we denote Mmax . Hence, in order
to approximate the function g(·) at every step we generate new

samples z j such that Mmin  z j  Mmax . After generating
the samples we compute the oracle values O(z j ) and add new
samples to the set Mmin if O(z j ) = 0, or to the set Mmax if
O(z j ) = 1. After new samples have been added there might
exist w, z in Mmin (respectively, Mmax ) such that w  z
(respectively, w  z). In these cases all such samples w must
be removed from Mmin (respectively, Mmax ). We call this
procedure pruning of the sets.
All that is left is to explain the procedure of generation of
new samples. First of all, let A = {z|Mmin  z  Mmax }
with the Lebesgue measure |A|. The measure of the whole
space, which is the defined above hypercube B, is its volume
|B|. The value |A| represents the current error of the algorithm.
In order to sample from the actual level set y = g(z) the
measure |A| has to be equal to zero.
We consider two ways of generating samples: random and
greedy. In “random” generation we randomly
generate points

in the set A. Let z = z1 . . . zn be a sample, which we
need to generate. First, we generate the scalars z1 , . . . , zn−1
on the edge the hypercube B with zn = b0n . We do it by using a
probability distribution δ2 with the support on the whole edge.
Then we compute the limits for the generation of zn such that
Mmin K z K Mmax . Next, we generate zn using the same
distribution δ2 , while adjusting its support. We generate Nrn
samples in this manner.
The “greedy” generation is meant to exploit more information about |A| and decrease the value |A|, which represents
the
S max
i
min
error of the algorithm. For every
sample
z
in
M
M
S
we find points wj from Mmin Mmax such that wj  z i and
such that there is no other ξ satisfying wj  ξ  z i . Then
we form hypercubes Dji with vertices wj and z i and choose
Dji with the largest volume. We denote these hypercubes
Di .
S
Having done it for every sample z i in Mmin Mmax we
compute Ngr hypercubes Di with largest volumes. After that
we generate Ngr samples in each of the hypercubes Di using a
probability distribution δ1 . If no information about the system
is known then we choose δ1 , δ2 as uniform distributions.
Otherwise, we can choose a Beta distribution so that the
median approximates the most probable location of the curve.
The procedure is summarized in Algorithm 1. Note that
our sample generation guarantees that the measure |A| is
nonincreasing with generation of new samples, while random
sampling guarantees that the value |A| eventually converges
to zero with the rate of convergence depending on specific
functions g. However, in our examples we observed exponential empirical convergence, which can be explained by the fact
that our algorithm can be seen as a bisection procedure using
random sampling.
1) Oracles for Computation of Isostables and Basins of
Attraction: Without loss of generality we assume that the
stable equilibrium lies at the origin. Since the sets ∂Bα (0)
and ∂B(0) are level sets of some increasing function y = g(x)
according to Proposition 4, we can use Algorithm 1 to compute
them. In order to do so, we can use the following oracle:
(
0 if s1 (x) < α and kφ(T, x) − x∗ k < ε,
O(x) =
(9)
1 otherwise,

Algorithm 1 Computation of a level set of incomparable
points of a function α = g(x)
1: Inputs: Oracle O, the number of greedy samples Ngr ,
the number of random samples Nrn , the total number of
samples Ntot , the boundary points b0 , b1
2: Outputs: The sets of points Mmin , Mmax .
3: Set N = Ngr + Nrn
4: Adjust the boundary points b0 and b1
5: Set Mmax = {ci }, Mmin = {di }
6: for i = 1, . . . , [Ntot /N ] do
7:
Compute the hypercubes Di , where i = 1, . . . , Ngr
8:
Generate Ngr samples z j in the hypercubes Di using
a probability distribution δ1 .
9:
Generate Nrn samples z j in the admissible set A using
a probability distribution δ2 .
10:
for j = 1, . . . , N do
11:
If O(z j ) = 0 add the sample z j to the set Mmin ,
otherwise add z j to the set Mmax .
12:
end for
13:
Prune the sets Mmin , Mmax for comparable samples
z j and update the admissible set A.
14: end for
Isostables and Boundaries of Basins of Attraction
100

We compute the isostables ∂Bα (x∗ (·)) with α = 0, 2 · 103 , ∞
for systems with parameters p, pmin , pmax , where ∂B∞ (x∗ (·))
is the boundary of the basin of attraction B(x∗ (·)), and the
matrix of parameters p is equal to:


2 1000 4 1
p=
.
1 1000 3 2

80

State

60
40
20
0

system for synthetic biology [31]. We consider the following
model of the toggle switch
p12
− p14 x1 ,
ẋ1 = p11 +
1 + xp213
(10)
p22
ẋ2 = p21 +
p23 − p24 x2 ,
1 + x1
where all pij ≥ 0. The states xi represent the concentration of
proteins LacI and TetR, whose mutual repression is modeled
via a rational function. The parameters p11 and p21 model
the basal synthesis rate of each protein. The parameters p14
and p24 are degradation rate constants, and p12 , p22 describe
the strength of mutual repression. The parameters p13 , p23
are called Hill coefficients. By means of Proposition 1 we
can readily check that the model is monotone on R2≥0 for
all nonnegative
 parameter values with respect to the orthant
diag 1 −1 R2≥0 . Moreover, the model is monotone with
respect to all parameters but p13 , p23 . The stable equilibrium
x∗ has the state x2 “switched on” (i.e., x2 is much larger than
x1 ), while x• has the state x1 “switched on” (i.e., x1 is much
larger than x2 ).
In order to evaluate Theorem 1 we consider the set of
admissible parameters P = {q|pmax p q p pmin }, where




1.8 950 4 1
2.2 1100 4 1
pmin =
pmax =
.
1.2 1050 3 2
0.7 900 3 2

0

5

15

10

20

25

State
Fig. 2. Isostables ∂B0 (x∗ (q)), ∂B2·103 (x∗ (q)), and boundaries of basins
of attraction ∂B(x∗ (q)) for q equal to p, pmin , or pmax . The dots represent
the equilibria x∗ (p) (blue), x∗ (pmin ) (red), x∗ (pmax ) (green) for different
parameter values

where the value of the eigenfunction s1 (x) is computed using
Laplace averages (6), and T is a large enough time. We need
the second condition kφ(T, x) − x∗ k < ε to make sure that
the points also lie in B(0). Even though it is theoretically
unlikely to have a point x with a finite (nonzero) s1 (x) not
lying in B(0), such situations can occur numerically. If we
need to compute B(0), then we can drop the first condition
s1 (x) < α, since in this case α = ∞.
We can also apply this algorithm to the computation of the
so-called switching separatrix defined in the context of pulsecontrolled bistable systems [26], since it is also the level set
of an increasing function.
V. N UMERICAL E XAMPLES
A. A Two-State LacI-TetR Switch
The genetic toggle switch is composed of two mutually
repressive genes LacI and TetR and was a pioneering genetic

The computational results are depicted in Figure 2. These
results suggest that for all parameter values p ∈ P =
{q|pmax Kp q Kp pmin } the manifold ∂B(x∗ (p)) will
lie between the manifolds ∂B(x∗ (pmin )) and ∂B(x∗ (pmax )).
It appears that ∂B0 (x∗ (p)) also lies between the manifolds
∂B0 (x∗ (pmin )) and ∂B0 (x∗ (pmax )), however, in a different
order. This change of order and continuity of s1 implies that
there exists an α such that at least two manifolds ∂Bα (x∗ (p)),
∂Bα (x∗ (pmax )), ∂Bα (x∗ (pmin )) intersect. This case is also
depicted with α = 2 · 103 . This observation implies that
s1 (x, p) is not an increasing function in p. This is consistent
with the linear case, where changes
in the drift matrix A, will

simply rotate the hyperplane x w1T x = 0 , where w1 is the
left dominant eigenvector of A.
B. A Four State Toggle Switch.
We can compute the basins of attraction and isostables for n
dimensional models, as well. Naturally, we can visualize only
3-D cross-sections of basins of attraction of n dimensional systems. We consider again a toggle switch, while modeling also
mRNA concentrations activating their corresponding proteins.
The resulting model is as follows
p21
− p23 x2 ,
ẋ1 = p11 xp212 − p13 x1 , ẋ2 =
1 + xp322
(11)
p41
ẋ3 = p31 xp432 − p33 x3 , ẋ4 =
p42 − p43 x4 ,
1 + x1

Isostables and the Boundary of the Basin of Attraction
12000

State

8000
4000
0
2000

2000
1000

State

1000
0

0

State

Fig. 3. Isostables and the basin of attraction of a four state toggle switch.
We depict cross section of isostables ∂B104 , ∂B2·104 and ∂B105 , and the
boundary of the basin of attraction ∂B for a four state system with x2 =
10−7 .

where the parameters have the following values
p11 = 1000, p12 = 1, p13 = 1, p21 = 2, p22 = 3, p23 = 1
p31 = 1000, p32 = 1, p33 = 2, p41 = 1, p42 = 3, p43 = 2,
Again, by means of Proposition 1, we can check that the
model is monotone on R4≥0 for all nonnegative parameter

values with respect to the orthant diag 1 1 −1 −1 R4≥0 .
We will plot a cross-section with x2 = 10−7 . In Figure 3, we
plot isostables ∂− B1 , ∂+ B104 , ∂+ B5·104 and ∂+ B105 . We also
plot the boundary of the domain of attraction.
It is noticeable that the sets ∂Bα do not depend on the state
x4 in a significant manner. Hence the state x4 can be reduced
using time-scale separation methods. Similarly, we can reduce
x2 . This is in an agreement with “a rule of thumb” in proteinmRNA interactions, which states that mRNA dynamics are
much faster than protein dynamics and hence can be reduced.
Moreover, for every fixed x2 , x4 the shapes of isostables are
similar to the shape of the isostables of a two state toggle
switch.
VI. C ONCLUSION
In this paper, we study geometric properties of monotone
systems such as properties of isostables and basins of attraction. Isostables contain the initial states of trajectories that
converge synchronously toward the equilibrium and can be
computed using the Koopman operator framework. Isostables
are boundaries of forward-invariant sets of dynamical systems,
and hence serve as convenient refinement of the basin of
attraction notion. We show that the isostables and the boundary
of basin of attraction of monotone systems share the same
geometric properties. For example, these manifolds do not
contain strictly comparable points in the order. Our derivations
are based on a positivity result for the Koopman semigroup,
which is related to the statement of the Perron-Frobenius
theorem. We then focus on properties of basins of attraction
of bistable systems. We show how to estimate basins of
attraction of monotone parameter-dependent systems subject to
parametric uncertainty, and illustrate our results on numerical
examples. Future work includes more practical implications
of our results in synthetic biology, and a further study of
the properties of our algorithm such as rate of convergence,
efficiency and scalability.

R EFERENCES
[1] W. W. Leontief, Input-output economics. Oxford University Press on
Demand, 1986.
[2] M. Hirsch, H. Smith, et al., “Monotone dynamical systems,” Handbook
of differential equations: ordinary differential equations, vol. 2, pp. 239–
357, 2005.
[3] P. G. Coxson and H. Shapiro, “Positive input reachability and controllability of positive systems,” Linear Algebra and its Applications, vol. 94,
pp. 35–53, 1987.
[4] L. Farina and S. Rinaldi, Positive linear systems: theory and applications. John Wiley & Sons, 2010.
[5] A. Sootla and A. Rantzer, “Scalable positivity preserving model reduction using linear energy functions,” in Proc IEEE Conf Decision Control,
Dec. 2012, pp. 4285–4290.
[6] C. Grussler and T. Damm, “A symmetry approach for balanced truncation of positive linear systems,” in 51st IEEE Conference on Decision
and Control, Maui, Hi, USA, Dec. 2012, pp. 4308–4313.
[7] A. Rantzer, “Scalable control of positive systems,” European Journal of
Control, vol. 24, pp. 72–80, 2015.
[8] T. Tanaka and C. Langbort, “The bounded real lemma for internally
positive systems and h-infinity structured static state feedback,” IEEE
transactions on automatic control, vol. 56, no. 9, pp. 2218–2223, 2011.
[9] A. Berman and R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences. SIAM, 1994, vol. 9.
[10] M. G. Krein and M. A. Rutman, “Linear operators leaving invariant a
cone in a Banach space,” Amer. Math. Soc. Translation, p. 128, 1950.
[11] I. Mezić, “Spectral properties of dynamical systems, model reduction
and decompositions,” Nonlinear Dynam, vol. 41, no. 1-3, pp. 309–325,
2005.
[12] M. Budišić, R. Mohr, and I. Mezić, “Applied koopmanism,” Chaos,
vol. 22, no. 4, p. 047510, 2012.
[13] I. Mezic, “Analysis of fluid flows via spectral properties of the Koopman
operator,” Annual Review of Fluid Mechanics, vol. 45, pp. 357–378,
2013.
[14] A. Mauroy, I. Mezić, and J. Moehlis, “Isostables, isochrons, and Koopman spectrum for the action–angle representation of stable fixed point
dynamics,” Physica D, vol. 261, pp. 19–30, 2013.
[15] A. Mauroy and I. Mezic, “Global stability analysis using the eigenfunctions of the Koopman operator,” IEEE Tran Autom Control, 2016, In
press.
[16] B. Lemmens and R. Nussbaum, Nonlinear Perron-Frobenius Theory.
Cambridge University Press, 2012, vol. 189.
[17] H. Smith and H. Thieme, “Stable coexistence and bi-stability for
competitive systems on ordered banach spaces,” Journal of Differential
Equations, vol. 176, no. 1, pp. 195–222, 2001.
[18] J. Jiang, X. Liang, and X.-Q. Zhao, “Saddle-point behavior for monotone
semiflows and reaction–diffusion models,” J. Differential Equations, vol.
203, no. 2, pp. 313–330, 2004.
[19] B. S. Rüffer, H. Ito, and P. Dower, “Computing asymptotic gains of
large-scale interconnections,” in IEEE Proc Conf Decision Control,
2010, pp. 7413–7418.
[20] D. Angeli and E. Sontag, “Monotone control systems,” IEEE Trans
Autom Control, vol. 48, no. 10, pp. 1684–1698, 2003.
[21] I. Mezic, “On applications of the spectral theory of the koopman operator
in dynamical systems and control theory,” in IEEE Conf Decision
Control, 2015, pp. 7034–7041.
[22] A. Sootla and A. Mauroy, “Operator-theoretic characterization of eventually monotone systems,” Provisionally accepted for publication in
Automatica, March 2016, http://arxiv.org/abs/1510.01149v2.
[23] I. Herbst and L. Pitt, “Diffusion equation techniques in stochastic
monotonicity and positive correlations,” Probab Theory Rel, vol. 87,
no. 3, pp. 275–312, 1991.
[24] S. Gaubert and J. Gunawardena, “The Perron-Frobenius theorem for
homogeneous, monotone functions,” Transactions of the American
Mathematical Society, vol. 356, no. 12, pp. 4931–4950, 2004.
[25] D. Aeyels and P. De Leenheer, “Extension of the Perron–Frobenius
theorem to homogeneous systems,” SIAM Journal on Control and
Optimization, vol. 41, no. 2, pp. 563–582, 2002.
[26] A. Sootla, D. Oyarzún, D. Angeli, and G.-B. Stan, “Shaping pulses to
control bistable systems: Analysis, computation and counterexamples,”
Automatica, vol. 63, pp. 254–264, Jan. 2016.
[27] D. Henrion and M. Korda, “Convex computation of the region of
attraction of polynomial control systems,” IEEE Tran Autom Control,,
vol. 59, no. 2, pp. 297–312, 2014.
[28] G. Valmorbida and J. Anderson, “Region of attraction analysis via
invariant sets,” in Proc Am Control Conf, 2014, pp. 3591–3596.

[29] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” Journal of Fluid Mechanics, vol. 656, pp. 5–28, 2010.
[30] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N.
Kutz, “On dynamic mode decomposition: Theory and applications,” J
Comput Dynamics, vol. 1, no. 2, pp. 391 – 421, December 2014.
[31] T. Gardner, C. R. Cantor, and J. J. Collins, “Construction of a genetic
toggle switch in escherichia coli,” Nature, vol. 403, pp. 339–342, 2000.

