Approximating the Frequency Response of Contractive Systems

arXiv:1702.06576v1 [] 21 Feb 2017

Michael Margaliot∗

Samuel Coogan†

Abstract
We consider contractive systems whose trajectories evolve on a compact and convex state-space. It is
well-known that if the time-varying vector field of the system is periodic then the system admits a unique
globally asymptotically stable periodic solution. Obtaining explicit information on this periodic solution
and its dependence on various parameters is important both theoretically and in numerous applications.
We develop an approach for approximating such a periodic trajectory using the periodic trajectory of a
simpler system (e.g. an LTI system). Our approximation includes an error bound that is based on the
input-to-state stability property of contractive systems. We show that in some cases this error bound
can be computed explicitly. We also use the bound to derive a new theoretical result, namely, that a
contractive system with an additive periodic input behaves like a low pass filter. We demonstrate our
results using several examples from systems biology.

1

Introduction

A dynamical system is called contractive if any two trajectories approach each other [1, 2]. This is a strong
property with many important implications. For example, if the trajectories evolve on a compact and convex
state-space Ω then the system admits an equilibrium point e ∈ Ω, and since every trajectory converges to
the trajectory emanating from e, e is globally asymptotically stable. Note that establishing this does not
require an explicit description of e.
More generally, contractive systems with a periodic excitation entrain, that is, their trajectories converge
to a periodic solution with the same period as the excitation. This property is very important in applications
ranging from entrainment of biological systems to periodic excitations (e.g., the 24h solar day or the periodic
cell-cycle division program) to the entrainment of synchronous generators to the frequency of the electric
grid. However, the proof of the entrainment property of contractive systems is based on implicit arguments
(see, e.g. [3]) and provides no explicit information on the periodic trajectory (except for its period).
Contraction theory has found numerous applications in systems and control theory, systems biology [4],
and more (see e.g. the recent survey [2]). A particularly interesting line of research is based on combining
contraction theory and graph theory in order to study various networks of multi-agent systems (see, e.g. [5,
6, 7, 8]).
As already noted by Desoer and Haneda [9], contractive systems satisfy a special case of the input-tostate stability (ISS) property (see the survey paper [10]). Desoer and Haneda used this to derive bounds on
the error between trajectories of a continuous-time contractive system and its time-discretized model. This
is important when computing solutions of contractive systems using numerical integration methods [11].
Sontag [12] has shown that contractive systems satisfy a “converging-input converging output” property. A
recent paper [13] used the ISS property to derive a bound on the error between trajectories of a continuoustime contractive system and those of some “simpler” continuous-time system (e.g. an LTI system). This
bound is particularly useful when the simpler model can be solved explicitly.
Here, we derive new bounds on the distance between the periodic trajectory of a contractive system
and the periodic trajectory of a “simpler” system, e.g. an LTI system with a periodic forcing. We show
∗ M. Margaliot is with the School of Electrical Engineering and the Sagol School of Neuroscience, Tel-Aviv. University,
Tel-Aviv 69978, Israel. E-mail: michaelm@eng.tau.ac.il. The research of MM is partially supported by research grants
from the Israeli Ministry of Science, Technology & Space, the US-Israel Binational Science Foundation, and the Israel Science
Foundation (ISF grant 410/15)
† S. Coogan is with the Department of Electrical Engineering, University of California, Los Angeles.
E-mail:
scoogan@ucla.edu.

1

several cases where the periodic trajectory of the simpler system is explicitly known and the bound is also
explicit, so this provides considerable information on the unknown periodic trajectory of the contractive
system. The explicit bounds also pave the way for new theoretical results. We demonstrate this by using
one of the bounds to prove that any contractive system with an additive sinusoidal forcing behaves like a
low-pass filter, that is, as the frequency of the sinusoidal signal goes to infinity the corresponding solution
of the system converges to an equilibrium state. This generalizes the well-known behavior of asymptotically
stable LTI systems.
The remainder of this paper is organized as follows. The next section reviews some properties of contractive systems and in particular their ISS property. For more details, including the historic development
of contraction theory, see e.g. [14, 15]. The next three sections describe our main results. Section 3 develops
a bound for the difference between the periodic trajectories of two systems: a contractive system and some
simpler “approximating” system. We show using an example that in general this bound cannot be improved.
Section 4 suggests two possible approximating systems for the case of a contractive system with a periodic
forcing. Section 5 shows how the explicit bounds can be used to derive a new theoretical result on the
frequency response of contractive systems. The final section concludes and describes possible directions for
further research.

2

Preliminaries

Consider the time-varying dynamical system
ẋ(t) = f (t, x(t)),

(1)
n

with the state x evolving on a positively invariant convex set Ω ⊆ R . We assume that f (t, x) is differentiable with respect to x, and that both f (t, x) and its Jacobian J(t, x) := ∂f
∂x (t, x) are continuous in (t, x).
Let x(t, t0 , x0 ) denote the solution of (1) at time t ≥ t0 for the initial condition x(t0 ) = x0 . For the sake of
simplicity, we assume from here on that x(t, t0 , x0 ) exists and is unique for all t ≥ t0 ≥ 0 and all x0 ∈ Ω.
The system (1) is said to be contractive on Ω with respect to a vector norm | · | : Rn → R+ if there
exists η > 0 such that
|x(t, t0 , a) − x(t, t0 , b)| ≤ e−(t−t0 )η |a − b|
(2)

for all t ≥ t0 ≥ 0 and all a, b ∈ Ω. This means that any two trajectories approach one another at an
exponential rate η. This implies in particular that the initial condition is “quickly forgotten”.
Note that contraction can be defined in a more general way, for example with respect to a time- and
space-varying norm [1] (see also [16]). We focus here on exponential contraction with respect to a fixed
vector norm because there exist easy to check sufficient conditions, based on matrix measures, guaranteeing
that (2) holds. A vector norm | · | : Rn → R+ induces a matrix measure µ : Rn×n → R defined by
1
µ(A) := lim (||I + εA|| − 1),
ε↓0 ε

where || · || : Rn×n → R+ is the matrix norm induced by | · |. For example, for the ℓ1 vector norm,
denoted | · |1 , the induced matrix norm is the maximum absolute column sum of the matrix, and the induced
matrix measure is
µ1 (A) = max{c1 (A), . . . , cn (A)},
where
cj (A) := Ajj +

X

1≤i≤n
i6=j

|Aij |,

i.e., the sum of the entries in column j of A, with non-diagonal elements replaced by their absolute values.
Matrix measures satisfy several useful properties (see, e.g. [17, 9]). We list here two properties that will be
used later on:
µ(A + B) ≤ µ(A) + µ(B),
µ(cA) = cµ(A) for all c ≥ 0,

(subadditivity),
(homogeneity).
2

If the Jacobian of f satisfies
µ(J(t, x)) ≤ −η,

for all x ∈ Ω and all t ≥ t0 ≥ 0,

(3)

then (2) holds (see [3] for a self-contained proof). This is in fact a particular case of using a LyapunovFinsler function to prove contraction [16]. We will focus on the case where η > 0, but some of our results
hold when η ≤ 0 as well. In this case, (2) provides a bound on how quickly can trajectories of (1) separate
from one another.
Often it is useful to work with scaled vector norms (see, e.g. [18, 19]). Let | · |∗ : Rn → R+ be some
vector norm, and let µ∗ : Rn×n → R denote its induced matrix measure. If D ∈ Rn×n is an invertible
matrix, and | · |∗,D : Rn → R+ is the vector norm defined by |z|∗,D := |Dz|∗ , then the induced matrix
measure is µ∗,D (A) = µ∗ (DAD−1 ). For example, the matrix measure induced by the Euclidean norm | · |2 is
µ2 (A) = max{λ : λ is an eigenvalue of (A + A′ )/2}, so
µ2,D (A) = µ2 (DAD−1 )

(4)

= max{λ : λ is an eigenvalue of (DAD

−1

+ (DAD

−1 ′

) )/2}.

The next result describes an ISS property of contractive systems with an additive input.
Theorem 1 [9] Consider the system
ẋ(t) = f (t, x(t)) + u(t),

(5)

where y → f (t, y) is C 1 for all t ≥ t0 . Fix some vector norm | · | : Rn → R+ and suppose that (3) holds for
the induced matrix measure µ(·). Then the solution of (5) with x(t0 ) = x0 satisfies
Z t
e−η(t−s) |u(s)| ds
|x(t, t0 , x0 )| ≤ e−η(t−t0 ) |x0 | +
t0

for all t ≥ t0 .
Ref. [13] has applied the ISS property to derive a bound on the error between trajectories of the contractive
system (1) and those of a “simpler” dynamical system ẏ = g(t, y(t)). For such a system, pick y0 ∈ Ω, and
let τ ≥ t0 be such that the solution y(t, t0 , y0 ) belongs to Ω for all t ∈ [t0 , τ ]. Then the difference between
the trajectories of the two systems d(t) := x(t, t0 , x0 ) − y(t, t0 , y0 ) satisfies
Z t
−η(t−t0 )
e−η(t−s) |f (s, y(s, t0 , y0 )) − g(s, y(s, t0 , y0 ))| ds
(6)
|d(t)| ≤ e
|x0 − y0 | +
t0

for all t ∈ [t0 , τ ]. The proof of this result is based on noting that
˙ = f (t, x(t)) − f (t, y(t)) + f (t, y(t)) − g(t, y(t))
d(t)
= M (t)d + u(t),

R1
where M (t) := 0 J(t, sx(t) + (1 − s)y(t)) ds, and u(t) := f (t, y(t)) − g(t, y(t)). Since y(t) ∈ Ω for all t ∈
[0, τ ] and Ω is convex, sx(t) + (1 − s)y(t) ∈ Ω for all t ∈ [0, τ ] and all s ∈ [0, 1]. Using (3) and the
subadditivity of matrix measures, which, by continuity, extends to integrals yields µ(M (t)) ≤ −η for all t ∈
˙ = M (t)d(t) + u(t) is a contractive system with an additive “disturbance” u and
[0, τ ]. Summarizing, d(t)
applying the ISS property of contractive systems yields (6).
Note that the integrand in (6) depends on the difference between the vector fields f and g evaluated along
the trajectory of the y system. This is useful, for example, when the trajectory of the y system is explicitly
known.
The applications studied in [13] were contractive systems with time-invariant vector fields approximated
by time-invariant LTI systems. Here, we consider a different case, namely, when the vector field f (t, x) is
time-varying and T -periodic for some T > 0, that is,
f (t, z) = f (t + T, z)
3

for all t ≥ t0 and all z ∈ Ω. It is well-known that in this case every trajectory of (1) converges to a unique
periodic solution γ(t) of (1) with period T (see [3] for a self-contained proof). This entrainment property
is very important in applications (see, e.g. [20, 3]). However, the proof of entrainment is based on implicit
arguments and provides no information on the properties of the period trajectory (except for its period).
Our goal here is to develop a suitable bound for the difference between γ(t) and the periodic solution κ(t) of
some simpler approximating y system, and to suggest suitable approximating systems. We also show that
these explicit bounds can be used to derive new theoretical results on the response of contractive systems to
a sinusoidal input. The next three sections present our main results.

3

Bounds on the difference between two periodic trajectories

In this section, we consider the T -periodic orbit of a T -periodic contractive system. Theorem 2 below is our
main result in this section, and provides a bound on the distance of this periodic orbit to a T -periodic orbit
of some approximating system.
Theorem 2 Consider the system
ẋ = f (t, x)

(7)

whose trajectories evolve on a compact and convex state-space Ω ⊆ Rn . Suppose that f (t, x) is T -periodic
and that f (t, x) and J(t, x) are continuous in (t, x). Let | · | be some vector norm on Rn and µ(·) its induced
matrix measure, and suppose that µ(J(t, x)) ≤ −η < 0 for all t ≥ 0 and all x ∈ Ω. Let γ(t) be the unique
periodic trajectory of (7) with period T . Consider another time-varying system
ẏ = g(t, y)

(8)

and suppose that g(t, y) is also T -periodic and that κ(t) is a T -periodic trajectory of (8) with κ(t) ∈ Ω for
all t ∈ [0, T ]. Define c : R+ → R+ by
Z α
c(α) :=
e−η(α−s) |f (s, κ(s)) − g(s, κ(s))| ds.
(9)
0

Then the difference between the two periodic trajectories satisfies
|γ(τ ) − κ(τ )| ≤

e−ητ
c(T ) + c(τ )
1 − e−ηT

(10)

for all τ ∈ [0, T ].
Note that the bound here depends on the difference between the vector fields f and g evaluated along
the periodic trajectory κ(s) of the “simpler” y system. This is useful for example when the y system is an
asymptotically stable LTI system with a sinusoidal forcing term, as then κ(t) is known explicitly.
Proof. Define the T -periodic function u(t) := f (t, κ(t)) − g(t, κ(t)). For any t ≥ 0, Theorem 1 gives
|γ(t) − κ(t)| ≤ e

−ηt

|γ(0) − κ(0)| +

Z

t

0

e−η(t−s) |u(s)| ds.

(11)

There exist τ ∈ [0, T ) and a non-negative integer k such that t = kT + τ . Observe that γ(t) − κ(t) =
γ(τ ) − κ(τ ). Write the integral on the right-hand side of (11) as
Z

0

t

e

−η(t−s)

|u(s)| ds = e

−ητ

=e

−ητ

Z

kT +τ

e−η(kT −s) |u(s)| ds

0
k−1
X Z (i+1)T
i=0

e

−η(kT −s)

iT

4

|u(s)| ds +

Z

0

τ

e−η(τ −p) |u(p)| dp.

Moreover,
k−1
X Z (i+1)T
i=0

iT

e−η(kT −s) |u(s)| ds = c(T )

1 − e−ηk
1 − e−ηT

so that
|γ(τ ) − κ(τ )| = |γ(t) − κ(t)| ≤ e−η(kT +τ ) |γ(0) − κ(0)| + e−ητ c(T )

1 − e−kηT
+ c(τ )
1 − e−ηT

for all τ ∈ [0, T ). Taking k → ∞ completes the proof.
Note that the bound is actually based on taking the time t → ∞. This is possible because we are
considering the difference between two periodic trajectories.
The next example is important, as it shows that in general the bound (10) cannot be improved.
Example 1 Consider the scalar system
ẋ = f (t, x) := −x + 1 + sin(2πt/T ),

(12)

with T > 0. Note that Ω := [0, 2] is an invariant set of this dynamics, and that f is T -periodic. The Jacobian
of f is J(x) = −1, so for any vector norm the induced matrix measure satisfies µ(J(x)) = −1. For any initial
condition the solution of (12) converges to the T -periodic trajectory:
γ(t) = 1 +

T 2 sin(2πt/T ) − 2πT cos(2πt/T )
.
4π 2 + T 2

(13)

Consider the approximating system ẏ = −y, which is (vacuously) T -periodic, and admits the T -periodic
solution κ(t) ≡ 0, that belongs to Ω for all t. In this case, (9) yields
Z α
e−(α−s) |1 + sin(2πs/T )| ds
c(α) =
0

= 1 − e−α +

2πT e−α − 2πT cos(2πα/T ) + T 2 sin(2πα/T )
.
4π 2 + T 2

Thus for large values of α,
c(α) ≈

−2πT cos(2πα/T ) + T 2 sin(2πα/T )
+ 1.
4π 2 + T 2

(14)

Now consider the case where T → ∞ and τ = T − ε, with ε > 0 and very small. Then (13) implies that the
term on the left-hand side of (10) is
|γ(τ )| ≈ 1 + sin(2π(T − ε)/T ),
whereas (14) implies that the term on the right-hand side of (10) is
e−ητ
c(T ) + c(τ ) ≈ c(τ )
1 − e−ηT
≈ 1 + sin(2π(T − ε)/T ).
Thus, this example shows that in general the bound (10) cannot be improved.
We now derive a simpler (and less tight) bound. By the definition of c(·),
c(α) ≤

1 − e−ηα
max |f (t, κ(t)) − g(t, κ(t))|.
η
t∈[0,α]

for all α ≥ 0, and combining this with (10) yields the following result.
5



Corollary 1 Under the hypotheses of Theorem 2,
|γ(τ ) − κ(τ )| ≤

1
max |f (t, κ(t)) − g(t, κ(t))|
η t∈[0,T ]

for all τ ≥ 0.
This bound is useful in cases where one can establish a bound on the difference between the vector fields f
and g along the periodic trajectory κ of the approximating system. Note that the bound here demonstrates
a clear tradeoff: if g is “close” to f then the error f − g will be small, yet κ may be an unknown complicated
trajectory (as we assume that f is a nonlinear vector filed). On the other hand, if g is relatively simple (e.g.,
the vector field of an LTI system) then κ may be known explicitly yet that difference |f − g| may be large.
To summarize, Theorem 2 and Corollary 1 provide a bound on the distance of the unique T -periodic
trajectory of a contractive system and some T -periodic trajectory of an approximating system. The next
step is to determine a suitable approximating system. We propose two natural approximating systems for
the case where the periodic vector field arises via a periodic forcing function. The first approximating system
considers the time-averaged periodic forcing function to arrive at an autonomous dynamical system with a
unique equilibrium. The second approximating system results from a linearization of the dynamics, keeping
the periodic excitation as is.

4

Approximating Systems

From hereon, we consider a special case of the contractive system (7) with the form
ẋ(t) = f (t, x(t)) = F (x(t), u(t))
where u(t) is a given m-dimensional, T -periodic excitation.

4.1

Averaging the input

Our first result is based on using a “simpler” y system derived by averaging the excitation u over a period.
The excitation in the y system is thus constant. We assume that the y system admits an equilibrium
point e ∈ Ω, and apply Theorem 2 to derive a bound on the distance between the periodic trajectory γ(t) of
the original x system and the point e.
Theorem 3 Consider the system
ẋ = F (x, u),

(15)

where u is an m-dimensional periodic excitation with period T ≥ 0. Suppose that the trajectories of (15)
evolve on a compact and convex state space Ω ⊂ Rn . Assume that for some vector norm | · | : Rn → R+ and
induced matrix measure µ : Rn×n → R,


∂F
µ
(x, u(t)) ≤ −η < 0
∂x
for all t ≥ 0 and all x ∈ Ω. Let γ(t) be the unique, attracting, T -periodic orbit of (15) in Ω. Then, for any
z ∈ Ω,
Z t
|x(t, 0, z) − z| ≤
e−η(t−s) |F (z, u(s))| ds
(16)
0

for all t ≥ 0. In particular, for all t ≥ 0,

|x(t, 0, z) − z| ≤ (1 − e−ηt )c/η,
where c := maxt∈[0,T ] |F (z, u(t))|. Moreover, for all τ ∈ [0, T ],
Z τ
Z T
e−ητ
−η(T −s)
e−η(τ −s) |F (z, u(s))| ds
e
|F
(z,
u(s))|
ds
+
|γ(τ ) − z| ≤
1 − e−ηT 0
0
≤ c/η.
6

(17)
(18)

Proof. Define the approximating system ẏ = G(y) ≡ 0 so that any z ∈ Ω is an equilibrium. Then (16)
follows from the bound (6), and the bounds (17) and (18) follow from Theorem 2 and Corollary 1.
The following simple example demonstrates a case where the bounds in Theorem 3 are tight.
Example 2 Consider the scalar system ẋ = F (x, u) := −ax + b with a > 0. Then γ(t) ≡ b/a =: e is a
periodic trajectory. This system is contracting with rate η = a. The bound (18) gives |e − z| ≤ |− az + b|/a =
|e − z| for any z ∈ R so that this bound is tight.

The next two examples demonstrate that a natural choice for z in Theorem 3 is the equilibrium point
induced by the average of the periodic excitation.
Example 3 Our focus here is on nonlinear dynamical systems, but it is still useful to begin by considering
the linear system
ẋ = Ax + Bu,
(19)
where A ∈ Rn×n is Hurwitz, B ∈ Rn×m , and u is an m-dimensional T -periodic control. It is well-known
that such a system is contractive [2, 21]. For the sake of completeness we repeat the argument here. We use
the notation Q > 0 to denote that a matrix Q is symmetric and positive-definite. Since A is Hurwitz, there
exist η > 0 and Q > 0 such that
QA + A′ Q ≤ −2ηQ.
(20)

Let P > 0 be a matrix such that P 2 = Q. Then multiplying (20) by P −1 on the left and on the right yields
P AP −1 + P −1 A′ P ≤ −2ηI.

(21)

This means that the Jacobian A of (19) satisfies µ2,P (A) ≤ −η, where µ2,P is the matrix measure induced
by the scaled Euclidean norm |z|2,P := |P z|2 (see (4)). Thus, (19) is contractive with respect to this scaled
norm with contraction rate η, and every solution of (19) converges to the unique T -periodic solution γ(t)
RT
of (19). Let ū := T1 0 u(s)ds and choose z = A−1 B ū =: e, the equilibrium of the time-invariant system
with input equal to ū. To apply the bound (18), note that
F (e, u(s)) = Ae + Bu(s) = B(u(s) − ū).
Thus, |F (e, u(s))|2,P = ((u(s) − ū)′ B ′ P ′ P B(u(s) − ū))
|γ(τ ) − e|2,P ≤

1/2

, and the bound (18) yields

1
1/2
max ((u(t) − ū)′ B ′ P ′ P B(u(t) − ū))
η t∈[0,T ]

(22)

for all τ ∈ [0, T ].
Of course, for linear systems the periodic solution corresponding to sinusoidal excitations is known explicitly in terms of the system’s frequency response. Nevertheless, (22) seems to be new and provides considerable
intuition: the bound on the distance between γ(t) and e decreases when: the contraction rate η increases;
the input channel B becomes “more orthogonal” to the matrix P in (21); or maxt∈[0,T ] |u(t) − ū| decreases,
that is, the periodic excitation becomes more similar to its mean.

The next example demonstrates an application of Theorem 3 for a nonlinear contractive system.
Example 4 The ribosome flow model (RFM) [22] is a nonlinear compartmental model describing the unidirectional flow of particles along a 1D chain of n sites using n non-linear first-order differential equations:
ẋ1 = λ0 (1 − x1 ) − λ1 x1 (1 − x2 ),
ẋ2 = λ1 x1 (1 − x2 ) − λ2 x2 (1 − x3 ),
ẋ3 = λ2 x2 (1 − x3 ) − λ3 x3 (1 − x4 ),
..
.

ẋn−1 = λn−2 xn−2 (1 − xn−1 ) − λn−1 xn−1 (1 − xn ),
ẋn = λn−1 xn−1 (1 − xn ) − λn xn .
7

(23)

Here xi (t) ∈ [0, 1] represents the level of occupancy of site i at time t, normalized such that xi (t) = 1
[xi (t) = 0] means that site i is completely full [empty]. The state-space is thus [0, 1]n , and this is an
invariant set of (23) (see [20]). The transition rate λi > 0 controls the flow from site i to site i + 1,
with λ0 [λn ] called the initiation [exit] rate. To understand these equations, note that they may be written
as ẋi = gi−1 (x) − gi (x), where gk (x) is the flow from site k to site k + 1 at time t. This flow increases
with xk and decreases with xk+1 . In other words, the flow satisfies a “soft” excluded volume principle: as
site k + 1 becomes fuller the flow from site k to site k + 1 decreases. This models the fact that the particles
have volume and thus cannot overtake one another. The rate at which particles leave the chain, that is,
R(t) := λn xn is called the production rate.
The RFM with n > 2 is actually not contractive in the sense defined above on [0, 1]n , as there exists p ∈
[0, 1]n such that J(p) is singular, but it is “weakly” contractive in a well-defined sense; see [23, 24].
Recently, the RFM has been used to model and analyze the flow of ribosomes (the particles) along groups
of codons (the sites) along the mRNA molecule during translation (see, e.g. [25, 26, 27, 20, 28, 29, 30, 31, 32]).
In this case, every ribosome that leaves the chain releases the produced protein, so R(t) is the protein
production rate at time t. The values of the transition rates depend on various biophysical properties, e.g.
the abundance of tRNA molecules that carry the corresponding amino-acids.
Consider the RFM with n = 2 and a time-varying initiation rate u0 (t), that is,
ẋ1 = (1 − x1 )u0 − λ1 x1 (1 − x2 ),
ẋ2 = λ1 x1 (1 − x2 ) − λ2 x2 ,

(24)

where λ1 , λ2 are positive constants. Suppose that u0 (t) = λ0 + sin(2πt/T ), with λ0 > 1, T > 0, i.e.
the initiation rate is a strictly positive periodic function with (minimal) period T . The state space here
is Ω := [0, 1]2 . The Jacobian of (24) is


−u0 (t) − λ1 (1 − x2 )
λ1 x1
.
J(t, x) =
λ1 (1 − x2 )
−λ1 x1 − λ2
The off-diagonal terms are non-negative for any x ∈ [0, 1]2 , so µ1 (J(t, x)) = max{−u0 (t), −λ2 } for
all t ≥ 0 and all x ∈ [0, 1]2 . Thus, the system is contractive with respect to the ℓ1 norm with contraction
rate η := min{λ0 − 1, λ2 } > 0. This means that it admits a unique periodic solution γ ∈ [0, 1]2 , with
period T , and that every solution converges to γ. Entrainment in mRNA translation is important as
biological organisms are often exposed to periodic excitations, for example the periodic cell-cycle division
process. Proper biological functioning requires entrainment to such excitations [20].
RT
Let ū0 = T1 0 u0 (s) ds = λ0 and consider the system
ẏ1 = λ0 (1 − y1 ) − λ1 y1 (1 − y2 ),

ẏ2 = λ1 y1 (1 − y2 ) − λ2 y2 .
This system admits an equilibrium point
h
√
λ2 −λ1 λ2 + d
e = λ0 λ1 −λ02λ
λ
0

1

√ i′
λ0 λ1 +λ0 λ2 +λ1 λ2 − d
2λ1 λ2

(25)

∈ (0, 1)2 ,

(26)

where d := 4λ20 λ1 λ2 + (λ0 λ1 − λ0 λ2 − λ1 λ2 )2 .
Here,


(λ0 + sin(2πs/T ))(1 − e1 ) − λ1 e1 (1 − e2 )
,
F (e, u(s)) =
λ1 e1 (1 − e2 ) − λ2 e2

′
and since e is an equilibrium point of (25), F (e, u(s)) = (1 − e1 ) sin(2πs/T ) 0 . Thus, the bound (16)
yields
Z t
|x(t, 0, e) − e|1 ≤ (1 − e1 )
e−η(t−s) | sin(2πs/T )| ds.
(27)
0

8

x(t, 0, e)
0.102

0.101
x2

e

0.1

0.099
0.87

0.88

0.89

0.9

0.91

0.92

x1

Figure 1: RFM in Example 4: Averaging the periodic excitation leads to an autonomous approximating
system with a unique equilibrium e. For the periodic excitation the trajectory x(t, 0, e) converges to the
unique periodic solution γ(t) of the RFM, and Theorem 1 provides a bound on the distance between x(t, 0, e)
and e for all t ≥ 0.
Likewise, (17) implies
e−ητ
|γ(τ ) − e|1 ≤ (1 − e1 )
1 − e−ηT

Z

T

e

0

−η(T −s)

| sin(2πs/T )| ds + (1 − e1 )

Z

τ

0

e−η(τ −s) | sin(2πs/T )| ds

(28)

for all τ . Furthermore,

|F (e, u(t))|1 = (1 − e1 )| sin(2πt/T )| ≤ 1 − e1
so (17) implies the simpler yet more conservative bound
|γ(τ ) − e|1 ≤ (1 − e1 )/η

(29)

for all τ .
Note that all the bounds above can be computed explicitly. For example, a tedious yet straightforward
calculation of the integrals in (28) yields
|γ(τ ) − e|1 ≤

2πT (1 − e1 ) coth(ηT /4)
T (1 − e1 )
r(τ ),
+ 2
eητ (4π 2 + η 2 T 2 )
4π + η 2 T 2

where
r(τ ) :=

(

2πe−ητ + ηT sin(2πτ /T ) − 2π cos(2πτ /T )
if 0 ≤ τ < T /2,
2πe−ητ (1 + 2eηT /2 ) − ηT sin(2πτ /T ) + 2π cos(2πτ /T ) if T /2 ≤ τ < T.

Summarizing, in this example, we have an analytical expression both for e and for the error bounds. Taken
together, this provides considerable explicit information on the periodic trajectory γ.
For the case λ0 = 4, λ1 = 1/2, λ2 = 4, and T = 2, Fig. 1 depicts the periodic trajectory of (24) and the
equilibrium point e of (25), and Fig. 2 depicts the error |x(t, 0, e) − e|1 and the bound (27). In this case, (26)

′
yields e ≈ 0.8990 0.1010 (all numerical values in this note are to four digit accuracy). Fig. 3 illustrates
the other bounds on the periodic trajectory. It may be seen that these bounds indeed provide a reasonable
approximation for the ℓ1 distance between the unknown periodic trajectory and the point e.

Thm. 3 is based on averaging the excitation over a period, thus obtaining a constant input. Such an
RT
approximation is not always suitable. For example, when u(t) = sin(2πt/T ) then ū := T1 0 u(t)dt = 0 for
all T . This may obscure the effect of the frequency of the excitation in the derived bounds. The approach
in the next subsection tries to overcome this using a different approximating system, namely, an LTI system
that is excited by the original periodic input.
9

0.03

0.02

0.01

0

0

2

4

6

t

Figure 2: The error |x(t, 0, e) − e|1 (solid line) and the bound provided by (27) (dashed line) as a function
of time for Example 4.
0.04
0.03
0.02
0.01
0

0

1
T
2

T

t

Figure 3: RFM in Example 4. The error |γ(t)−e|1 (solid line) and the bounds in Theorem 3: the bound (28)
(dashed line) and the bound (29) (dotted line). These bounds can be obtained analytically for this example.

4.2

An LTI approximation

Theorem 4 Consider the system
ẋ = F (x, u),

(30)

where u is an m-dimensional periodic excitation with period T > 0. Suppose that the trajectories of (30)
evolve on a compact and convex state space Ω ⊂ Rn . Assume that for some vector norm | · | : Rn → R+ and
the induced matrix measure µ : Rn×n → R,


∂F
(x, u(t)) ≤ −η < 0
µ
∂x
for all t ≥ 0, all x ∈ Ω. Let γ(t) be the unique, attracting, T -periodic orbit of (30) in Ω.
Suppose also that for the unforced dynamics, i.e. ẋ = F (x, 0), there exists a locally stable equilibrium
point e ∈ Ω, and without loss of generality, that e = 0. Let
A :=

∂F
(0, 0),
∂x

B :=

∂F
(0, 0),
∂u

and consider the LTI approximating system
ẏ = Ay + Bu := G(y, u).

(31)

Pick x0 , y0 ∈ Ω and let τ ≥ 0 be such that y(t) ∈ Ω for all t ∈ [0, τ ] where y(t) is the solution of (31) with
10

y(0) = y0 . Then
|x(t) − y(t)| ≤ e−ηt |x0 − y0 | +

Z

0

t

e−η(t−s) |F (y(s), u(s)) − G(y(s), u(s))| ds

for all t ∈ [0, τ ]. Moreover, let κ(t) be the unique T -periodic trajectory of (31) and assume that κ(t) ∈ Ω for
all t. Then, for all τ ∈ [0, T ],
|γ(τ ) − κ(τ )| ≤

Z T
e−ητ
e−η(T −s) |F (κ(s), u(s)) − G(κ(s), u(s))| ds
1 − e−ηT 0
Z τ
e−η(τ −s) |F (κ(s), u(s)) − G(κ(s), u(s))| ds
+

(32)

0

≤

1
max |F (κ(t), u(t)) − G(κ(t), u(t))|.
η t∈[0,T ]

We emphasize again that the advantage of the bounds here is that the integrand depends on the difference
between the vector fields F and G evaluated along the solution κ of the LTI system (31). Note that our
assumptions imply that A is Hurwitz and thus, for any initial condition, y(t) converges to the periodic
trajectory κ(t). In some cases, this solution can be written explicitly, and the integral can be computed
explicitly. For example, if u(t) is a complex exponential, then κ(t) is also a complex exponential and can
be easily computed using a Fourier transform. Then a bound on |F (κ(t), u(t)) − G(κ(t), u(t))|, t ∈ [0, T ],
may be straightforward to establish. This leads to the following corollary of Theorem 4. For the sake of
simplicity, we state this for the case of a scalar control.
Pp
Corollary 2 Consider the system (30) where u(t) = i=1 ai cos(ωi t) with ai ∈ R and every ωi has the
form ωi = 2πki /T , with ki a non-negative integer. Suppose that the trajectories of (30) evolve on a compact
and convex state space Ω ⊂ Rn . Assume that for some vector norm | · | : Rn → R+ the induced matrix
measure µ : Rn×n → R satisfies
!
p
X
∂F
ai cos(ωi t) ≤ −η < 0
x,
∂x
i=1
for all t ≥ 0 and all x ∈ Ω. Let γ(t) be the unique, attracting, T -periodic orbit of (30) in Ω.
Suppose also that the unforced dynamics, i.e. ẋ = F (x, 0) admits a locally stable equilibrium point e ∈ Ω,
and without loss of generality, that e = 0. Let
A :=

∂F
(0, 0),
∂x

B :=

∂F
(0, 0),
∂u

and consider the approximating system
ẏ = Ay + Bu := G(y, u).

(33)

Let ĝ(s) := (sI − A)−1 b and let κ(t) be the unique T -periodic trajectory of (33), that is,
κr (t) =

p
X
i=1

ai |ĝr (jωi )| cos(ωi t + ∠ĝr (jωi )),

r = 1, . . . , n,

and assume that κ(t) ∈ Ω for all t ∈ [0, T ]. Then for all τ ∈ [0, T ],
1
|γ(τ ) − κ(τ )| ≤
max H
η t∈[0,T ]
where H(z, v) := F (z, v) − G(z, v).
11

κ(t),

p
X
i=1

!

ai cos(ωi t)

(34)

Remark 1 Our focus here is on cases where κ(t) is explicitly known. However, the derived bounds are useful
even when this is not the case. For example, suppose that κ(t) is not known, yet the bounds
Cr := max |κr (t)|,
t∈[0,T ]

r = 1, . . . , n,

Cv := max |u(t)|,
t∈[0,T ]

are known. Then (34) implies that
|γ(τ ) − κ(τ )|≤

1
max |H(z, v)|
η z,v
subject to |zr | ≤ Cr , r = 1, . . . , n
|v| ≤ Cv .

This nonlinear optimization program is useful because the feasible set is a box constraint and thus is a convex
set.
Example 5 We again consider the RFM with n = 2 and the periodic initiation rate u0 (t) := λ0 + u(t),
with λ0 > 1 and u(t) = sin(2πt/T ). Again, let e be the unique equilibrium of the system when the initiation
˙ = Aδx + bu, where
rate is λ0 (see (26)). Let δx := x − e. Then the linearized system is δx




1 − e1
−λ0 − λ1 (1 − e2 )
λ1 e1
.
, b=
A=
0
λ1 (1 − e2 )
−λ1 e1 − λ2
Note that µ1 (A) = max{−λ0 , −λ2 } < 0, so, in particular, A is Hurwitz. Thus, the approximating system is
ẏ = A(y − e) + bu =: G(y, u),

u(t) = sin(2πt/T ).

(35)

The difference between the vector fields evaluated along a solution of the y system is


λ (y − e1 )(y2 − e2 ) − (y1 − e1 ) sin(2πt/T )
F (y, sin(2πt/T )) − G(y, sin(2πt/T )) = 1 1
.
−λ1 (y1 − e1 )(y2 − e2 )


ĝ1 (s)
Let ĝ(s) :=
= (sI − A)−1 b, and let κ(t) : R → R2 be the unique periodic trajectory of (35) defined
ĝ2 (s)
for all −∞ < t < ∞. Then


|ĝ1 (jω)| sin(ωt + ∠ĝ1 (jω))
κ(t) − e =
,
|ĝ2 (jω)| sin(ωt + ∠ĝ2 (jω))
with ω := 2π/T . By Remark 1,
|γ(t) − κ(t)|1 ≤

1
(|λ1 z1 z2 − z1 v| + |λ1 z1 z2 |)
η
subject to |z1 | ≤ |ĝ1 (jω)|
max

z1 ,z2 ,v

|z2 | ≤ |ĝ2 (jω)|

|v| ≤ 1
1
= (2λ1 |ĝ1 (jω)||ĝ2 (jω)| + |ĝ1 (jω)|)
η

(36)

where η := min{λ0 − 1, λ2 } as before. Note that the bound here depends on the frequency of the periodic
excitation. The more exact bound in (32) can be computed numerically.
For the parameters λ1 = 1/2, λ2 = 4, and T = 2, Figure 4 shows the equilibrium point when λ0 = 4, the
periodic trajectory for the case when the initiation rate is u0 (t) = 4 + sin(2πt/T ), and the periodic trajectory
of the linearized system. Figure 5 illustrates the bounds from Theorem 4. It may be observed that these
bounds provide a reasonable estimate of the error.

12

0.102

x2

0.101

0.1

0.099
0.87

0.88

0.89

0.9

0.91

0.92

x1

Figure 4: RFM in Example 5. The equilibrium e for λ0 = 4 is marked by a dot. The periodic trajectory γ(t)
of the RFM (solid line) and the periodic trajectory κ(t) of the linearized system (dashed line) when u0 (t) =
4 + sin(2πt/T ).
0.006

0.004

0.002

0

0

1
T
2

T

t

Figure 5: RFM in Example 5. The error |γ(t) − κ(t)|1 (solid line) and the bounds (32) (dashed lines) and
(36) (dotted line).
The bound (36) has some interesting implications. For example, if ĝ1 (jω) = 0 for some ω then (36) implies
that γ(t) ≡ κ(t) for a sinusoidal excitation with frequency ω. Similarly, if limω→∞ ĝ1 (jω) = 0 then (36)
implies that for a high frequency sinusoidal forcing term, γ will approach κ. Note that the conclusions on γ
here are based on properties of the LTI system. In the next section, we use this idea to derive a theoretical
result on the response of contractive systems to a sinusoidal input.

5

Contractive Systems as Low-Pass Filters

We consider a contractive systems with an additive input and show that for a high-frequency sinusoidal
input, the periodic trajectory of the contractive system is very similar to that of a suitable LTI system. For
the sake of simplicity, we state this for the case of a scalar control.
Theorem 5 Consider the system
ẋ = f (x) + bu
where
u(t) = a cos(ωt + φ).

13

(37)

Suppose that the trajectories of (37) evolve on a compact and convex state space Ω ⊂ Rn . Assume that for
some vector norm | · | : Rn → R+ , and induced matrix measure µ : Rn×n → R,


∂f
µ
(x) ≤ −η < 0
∂x
for all x ∈ Ω. Denote T := 2π/ω, and let γ(t) be the unique, attracting, T -periodic orbit of (37) in Ω.
Suppose also that for the unforced dynamics, i.e. ẋ = f (x) there exists a locally stable equilibrium point e ∈
Ω, and without loss of generality, that e = 0. Let A := ∂f
∂x (0), and consider the approximating system
ẏ = Ay + bu := G(y, u).

(38)

Let ĝ(s) := (sI − A)−1 b and let κ(t) be the unique T -periodic trajectory of (38), that is,
κr (t) = a|ĝr (jω)| cos(ωt + φ + ∠ĝr (jω)),

r = 1, . . . , n.

(39)

Then
max |γ(t) − κ(t)| = o(1/ω).

t∈[0,T ]

(40)

Proof. It follows from Corollary 2 that for any τ ,
1
max |f (κ(t)) − Aκ(t)|
η t∈[0,T ]
1
≤
max o(|κ(t)|).
η t∈[0,T ]

|γ(τ ) − κ(τ )| ≤

(41)

adj(sI−A)
b, where adj denotes the adjugate, (39) implies that |κr (t)| = O(1/ω) for all r and
Since ĝ(s) = det(sI−A)
all t ∈ [0, T ]. Combining this with (41) completes the proof.
The next two examples demonstrate Theorem 5.

Example 6 We consider a basic model for an externally driven transcriptional module that is ubiquitous
in both biology and synthetic biology (see, e.g., [33, 3]):
ẋ1 =u − δx1 + k1 x2 − k2 (eT − x2 )x1 ,
ẋ2 = − k1 x2 + k2 (eT − x2 )x1 ,

(42)

where δ, k1 , k2 , eT are strictly positive parameters. Here x1 (t) is the concentration at time t of a transcriptional factor X that regulates a downstream transcriptional module by binding to a promoter with
concentration e(t) yielding a protein-promoter complex Y with concentration x2 (t). The binding reaction is
reversible with binding and dissociation rates k2 and k1 , respectively. The linear degradation rate of X is δ,
and as the promoter is not subject to decay, its total concentration, eT , is conserved, so e(t) = eT − x2 (t)
for all t ≥ 0. The input u(t) might represent for example the concentration of an enzyme or of a second
messenger that activates X, so we assume that u(t) ≥ 0 for all t ≥ 0.
Trajectories of (42) evolve on [0, ∞) × [0, eT ]. For an input satisfying 0 ≤ u(t) ≤ c for all t ≥ 0, the set
Ω := [0, (c + k1 eT )/δ] × [0, eT ] is a convex and compact invariant set.
Ref. [3] has shown that (42) is contractive with respect to a certain weighted L1 norm. Indeed, the
Jacobian of (42) is


−δ − k2 (eT − x2 ) k1 + k2 x1
,
J(x) =
k2 (eT − x2 )
−k1 − k2 x1
so for D := diag(d, 1), with d > 0,
DJ(x)D

−1


−δ − k2 (eT − x2 )
=
k2 (eT − x2 )/d
14


(k1 + k2 x1 )d
.
−k1 − k2 x1

(43)

0.2
1

x2

x2

0.1
0

0
−0.1

−1
−0.1

0
x1

−0.2
−0.1

0.1

0
x1

0.1

Figure 6: Trajectories γ (solid line) and κ (dashed line) for the system in Example 6 for ω = 1 (top) and ω = 5
(bottom). Note the different scales in the figures.

The off-diagonal terms here are non-negative, and this means that for any d ∈ ( k2ke2TeT+δ , 1),

µ1,D (J(x)) ≤ −η, for all x1

x2

′

∈ Ω,

where η := min{k1 (1 − d), δ + k2 eT (1 − d−1 )} > 0. Thus, (42) is contractive with respect to the scaled
norm | · |1,D with contraction rate η.
Linearizing (42) yields ẏ = G(y, u) = Ay + bu, with

 

1
−δ − k2 eT
k1
, b :=
A :=
,
0
k2 eT
−k1
and
−1

ĝ(s) = (sI − A)



1
s + k1
.
b= 2
k1
s + (δ + k1 + k2 eT )s + δk1


′
Since f (y) − Ay = k2 y1 y2 1 −1 , the bound (41) yields


′
k2
max |κ1 (t)κ2 (t) 1 −1 |1,D
η t∈[0,T ]
k2 (d + 1)
max |κ1 (t)κ2 (t)|.
≤
η
t∈[0,T ]

|D(γ(τ ) − κ(τ ))|1 ≤

Note that for any input in the form
u(t) =

p
X

(44)

ai cos(ωi t),

i=1

the periodic trajectory κ(t) is explicitly known and thus the bound (44) is explicit. Figure 6 depicts the
trajectories of both the contractive system (42) and of the LTI system for the parameters k1 = 1, k2 = 5,
δ = 1, eT = 2, and the excitation u(t) = cos(ωt) for two different values of ω.1 It may be seen that for a
larger value of ω the difference between γ and κ decreases, as anticipated by (40).

The next example demonstrates the result in Theorem 5 using a nonlinear system for which the frequency
response has been computed explicitly in [34].
1 This control is not positive for all times, yet for the initial conditions in the simulations the trajectory remains in a convex
and compact region in which the off-diagonal terms in (43) are non-negative and contraction holds.

15

Example 7 Consider the system:
ẋ1 = −x1 + x22 ,

ẋ2 = −x2 + u,
where the excitation is u(t) = a sin(ωt), with
 a, ω > 0.
−1 2x2
set of x2 . The Jacobian of (45) is J(x) =
.
0
−1


−1 2x2 /c
DJ(x)D−1 =
, so µ1 (DJ(x)D−1 ) ≤ −1 + 2a
c
0
−1
contractive with respect to the scaled norm |z|1,D := |Dz|1
η =1−

(45)

It is clear that Ω2 := [−a, a] is an invariant
For any c > 0 and D := diag(1, c), we have
for all x2 ∈ Ω2 . If c > 2a, then this systems is
with contraction rate

2a
.
c

(46)

Note that, by taking c arbitrarily large, we may obtain a contraction rate arbitrarily close to 1. The periodic
trajectory γ(t) can be computed explicitly as follows. First, it is clear that
a
sin(ωt − tan−1 (ω)),
γ2 (t) = √
2
1+ω
and substituting this in the first equation of (45) yields


γ1 (t) = M 1 + 5ω 2 + 4ω 4 + (5ω 2 − 1) cos(2ωt) + 2ω(ω 2 − 2) sin(2ωt) ,

(47)

2

where M := 2(1+ω2 )a2 (1+4ω2 ) .
Note that the unforced dynamics
system is ẏ =

′ admits an equilibrium e = 0. The approximating

′
G(y, u) = −y + bu, with b := 0 1 and ĝ(s) = (s + 1)−1 b. Thus, κ(t) = 0 γ2 (t) , so γ(t) − κ(t) =

′
γ1 (t) 0 , and |γ(t) − κ(t)|1,D = |γ1 (t)|. To apply Corollary 2 note that H(z, v) := F (z, v) − G(z, v) =
 2
′
z2 0 , so applying the bound (34) gives
max |γ(t) − κ(t)|1,D = max |γ1 (t)|

t∈[0,T ]

t∈[0,T ]

≤ η −1 max |κ22 (t)|
t∈[0,T ]

= η −1

a2
1 + ω2

(48)

where η is given in (46) with c > 2a, and T = 2π/ω. Taking c → ∞ gives the explicit bound maxt∈[0,T ] |γ(t)−
κ(t)|1,D ≤ a2 /(1 + ω 2 ). In fact, it follows from (47), after some calculation, that
√
a2 (1 + 4ω 2 + 1)
√
.
(49)
max |γ1 (t)| =
t∈[0,T ]
2(1 + ω 2 ) 4ω 2 + 1
Fig. 7 depicts the exact difference maxt∈[0,T ] |γ(t) − κ(t)|1,D given in (49) and the bound a2 /(1 + ω 2 )
implied by (48), as a function of ω for a = 1. Theorem 5 guarantees maxt∈[0,T ] |γ(t) − κ(t)| = o(1/ω), as
seen in the figure.


6

Discussion

Contractive systems entrain to periodic excitations. Analyzing the corresponding periodic solution of the
contractive system and its dependence on various parameters is an important theoretical question with
many potential applications. We developed approximation schemes for this periodic solution using LTI
systems and, using the ISS property of contractive systems, provided bounds on the approximation error.
An important advantage of these bounds is that in some cases they can be computed explicitly. This also
led to a new theoretical result on the behavior of contractive systems for a high frequency excitation.
16

maxt |γ(t) − κ(t)|

100

Bound
Exact

10−1
10−2
10−3
10−2

100
ω

10−1

101

102

Figure 7: Maximum distance between the exact periodic trajectory γ(t) and the approximate periodic
trajectory κ(t) as a function of the excitation frequency ω for Example 7 with a = 1. The solid line is the
exact difference, and the dashed line is the bound on the difference determined by Corollary 2.
More generally, it is well-known that contractive systems whose solutions evolve on a compact statespace have a well-defined frequency response [34, 35]. For the contractive system ẋ = F (x, u), with x ∈ Rn
and u ∈ R, this means that there exists a continuous function α : R3 → Rn such that the following property
holds. For the sinusoidal input u(t) = a sin(ωt), with frequency ω := 2π/T and amplitude a ≥ 0, the solution
of the contractive system converges to a periodic solution γaω satisfying
γaω (t) = α(a sin(ωt), a cos(ωt), ω)
(see [34, Theorem 3]). The function α(v1 , v2 , w) is called the state frequency response. For the special case
of a linear system, i.e. F (x, u) = Ax + bu the state frequency response is known explicitly:
′

α(v1 , v2 , ω) = Π(ω) v1 v2


with Π(ω) := Re(ĝ(jω)) Im(ĝ(jω)) ∈ Rn×2 , and ĝ(s) := (sI − A)−1 b. That is, for linear systems, the
state frequency response recovers the standard notion of frequency response.
Of course, for nonlinear systems it is typically not possible to compute the frequency response analytically.
Our results may be interpreted in this context as follows. Considering Theorem 3, we have that ū =
R
1 T
T 0 u(t)dt = 0 for any a and ω and that ẏ = F (y, 0) admits an equilibrium point e. Thus, e = α(0, 0, ω)
(where e is in fact independent of ω), and (17) may be interpreted as providing bounds on
|α(a sin(ωt), a cos(ωt), ω) − e|.
On the other-hand, the results in Theorem 4 may be interpreted as bounds on the difference
|α(a sin(ωt), a cos(ωt), ω) − ᾱ(a sin(ωt), a cos(ωt), ω)| ,
where ᾱ is the state frequency response of the linearized system ẏ = Ay + Bu, with A = ∂F
∂x (0, 0) and B =
∂F
(0,
0).
∂u
An interesting topic for further research is deriving more theoretical results using the explicit bounds
described here. Other possible topics include the design of an excitation signal that yields a pre-specified
periodic trajectory for a contractive system. This issue is important for example in synthetic biology, where
an important goal is to design programmable biochemical oscillators (see e.g., [36, 37, 38, 39]). Another
possible research topic is the extension of the results presented here to more general classes of dynamical
systems (see, e.g., [40] for a special class of infinite dimensional systems that admit a frequency response).

Acknowledgments
We are grateful to Eduardo D. Sontag for reading an earlier version of this paper and providing us with
many useful comments.
17

References
[1] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear systems,” Automatica, vol. 34,
pp. 683–696, 1998.
[2] Z. Aminzare and E. D. Sontag, “Contraction methods for nonlinear systems: A brief introduction and
some open problems,” in Proc. 53rd IEEE Conf. on Decision and Control, Los Angeles, CA, 2014, pp.
3835–3847.
[3] G. Russo, M. di Bernardo, and E. D. Sontag, “Global entrainment of transcriptional systems to periodic
inputs,” PLOS Computational Biology, vol. 6, p. e1000739, 2010.
[4] G. Russo, M. di Bernardo, and J. J. Slotine, “Contraction theory for systems biology,” in Design
and Analysis of Biomolecular Circuits: Engineering Approaches to Systems and Synthetic Biology,
H. Koeppl, G. Setti, M. di Bernardo, and D. Densmore, Eds. New York, NY: Springer, 2011, pp.
93–114.
[5] G. Russo, M. di Bernardo, and E. D. Sontag, “A contraction approach to the hierarchical analysis and
design of networked systems,” IEEE Trans. Automat. Control, vol. 58, pp. 1328–1331, 2013.
[6] M. Arcak, “Certifying spatially uniform behavior in reaction-diffusion PDE and compartmental ODE
systems,” Automatica, vol. 47, no. 6, pp. 1219–1229, 2011.
[7] G. Russo, M. di Bernardo, and J. J. E. Slotine, “A graphical approach to prove contraction of nonlinear
circuits and systems,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 58, no. 2, pp.
336–348, 2011.
[8] S. Coogan and M. Arcak, “A compartmental model for traffic networks and its dynamical behavior,”
IEEE Trans. Automat. Control, vol. 60, no. 10, pp. 2698–2703, 2015.
[9] C. Desoer and H. Haneda, “The measure of a matrix as a tool to analyze computer algorithms for circuit
analysis,” IEEE Trans. Circuit Theory, vol. 19, pp. 480–486, 1972.
[10] E. D. Sontag, “Input to state stability: Basic concepts and results,” in Nonlinear and Optimal Control
Theory, P. Nistri and G. Stefani, Eds. Berlin, Heidelberg: Springer, 2008, pp. 163–220.
[11] J. Maidens and M. Arcak, “Reachability analysis of nonlinear systems using matrix measures,” IEEE
Trans. Automat. Control, vol. 60, no. 1, pp. 265–270, 2015.
[12] E. D. Sontag, “Contractive systems with inputs,” in Perspectives in Mathematical System Theory,
Control, and Signal Processing, J. Willems, S. Hara, Y. Ohta, and H. Fujioka, Eds. Berlin Heidelberg:
Springer-Verlag, 2010, pp. 217–228.
[13] M. Botner, Y. Zarai, M. Margaliot, and L. Grüne, “On approximating contractive systems,” IEEE Trans.
Automat. Control, 2017, To appear. [Online]. Available: http://ieeexplore.ieee.org/document/7814289/
[14] G. Soderlind, “The logarithmic norm. History and modern theory,” BIT Numerical Mathematics, vol. 46,
pp. 631–652, 2006.
[15] J. Jouffroy, “Some ancestors of contraction analysis,” in Proc. 44th IEEE Conf. on Decision and Control,
Seville, Spain, 2005, pp. 5450–5455.
[16] F. Forni and R. Sepulchre, “A differential Lyapunov framework for contraction analysis,” IEEE Trans.
Automat. Control, vol. 59, no. 3, pp. 614–628, 2014.
[17] M. Vidyasagar, Nonlinear Systems Analysis.

Englewood Cliffs, NJ: Prentice Hall, 1978.

[18] I. W. Sandberg, “On the mathematical foundations of compartmental analysis in biology, medicine, and
ecology,” IEEE Trans. Circuits and Systems, vol. 25, no. 5, pp. 273–279, 1978.

18

[19] S. Coogan, “Separability of Lyapunov functions for contractive monotone systems,” in Proc. 55th IEEE
Conf. on Decision and Control, Las Vegas, NV, 2016, pp. 2184–2189.
[20] M. Margaliot, E. D. Sontag, and T. Tuller, “Entrainment to periodic initiation and transition rates in
a computational model for gene translation,” PLOS ONE, vol. 9, no. 5, p. e96039, 2014.
[21] T. Strom, “On logarithmic norms,” SIAM J . Numerical Analysis, vol. 12, pp. 741–753, 1975.
[22] S. Reuveni, I. Meilijson, M. Kupiec, E. Ruppin, and T. Tuller, “Genome-scale analysis of translation
elongation with a ribosome flow model,” PLOS Computational Biology, vol. 7, p. e1002127, 2011.
[23] M. Margaliot, E. D. Sontag, and T. Tuller, “Checkable conditions for contraction after small transients
in time and amplitude,” in Feedback Stabilization of Controlled Dynamical Systems-In Honor of Laurent
Praly, ser. Lecture Notes in Control and Information Sciences, N. Petit, Ed. Springer-Verlag, 2017,
vol. 466.
[24] ——, “Contraction after small transients,” Automatica, vol. 67, pp. 178–184, 2016.
[25] Y. Zarai, M. Margaliot, and T. Tuller, “Explicit expression for the steady-state translation rate in the
infinite-dimensional homogeneous ribosome flow model,” IEEE/ACM Trans. Computational Biology
and Bioinformatics, vol. 10, pp. 1322–1328, 2013.
[26] Margaliot, M. and Tuller, T., “Ribosome flow model with positive feedback,” J. Royal Society Interface,
vol. 10, p. 20130267, 2013.
[27] M. Margaliot and T. Tuller, “Stability analysis of the ribosome flow model,” IEEE/ACM Trans. Computational Biology and Bioinformatics, vol. 9, pp. 1545–1552, 2012.
[28] A. Raveh, Y. Zarai, M. Margaliot, and T. Tuller, “Ribosome flow model on a ring,” IEEE/ACM Trans.
Computational Biology and Bioinformatics, vol. 12, no. 6, pp. 1429–1439, 2015.
[29] A. Raveh, M. Margaliot, E. D. Sontag, and T. Tuller, “A model for competition for ribosomes in the
cell,” J. Royal Society Interface, vol. 13, no. 116, 2016.
[30] G. Poker, M. Margaliot, and T. Tuller, “Sensitivity of mRNA translation,” Sci. Rep., vol. 5, p. 12795,
2015.
[31] Y. Zarai, M. Margaliot, and T. Tuller, “On the ribosomal density that maximizes protein translation
rate,” PLOS ONE, vol. 11, no. 11, pp. 1–26, 2016.
[32] ——, “Optimal down regulation of mRNA translation,” Sci. Rep., vol. 7, no. 41243, 2017.
[33] D. Del Vecchio, A. J. Ninfa, and E. D. Sontag, “Modular cell biology: Retroactivity and insulation,”
Molecular Systems Biology, vol. 4, no. 1, p. 161, 2008.
[34] A. Pavlov, N. van de Wouw, and H. Nijmeijer, “Frequency response functions for nonlinear convergent
systems,” IEEE Trans. Automat. Control, vol. 52, no. 6, pp. 1159–1165, 2007.
[35] B. S. Ruffer, N. van de Wouw, and M. Mueller, “Convergent systems vs. incremental stability,” Systems
Control Lett., vol. 62, no. 3, pp. 277–285, 2013.
[36] M. B. Elowitz and S. Leibler, “A synthetic oscillatory network of transcriptional regulators,” Nature,
vol. 403, pp. 335–338, 2000.
[37] E. Fung, W. W. Wong, J. K. Suen, T. Bulter, S.-g. Lee, and J. C. Liao, “A synthetic gene-metabolic
oscillator,” Nature, vol. 435, pp. 118–122, 2005.
[38] J. Stricker, S. Cookson, M. R. Bennett, W. H. Mather, L. S. Tsimring, and J. Hasty, “A fast, robust
and tunable synthetic gene oscillator,” Nature, vol. 456, pp. 516–539, 2008.

19

[39] M. Weitz, J. Kim, K. Kapsner, E. Winfree, E. Franco, and F. C. Simmel, “Diversity in the dynamical
behaviour of a compartmentalized programmable biochemical oscillator,” Nature Chemistry, vol. 6, pp.
295–302, 2014.
[40] V. Natarajan and G. Weiss, “Behavior of a stable nonlinear infinite-dimensional system under the
influence of a nonlinear exosystem,” in Proc. 1st IFAC Workshop on Control of Systems Governed by
Partial Differential Equations, Paris, France, 2013, pp. 155–160.

20

