Background
Model of an oscillating biochemical network
In [ 8 ] , Laub and Loomis propose a model of the
molecular network underlying adenosine 3',5'-cyclic
monophosphate (cAMP) oscillations observed in fields of
chemotactic
Dictyostelium discoideum cells. The
model, based on the network depicted in Fig. 1, induces
the spontaneous oscillations in cAMP observed during the
early development of
D. discoideum.
In this model, changes in the enzymatic activities of
these proteins are described by a system of seven
non-linear differential equations:
where the state variable
x = [
x
1 ,...,
x
7 ] represents the concentrations of
the seven proteins:
x
1 = [ACA],
x
2 = [PKA],
x
3 = [ERK2],
x
4 = [REG A],
x
5 = [Internal cAMP],
x
6 = [External cAMP] and
x
7 = [CAR1] and the fourteen different
k
i
represent system parameter values. It was shown
numerically in [ 8 ] that spontaneous oscillations appear
at the nominal parameter values found in Table 1. Note
that because there are typographical errors in the
original paper, the values being used here for the
nominal parameters were obtained directly from the
authors of [ 8 ] .
This particular model is primarily concerned with
describing self-sustaining oscillations in the
biochemical system. From Fig. 2, it is clear that at the
nominal parameter values of the model, this is achieved.
We seek to determine whether the system is robust - that
is, if we change these kinetic parameters, will the
systems oscillatory behaviour persist? We next present
two possible means, based on whether parameters are
changed one at a time or in groups.
Methods
Single parameter robustness: Bifurcation
analysis
Self-sustained oscillations such as those being
modelled here appear as stable limit cycles in trajectory
of the underlying dynamical system [ 11 ] . The existence
and stability of these limit cycles may change under
parametric perturbations. Whenever a stable periodic
solution loses stability as we vary the underlying
parameters of the system and this solution transitions to
another qualitative solution - for example, a
steady-state equilibrium - we say that the system
undergoes a
Hopf bifurcation. It is therefore
possible to use bifurcation theory as a means of
quantifying the robustness of this oscillatory network
model [ 12 13 ] .
Using the bifurcation analysis package AUTO [ 14 ] ,
it is possible to produce one-parameter bifurcation
diagrams for each of the model parameters
k
i
. These diagrams illustrate the steady-state
behaviour of the systems as the parameter values are
changed. Suppose that Hopf bifurcations occur at
k
i
and
i
so that (stable) limit cycles occur for the range (
k
i
,
i
). Both the size of this interval as well as the
proximity of the nominal parameter value to either
boundary are measures of the robustness of the system. To
compare the robustness of the system to the different
parameters, we can define a series of parametric
robustness measures. We define the
degree of robustness (DOR) for each
parameter
k
i
as:
It is straightforward to see that this value is always
between zero and one. The former indicates extreme
parameter sensitivity whereas the latter implies large
insensitivity.
Bifurcation diagrams provide an excellent means of
determining the robustness of systems to single parameter
perturbations. We next describe a method for analysing
and quantifying robustness to simultaneous changes in
several parameters.
Multiparametric robustness: Structural singular
values
Determining the limit cycle: Harmonic balance
method
The first step is to obtain a mathematical
expression for the limit cycle oscillation. The
harmonic balance method can be used [ 17 ] . The basic
idea is to represent the limit cycle by a Fourier
series with unknown coefficients (
a
n ,
i , φ
n ,
i ) and period
T :
The non-linear differential equation can be used to
set up a series of algebraic equations that the
coefficients must satisfy. These equations can be
solved using numerical packages such as Mathematica or
Maple. Depending on the particular form of the limit
cycle, a small finite number of coefficients can be
used. We can denote this periodic solution as
x *(
t ).
Linearization
The non-linear differential equation must now be
linearized about this periodic orbit [ 17 ] . Writing
the state vector
x (
t ) as
x (
t ) =
x *(
t ) +
x
δ (
t )
then the local behaviour of the non-linear system is
governed by that of the linearized system:
δ (
t ) ≅
J (
x *(
t ))
x
≈ (
t )
where
J is the Jacobian matrix of the
vector field
f . Note that since the
linearization is performed about a periodic orbit, the
linear system is periodic.
Restructuring into nominal/uncertainty
systems
The Jacobian matrix includes all uncertain
parameters. At this point we need to separate the
system into a nominal model and a feedback
interconnection that involves all parametric
uncertainty. We first write each parameter as
where
k
i
is the nominal value and δ
i
is the relative amount of perturbation in the
i thparameter. We now separate
the Jacobian matrix as
J (
x *(
t )) =
A
0 (
t ) +
B
0 (
t ) Δ
C
0 (
t ) (4)
where
A
0 (
t ) is the Jacobian matrix with
all parameters at their nominal value, and Δ is a
diagonal matrix containing all the uncertainties δ
i
. Let
y (
t ) =
C
0 (
t )
x
δ (
t ) and
u (
t ) = Δ
y (
t ), the system is now of the
form of Eqn. (2) (with
x (
t ) replaced by
x
δ (
t )).
Discretization
The system can be discretized by sampling the state
and output with sampling period
h =
T/n , where
n is a positive integer and
assuming that the inputs are piecewise constant; this
is also a standard technique in control engineering [
18 ] . In particular, a linear continuous-time system
governed by Eqn. (2) gives rise to the discrete-time,
linear system
ξ (
k + 1) =
A
d
(
k ) ξ(
k ) +
B
d
(
k )
v (
k )
η (
k ) =
C
d
(
k ) ξ(
k )
where
A
d
(
k ) = Φ (
kh +
h ,
kh ),
B
d
(
k ) = Φ (
kh +
h ,
s )
B
0 (
s )
ds ,
C
d
(
k ) =
C
0 (
kh ), and Φ (
t , τ) is the transition matrix
of
A
0 (
t ) [ 19 ] . The discretized
signals are
v (
k ) =
u (
kh ), η(
k ) =
y (
kh ), and ξ(
k ) =
x (
kh ). Periodicity of
A
d
and
B
d
is preserved due to the periodicity of the
transition matrix. Moreover, it is not difficult to
confirm that
A
d
,
B
d
and
C
d
are periodic with period
n . The uncertainty matrix after
discretization is now Δ
d
. The discretization step is illustrated in Fig.
4A.
Lifting
The final step in preparing the system for SSV
analysis is to transform the periodic, linearized
system into an equivalent time-invariant one. The
technique for this is known as
lifting [ 18 ] . Rather than
giving the general formulae, it is easier to illustrate
the general principle with an example.
Suppose that a discrete-time system with state
variable ξ, input
v , and output η obeys the
difference equation
ξ (
k + 1) =
a (
k )ξ(
k ) +
b (
k )
v (
k )
η (
k ) =
c (
k )ξ(
k )
where the time varying coefficients
a (
k ),
b (
k ) and
c (
k ) are all periodic with period
two. Calculating the state variable and output
step-by-step leads to:
for any integer
p . By defining "lifted" inputs
and outputs
and considering the system state only at the even
time points ( (
p ) = ξ(2
p )) we arrive at an equivalent
time-invariant system.
The lifting technique has been illustrated above for
a discrete-time system with period two; however, it can
be applied to systems with arbitrary period - though
the corresponding formulae are considerably more
complicated; see [ 18 ] .
Computation of SSV
There is considerable literature in control theory
on the computation of the SSV; see for example [ 20 21
22 ] . For general classes of uncertainty, computing μ
Δ is known to be NP-hard [ 21 ] .
Typically, given the feedback loop consisting of
G and Δ we compute upper and
lower bounds for the SSV [ 15 ] . The lower bound is
exactly equal to μ
Δ [ 15 ] ; unfortunately computing
this lower bound involves a search over a non-convex
set and therefore may converge to local optimums that
are not global. In contrast, the upper bound can be
rewritten in terms of a convex optimisation problem, so
that a global minimum can be obtained. However, this
upper bound is, in general not tight. A software
package is commercially available that can compute μ
upper and uses a power algorithm
to attempt to compute μ
lower [ 22 ] .
Results
Single parameter robustness
The robustness of Laub and Loomis's oscillatory model
was first analysed by means of single-parameter
bifurcation diagrams. Four typical diagrams are shown in
Fig. 5. The activity of internal cAMP (
x
5 ) is plotted as a function of the
variation of each parameter. We use internal cAMP in the
diagram as it is the element that is usually observed
experimentally [ 23 ] . In each of the diagrams, there
are three types of solutions: stable steady state,
unstable steady state and limit cycle oscillations.
These diagrams illustrate that Hopf bifurcations occur
for each parameter; that is, the oscillatory behaviour
exists only in a limited range of parameters around the
nominal value. For each of these parameters, the
respective intervals and values for degree-of-robustness
are found in Table 1.
Structural singular value
From the numerical simulation (Fig. 2) of the
non-linear model, we observed that the oscillatory curves
did not deviate greatly from a simple harmonic oscillator
plus a constant offset. Thus, to obtain an analytic
expression for the periodic orbits we assume that the
state variables are expanded into Fourier series
containing only the fundamental and constant terms:
for each of the seven states. Since it is the relative
phase shift between each state variable that is relevant,
we assume that θ
1 = 0. The substitution of the Fourier
series into the original equations leads to a series of
real algebraic equations for the coefficients (not shown)
whose solution was obtained using Mathematica. This leads
us to obtain the corresponding periodic solutions where
the values of
A
0,
i ,
A
1,
i and θ
i
are found in Table 2. The period
T is approximately 7.31 minutes.
This analytic solution matches well with the numerical
simulation except for an arbitrary phase shift, which
does not affect the shape and location of the limit cycle
in the phase space and can thereby be ignored (not
shown).
Following our prescribed methods, we next linearized
the system. The Jacobian matrix is obtained and was
decomposed as in Eqn. (4) to obtain:
The matrix
B
0 (
t ) = {
B
i ,
j } where
Similarly, the matrix
C
0 (
t ) = {
C
i ,
j } where all coefficients are
zero except for the following:
Finally,
D
0 = 0 and the perturbation structure
is given by
Δ = diag {δ
1 ,δ
2 ,δ
2 ,δ
3 ,δ
4 ,δ
5 ,δ
6 ,δ
6 ,δ
8 ,δ
8 ,δ
9 ,δ
10 ,δ
10 ,δ
11 ,δ
12 ,δ
13 ,δ
14 }
Note that since the nominal trajectory is periodic,
the matrix functions in the nominal description are also
periodic. Note also that in the uncertainty matrix, Δ,
some uncertainties are repeated (δ
2 ,δ
6 ,δ
8 and δ
10 ) while δ
7 is missing.
The system was then discretized and lifted following
the procedure outlined above. A comparison of the system
response for each of the approximations at the nominal
parameter values is given in Fig. 4B.
For the sampling time we tried various values of
n but found negligible differences
for values above eight. Finally, we computed the bounds
on the SSV. Once again, we found that the values of these
two bounds were not affected much by the sampling
frequency provided that
n is greater than eight. The upper
bound was successfully computed using [ 22 ] . The
maximum over all frequencies is approximately 12.06.
However, the high dimension of system causes convergence
problem during the computation of μ
lower using this package. To obtain
an acceptable lower bound, we calculate the spectral
radius at each frequency. This gives us a lower bound for
μ
lower . The plot of the bounds for
the SSV when
n = 16 is shown in Fig. 6. We can
use μ
lower to obtain a conservative
region for robust stability. The highest value over all
frequencies for μ
lower is approximately 2.636.
Discussion
Recent years has seen an appreciation that key cellular
properties are robust to variations in individual parameter
values. Based on the topology of many of these networks,
this should not be surprising. Feedback - both negative and
positive - control systems are ubiquitous in most
biological networks [ 24 ] and one of the reasons for using
feedback is that it reduces sensitivity of a system's
behaviour to its parameter values.
In modelling biological networks, it is important that
this robustness also be in evidence. The particular
behaviour being characterized by the model should not rely
on precise values of the model's parameters - for example,
reaction rate constants or protein concentrations. In
particular, a precise measurement of these constants is
difficult whereas protein concentrations will vary from one
cell to another or throughout the lifetime of any
individual. Deviations from the nominal model parameter
values should not result in a loss of the network's
performance; thus, parameter sensitivity can be used to
validate mathematical models of biochemical system. That
is, the more insensitive the system response is to the
accuracy of the parameter, the more faith we should have in
the model [ 25 ] .
In looking at certain classes of behaviour, where
qualitative changes in the stability of the system are
possible, bifurcation diagrams provide an elegant means of
evaluating robustness. For example, in evaluating the
robustness of the model of Laub and Loomis, of primary
importance is determining whether the oscillatory behaviour
will persist if the parameter values are altered. This
qualitative difference in performance - from limit cycle
oscillations to constant steady states - can be quantified
and compared across parameters or from one model to
another. Once the robustness of the oscillatory behaviour
is established, further investigations of the robustness of
some of the oscillatory features, for example frequency and
amplitude can further be evaluated.
From the bifurcation diagrams obtained for each of the
fourteen parameters, we know that oscillations exist only
in a limited range around the nominal value. We find the
system to be quite sensitive to variations in
k
2 ,
k
4 ,
k
10 and
k
14 and mostly insensitive to the others.
Single-parameter bifurcation analysis also shows that the
amplitude of the oscillation is greatly affected by the
variation of 9 parameters (
k
1 ,
k
2 ,
k
4 ,
k
6 ,
k
7 ,
k
10 ,
k
11 ,
k
12 , and
k
14 ).
Based on the SSV stability to interpret multiple
parameter sensitivity, we can conclude that robust
stability of the periodic orbits will be maintained,
provided that
Since the uncertainty matrix consists only of diagonal
entries, this bound applies to each of the individual
parameters. Thus, we can guarantee that the system will be
robustly stable provided that no single parameter differs
more than 8.3% from its nominal value.
In our analysis we found a large gap between μ
upper and our lower bound for μ. As
we later show, for this system the upper bound is fairly
tight, as we are able to obtain a destabilizing
perturbation of size 9%. For general biological models, a
robustness measure based on the upper bound μ
upper may also be more appropriate.
Robustness bounds for systems in which arbitrarily
slowly-time-varying parameter values are allowed are known
[ 26 ] . For these systems it has been shown that the
bounds converge as the time-variations approach zero to the
upper bound μ
upper [ 26 ] . Since many of the
parameters in models of biochemical networks represent
features that will vary over time, such as enzyme
concentrations, this number may therefore be more
indicative of the model's true robustness.
The ability to consider the effect of time-variations on
the robustness of the system is one great advantage of the
SSV over other methodologies. One drawback of the SSV
approach compared to the bifurcation theory is that it does
not provide the precise combination of parameters that
destabilizes the system - only its size. Also, since the
upper bound is only sufficient to guarantee robustness,
this number may, in general, give an overly conservative
notion of robustness.
It must be emphasized that the SSV approach denoted here
is based on the linearized model of the system. For some
classes of systems this linearization may not be possible -
in this case, the linear SSV approach documented here will
not be applicable. However, for most models used to
describe biochemical reactions, this should not be a
problem.
Because we are concentrating exclusively on the local
stability of the linearized model, important parameters of
the oscillatory behaviour such as robustness of the
frequency and amplitudes of oscillation are not evaluated.
Also, the effect of parameter variations on the equilibrium
orbit are omitted. In particular, varying the kinetic
parameters will change the behaviour of the system in two
different ways: the equilibrium periodic orbit will change
and the stability of deviations about this orbit will also
change. The SSV allows one to quantify the robustness of
the second of these two effects. It does not say anything
directly regarding the effect of parameter variations on
the equilibrium periodic orbit. One way of bounding the
effect of these parameter changes is to write the original
differential equation as
(
t ) =
f (
x ,
k )
where
k =
k
0 + δ is the set of kinetic parameters
with nominal values
k
0 . If the nominal periodic orbit (when
δ = 0) is given by
x *(
t ) =
x (
t ) -
x
δ (
t )
then, linearizing about this orbit yields
δ (
t ) =
A (
t )
x
δ (
t ) +
v
where δ is a constant vector that includes the effect of
this parametric uncertainty. Thus, the system can be
considered as being perturbed by a constant input signal
v . Provided that the homogeneous
system is exponentially stable (and this is guaranteed by
the existence of a stable limit cycle) and that
v is not "too large", the perturbed
system's state will remain in a neighbourhood of the origin
if the
f (
x ,
p ) in the original equation is
reasonably well behaved in
k . Detailed bounds and conditions on
f are given in Theorem 5.1 of [ 17 ]
, though it should be emphasized that these bounds tend to
be overly conservative in practice.
To illustrate the local nature of the SSV analysis for
this system, we perturbed the system parameters by varying
amounts. The particular parameters were either increased or
decreased so as to bring them closer to the Hopf
bifurcation. For example, the nominal value of
k
1 is closer to
k
1 than to
1 so that we reduced
k
1 whereas
k
4 is closer to
4 than to
k
4 so we increased
k
4 . In Fig. 7we show the response of
these systems to changes of 7% and 9% both for the
linearized system - where the linearized response has been
superimposed on the nominal limit cycle (Fig. 7A) and the
original non-linear system (Fig. 7B). For the smaller
value, the linearized response is stable and we see that,
after a transient, the response settles to the nominal
limit cycle. We also see this same behaviour in the
response of the non-linear system with this level of
parameter perturbation. For a 9% change in the parameters,
however, the linearized system is unstable. We see this as
a deviation from the nominal limit cycle. In the non-linear
system's response, this translates into an end to the
stable limit cycle. The response does not "blow up" but
instead settles into a fixed point.
This example illustrates how a robustness analysis of
the linearized system can be used to deduce the robustness
of the original non-linear system, as it shows that when
the linearized system is unstable, the desired behaviour of
the non-linear system will no longer be present. This
example also points out the fact that the upper bound μ
upper ≅ 8.3% is not overly
conservative for this system as we were able to produce a
destabilizing perturbation of size 9%.
Finally, we note that multiparametric robustness
analysis considered here is based on local properties of
the dynamical system, since we are evaluating the
robustness of the linearized model. Extensions to the
non-linear model are the subject of active investigation [
27 ] .
However, it is by combining the robustness analysis of
both single and multiple parameters, we can obtain a more
thorough understanding of the region of stability of the
periodic solution in the high dimensional parameter space
and use this to improve upon the model. In this particular
example, we find that the system's robustness is governed
by several "robustness limiting" parameters,
k
2 ,
k
4 ,
k
10 and
k
14 .
Conclusions
Determining the robustness of mathematical models of
biological systems is important for several reasons. First,
there is growing evidence that many aspects of the networks
being modelled have evolved in such a way so that they are
robust as this allows them to tolerate natural variations
in the environment. Thus, faithful models should replicate
this robustness. Second, robustness of the models provide a
means of validating model quality since the performance of
the models should not rely on precisely tuned parameter
values that are impossible - or at best - difficult to
measure exactly.
In this paper, we illustrated the use of two tools
developed in dynamical systems theory and control
engineering to assess robustness quantitatively. For an
example, we considered an oscillatory molecular network
model due to Laub and Loomis that aims at describing
oscillatory behaviour in cAMP signalling observed in the
social amoeba
D. discoideum. This behaviour appears
as a stable limit cycle of the equations describing the
model. We have evaluated the degree to which this limit
cycle is robust to variations in all the system
parameters.
The robustness of the oscillatory behaviour to single
parameter variations was quantified using bifurcation
analysis. Using the bifurcation analysis software tool AUTO
we determined that single parameter changes as small as 20%
from the nominal value can cause the limit cycle to
disappear and a stable equilibrium to appear. In addition
to the stability robustness, AUTO is also able to evaluate
the sensitivity of the amplitude of the oscillation to
these parameter changes.
To investigate the robustness of the model to
simultaneous changes in parameter values, the structured
singular value (SSV) analysis tool was used. Once the
system was in the correct framework for SSV analysis, we
were able to determine that the system can only tolerate
very small changes in the parameter values - in the order
of 8% - if we allow these parameters to vary with time
arbitrarily slowly.
Finally, it is important to note that to understand
completely the robustness properties of a model, it is
appropriate to combine single and multiple parameter
sensitivity analyses.
Authors' contributions
LM carried out the computational studies and analysis.
PAI conceived of the study and participated in its design
and coordination. All authors read and approved the final
manuscript.
