BMC Systems Biology BioMed Central

Research article Open Access

Dynamical pathway analysis
Hao Xiong* and Yoonsuck Choe

Address: Department of Computer Science, Texas A&M University, College Station, TX 77843, USA
Email: Hao Xiong* - hxiong@cs.tamu.edu; Yoonsuck Choe - choe@cs.tamu.edu
* Corresponding author Equal contributors

Published: 27 January 2008 Received: 28 August 2007

Accepted: 27 January 2008
BMC Systems Biology 2008, 2:9 doi:10.1186/1752-0509-2-9
This article is available from: http://www.biomedcentral.com/1752-0509/2/9
2008 Xiong and Choe; licensee BioMed Central Ltd.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract
Background: Although a great deal is known about one gene or protein and its functions under
different environmental conditions, little information is available about the complex behaviour of
biological networks subject to different environmental perturbations. Observing differential
expressions of one or more genes between normal and abnormal cells has been a mainstream
method of discovering pertinent genes in diseases and therefore valuable drug targets. However,
to date, no such method exists for elucidating and quantifying the differential dynamical behaviour
of genetic regulatory networks, which can have greater impact on phenotypes than individual genes.
Results: We propose to redress the deficiency by formulating the functional study of biological
networks as a control problem of dynamical systems. We developed mathematical methods to
study the stability, the controllability, and the steady-state behaviour, as well as the transient
responses of biological networks under different environmental perturbations. We applied our
framework to three real-world datasets: the SOS DNA repair network in E. coli under different
dosages of radiation, the GSH redox cycle in mice lung exposed to either poisonous air or normal
air, and the MAPK pathway in mammalian cell lines exposed to three types of HIV type I Vpr, a wild
type and two mutant types; and we found that the three genetic networks exhibited fundamentally
different dynamical properties in normal and abnormal cells.
Conclusion: Difference in stability, relative stability, degrees of controllability, and transient
responses between normal and abnormal cells means considerable difference in dynamical
behaviours and different functioning of cells. Therefore differential dynamical properties can be a
valuable tool in biomedical research.

Background genes in disease pathogenesis [2], but these methods offer

Cell functions are complex temporal processes and only static views and steady-state explanations and thus
should be studied as complex dynamical processes rather fail to account for the transient behaviours that influence
than only in their individual steady states. It is increas- phenotypes. Genetic regulatory networks seek to model
ingly recognized that it is the dynamics and the internal complex interactions and dynamics of gene regulations.
structures of the biological systems that give rise to the Genetic networks should behave differently in sick cells
functioning of cells [1]. Currently, uncovering co- vs. healthy cells because genes that cause diseases behave
expressed genes and discovering differentially expressed fundamentally differently, and that difference should be
genes are the primary methods for discovering the role of reflected in their dynamical properties. Dynamical prop-

Page 1 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

erties of genetic networks such as their response time have the inputs and the outputs can be environmental stimuli
been studied mostly in the context of network motifs or expression levels. Because genetic networks have many
[3,4], but now we propose that they be investigated for unknown quantities, state-space models can serve as a
their difference in normal vs. abnormal cells. good modelling framework.

In this report we studied four dynamical properties: stabil- In this paper, the parameters in state-space models were
ity, relative stability, controllability, and transient behav- estimated from the time course of gene expressions using
iours (overshoot, settling time, and rise time). Stability Kalman filter and the constrained expectation-maximiza-
governs how a system responds to internal noise and tion (EM) algorithm (a modified EM algorithm that incor-
external perturbation and determines whether the system porates prior knowledge about the structure of genetic
returns to steady states and whether the effect of noise and networks). The regular EM algorithm is commonly used
perturbation diminishes over time. Biologically, an unsta- to estimate parameters in the presence of hidden quanti-
ble cellular system is very brittle and the slightest distur- ties, and they comprise two steps, E-step (expectation)
bance can drive the system beyond tolerance and possibly and M-step (maximization), where the E-step estimates
result in cell death. Prill et al. [5] used stability as a crite- the hidden states, and the M-step the parameters [23]. We
rion to discern network motifs and their organizing prin- applied EM algorithm to three sets of time course data and
ciples, and synthetic biologists are beginning to pay close estimated three genetic networks for analysis.
attention to the stability of their artificial networks [6].
Furthermore, the stability of the system under pure gain The first network we used is the SOS DNA repair system.
feedback control can be analyzed by the root-locus The SOS network is a highly conserved system [8,24] and
method and the result can be interpreted as a measure of consists of about 30 genes, the master regulator being
relative stability. In control theory, the root-locus method gene lexA. The lexA gene inhibits the rest of the SOS net-
is a design tool but it is also used as an analytic tool, to see work's genes under normal conditions, but when DNA
how large a gain can drive the system unstable with feed- damage is sensed, protein LexA is cleaved and the genes
back loops: the larger margins of stabilizing gains, the bet- normally suppressed are activated. A diagram of the SOS
ter. Related to feedback control, controllability is another network with 8 essential genes is shown in Fig. 1. Shown
pivotal concept in control theory. It and its dual property, in Fig. 2 is the second system we modelled, the glutath-
observability, were originally conceived as solutions to ione (GSH) redox cycle with one gene from the urea cycle
existence and uniqueness problems of optimal control that interacts with the redox cycle [25]. The data are from
[7], and the controllability of a dynamical system roughly Sciuto et al. [26] who investigated the differential gene
refers to the ability to move the states of the system expressions in mice lung cells exposed to either carbonyl
around the state space with reasonable efforts. Although chloride (phosgene) or normal air. They found elements
controllability is a binary question, there is a measure of of the GSH redox cycle differentially expressed, which is
the degree of controllability, the idea being that the more not surprising given that the redox cycle is heavily
controllable a system is the less effort is needed to move involved in protecting organisms from reactive oxygen
the system. Less theoretical than stability and controllabil- species, that it is heavily present in the lung, and that
ity are transient behaviours like settling time and over- phosgene causes massive lung damages. The third system
shoots, which have also received attention from systems we investigated is the mitogen-activated protein kinase
biologists [3,4,8]. These four dynamical properties are
determined by the parameters of the dynamical system
and the unknown parameters of biological systems need
to be estimated.

Parameter estimation must be done under a particular

modelling framework. Several modelling frameworks
have been proposed: Boolean networks [9-12], differen-
tial equations [13], S-system [14,15], and dynamical
Bayesian networks [16,17]. A special case of dynamical
Bayesian networks is the state-space model, which has
been used to model genetic regulatory networks [18-22].
A state-space model has states, inputs, and outputs, where
hidden states contain complete information of the system 

driven by the inputs, and the outputs are the measure- Figure
The diagram
1 of SOS DNA repair network
ments made by scientists. In the state-space models of The diagram of SOS DNA repair network.
genetic networks, states are the regulatory elements, and

Page 2 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

Figure
The diagram
2 of GSH redox cycle
The diagram of GSH redox cycle.

(MAPK) network in cell lines disturbed by either the wild Figure

The diagram
3 of MAPK network
The diagram of MAPK network.
type HIV type I Vpr or the mutant type R73A or the mutant
type R80A. HIV-1 Vpr is an important protein in promot-
ing the pathogenesis of AIDS by facilitating apoptosis and
cell cycle stall at G2. Yoshizuka et al. [27] studied the where x(t) is the state vector, y(t) the output vector, and
effects of Vpr on MAPK-network-related genes in stalling u(t) the input vector, all at time t; w and v are independent
cell cycle, so they obtained cell lines that can express wild noise terms assumed to be white Gaussian with zero
type or mutant Vpr under a tetracycline-inducible pro- means and covariance Q and R, respectively. Matrix A is
moter. They found that many genes related to the MAPK called the state transition matrix, B the input matrix, C the
network differentially expressed when subjected to differ- output matrix, and D the feed-forward matrix. Matrices A,
ent types of Vpr. The MAPK network used for this report is B, C, D and covariance matrices Q and R together make up
shown in Fig. 3. All those data sets compare the organ- the parameters of the dynamical system.
ism's reactions to different environmental perturbations,
and from estimated genetic networks we hope to discover The states represent the biological forces that regulate
the differential dynamical properties of genetic networks gene regulation; they describe the behaviours of gene tran-
under stress. scription but are hidden. The outputs denote the gene
expression levels and are measured, and it is assumed that
We applied our framework to three real-world time series the expression level of a gene is determined by the state of
datasets above and found differential stability, transient the regulated gene. The inputs can be any external stimuli
responses, and controllability of genetic networks in nor- that influence gene regulation: substances like drugs, pro-
mal vs. abnormal cells. teins, RNAs, or expression levels of other genes.

Results Estimated system

Models of genetic networks and their application to real For the SOS system, x2 is the discretized first derivative of
data sets x1, whereas x1 is the expression level of gene lexA, x3 gene
We modeled genetic networks as dynamical systems, polB, x4 gene umuD, x5 gene uvrD, x6 gene uvrA, x7 gene
more specifically as linear state-space systems. A linear uvrY, and x8 gene ruvA. The outputs are the measured
state-space model of dynamical systems can be written as expression levels of the seven genes listed above, and the
input is gene recA. In Fig. 4 and Fig. 5, we included the
x(t + 1) = Ax(t ) + Bu(t ) + w estimated outputs and the measured outputs superim-
(1) posed into one plot, as well as estimation errors in a sep-
y(t ) = Cx(t ) + Du(t ) + v
arate panel for each gene. From the plots we can see that
the estimated trajectory largely follows measured values.

Page 3 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

a b c

d e f

Figure 4 expression levels and error in estimation

Estimated
Estimated expression levels and error in estimation. For each gene of the SOS system, we have superimposed the esti-
mated expression levels on measured expression levels and plotted the error in estimation in a separate panel. We have done
this for the low radiation level data set in Figure 4. The estimations generally show good behaviors. Figs. 4a, 4b, 4c, 4d, 4e, 4f,
4g are for the low radiation level data set, and plotted genes lexA, polB, umuD, uvrD, uvrA, uvrY, and ruvA, respectively. Each
gene has two plots; the bottom panel shows estimated expression level superimposed on measured expression level, while the
top panel is the estimation error.

The estimated system parameters are listed below for the

x1(t + 1) = x1(t ) + x 2(t ) y1(t ) = x1(t )
low level of radiation:
x 2(t + 1) = 0.17 x1(t ) + 0.59x 2(t ) + 0.084u(t ) y 2(t ) = x 3(t )
x 3(t + 1) = 0.009x(t )1 + 0.81x 3(t ) y 3(t ) = x 4 (t )
x 4 (t + 1) = 0.037 x1(t ) + 0.74 x 4 (t ) y 4 (t ) = x 5(t )
x 5(t + 1) = 0.007 x1(t ) + 0.964 x 5(t ) y 5(t ) = x 6(t )
x 6(t + 1) = 0.037 x1(t ) + 0.965x 6(t ) y 6(t ) = x 7(t )
x 7(t + 1) = 0.008 x1(t ) + 0.697 x 7(t ) y 7(t ) = x 8(t )
x 8(t + 1) = 0.009x1(t ) + 0.621x 8(t ).
(2)

For the high level of radiation, the estimated system is

Page 4 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

a b c

d e f



Figure 5 expression levels and error in estimation

Estimated
Estimated expression levels and error in estimation. We have superimposed the estimated expression levels on meas-
ured expression levels and plotted the error in estimation in a separate panel for the high radiation level data set in Figure 5.
Figs. 5a, 5b, 5c, 5d, 5e, 5f, 5g are for the high radiation level data set plotted genes lexA, polB, umuD, uvrD, uvrA, uvrY, and
ruvA, respectively. Each gene has two plots; the bottom panel shows estimated expression level superimposed on measured
expression level, while the top panel is the estimation error.

elled with second order dynamics so the last four states x5,
x1(t + 1) = x1(t ) + x 2(t ) y1(t ) = x1(t ) x6, x7 and x8 are the discretized first derivatives of x1, x2, x3
x 2(t + 1) = 0.242 x1(t ) + 0.329x 2(t ) 0.014u(t ) y 2(t ) = x 3(t ) and x4, respectively. Here, gene GCLC is x1, gene GCLM x2,
x 3(t + 1) = 0.008 x1(t ) + 0.832 x 3(t ) y 3(t ) = x 4 (t ) gene GSS x3, and gene IDH2 x4. The estimated system for
exposure to normal air is
x 4 (t + 1) = 0.051x1(t ) + 0.653x 4 (t ) y 4 (t ) = x 5(t )
x 5(t + 1) = 0.01x1(t ) + 0.889x 5(t ) y 5(t ) = x 6(t )
x 6(t + 1) = 0.366 x1(t ) + 1.22 x 6(t ) y 6(t ) = x 7(t )
x 7(t + 1) = 0.002 x1(t ) + 0.906 x 7(t ) y 7(t ) = x 8(t )
x 7(t + 1) = 0.002 x1(t ) + 0.906 x 7(t ).
(3)

For the GSH redox cycle there are two inputs, gene
ALD2A1 as u1 and GPX4 as u2. All the states were mod-

Page 5 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

x1(t + 1) = x1(t ) + x 5(t ) x1(t + 1) = x1(t ) + x 5(t )

x 2(t + 1) = x 2(t ) + x 6(t ) x 2(t + 1) = x 2(t ) + x 6(t )
x 3(t + 1) = x 3(t ) + x 7(t ) y1(t ) = x1(t ) x 3(t + 1) = x 3(t ) + x 7(t ) y1(t ) = x1(t )
x 4 (t + 1) = x 4 (t ) + x 8(t ) y 2(t ) = x 2(t ) x 4 (t + 1) = x 4 (t ) + x 8(t ) y 2(t ) = x 2(t )
x 5(t + 1) = 0.37 x1(t ) 0.39x 5(t ) + 0.814u1(t ) y 3(t ) = x 3(t ) x 5(t + 1) = 0.582 x1(t ) 0.821x 5(t ) + 0.082u1(t ) 0.085u 2(t ) y 3(t ) = x 3(t )
x 6(t + 1) = 0.429x 2(t ) 0.006 x 6(t ) + 0.632u 2(t ) y 4 (t ) = x 4 (t ) x 6(t + 1) = 0.28 x 2(t ) 0.836 x 6(t ) 0.149u1(t ) 0.009u 2(t ) y 4 (t ) = x 4 (t )
x 7(t + 1) = 0.095x1(t ) 0.015x 2(t ) 0.217 x 3(t ) 0.128 x 7(t ) x 7(t + 1) = 0.019x1(t ) + 0.273x 2(t ) 0.056 x 3(t ) 0.249x 7(t )
x 8(t + 1) = 0.753x 3(t ) 0.409x 4 (t ) 0.867 x 8(t ) + 0.017u 2(t ) x 8(t + 1) = 0.112 x 3(t ) 0.467 x 4 (t ) 1.248 x 8(t ).

(4) (8)
For exposure to phosgene, the estimated model is Although the number of parameters is small compared
with the number of states, which agrees with the knowl-
x1(t + 1) = x1(t ) + x 5(t ) edge that genetic networks are sparse [28], it is still hard to
x 2(t + 1) = x 2(t ) + x 6(t ) see at a glance whether they differ in any fundamental
x 3(t + 1) = x 3(t ) + x 7(t ) y1(t ) = x1(t ) way. For that, we must apply systematic analysis to the
x 4 (t + 1) = x 4 (t ) + x 8(t ) y 2(t ) = x 2(t ) estimated systems.
x 5(t + 1) = 0.141x1(t ) + 1.95x 5(t ) + 0.77u1(t ) y 3(t ) = x 3(t )
x 6(t + 1) = 0.076 x 2(t ) 1.09x 6(t ) 0.13u 2(t ) y 4 (t ) = x 4 (t ) Differential stability of systems under different
x 7(t + 1) = 0.05x1(t ) 0.161x 2(t ) 0.09x 3(t ) 0.705x 7(t )
perturbations
x 8(t + 1) = 0.336 x 3(t ) 0.179x 4 (t ) 1.126 x 8(t ) + 0.103u 2(t ).
Stability is a very important property of a biological sys-
(5) tem, for an unstable system puts great stress on neigh-
As for the MAPK system, the inputs are gene BRAF as bouring systems and may even lead to cell death. A system
u1and gene RAF1 as u2. The states x1, x2, x3 and x4 are genes is stable if it will converge to steady states after distur-
MAP2K1, MAP2K2, MAPK1, and MKNK2, respectively; bance; it is unstable otherwise. The stability of a discrete
the other four states are the discretized first derivatives as linear system can be determined by the eigenvalues of its
in the system for the GSH redox cycle. The estimated sys- state transition matrix A: if all the eigenvalues are within
tem for the wild type Vpr is the unit circle in the complex plane, then the discrete sys-
tem is stable. The eigenvalues of the three analyzed net-
x1(t + 1) = x1(t ) + x 5(t )
works are listed in Table 1, 2, and 3, and their
x 2(t + 1) = x 2(t ) + x 6(t ) implications discussed below.
x 3(t + 1) = x 3(t ) + x 7(t ) y1(t ) = x1(t )
x 4 (t + 1) = x 4 (t ) + x 8(t ) y 2(t ) = x 2(t ) We analyzed the SOS DNA network under low and high
x 5(t + 1) = 1.48 x1(t ) 1.32 x 5(t ) + 0.14u1(t ) + 0.2u 2(t ) y 3(t ) = x 3(t ) dosage of radiation and discovered that the network was
x 6(t + 1) = 0.098 x 2(t ) 0.52 x 6(t ) 0.079u1 0.314u 2(t ) y 4 (t ) = x 4 (t ) stable for low dosage and unstable for high dosage. We
x 7(t + 1) = 0.498 x1(t ) + 0.052 x 2(t ) 0.215x 3 0.618 x 7(t ) found that the eigenvalues of SOS network under low dos-
x 8(t + 1) = 0.123x 3(t ) 0.169x 4 (t ) 0.602 x 8(t ), age of radiation to have the eigenvalues' norm all less than
(6) one, and therefore the network was stable. On the other
hand, the SOS network was unstable under high dosage of
and for the R73A mutant
radiation, as the norm of one of its eigenvalues was greater
than one.
x1(t + 1) = x1(t ) + x 5(t )
x 2(t + 1) = x 2(t ) + x 6(t )
We also analyzed the redox cycle in mice lung cells that
x 3(t + 1) = x 3(t ) + x 7(t ) y1(t ) = x1(t )
were exposed to either carbonyl chloride (phosgene), an
x 4 (t + 1) = x 4 (t ) + x 8(t ) y 2(t ) = x 2(t )
x 5(t + 1) = 1.098 x1(t ) 1.087 x 5(t ) + 0.068u1(t ) 0.23u 2(t ) y 3(t ) = x 3(t ) Table 1: Differential stability of the SOS network
x 6(t + 1) = 0.76 x 2(t ) x 6(t ) 0.06u1(t ) 0.03u 2(t ) y 4 (t ) = x 4 (t )
x 7(t + 1) = 0.073x1(t ) + 0.5x 2(t ) 0.355x 3(t ) 0.647 x 7(t ) Low Dosage High Dosage
x 8(t + 1) = 0.79x 3(t ) 0.876 x 4 (t ) 1.179x 8(t ),
(7) 0.8117 0.8321
0.7367 0.6530
and for R80A mutant 0.9637 0.8893
0.9652 1.2216 (unstable)
0.6969 0.9062
0.6219 0.6291
0.7952 + 0.3630i 0.6647 + 0.3597i
0.7952 - 0.3630i 0.6647 - 0.3597i

Page 6 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

Table 2: Differential stability of GSH redox cycle by and poles; the zeros and poles of the overall system are
called the closed-loop zeros and poles. A dynamical sys-
Normal Air Phosgene
tem's zeroes are the roots of the numerator of the transfer
-0.6141 2.0830 (unstable) function (for an explanation of the transfer function, see
0.1177 -1.0383 (unstable) Methods), and the poles are the roots of the denominator.
0.4803 + 0.3199i -0.7561 The stability of closed-loop systems depend on the closed-
0.4803 - 0.3199i 1.0512 (unstable) loop poles. The root-locus method generates a plot that
0.7467 0.9120 traces the closed-loop poles as the gain of the controller is
0.7542 0.8696 varied, and the portion of gains that make the closed-loop
0.4972 + 0.4196i 1.0470 + 0.2711i
stable is called the stability margin. In the root locus plot,
0.4972 - 0.4196i 1.0470 - 0.2711i
the open-loop zero is represented by a circle ( ), the open-
loop pole by a cross (), and if there is a zero-pole cancel-
industrial toxin, or normal air; and we found that GSH lation we will see a circle and a cross on top of each other
redox system in normal lung cells was stable all eigen- (). The root-locus method can only study systems with
values were within the unit circle, and that the network single input and single output (SISO), but the dynamical
exposed to toxin was unstable some eigenvalues were properties of SISO systems is a reflection of the overall sys-
outside the unit circle. Whether the unstable detoxifica- tem's dynamical properties, so that the performance of the
tion system contributed to the death of mice exposed to SISO system will manifest itself in the overall system's per-
phosgene is not yet known, but Sciuto et al. [26] specu- formance.
lated that the poison might have overwhelmed the detoxi-
fication system. In the SOS DNA repair network, the recA to uvrA SISO sys-
tem showed differential root-locus plots, depending on
We also analyzed the activity data from the MAPK net- radiation levels. Their respective root-locus plots, for both
work in mammalian cells that expressed either wild type negative and positive feedbacks, are shown in Fig. 6a, 6b,
Vpr, mutant R73A Vpr, or mutant R80A Vpr; and we 6c, and 6d. Under low level of radiation, we found that
found that both the wild type and R73A produced stable the SISO system was comfortable with positive feedback,
behaviours, and R80A caused the network to become which had larger margin of stabilizing gains, whereas neg-
unstable. A stable MAPK network helps the virus most, for ative feedback allowed far narrower choices. Under high
Yoshizuka et al. [27] found the HIV virus uses MAPK net- level of radiation, the opposite was true: positive feedback
work to cause cell cycle G2 arrest, and over-expression of had no stabilizing gain whereas negative feedback had a
MAP2K2 reversed the arrest. large margin. The need for positive feedback loop in low
radiation level is an interesting discovery from our root-
Differential relative stability analyzed via root locus locus analysis, because it runs counter to the common
The relative stability of genetic networks is also important; perception that negative feedback loop promotes stability
it is a measure of robustness. We studied relative stability and positive feedback loop leads to instability. Perhaps
by examining the stability margins of pure gain feedback under low radiation level, the SOS network is not suffi-
loops through root-locus plots. Given a dynamical sys- ciently stimulated and positive feedback fully activates the
tem, one forms a feedback loop from the output to the network which then leads to overall stability.
input through only a pure gain controller. Depending on
whether the control signal is negated as it is fed into the In the GSH redox network, we discovered that the
input, the feedback can be positive (not negated) or nega- ALDH2A1 to IDH2 SISO system showed a simpler but
tive (negated). The original system is called the open-loop more striking difference under different environmental
system, and its zeros and poles are the open-loop zeros conditions. When exposed to normal air, the SISO system
was stable and the root-locus plot in Fig. 7a shows that
Table 3: Differential stability of the MAPK network sizeable gain values do not destabilize the closed-loop
system, which represents a nice scenario, because the sub-
Wild type R73A R80A
system can sustain a lot of stress. But, as we can see in Fig.
7b and Fig. 7c, the same SISO system, when exposed to
0.8527 0.7448 -1.0169 (unstable)
0.8862 -0.0437 + 0.0925i -0.4078 toxin, not only had an unstable open-loop system, but the
-0.1615 + 0.3646i -0.0437 - 0.0925i -0.2023 closed-loop system also remained unstable no matter
-0.1615 - 0.3646i -0.6472 0.5867 what value of the gain was, positive or negative. This
-0.4884 -0.3913 0.9534 means that not only the ALDH2A1 to IDH2 SISO system
-0.4601 0.4680 0.7685 was very unstable, but that a higher order controller must
0.9324 0.4877 -0.6676 be used to produce a stable closed-loop system, a sign of
-0.4477 -0.4916 0.8315
very serious damage.

Page 7 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

a b

c      d

 
Figure
Root locus
6 plots of recA to uvrA SISO system
Root locus plots of recA to uvrA SISO system. These plots trace the poles of the closed-loop system as the gain K is var-
ied from zero to infinity. The trajectories start at the open-loop poles which are represented by the cross, and could end at
the open-loop zeros which are represented by an open circle, or they could go on infinitely in some direction. The different
colours represent distinct trajectories of different closed-loop poles. Fig. 6a This is the root locus plot of recA to uvrA system
under low level of radiation with negative feedback, where the locus on the real axis goes out of the unit circle quickly and
therefore shows small stability margins. (The dotted circle is the unit circle.). Fig. 6b This is the root locus plot of recA to uvrA
system under low level of radiation with positive feedback, with some stability margins. Fig. 6c This is the root locus plot of
recA to uvrA system under high level of radiation with negative feedback, where a good portion of all three loci stays within
the unit circle and therefore exhibits large stability margins. Fig. 6d This is the root locus plot of recA to uvrA system under
high level of radiation with positive feedback, which has no stability margin.

We also found that MAPK network in mammalian cell margin of gain values for which the closed-loop system
lines subject to different versions of Vpr of HIV type I virus was stable. The SISO system under R80A mutant protein
had similarly different root locus plots, which are shown exhibited a stable closed-loop system with only a small
in Fig. 8a, 8b, 8c, and 8d. The RAF1 to MKNK2 SISO sys- margin of gain with positive feedback and none with neg-
tem was stable under both the wild type and the R73A ative feedback. If that small margin does not include a
mutant Vpr perturbation, and both showed comfortable

Page 8 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

a b

c

Figure
Root locus
7 plots of recA to uvrA SISO system
Root locus plots of recA to uvrA SISO system. Fig. 7a Root locus plot of ALDH2A1 to IDH2 system exposed to normal
air with negative feedback is shown, with large stability margins. Fig. 7b Root locus plot of ALDH2A1 to IDH2 system exposed
to poisonous air with negative feedback, where the locus on the positive real axis is entirely outside of the unit circle and
therefore it has no stability margin. Fig. 7c Root locus plot of ALDH2A1 to IDH2 system exposed to poisonous air with posi-
tive feedback, showing no stability margin because of the locus at the right.

gain that can produce a closed system with satisfactory inputs is called controllability, which is a pivotal concept
performance, then a higher order controller is called for. in linear time systems theory [7]. Controllability can be
tested by the rank of controllability matrix; if the control-
Differential degree of controllability lability matrix is of full rank, then the system is controlla-
Since one goal of systems biology is to aid the develop- ble, otherwise uncontrollable. Beyond the binary test
ment of therapeutic treatments, which in the context of (controllable or not) there are also degrees of controllabil-
genetic networks is to bring the network from undesirable ity. The condition number of the controllability matrix
states to healthy states by manipulating inputs, the rela- can be considered as a measure of the degree of controlla-
tive ease of moving around in the state space is an impor- bility, the bigger the number the less the controllability. A
tant issue. The ability to move a system from one point in system with less controllability may require much greater
the state space to another in finite time with only finite inputs to achieve the desired final state, which could be a

Page 9 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

a b

c d

Figure
Root locus
8 plots of RAF1 to MKNK2 SISO system
Root locus plots of RAF1 to MKNK2 SISO system. Fig. 8a Root locus plot of RAF1 to MKNK2 system perturbed by wild
type Vpr with negative feedback where a large portion of the locus can be seen within the unit circle. Fig. 8b Root locus plot of
RAF1 to MKNK2 system perturbed by R73A mutant Vpr with negative feedback showing very good stability margins. Fig. 8c
Root locus plot of RAF1 to MKNK2 system perturbed by R80A mutant Vpr with negative feedback, where there is basically no
stability margin due to the two loci on the real axis. Fig. 8d Root locus plot of RAF1 to MKNK2 system perturbed by R80A
mutant Vpr with positive feedback and small stability margins.

problem as the inputs for biological systems are drugs, tion numbers differed, for one significantly. We
radiation therapy, things in limited supply and subject to discovered that the SOS DNA repair system under high
cost factors. As we will see, different systems could have dose of radiation had a condition number of 2.8109 for
radically different degrees of controllability. its controllability matrix, and that under low dosage the
condition number was 2.6109. The similarly large condi-
Although we found all the three genetic regulatory net- tion numbers suggest the SOS system under study is diffi-
works controllable under all circumstances, their condi- cult to control; whether this is due to radiation is not

Page 10 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

known. On the other hand, in mice lung exposed to nor- differential rise time and settling time, the recA to uvrD
mal air we saw that the redox system had a condition system. The SISO system, when exposed to high radiation
number of 567 for its controllability matrix, and that dosage, was almost twice as fast as the system exposed to
those exposed to toxin had 70267. The different condition low dosage of radiation, to reach their respective steady
numbers peg the redox system as much more difficult to states. This suggests that the SOS system needs uvrD to
control after exposure to poison, perhaps due to damages respond faster to, and therefore has faster dynamics
or the fact that the network was being overwhelmed by the under, higher levels of radiation. With no overshoot in
effects of the toxin. The third network, the MAPK network both cases and a smaller settling time for a higher dosage,
in mammalian cell lines, was found to have a condition the recA to uvrD system under high radiation level stayed
number of 62.15 when exposed to the wild type Vpr, 88.5 closer to the steady states. The rise time and settling time
for those exposed to the R80A mutant, and 285.4 for the are listed in Table 4.
R73A mutant. It is obvious that R73A mutant results in a
stodgier MAPK network than other variants, but overall The MAPK network in mammalian cells exhibited differ-
the MAPK system retains good controllability, making it a ential transient responses to three types of Vpr of HIV type
good target for treatment. I virus. The BRAF to MAP2K2 SISO system displayed
slower dynamics and were more distant from the steady
Differential transient responses state under the wild type than both mutants, and among
To study cell functions as temporal processes means we the mutants, R73A had faster dynamics and better ability
must take stock of transient behaviors in addition to to stay close to the steady state. On the other hand, the
steady states. One way to characterize transient behaviors BRAF to MAPK1 system's transient behaviors in response
is through the transient response of the dynamical system to the wild type Vpr were dominated by a 440% over-
to inputs like step input and impulse input, but because shoot, and with its long settling time the system's tran-
the step responses and impulses responses give same sient responses were far from the steady state. The system
information for linear systems, we will concentrate on the perturbed by the wild type protein also had faster dynam-
step input responses. A step input is a constant input, a ics due to its smaller rise time, and the R73A mutant pro-
unit step, a constant unity. The rise time is a measure of duced a system that had relatively fast dynamics and
the speed of the dynamics, and the settling time and the transient response closer to the steady state. The R80A
overshoot gauge how close to the steady state the transient mutant resulted in a system with slow dynamics and tran-
behaviors are. Of course, systems that exhibit differential sient responses distant from the steady state with its rela-
stability will have different transient responses, but tive large rise time and settling time and no overshoot.
because differential stability is addressed earlier, we will The respective rise time, settling time, and overshoots are
disregard any difference in transient responses due to sta- in Table 5.
bility difference.
Although overshoot is generally considered undesirable
The transient responses are by their nature studied as in engineering (whether fast dynamics or staying close to
input-output pairs, also called a single-input-single-out- the steady states are good or bad depends on the circum-
put (SISO) system. Although we will look at individual stances and cannot be determined a priori;) a large over-
SISO system extracted from multiple-input-multiple-out- shoot can be a fast way of signalling, or it can be an
put (MIMO) systems, the transient responses are still the unbearable disturbance to cells. But being aware of the
intrinsic properties of the original system, and differential difference in transient responses is the first step toward
transient responses suggest fundamentally different devising treatment strategies that shape biological sys-
dynamical behaviors of the original system in response to tems' dynamics to our liking.
external perturbations.
Table 5: Different transient responses of the MAPK network
The SOS DNA repair network has only one SISO system,
besides those due to differential stability, that exhibited BRAF to MAP2K2 BRAF to MAPK1

Wild type RiseTime: 31.5 RiseTime: 0.26

Table 4: Different transient responses of the SOS network SettlingTime: 56.5 SettlingTime: 76.1
Overshoot: 0 Overshoot: 440.4
recA to uvrD R73A Mutant RiseTime: 2.6 RiseTime: 8.2
SettlingTime: 5.8 SettlingTime: 17.0
Low Radiation Dosage RiseTime: 58.1993 Overshoot: 0 Overshoot: 0
SettlingTime: 108.1459 R80A Mutant RiseTime: 11.7 RiseTime: 48.5
High Radiation Dosage RiseTime: 17.9888 SettlingTime: 21.7 SettlingTime: 90.2
SettlingTime: 35.8853 Overshoot: 0 Overshoot: 0

Page 11 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

Discussion closed-loop system was stable for all possible gains, which
Discovering differentially expressed genes and clustering makes this SISO system robust in normal tissues. But
co-expressed genes into functional groups have given when exposed to toxin, not only was it unstable in itself,
researchers hints about the role of genes in pathogenesis. but no gain value could make the closed-loop system sta-
However, with increasing recognition that cell functions ble, which makes the system brittle. Systems that change
are temporal processes and that the dynamics of gene from high relative stability to low relative stability can be
expression levels and gene interactions play a vital role in considered for association with diseases, because they
determining the health of the organism [1,29], there is a impact the outer loop systems and could make the overall
need to distinguish peculiar dynamical behaviors that system unstable. However, relative stability is not the only
result in sickness from those that do not. Dynamical prop- thing root-locus plots can show. In the recA to uvrA SISO
erties succinctly characterize dynamical behaviors, and system of the SOS network, positive feedback enabled a
differential dynamical properties of gene networks can be lot of stabilizing gains for the SISO system exposed to low
seen as a natural extension of differentially expressed level of radiation, as opposed to the same system exposed
genes. to high level of radiation, which needed negative feedback
for large margins of stabilizing gains. This may portend
In this report we analyzed the dynamical properties of drastic changes in the outer loops, as changing from pro-
genetic networks, such as their stability, their closed-loop motion to inhibition is not easy for biochemical reac-
stability embodied in the root-locus plot, their step tions, and it could be a major sign that this system is
responses, and their controllability. First, we estimated associated with unhealthy conditions.
the state-space models of three genetic networks: the SOS
DNA repair network, the GSH redox cycle system, and the The last dynamical property we looked at was controlla-
MAPK network; then we performed analysis on the esti- bility. Therapeutic treatments can be seen as pushing gene
mated models. From the preliminary results, we found expression levels from unhealthy states to healthy states,
that significant differences in dynamical properties exist and controllability is a theoretical guarantee that there are
in all three networks. possible inputs that can achieve healthy states. Although
we found all systems to be controllable, we did find differ-
All three genetic networks exhibited differential stability. ent degrees of controllability. The condition number of
Stability is a fundamental dynamical property in any the controllability matrix was taken as a measure of degree
dynamical system. A dynamical system is unstable if it of controllability and the redox cycle system in mice lung
diverges or oscillates after being subjected to perturba- exhibited over 100 times difference in its condition
tions. An unstable system is sensitive to perturbation or number, suggesting a much higher possibility that unac-
noise, and it will have erratic behaviors, possibly causing ceptably large inputs are required to move the system into
irreparable cell damage, leading to impairment of cell desired states.
functions and maybe even cell death. A stable genetic net-
work on the other hand confers a degree of robustness Of course much work remains. So far in this report we
against noise on the overall organism. Recently Hornstein have only analyzed a small number of dynamical proper-
and Shomron [30] proposed that miRNAs play a stabiliz- ties while many more remain. Robustness is an important
ing role for a number of genetic networks and the stability property that some consider an organizing principle of
was necessary for the proper functioning of organisms. It complex biological system [31,32], yet we have not inves-
would be interesting to see whether restoring stability to tigated it. There is also the issue of the robustness of esti-
an organism's genetic networks restores the organism's mation. Due to inherent noise in measurements, there are
health. inevitable uncertainties in any parameter estimation. In
general, increasing the sample size will increase the relia-
Besides stability, we also studied relative stability. Root- bility of the results for dynamic properties. Another way to
locus plots track the stability of the closed-loop system deal with this is to obtain confidence intervals for esti-
under the influence of a pure gain controller for single- mated parameter values. However, confidence intervals
input-single-output (SISO) systems, and they can be seen on individual parameters do not directly translate into
as a measure of the relative stability of the SISO system. As confidence intervals for dynamical properties, especially
biological systems are often under control of other, bigger because we have imposed constraints on the parameter
systems, wide margins of stabilizing gains give more lee- space. This should be a topic for further study.
way to, and can sustain some stress from, the controlling
systems, and therefore they are more relatively stable than On the issue of scalability, it is known that the number of
those with narrow margins. The redox cycle system in floating point operations roughly grows to the cubic
mice lungs is the clearest example. Exposed to normal air, power of the number of states [23], assuming that the
the ALDH2A1 to IDH2 system was itself stable and the number of states is larger than either the number of inputs

Page 12 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

or that of outputs. We have implemented our method in Yoshizuka et al. [27] observed the effect of viral protein R
Matlab and for the systems studied in this report compu- (Vpr) on cell cycles. They transfected plasmids that
tation time is around ten minutes on a 1.6 GHz Core Duo expressed wild type Vpr and mutated Vpr (R73A and
laptop, so we expect our program to have no difficulty R80A) into mammalian cells. The microarrays (Hs
with a network with dozens of genes. For large systems, we Operon V2) containing 22,434 oligonucleotide (60- to
should investigate hierarchical system identification 70-mer) spots on a glass slide were used to generate the
method [33]. data. There were three replications for each time point.

Conclusion Our analysis in this paper was done exclusively on the

Dynamical properties are considered to be pivotal in three data sets above.
determining cellular functions such as apoptosis, cell divi-
sion, proliferation, etc. [34], and it follows that differen- Transfer functions and dynamical properties
tial dynamical properties can serve as important A transfer function is a Laplace transform of a linear ordi-
indicators for discovering the role of specific biological nary differential equation of constant coefficients with
processes in causing the malfunction of cells. Only by zero initial conditions. A single transfer function repre-
comparing fundamentally different dynamical behav- sents a single-input-single-output (SISO) system and one
iours between normal and abnormal cells can we begin to can obtain a series of transfer functions from a state-space
untangle the complex interactions and roles of genes in representation of a dynamical system and vice versa [36].
pathogenesis. This will not only add to our understanding The zeroes are roots of the numerator. The characteristic
of diseases but could also be a step toward effective treat- equation of the transfer function is the denominator
ments. equal to zero, and it determines a lot of the dynamical
properties of the system. In particular, the roots of the
Methods characteristic equation are the poles of the system, which
Data sources determine the stability of the system and have great influ-
To test our method on real-world data, we obtained three ence over other dynamical properties.
data sets: E. coli under radiation, mice lung cells exposed
to the normal air and a toxin, and mammalian cell lines Stability analysis
under the influence of various types of Vpr. They were For discrete linear time-invariant systems, the system is
chosen because they all have time course data of organ- stable (its steady states do not diverge) if and only if all of
isms reacting to different perturbations and therefore the eigenvalues of the state transition matrix or all of the
could embody differential dynamical properties. poles of all the transfer functions have magnitude less
than 1 [7]. For continuous systems the requirement is that
Ronen et al. [35] irradiated E. coli and used green fluores- all eigenvalues or poles have negative real part. The sim-
cent protein (GFP) to obtain the rate of transcription of plicity of determining stability belies its importance, for it
various genes in the SOS network. Ronen et al. tracked 8 is one of the most important, best analyzed, and best
genes of the SOS network as they reacted to different irra- known dynamical property. Feedback control's first task is
diation levels, 5 Jm-2and 20 Jm-2. Each level had two sam- to ensure stability and robust control spends a great deal
ples and each sample had 50 time points. They monitored of efforts to ensure stability for uncertain models [37,38].
eight genes: uvrD, lexA, umuD, recA, uvrA, uvrY, ruvA,
and polB. They performed extensive data pre-processing Root-locus plots
on the raw data using hybrid Gaussian median filter and The root-locus method graphically illustrates how the
polynomial fit for smoothing. They also assumed that the poles of the closed-loop system change as the gain of a
rate of accumulation of GFP was proportional to tran- pure gain controller is varied. Later it is generalized to
script production, so we shall make the same assumption. show how the roots change as any parameter of the char-
acteristic equation varies. The parameters must be in the
Sciuto et al. [26] measured the effects of carbonyl chloride form of 1 + KG(s) = 0 where K is the gain (or the parame-
(phosgene) on mice lung. They exposed the mice to either ter), G(s) is a transfer function, and s is a complex variable.
normal air or phosgene for 20 minutes at a concentration The gain is required to be non-negative but this is not a
of 32 42 mg/m3 and sacrificed some of the mice at each problem because we could just make-G(s) the new nomi-
time point. Each time point had 3 samples for air or phos- nal system. We only need two criteria to determine the tra-
gene and two replications. All experimental data were col- jectory
lected using Affymetrix Mouse 430A oligonucleotide
arrays. The raw data were transformed by adding a con- |KG(s)| = 1 (9)
stant first, and then they performed a log transformation.
KG(s) = 180 + k360 (10)

Page 13 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

where k is some integer. The root-locus plot lies in the root-locus technique is one way to use feedbacks to design
complex plane. The path of roots starts at the open-loop a closed-loop system with better rise time, better settling
poles and ends at the open-loop zeros, and if part of the time, and better overshoot.
path lies on the real axis, then it lies to the left of an odd
number of poles [36]. Parameter estimation
We used expectation-maximization (EM) to estimate the
Controllability parameters of genetic networks. EM is a well known and
Controllability is a concept central in systems theory. It is well studied method [23,39], and its application to esti-
about the ability of a system to move from any initial state mating parameters of linear dynamical system has
received attention recently from Gibson et al. [23], whose
to any final state with final control in finite time. The con-
notations we shall follow. We first rewrite
trollability matrix is defined as H = [B AB A2B] for a lin-
ear time invariant system (LTI) of x = Ax + Bu where u is
x t +1 = Ax t + Bu t + w
an m 1 vector, x an n 1 vector, A an n n matrix, and (11)
y t = Cx t + Du t + v
B an n m matrix. If the controllability matrix has full
rank, then LTI is controllable; otherwise it is uncontrolla- as
ble. Another way of saying that a matrix is not full rank is
that it is singular, and due to numerical inaccuracy of dig- x t +1 A B x t w
ital computers and model uncertainty, condition number y = C D u + v , (12)
t t
is used to measure how close to singularity a matrix is. The
and make the following definition for the sake of conven-
condition number of a matrix is defined to be ||H||||H-
1|| where |||| is any matrix norm. We used 2-norm in this
ience:

report. The condition number of the controllability matrix

xt x t +1 A B Q 0
can be seen as a measure of the degree of controllability. z t = xt = , = = .
The larger the condition number is, the greater the inputs ut yt C D 0 R
are needed to reach a target state, even though reaching (13)
nearby states requires no great efforts. Equation (12) then becomes

Unit-step signal and step-response plots

w w 0
A unit-step signal is a constant signal of strength one. The xt = z t + , , . (14)
step response is the output of a dynamical system in v
v 0
response to a unit-step input. The step-response plot We shall also denote all the observations (or outputs) as
graphically gives much information about the dynamical Y, all the inputs as U, and all the states as X. We will add
properties of a system. The most important property the appropriate subscripts if we mean they are up to a certain
step response manifests is stability. A stable system's plot time step for a particular time course expressed as a super-
will converge to a steady state while an unstable system script. The E-step needs to estimate the conditional expec-
will diverge or oscillate. Step-response plots also show set- tation
tling time, rise time, and percent overshoot. Settling time
measures how fast the system achieves the steady state and Q(, ') = E'[log P(X, Y|U)|Y, U] (15)
rise time how quickly the system responds to perturba-
tion. Rise time is defined to be the time for the output to
where is a vector of model parameters and ' is the cur-
go from 10% to 90% of the steady state. Settling time is
rent estimate of the parameters. First, the likelihood func-
defined to be the time for the output to reach and stay
tion of the nth time series, whose length is n, is
within a 2% neighborhood of the steady-state value. Per-
cent overshoot or undershoot is the percentage of the
tn

P (x
maximum or minimum minus the steady state and
divided by the steady state. Rise time is generally associ- Pq ( Ytn , X tn +1 | Utn ) = Pq ( x1) q t +1 , y t | x t , u t ),
ated with the speed of the dynamics, that is, how fast the t =1

system responds to inputs, while overshoot and settling (16)

time measure how close the transient responses stay
within the vicinity of the steady states. They are also
inversely related in nature, that is, both rise time and set-
tling time cannot be kept small: decrease in one necessi-
tates increase in the other if nothing else changes. The

Page 14 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

where Pq ( x1) ~ (m ,1) and where x t|t , Pt/t, Pt|t-1 are calculated from the Kalman filter
x t +1
Pq x t , u t ~ (z t , ) . Expanding equation
yt Pt|t 1 = APt 1|t 1 A T + Q

(16) and taking its logarithm gives G t = Pt|t 1C T (CPt t 1C T + R) 1
Pt|t = Pt|t 1 G t CPt|t 1
2 log Pq ( Ytn , X tn +1 | Utn ) = log |1| +( x1 m) T 11( x1 m)
x t|t 1 = Ax t 1|t 1 + BU t 1
tn
+tn log | | + (x z )
t =1
t t
T 1
(xt z t ).
x t|t = x t|t 1 + G t (y t Cx t|t 1 DU t ), t = 1,..., tn
Mtn |tn = (I Gtn C)APtn 1|tn 1.
(17) (21)
Then we define the following notations for N time series: That constituted the E-step. The M-step, due to the con-
straints imposed by the network structure, is
N N n

1
E (x ) [new - ] M = 0
T
= tn , = n n
| Ytnn , U tnn }, (22)
q {xt t

n =1 n =1 t =1
N n new = { - newT - newT + newnewT} I
1
E ( )
T
n
= q {xt z tn | Ytnn , U tnn }, (23)

n =1 t =1
N n where new and new are the updated parameters, and M is
1
( )
T
= E q {z tn z tn | Ytnn , U tnn }, a constraint matrix of the same size as , so that if an entry
of is constrained then it is zero and otherwise one. We
n =1 t =1

(18) also assume all noise covariance matrices are diagonal.

where n is the number of time points in the nth time Higher order dynamics
series. Taking expectation of equation (17) gives If we stick with one gene for one state, then the system will
only have first order dynamics, which are either exponen-
2 Q(q , q ) = log |1| + trace{11 E q {( x1 m) T ( x1 m)} + log | |
tial decay or exponential growth, associated with all the
genes, but because oscillation is widely observed in biol-
{
+ trace 1 T T + T .
} ogy, at least second order dynamics should be available to
(19) models of genetic networks. We will give a simple deriva-
We need the following quantities for equation (19): tion of how to add second order dynamics for the individ-
ual nodes of the genetic networks using the principle of
continuous to discrete conversion. This is similar to
E q {y t x tT | Ytn , Utn } = y t x tT|tn d'Alch-Buc's method [28]. The third or higher order
dynamics can be similarly added but we do not make use
E q {x t x tT | Ytn , Utn } = x t|tn x tT|tn + Pt|tn
of it in this report.
E q {x t x tT1 | Ytn , Utn } = x t|tn x tT1|tn + Mt|tn ,
We shall focus on one node in genetic networks but the
where x t|tn = E x t | Ytn , Utn , results are easily extrapolated to the entire network. Sup-
pose we have a second order linear differential equation
Pt|tn = var x t | Ytn , Utn , and describing the dynamics of a node:

Mt|tn = cov x t , x t 1 | Ytn , Utn ; and they are obtained

from the Kalman smoother
x + l1x + l2 x = w z ,
j
j j (24)

where x is the node we are interested in, zj is the jth nodes'

J t = Pt|t A T Pt+11|t expression levels, wj its corresponding weights, and 1 and
x t|tn = x t|t + J t [x t +1|tn Ax t|t Bu t R 1y t ]
2 parameters. Let
(20)
Pt|tn = Pt|t + J t [Pt +1|N Pt +1|t ]J tt n x1 = x , x 2 = x. (25)
Mt|tn = Pt|t J tT1 + J t [Mt +1|tn APt t ]J tT1
Then we get

Page 15 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

6. de Lorenzo V, Serrano L, Valencia A: Synthetic biology: chal-

lenges ahead. Bioinformatics (Oxford, England) 2006, 22(2):127-128.
x2 7. Kailath T: Linear Systems. Englewood Cliffs , Prentice-Hall, Inc.;
x1 0 1 x1 0 z1
x
2
=
w j z j l1x 2 l2 x1 = l
2
+
l1 x 2 w1 .
8.
1980.
Camas FM, Blazquez J, Poyatos JF: Autogenous and nonautoge-
j
nous control of response in a genetic network. Proceedings of
(26) the National Academy of Sciences of the United States of America 2006,
103(34):12718-12723.
9. Akutsu T, Miyano S, Kuhara S: Identification of genetic networks
If the steps are uniform, then we can represent the deriva- from a small number of gene expression patterns under the
tives of x as Boolean network model. Pacific Symposium on Biocomputing
1999:17-28.
10. Martin S, Zhang Z, Martino A, Faulon JL: Boolean Dynamics of
dx x Genetic Regulatory Networks Inferred from Microarray
, which becomes x = x(k + 1) x(k), Time Series Data. Bioinformatics (Oxford, England) 2007,
dt t 23(7):86-874.
(27) 11. Szallasi Z, Liang S: Modeling the normal and neoplastic cell
cycle with "realistic Boolean genetic networks": their appli-
cation for understanding carcinogenesis and assessing thera-
where k is the time step. The equation (26) then becomes peutic strategies. Pacific Symposium on Biocomputing 1998:66-76.
12. Kauffman S, Peterson C, Samuelsson B, Troein C: Random Boolean
network models and the yeast transcriptional network. Pro-
x1(k) + x 2(k) ceedings of the National Academy of Sciences of the United States of Amer-
x1(k + 1) ica 2003, 100(25):14796-14799.

2
=
x (k + 1)

w j z j (k) (l1 1)x 2(k) l2 x1(k)

13. Chen T, He HL, Church GM: Modeling gene expression with dif-
ferential equations. Pacific Symposium on Biocomputing 1999:29-40.
j 14. Kikuchi S, Tominaga D, Arita M, Takahashi K, Tomita M: Dynamic
modeling of genetic networks using genetic algorithm and S-
1 1 x1 0 z1
= + .
system. Bioinformatics (Oxford, England) 2003, 19(5):643-650.
l2 1 l1 x 2 w1
15. Kimura S, Ide K, Kashihara A, Kano M, Hatakeyama M, Masui R, Nak-
agawa N, Yokoyama S, Kuramitsu S, Konagaya A: Inference of S-
(28) system models of genetic networks using a cooperative
coevolutionary algorithm. Bioinformatics (Oxford, England) 2005,
The ones and zeros in equation (28) are fixed except in 1 21(7):1154-1163.
- 1 where the whole term is variable. An interesting obser- 16. Dojer N, Gambin A, Mizera A, Wilczynski B, Tiuryn J: Applying
dynamic Bayesian networks to perturbed gene expression
vation is that all interactions and inputs should be on the data. BMC bioinformatics 2006, 7:249.
second order term x2. 17. Husmeier D: Reverse engineering of genetic networks with
Bayesian networks. Biochemical Society transactions 2003, 31(Pt
6):1516-1518.
Authors' contributions 18. Wu FX, Zhang WJ, Kusalik AJ: Modeling gene expression from
Except that HX designed the algorithm and programmed microarray expression data with state-space equations.
Pacific Symposium on Biocomputing 2004:581-592.
the Matlab code, HX and YC both contributed to the man- 19. Rangel C, Angus J, Ghahramani Z, Lioumi M, Sotheran E, Gaiba A,
uscript equally. Wild DL, Falciani F: Modeling T-cell activation using gene
expression profiling and state-space models. Bioinformatics
(Oxford, England) 2004, 20(9):1361-1372.
Acknowledgements 20. Rangel C, Angus J, Ghahramani Z, Wild DL: Modeling genetic reg-
We wish to thank Sciuto AM, Pilips CS, Orzokek LD, Hege AI, Moran TS, ulatory networks using gene expression profiling and state
Dillman JF, Yoshizuka N, Yoshizuka-Chadani Y, Krishnan V, Zeichner SL, space models. In Applications of Probabilistic Modelling in Medical
Informatics and Bioinformatics Edited by: Husmeier D, Robert S,
Rosenfild N, and Alon U for making their data available and Gibson S and Dybowski R. Springer-Verlag; 2004:269-293.
Ninness B for making available the Kalman smoother code. 21. Yamaguchi R, Higuchi T: State-space aproach with the maxi-
mum likelihood principle to identify the system generating
This research was funded in part by Texas A&M University and the Texas time-course gene expression data of yeast. International Journal
of Data Mining and Bioinformatics 2006, 1(1):77-87.
Engineering Experiment Station. 22. Li Z, Shaw SM, Yedwabnick MJ, Chan C: Using a state-space
model with hidden variables to infer transcription factor
References activities. Bioinformatics (Oxford, England) 2006, 22(6):747-754.
1. Wolkenhauer O, Mesarovic M: Feedback dynamics and cell func- 23. Gibson S, Ninness B: Robust Maximum-Likelihood Estimation
tion: Why systems biology is called Systems Biology. Molecu- of Multivariable Dynamic Systems. Automatica 2005,
lar bioSystems 2005, 1(1):14-16. 41(10):1667-1682.
2. Herrero J, Vaquerizas JM, Al-Shahrour F, Conde L, Mateos A, Diaz- 24. Friedman N, Vardi S, Ronen M, Alon U, Stavans J: Precise temporal
Uriarte JS, Dopazo J: New challenges in gene expression data modulation in the response of the SOS DNA repair network
analysis and the extended GEPAS. Nucleic acids research 2004, in individual bacteria. PLoS biology 2005, 3(7):e238.
32(Web Server issue):W485-91. 25. B.M. Broom TJMD D. Subramanian: Predicting altered pathways
3. Rosenfeld N, Alon U: Response delays and the structure of using extendable scaffolds. International Journal of Bioinformatics
transcription networks. Journal of molecular biology 2003, Research and Applications 2006, 2(1):3-18.
329(4):645-654. 26. Sciuto AM, Phillips CS, Orzolek LD, Hege AI, Moran TS, Dillman JF
4. Rosenfeld N, Elowitz MB, Alon U: Negative autoregulation 3rd: Genomic analysis of murine pulmonary tissue following
speeds the response times of transcription networks. Journal carbonyl chloride inhalation. Chemical research in toxicology 2005,
of molecular biology 2002, 323(5):785-793. 18(11):1654-1660.
5. Prill RJ, Iglesias PA, Levchenko A: Dynamic properties of network 27. Yoshizuka N, Yoshizuka-Chadani Y, Krishnan V, Zeichner SL:
motifs contribute to biological network organization. PLoS Human immunodeficiency virus type 1 Vpr-dependent cell
biology 2005, 3(11):e343. cycle arrest through a mitogen-activated protein kinase sig-

Page 16 of 17
(page number not for citation purposes)
BMC Systems Biology 2008, 2:9 http://www.biomedcentral.com/1752-0509/2/9

nal transduction pathway. Journal of virology 2005,

79(17):11366-11381.
28. d'Alch-Buc F, Lahaye PJ, Perrin BE, Ralaviola L, Vujasinovic T, Mazu-
rie A, Bottani S: A Dynamic Model of Gene Regulatory Net-
works Based on Inertia Principle: Berlin. Edited by: Seiffert U,
Jain LC, Schweizer P. Springer; 2005:93-118.
29. Kitano H: Opinon - Cancer as a robust system: implications
for anticancer therapy. Nature Reviews Cancer 2004,
4(3):227-235.
30. Hornstein E, Shomron N: Canalization of development by
microRNAs. Nature genetics 2006, 38 Suppl:S20-4.
31. Kitano H: Biological robustness. Nature Reviews Genetics 2004,
5(11):826-837.
32. Kitano H: Biological robustness in complex host-pathogen
systems. Prog Drug Res 2007, 64:239, 241-263.
33. Ding F, Chen T: Hierarchical least squares identification meth-
ods for multivariable systems. IEEE Transactions on Automatic
Control 2005, 50(3):397-402.
34. Wolkenhauer O: Defining systems biology: an engineering per-
spective. IET systems biology 2007, 1(4):204-206.
35. Ronen M, Rosenberg R, Shraiman BI, Alon U: Assigning numbers
to the arrows: parameterizing a gene regulation network by
using accurate expression kinetics. Proceedings of the National
Academy of Sciences of the United States of America 2002,
99(16):10555-10560.
36. Dorf RC, Bishop RH: Modern control systems. 10th edition.
Upper Saddle River, NJ , Pearson Prentice Hall; 2005:xxv, 881 p..
37. Zhou K, Doyle JC: Essentials of robust control. Upper Saddle
River, N.J. , Prentice Hall; 1998:xvii, 411 p..
38. Bhattacharyya SP, Keel LH: Control of uncertain dynamic sys-
tems : a collection of papers presented at the International
Workshop on Robust Control, San Antonio, Texas, March
1991. Boca Raton , CRC Press; 1991:524 p..
39. Kailath T, Sayed AH, Hassibi B: Linear Estimation. In Prentice Hall
Information and System Sciences Series Prentice Hall; 2000.

Publish with Bio Med Central and every

scientist can read your work free of charge
"BioMed Central will be the most significant development for
disseminating the results of biomedical researc h in our lifetime."
Sir Paul Nurse, Cancer Research UK

Your research papers will be:

available free of charge to the entire biomedical community
peer reviewed and published immediately upon acceptance
cited in PubMed and archived on PubMed Central
yours you keep the copyright

Submit your manuscript here: BioMedcentral

http://www.biomedcentral.com/info/publishing_adv.asp

Page 17 of 17
(page number not for citation purposes)