Sensitivity Analysis and
Robust Experimental
Design of a Signal
Transduction Pathway
System
HONG YUE,
1,4
MARTIN BROWN,
2
FEI HE,
2
JIANFANG JIA,
3
DOUGLAS B. KELL
4
1
Department of Electronic and Electrical Engineering, University of Strathclyde, Glasgow G1 1QD, UK
2
School of Electrical and Electronic Engineering, University of Manchester, Manchester M60 1QD, UK
3
Department of Automation, North University of China, Taiyuan 030051, People’s Republic of China
4
School of Chemistry and Manchester Interdisciplinary Biocentre, University of Manchester, 131 Princess St.,
Manchester M1 7ND, UK
Received 29 October 2007; revised 2 April 2008, 20 April 2008; accepted 20 April 2008
DOI 10.1002/kin.20369
Published online in Wiley I nterScience (www.interscience.wiley.com).
ABSTRACT: Experimental design for cellular networks based on sensitivity analysis is studied
in this work. Both optimal and robust experimental design strategies are developed for the IκB-
NF-κB signal transduction model. Based on local sensitivity analysis, the initial IKK intensity is
calculated using an optimal experimental design process, and several scalarization measures of
the Fisher information matrix are compared. Global sensitivity analysis and robust experimental
design techniques are then developed to consider parametric uncertainties in the model.
The modified Morris method is employed in global sensitivity analysis, and a semidefinite
programming method is exploited to implement the robust experimental design for the problem
of measurement set selection. The parametric impacts on the oscillatory behavior of NF-κBin
the nucleus are also discussed.
C
2008 Wiley Periodicals, Inc. Int J Chem Kinet 40: 730–741, 2008
Correspondence to: Hong Yue; e-mail: hong.yue@eee.strath.
ac.uk.
Contract grant sponsor: National Natural Science Foundation of
China.
Contract grant number: 30770560.
Contract grant sponsor: UK Biotechnology and Biological Sci-
ences Research Council.
Contract grant number: BB/C007158/1.
Contract grant sponsor: Hong Kong Research Grants Council
(to Fei He).
Contract grant numbers: CityU SRG 7001821 and CityU
122305.
c
2008 Wiley Periodicals, Inc.
INTRODUCTION
Sensitivity analysis is used to understand how a
model’s output depends on variations in parameter val-
ues or initial conditions, and is perhaps best known
in metabolic systems biology via metabolic control
analysis [1–4]. It is particularly useful for complex
biological networks that involve a large number of
variables and parameters in which it is crucial to iden-
tify either the most important or the least relevant
parameters.
ANALYSIS AND DESIGN OF SIGNAL TRANSDUCTION PATHWAY SYSTEM 731
Based on the nominal parameter values, local sen-
sitivity analysis (LSA) measures the effects that small
changes in the parameters have on the output. It is
widely used in modeling and analysis of biological
systems, in which the nominal parameter values are es-
timated using experimental data or computation [4–7].
For continuous dynamic systems, the local sensitivi-
ties are defined as the first-order partial derivatives of
the system output with respect to the input parameters.
Such information reveals the gradient of a mathemati-
cal model’s output in parameter space at a given set of
parameter values, and therefore plays a central role in
many system identification problems.
LSA has a wide spectrum of applications in sys-
tems biology. However, for a complex and/or uncer-
tain model in which some parameter estimates are
most likely far from the true values, or for a signifi-
cantly nonlinear and interactive system, it is more rel-
evant to study global sensitivities. Global sensitivity
analysis (GSA) examines the effects of simultaneous
“arbitrary” variations of multiple parameters on the de-
pendent variables under conditions in which the varia-
tions are not local [8–10]. There are different ways to
perform GSA, such as screening techniques, variance-
based methods, Monte Carlo filtering approaches and
regression methods, and so on. In principle, GSA is
valid in a bounded region around the nominal value
for each parameter, and the effect of each parameter
is either aggregated [11] or some worst case measures
are taken for evaluation. It is not simply the result from
weighted local sensitivities, but a multidimensional av-
eraging over the whole parameter space, since when
one parameter is evaluated over its interval, all the
other parameters are also varying instead of keeping
their nominal values. It, therefore, reveals interactions
between parameters from simultaneous parameter vari-
ations. GSA has been applied in modeling, analysis,
and experimental design for a range of biological sys-
tems [12–16].
Optimal experimental design (OED) is one of the
techniques developed from local sensitivities, whose
purpose is to devise the necessary dynamic exper-
iments in such a way that the parameters are esti-
mated from the resulting experimental data with the
best possible statistical quality. The tasks of experi-
mental design include input signal design, sampling
rate optimization, measurement set selection, and so
on. Under the assumption of uncorrelated measure-
ment noise with zero-mean Gaussian distribution, the
information content of measurements can be quanti-
fied by the Fisher information matrix (FIM) [17,18].
In general, the smaller the joint confidence interval is
for the estimated parameters, the more information is
contained in the measurements. In many recent works
on modeling of biochemical networks, the FIM was
used to design the experiments to optimize the quality
of parameter estimation in a certain statistical sense
[18–23]. Several strategies for solving OED problems
in the context of parameter estimation for biochemical
models are discussed in [24].
The quality of optimal experimental design is de-
pendent on the accuracy of the mathematical models.
The true model parameters are in most cases rarely
known, and nominal or estimated parameter values are
used instead. These nominal parameters may be ob-
tained from preliminary experiments, the literature, or
from previous parameter estimation. When the qual-
ity of nominal parameters is poor, the experimental
design results may be overoptimistic or even mislead-
ing. In inverse modeling of complex biochemical net-
works, the normal way to surmount this problem is
to go through an iterative/sequential process for pa-
rameter estimation and experimental design. In each
iteration, the OED is implemented to provide “rich”
information for a better parameter estimation in the
subsequent iteration [19,25]. Using this approach, the
costs associated with experiments for several iterations
are nontrivial, especially for the expensive and time-
consuming data collection in cellular experiments. An
alternative method is the minimax experimental design,
which makes the OED at the worst case in a region
around the nominal parameter values. However, for a
nonlinear model with a large amount of parameters,
identifying the worst case for optimization of the FIM
can be computationally challenging and impractical.
In this work, the problem of experimental design
based on local sensitivities is addressed for biochem-
ical systems with particular interest in models with
parametric uncertainties. An optimal design on the in-
put activation intensity is studied first using the FIM,
so as to illustrate the principle of optimal experimental
design. The endeavors are then focused on the robust
experimental design and global sensitivity analysis, so
as to take into account model uncertainties. With the re-
sults given by local and global sensitivity analyses, the
influence of some important parameters on the oscilla-
tory behavior of NF-κB in the nucleus is investigated.
LOCAL SENSITIVITIES AND OPTIMAL
EXPERIMENTAL DESIGN
System Model: An Example of 1
κB-NF-κB
Signal Pathway
For a biochemical model with n reaction species and
m parameters, denote X = [x
1
x
2
··· x
n
]
T
as the
state vector, θ = [k
1
k
2
··· k
m
]
T
as the vector of
International Journal of Chemical Kinetics DOI 10.1002/kin
732 YUE ET AL.
parameters. The system model can be represented by
ordinary differential equations (ODEs) as
˙
X = f (X, θ, t ),X(t
0
) = X
0
Y = g(X, θ, t) + ω (1)
where f (·) is the nonlinear state transition function,
X
0
is the initial states vector at t
0
, g(·) is the mea-
surement function, Y ∈ ⺢
q
is the measurement output
vector, and ω is assumed to be a zero-mean, Gaussian
noise vector. There are two important assumptions for
using ODEs to represent a cellular system. First, the
system is assumed to be a homogeneous “well-stirred”
reaction system and thus spatial distributions within
each compartment are not considered, otherwise par-
tial differential equations should be used [26]. Second,
the reaction variables (molecular concentrations) are
continuous functions of time. This is valid only if the
number of molecules of each species in the reaction
volume is sufficiently large, otherwise the discrete na-
ture of molecules cannot be neglected and stochastic
differential equations or discrete stochastic models are
more appropriate [27].
The system under investigation is the NF-κB signal-
ing pathway. The NF-κB proteins regulate numerous
genes that play important roles in inter- and intracel-
lular signaling, cellular stress responses, cell growth,
survival, and apoptosis. As such, its specificity and its
role in the temporal control of gene expression are of
crucial physiological interest.
Models of the mainly downstream elements of an
IκB-NF-κB signal pathway have been described by
Hoffmann et al. [28], Nelson et al. [29], and Lipniacki
et al. [30]. While Hoffmann’s model measures the be-
havior of cell populations [28], Nelson’s work analy-
ses experimental data obtained from individual cells
[29], and Lipniacki’s model focuses on the role of
A20, which is regulated by NF-κB and acts as an
inhibitor of IKK [30]. The model used in this work
is the Hoffmann’s model, as described in [28,31–33].
It includes 26 reaction species participating in 64 re-
actions, out of which the concentrations of 24 species
are time varying and are defined as state variables in
the model. See Table I for the states definition. The
reactions and the parameter values can be found in our
previous paper [33].
Local Sensitivity Analysis and Optimal
Experimental Design
Under standard assumptions, parameter estimation
is obtained by solving the following optimization
Table I IκB-NF-κB R eaction Species and States
Species, States Species, States
IκBα, x
1
IKKIκBε-NF-κB, x
14
NF-κB, x
2
NF-κB
n
, x
15
IκBα-NF-κB, x
3
IκBα
n
, x
16
IκBβ, x
4
IκBα
n
-NF-κB
n
, X
17
IκBβ-NF-κB, x
5
IκBβ
n
, x
18
IκBε, x
6
IκBβ
n
-NF-κB
n
, x
19
IκBε-NF-κB, x
7
IκBε
n
, x
20
IKKIκBα, x
8
IκBε
n
-NF-κB
n
, x
21
IKKIκBα-NF-κB, x
9
Source (S = 1)
IKK, x
10
IκBα
−t
, x
22
IKKIκBβ, x
11
Sink (sink = 0)
IKKIκBβ-NF-κB, x
12
IκBβ
−t
, x
23
IKKIκBε, x
13
IκBε
−t
, x
24
problem:
ˆ
θ = arg min
θ
q
i=1
N
l=1
(y
i
(θ,t
l
) −
˜
y
i
(t
l
))
T
×W
l
(y
i
(θ,t
l
) −
˜
y
i
(t
l
)) (2)
where y
i
(θ,t
l
) and
˜
y
i
(t
l
) are model predictions and
measured values of the ith measurable state at time
t,(l = 1,...,N). W
l
is a square matrix with specified
weighting coefficients, which is often chosen as the
inverse of the measurement error covariance matrix Q,
that is, W
l
= Q
−1
. Here, measurement errors with con-
stant variances are considered and under the assump-
tions that the errors are additive, zero-mean Gaussian
processes,
ˆ
θ is the optimal, unbiased estimator.
At time t, the local sensitivity coefficient s
i,j
is de-
fined as the partial derivative of the ith state to the jth
parameter:
s
i,j
(t) =
∂x
i
(t)
∂θ
j
=
x
i
(θ
j
+ θ
j
,t) −x
i
(θ
j
,t)
θ
j
(3)
The sensitivity matrix S = ∂X/∂θ is composed of el-
ements of s
i,j
for i = 1,...,nand j = 1,...,m.The
sensitivity matrix is calculated using the direct differ-
ential method (DDM) where the column sensitivity
vector S
j
= ∂X/∂θ
j
with respect to the jth parameter
(θ
j
) and the system states are obtained simultaneously
by solving the joint state and sensitivity profiles:
˙
X = f (X, θ, t),X(t
0
) = X
0
˙
S
j
= J · S
j
+ F
j
,S
j
(t
0
) = 0
(4)
where J = ∂f /∂ X is the Jacobian matrix and F
j
=
∂f /∂ θ
j
is the parametric Jacobian matrix.
International Journal of Chemical Kinetics DOI 10.1002/kin
ANALYSIS AND DESIGN OF SIGNAL TRANSDUCTION PATHWAY SYSTEM 733
The FIM is a function of the local sensitivity matrix:
FIM =
N
l=1
S(t
l
)
T
Q
−1
S(t
l
)(5)
It is an approximation of the inverse of the parameter
estimation error covariance matrix (), that is, =
FIM
−1
. The lower bound of the variance for the ith
identifiable parameter is given by σ
2
i
≥ (FIM
−1
)
ii
.The
95% confidence interval (CI) for the ith parameter is
represented as CI =
ˆ
θ
i
± 1.96σ
i
.
The covariance matrix is a measure of attain-
able parameter estimation errors for a given set of
data/experimental conditions and therefore is used as
a basis for optimal experimental design. The perfor-
mance index in optimal experimental design is nor-
mally a scalar function of FIM, or equivalently, a
function of the error covariance matrix , and it
should be noted that the design depends on the es-
timated/nominal parameter values. Some commonly
used optimal, scalarization criteria for experimental
design are as follows:
A-optimal : max {trace(FIM)}or min
{
trace
(
)
}
D-optimal : max {det(FIM)}or min
{
det
(
)
}
E-optimal : max {λ
min
(FIM(θ))}or min
{
λ
max
(
)
}
Modified E-optimal : min
λ
max
(FIM)
λ
min
(FIM)
where λ
min
and λ
max
are the minimum and the maxi-
mum eigenvalues and det indicates the determinant of
a matrix. The A-optimal design minimizes the trace of
the covariance matrix (sum of eigenvalues). However,
this criterion is rarely used since it can lead to nonin-
formative experiments when is not positive definite
[34]. The D-optimal design minimizes the determinant
of the covariance matrix (product of eigenvalues) and
can thus be interpreted as minimizing the geometric
mean of the errors in the parameters. The largest error
is minimized by the E-optimal design, which corre-
sponds to a minimization of the largest eigenvalue of
. The modified E-optimal design minimizes the ra-
tio of the largest to the smallest eigenvector and thus
optimizes the functional shape of the confidence in-
tervals. The relationship between these various criteria
has been well studied [35]. Note that all these scalar-
ization criteria result in a convex optimization problem
when the FIM is a linear function of the experimental
design parameters.
Optimal Experimental Design
of IKK Activation Level
To illustrate the principle of optimal experimental de-
sign for the IκB-NF-κB signal pathway, the problem of
selecting the IKK activation intensity using dynamic
(local) sensitivity analysis will be addressed. The FIM
is used along with the Cramer–Rao theorem to deter-
mine the optimal step input signal so that the estimated
parameters have a minimum variance. Taking IKK as
the step input and nuclear NF-κB(x
15
)asthesystem
output, the previously described four criteria are ap-
plied in the design and the results are compared.
To make a better illustration, only three parameters
(k
9
, k
28
, and k
36
) are estimated in the example. These
three parameters are all identified to be influential in
our previous study of LSA [31,33]. The range of the
initial IKK concentration after equilibration is set to
[0.01, 1] μM. The simulation shows that the optimal
IKK input intensity is 0.06 μM with the A-, D-, and
E-optimal designs, and it is 0.01 μM for the modi-
fied E-optimal design. The 95% parameter confidence
intervals for these two results are given in Fig. 1. It
can be clearly seen that while the parametric uncer-
tainty region for the E-optimal design has a smaller
volume for all parameters, the modified E-optimal de-
sign produces a more circular uncertainty region. The
corresponding 95% confidence intervals for each pair
of parameters and their percentage estimation errors
are listed in Table II, for the E-optimal and the mod-
ified E-optimal design. The results in Fig. l and those
in Table II are consistent.
The above design shows that the optimal experi-
mental design is quite straightforward to be imple-
mented once the local sensitivity matrices are estab-
lished. Different optimization criteria can be employed
to determine the best experimental conditions for pa-
rameter estimation. Results from different criteria may
sometimes be different since these criteria represent
scalarized measures of the covariance matrix. Con-
sequently, preference information about the parameter
importance may be used to rank the different solutions.
It can be clearly seen that the OED results depend on
the quality of the model/sensitivity used and the re-
sults can be poor when the model uncertainty cannot
be ignored. In this case, robust strategies need to be
developed for the experimental design.
ROBUST EXPERIMENTAL DESIGN OF
MEASUREMENT SET SELECTION
The objective of robust experimental design (RED)
is to optimally design experiments when there exists
model uncertainty. As with OED, the aim is to de-
sign the experiments so that the uncertainties in the
estimated parameters are as small as possible. How-
ever, when the sensitivity matrix, and hence the FIM,
is a function of the parameters, the design results will
International Journal of Chemical Kinetics DOI 10.1002/kin
734 YUE ET AL.
Table II Ninety-five Percent Confidence Intervals and Estimation Errors
E-Optimal Design Modified E-Optimal Design
Parameters 95% CI Error% 95% CI Error%
(k
9
,k
28
)
ˆ
k
9
[0.02002, 0.02078] 1.87 [0.01491, 0.02589] 26.93
ˆ
k
28
[0.01507, 0.01793] 8.65 [0.01169, 0.02130] 29.11
(k
28
,k
36
)
ˆ
k
28
[0.01567, 0.01732] 5.01 [0.01416, 0.01887] 14.27
ˆ
k
36
[0.00377, 0.00439] 7.72 [0.00315, 0.00501] 22.84
(k
9
,k
36
)
ˆ
k
9
[0.01907, 0.02173] 6.52 [0.01704, 0.02379] 16.54
ˆ
k
36
[0.00389, 0.00428] 4.78 [0.00313, 0.00503] 23.28
Figure 1 Confidence intervals with respect to two different IKK activation intensities calculated by the E-optimal (0.06 μM)
and the modified E-optimal (0.01 μM) design. [Color figure can be viewed in the online issue, which is available at
www.interscience.wiley.com.]
be biased because the parameters are unknown. Se-
quential experimental design strategies that repeatedly
estimate the parameters and redesign the experiments
could be used; however, this is expensive in terms of
the number of experiments. RED is an alternative that
produces an experimental design that is valid across a
prespecified parameter range around the nominal val-
ues.
The problem of measurement set selection is ap-
plicable when there are a large number of states that
could be measured, but experimental constraints mean
that only a small number can actually be measured.
This is represented as the selection problem:
ξ =
x
1
··· x
n
ω
1
··· ω
n
,
n
i=1
ω
i
= 1,ω
i
≥ 0, ∀i (6)
where ω
i
is the nonnegative weight relating to the ith
state x
i
. The selection problem is originally an integer-
programming problem as the weights should be binary
variables; however, this is relaxed to produce an ap-
proximate, continuous optimization problem. RED es-
timates the weights and consequently focuses on those
measurements with larger weight values, that is, state
measurements that are more informative for parame-
ter estimation. This is a vector optimization problem
over the positive semidefinite cone, for which several
scalarizations have been proposed, such as D-optimal,
A-optimal, E-optimal, and so on.
Semidefinite Programming and Robust
Experimental Design
Ideally the design should start from parametric un-
certainties and map these uncertainty information into
uncertainty bound in FIM. This mapping is difficult
to be established explicitly, but the link via Taylor ex-
pansion can be managed. Considering a local model
with additive uncertainties to the model parameters,
the corresponding parametric sensitivity matrix can be
represented by a simple truncated first-order Taylor
expansion:
S(θ + θ) ≈ S(θ) +
dS(θ)
dθ
θ (7)
International Journal of Chemical Kinetics DOI 10.1002/kin
![Figure 4 Distribution of weight values at three uncertainty levels, ρ = 1e − 6 corresponds to small uncertainty, ρ = 10 corresponds to large uncertainty, andρ = 1e − 2 corresponds to the uncertainty level in between the other two. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]](/figures/figure-4-distribution-of-weight-values-at-three-uncertainty-cbowr3tz.webp)
![Figure 3 Change of weights with respect to ρ. These weights are ω2, ω5, ω11, ω8, ω24, ω12, ω7, ω13, ω14, ω3, ω18, ω20, ω6, ω4, ω2, ω15, ω19, ω23. Compared with the six weights in Fig. 2, they have smaller values when ρ is small. When ρ is large, all the weights also have the same value of 0.417. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]](/figures/figure-3-change-of-weights-with-respect-to-r-these-weights-2biib9w2.webp)
![Figure 2 Change of weights with respect to ρ. These six weight variables have larger values when ρ is small compared with the rest of the weight variables. When ρ is large, all the weights have the same value of 0.417. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]](/figures/figure-2-change-of-weights-with-respect-to-r-these-six-y8pflde8.webp)
![Figure 8 NF-κBn oscillation with respect to variations in k28. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]](/figures/figure-8-nf-kbn-oscillation-with-respect-to-variations-in-1vffst3f.webp)
![Figure 7 NF-κBn oscillation with respect to variations in k36. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]](/figures/figure-7-nf-kbn-oscillation-with-respect-to-variations-in-1bhkmfwa.webp)
![Figure 9 NF-κBn oscillations with respect to variations in k61. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]](/figures/figure-9-nf-kbn-oscillations-with-respect-to-variations-in-2v2cl4fw.webp)