- Research article
- Open Access
Dynamical pathway analysis
- Hao Xiong†^{1}Email author and
- Yoonsuck Choe†^{1}
https://doi.org/10.1186/1752-0509-2-9
© Xiong and Choe; licensee BioMed Central Ltd. 2008
Received: 28 August 2007
Accepted: 27 January 2008
Published: 27 January 2008
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.
Keywords
Background
Cell functions are complex temporal processes and should be studied as complex dynamical processes rather than only in their individual steady states. It is increasingly recognized that it is the dynamics and the internal structures of the biological systems that give rise to the functioning of cells [1]. Currently, uncovering co-expressed genes and discovering differentially expressed genes are the primary methods for discovering the role of genes in disease pathogenesis [2], but these methods offer only static views and steady-state explanations and thus fail to account for the transient behaviours that influence phenotypes. Genetic regulatory networks seek to model complex interactions and dynamics of gene regulations. Genetic networks should behave differently in sick cells vs. healthy cells because genes that cause diseases behave fundamentally differently, and that difference should be reflected in their dynamical properties. Dynamical properties of genetic networks such as their response time have been studied mostly in the context of network motifs [3, 4], but now we propose that they be investigated for their difference in normal vs. abnormal cells.
In this report we studied four dynamical properties: stability, relative stability, controllability, and transient behaviours (overshoot, settling time, and rise time). Stability governs how a system responds to internal noise and external perturbation and determines whether the system returns to steady states and whether the effect of noise and perturbation diminishes over time. Biologically, an unstable cellular system is very brittle and the slightest disturbance can drive the system beyond tolerance and possibly result in cell death. Prill et al. [5] used stability as a criterion to discern network motifs and their organizing principles, and synthetic biologists are beginning to pay close attention to the stability of their artificial networks [6]. Furthermore, the stability of the system under pure gain feedback control can be analyzed by the root-locus method and the result can be interpreted as a measure of relative stability. In control theory, the root-locus method is a design tool but it is also used as an analytic tool, to see how large a gain can drive the system unstable with feedback loops: the larger margins of stabilizing gains, the better. Related to feedback control, controllability is another pivotal concept in control theory. It and its dual property, observability, were originally conceived as solutions to existence and uniqueness problems of optimal control [7], and the controllability of a dynamical system roughly refers to the ability to move the states of the system around the state space with reasonable efforts. Although controllability is a binary question, there is a measure of the degree of controllability, the idea being that the more controllable a system is the less effort is needed to move the system. Less theoretical than stability and controllability are transient behaviours like settling time and overshoots, which have also received attention from systems 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], differential 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 measurements made by scientists. In the state-space models of genetic networks, states are the regulatory elements, and the inputs and the outputs can be environmental stimuli or expression levels. Because genetic networks have many unknown quantities, state-space models can serve as a good modelling framework.
In this paper, the parameters in state-space models were estimated from the time course of gene expressions using Kalman filter and the constrained expectation-maximization (EM) algorithm (a modified EM algorithm that incorporates prior knowledge about the structure of genetic networks). The regular EM algorithm is commonly used to estimate parameters in the presence of hidden quantities, and they comprise two steps, E-step (expectation) and M-step (maximization), where the E-step estimates the hidden states, and the M-step the parameters [23]. We applied EM algorithm to three sets of time course data and estimated three genetic networks for analysis.
We applied our framework to three real-world time series datasets above and found differential stability, transient responses, and controllability of genetic networks in normal vs. abnormal cells.
Results
Models of genetic networks and their application to real data sets
where x(t) is the state vector, y(t) the output vector, and u(t) the input vector, all at time t; w and v are independent noise terms assumed to be white Gaussian with zero means and covariance Q and R, respectively. Matrix A is called the state transition matrix, B the input matrix, C the output matrix, and D the feed-forward matrix. Matrices A, B, C, D and covariance matrices Q and R together make up the parameters of the dynamical system.
The states represent the biological forces that regulate gene regulation; they describe the behaviours of gene transcription but are hidden. The outputs denote the gene expression levels and are measured, and it is assumed that the expression level of a gene is determined by the state of the regulated gene. The inputs can be any external stimuli that influence gene regulation: substances like drugs, proteins, RNAs, or expression levels of other genes.
Estimated system
Although the number of parameters is small compared with the number of states, which agrees with the knowledge that genetic networks are sparse [28], it is still hard to see at a glance whether they differ in any fundamental way. For that, we must apply systematic analysis to the estimated systems.
Differential stability of systems under different perturbations
Differential stability of the SOS network
Low Dosage | High Dosage |
---|---|
0.8117 | 0.8321 |
0.7367 | 0.6530 |
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 |
Differential stability of GSH redox cycle by
Normal Air | Phosgene |
---|---|
-0.6141 | 2.0830 (unstable) |
0.1177 | -1.0383 (unstable) |
0.4803 + 0.3199i | -0.7561 |
0.4803 - 0.3199i | 1.0512 (unstable) |
0.7467 | 0.9120 |
0.7542 | 0.8696 |
0.4972 + 0.4196i | 1.0470 + 0.2711i |
0.4972 - 0.4196i | 1.0470 - 0.2711i |
Differential stability of the MAPK network
Wild type | R73A | R80A |
---|---|---|
0.8527 | 0.7448 | -1.0169 (unstable) |
0.8862 | -0.0437 + 0.0925i | -0.4078 |
-0.1615 + 0.3646i | -0.0437 - 0.0925i | -0.2023 |
-0.1615 - 0.3646i | -0.6472 | 0.5867 |
-0.4884 | -0.3913 | 0.9534 |
-0.4601 | 0.4680 | 0.7685 |
0.9324 | 0.4877 | -0.6676 |
-0.4477 | -0.4916 | 0.8315 |
We analyzed the SOS DNA network under low and high dosage of radiation and discovered that the network was stable for low dosage and unstable for high dosage. We found that the eigenvalues of SOS network under low dosage of radiation to have the eigenvalues' norm all less than one, and therefore the network was stable. On the other hand, the SOS network was unstable under high dosage of radiation, as the norm of one of its eigenvalues was greater than one.
We also analyzed the redox cycle in mice lung cells that were exposed to either carbonyl chloride (phosgene), an industrial toxin, or normal air; and we found that GSH redox system in normal lung cells was stable – all eigenvalues were within the unit circle, and that the network exposed to toxin was unstable – some eigenvalues were outside the unit circle. Whether the unstable detoxification system contributed to the death of mice exposed to phosgene is not yet known, but Sciuto et al. [26] speculated that the poison might have overwhelmed the detoxification system.
We also analyzed the activity data from the MAPK network in mammalian cells that expressed either wild type Vpr, mutant R73A Vpr, or mutant R80A Vpr; and we found that both the wild type and R73A produced stable behaviours, and R80A caused the network to become unstable. A stable MAPK network helps the virus most, for Yoshizuka et al. [27] found the HIV virus uses MAPK network to cause cell cycle G2 arrest, and over-expression of MAP2K2 reversed the arrest.
Differential relative stability analyzed via root locus
The relative stability of genetic networks is also important; it is a measure of robustness. We studied relative stability by examining the stability margins of pure gain feedback loops through root-locus plots. Given a dynamical system, one forms a feedback loop from the output to the input through only a pure gain controller. Depending on whether the control signal is negated as it is fed into the input, the feedback can be positive (not negated) or negative (negated). The original system is called the open-loop system, and its zeros and poles are the open-loop zeros and poles; the zeros and poles of the overall system are called the closed-loop zeros and poles. A dynamical system's zeroes are the roots of the numerator of the transfer function (for an explanation of the transfer function, see Methods), and the poles are the roots of the denominator. The stability of closed-loop systems depend on the closed-loop poles. The root-locus method generates a plot that traces the closed-loop poles as the gain of the controller is varied, and the portion of gains that make the closed-loop stable is called the stability margin. In the root locus plot, the open-loop zero is represented by a circle (○), the open-loop pole by a cross (×), and if there is a zero-pole cancellation we will see a circle and a cross on top of each other (⊗). The root-locus method can only study systems with single input and single output (SISO), but the dynamical properties of SISO systems is a reflection of the overall system's dynamical properties, so that the performance of the SISO system will manifest itself in the overall system's performance.
Differential degree of controllability
Since one goal of systems biology is to aid the development of therapeutic treatments, which in the context of genetic networks is to bring the network from undesirable states to healthy states by manipulating inputs, the relative ease of moving around in the state space is an important issue. The ability to move a system from one point in the state space to another in finite time with only finite inputs is called controllability, which is a pivotal concept in linear time systems theory [7]. Controllability can be tested by the rank of controllability matrix; if the controllability matrix is of full rank, then the system is controllable, otherwise uncontrollable. Beyond the binary test (controllable or not) there are also degrees of controllability. The condition number of the controllability matrix can be considered as a measure of the degree of controllability, the bigger the number the less the controllability. A system with less controllability may require much greater inputs to achieve the desired final state, which could be a problem as the inputs for biological systems are drugs, radiation therapy, things in limited supply and subject to cost factors. As we will see, different systems could have radically different degrees of controllability.
Although we found all the three genetic regulatory networks controllable under all circumstances, their condition numbers differed, for one significantly. We discovered that the SOS DNA repair system under high dose of radiation had a condition number of 2.8·10^{9} for its controllability matrix, and that under low dosage the condition number was 2.6·10^{9}. The similarly large condition numbers suggest the SOS system under study is difficult to control; whether this is due to radiation is not known. On the other hand, in mice lung exposed to normal air we saw that the redox system had a condition number of 567 for its controllability matrix, and that those exposed to toxin had 70267. The different condition numbers peg the redox system as much more difficult to control after exposure to poison, perhaps due to damages or the fact that the network was being overwhelmed by the effects of the toxin. The third network, the MAPK network in mammalian cell lines, was found to have a condition number of 62.15 when exposed to the wild type Vpr, 88.5 for those exposed to the R80A mutant, and 285.4 for the R73A mutant. It is obvious that R73A mutant results in a stodgier MAPK network than other variants, but overall the MAPK system retains good controllability, making it a good target for treatment.
Differential transient responses
To study cell functions as temporal processes means we must take stock of transient behaviors in addition to steady states. One way to characterize transient behaviors is through the transient response of the dynamical system to inputs like step input and impulse input, but because the step responses and impulses responses give same information for linear systems, we will concentrate on the step input responses. A step input is a constant input, a unit step, a constant unity. The rise time is a measure of the speed of the dynamics, and the settling time and the overshoot gauge how close to the steady state the transient behaviors are. Of course, systems that exhibit differential stability will have different transient responses, but because differential stability is addressed earlier, we will disregard any difference in transient responses due to stability difference.
The transient responses are by their nature studied as input-output pairs, also called a single-input-single-output (SISO) system. Although we will look at individual SISO system extracted from multiple-input-multiple-output (MIMO) systems, the transient responses are still the intrinsic properties of the original system, and differential transient responses suggest fundamentally different dynamical behaviors of the original system in response to external perturbations.
Different transient responses of the SOS network
recA to uvrD | |
---|---|
Low Radiation Dosage | RiseTime: 58.1993 |
SettlingTime: 108.1459 | |
High Radiation Dosage | RiseTime: 17.9888 |
SettlingTime: 35.8853 |
Different transient responses of the MAPK network
BRAF to MAP2K2 | BRAF to MAPK1 | |
---|---|---|
Wild type | RiseTime: 31.5 | RiseTime: 0.26 |
SettlingTime: 56.5 | SettlingTime: 76.1 | |
Overshoot: 0 | Overshoot: 440.4 | |
R73A Mutant | RiseTime: 2.6 | RiseTime: 8.2 |
SettlingTime: 5.8 | SettlingTime: 17.0 | |
Overshoot: 0 | Overshoot: 0 | |
R80A Mutant | RiseTime: 11.7 | RiseTime: 48.5 |
SettlingTime: 21.7 | SettlingTime: 90.2 | |
Overshoot: 0 | Overshoot: 0 |
Although overshoot is generally considered undesirable in engineering (whether fast dynamics or staying close to the steady states are good or bad depends on the circumstances and cannot be determined a priori;) a large overshoot can be a fast way of signalling, or it can be an unbearable disturbance to cells. But being aware of the difference in transient responses is the first step toward devising treatment strategies that shape biological systems' dynamics to our liking.
Discussion
Discovering differentially expressed genes and clustering co-expressed genes into functional groups have given researchers hints about the role of genes in pathogenesis. However, with increasing recognition that cell functions are temporal processes and that the dynamics of gene expression levels and gene interactions play a vital role in determining the health of the organism [1, 29], there is a need to distinguish peculiar dynamical behaviors that result in sickness from those that do not. Dynamical properties succinctly characterize dynamical behaviors, and differential dynamical properties of gene networks can be seen as a natural extension of differentially expressed genes.
In this report we analyzed the dynamical properties of genetic networks, such as their stability, their closed-loop stability embodied in the root-locus plot, their step responses, and their controllability. First, we estimated the state-space models of three genetic networks: the SOS DNA repair network, the GSH redox cycle system, and the MAPK network; then we performed analysis on the estimated models. From the preliminary results, we found that significant differences in dynamical properties exist in all three networks.
All three genetic networks exhibited differential stability. Stability is a fundamental dynamical property in any dynamical system. A dynamical system is unstable if it diverges or oscillates after being subjected to perturbations. An unstable system is sensitive to perturbation or noise, and it will have erratic behaviors, possibly causing irreparable cell damage, leading to impairment of cell functions and maybe even cell death. A stable genetic network on the other hand confers a degree of robustness against noise on the overall organism. Recently Hornstein and Shomron [30] proposed that miRNAs play a stabilizing role for a number of genetic networks and the stability was necessary for the proper functioning of organisms. It would be interesting to see whether restoring stability to an organism's genetic networks restores the organism's health.
Besides stability, we also studied relative stability. Root-locus plots track the stability of the closed-loop system under the influence of a pure gain controller for single-input-single-output (SISO) systems, and they can be seen as a measure of the relative stability of the SISO system. As biological systems are often under control of other, bigger systems, wide margins of stabilizing gains give more leeway to, and can sustain some stress from, the controlling systems, and therefore they are more relatively stable than those with narrow margins. The redox cycle system in mice lungs is the clearest example. Exposed to normal air, the ALDH2A1 to IDH2 system was itself stable and the closed-loop system was stable for all possible gains, which makes this SISO system robust in normal tissues. But when exposed to toxin, not only was it unstable in itself, but no gain value could make the closed-loop system stable, which makes the system brittle. Systems that change from high relative stability to low relative stability can be considered for association with diseases, because they impact the outer loop systems and could make the overall system unstable. However, relative stability is not the only thing root-locus plots can show. In the recA to uvrA SISO system of the SOS network, positive feedback enabled a lot of stabilizing gains for the SISO system exposed to low level of radiation, as opposed to the same system exposed to high level of radiation, which needed negative feedback for large margins of stabilizing gains. This may portend drastic changes in the outer loops, as changing from promotion to inhibition is not easy for biochemical reactions, and it could be a major sign that this system is associated with unhealthy conditions.
The last dynamical property we looked at was controllability. Therapeutic treatments can be seen as pushing gene expression levels from unhealthy states to healthy states, and controllability is a theoretical guarantee that there are possible inputs that can achieve healthy states. Although we found all systems to be controllable, we did find different degrees of controllability. The condition number of the controllability matrix was taken as a measure of degree of controllability and the redox cycle system in mice lung exhibited over 100 times difference in its condition number, suggesting a much higher possibility that unacceptably large inputs are required to move the system into desired states.
Of course much work remains. So far in this report we have only analyzed a small number of dynamical properties while many more remain. Robustness is an important property that some consider an organizing principle of complex biological system [31, 32], yet we have not investigated it. There is also the issue of the robustness of estimation. Due to inherent noise in measurements, there are inevitable uncertainties in any parameter estimation. In general, increasing the sample size will increase the reliability of the results for dynamic properties. Another way to deal with this is to obtain confidence intervals for estimated parameter values. However, confidence intervals on individual parameters do not directly translate into confidence intervals for dynamical properties, especially because we have imposed constraints on the parameter space. This should be a topic for further study.
On the issue of scalability, it is known that the number of floating point operations roughly grows to the cubic power of the number of states [23], assuming that the number of states is larger than either the number of inputs or that of outputs. We have implemented our method in Matlab and for the systems studied in this report computation time is around ten minutes on a 1.6 GHz Core Duo laptop, so we expect our program to have no difficulty with a network with dozens of genes. For large systems, we should investigate hierarchical system identification method [33].
Conclusion
Dynamical properties are considered to be pivotal in determining cellular functions such as apoptosis, cell division, proliferation, etc. [34], and it follows that differential dynamical properties can serve as important indicators for discovering the role of specific biological processes in causing the malfunction of cells. Only by comparing fundamentally different dynamical behaviours between normal and abnormal cells can we begin to untangle the complex interactions and roles of genes in pathogenesis. This will not only add to our understanding of diseases but could also be a step toward effective treatments.
Methods
Data sources
To test our method on real-world data, we obtained three data sets: E. coli under radiation, mice lung cells exposed to the normal air and a toxin, and mammalian cell lines under the influence of various types of Vpr. They were chosen because they all have time course data of organisms reacting to different perturbations and therefore could embody differential dynamical properties.
Ronen et al. [35] irradiated E. coli and used green fluorescent protein (GFP) to obtain the rate of transcription of various genes in the SOS network. Ronen et al. tracked 8 genes of the SOS network as they reacted to different irradiation levels, 5 Jm^{-2}and 20 Jm^{-2}. Each level had two samples and each sample had 50 time points. They monitored eight genes: uvrD, lexA, umuD, recA, uvrA, uvrY, ruvA, and polB. They performed extensive data pre-processing on the raw data using hybrid Gaussian median filter and polynomial fit for smoothing. They also assumed that the rate of accumulation of GFP was proportional to transcript production, so we shall make the same assumption.
Sciuto et al. [26] measured the effects of carbonyl chloride (phosgene) on mice lung. They exposed the mice to either normal air or phosgene for 20 minutes at a concentration of 32 – 42 mg/m^{3} and sacrificed some of the mice at each time point. Each time point had 3 samples for air or phosgene and two replications. All experimental data were collected using Affymetrix Mouse 430A oligonucleotide arrays. The raw data were transformed by adding a constant first, and then they performed a log transformation.
Yoshizuka et al. [27] observed the effect of viral protein R (Vpr) on cell cycles. They transfected plasmids that expressed wild type Vpr and mutated Vpr (R73A and R80A) into mammalian cells. The microarrays (Hs Operon V2) containing 22,434 oligonucleotide (60- to 70-mer) spots on a glass slide were used to generate the data. There were three replications for each time point.
Our analysis in this paper was done exclusively on the three data sets above.
Transfer functions and dynamical properties
A transfer function is a Laplace transform of a linear ordinary differential equation of constant coefficients with zero initial conditions. A single transfer function represents a single-input-single-output (SISO) system and one can obtain a series of transfer functions from a state-space representation of a dynamical system and vice versa [36]. The zeroes are roots of the numerator. The characteristic equation of the transfer function is the denominator equal to zero, and it determines a lot of the dynamical properties of the system. In particular, the roots of the characteristic equation are the poles of the system, which determine the stability of the system and have great influence over other dynamical properties.
Stability analysis
For discrete linear time-invariant systems, the system is stable (its steady states do not diverge) if and only if all of the eigenvalues of the state transition matrix or all of the poles of all the transfer functions have magnitude less than 1 [7]. For continuous systems the requirement is that all eigenvalues or poles have negative real part. The simplicity of determining stability belies its importance, for it is one of the most important, best analyzed, and best known dynamical property. Feedback control's first task is to ensure stability and robust control spends a great deal of efforts to ensure stability for uncertain models [37, 38].
Root-locus plots
The root-locus method graphically illustrates how the poles of the closed-loop system change as the gain of a pure gain controller is varied. Later it is generalized to show how the roots change as any parameter of the characteristic equation varies. The parameters must be in the form of 1 + KG(s) = 0 where K is the gain (or the parameter), G(s) is a transfer function, and s is a complex variable. The gain is required to be non-negative but this is not a problem because we could just make-G(s) the new nominal system. We only need two criteria to determine the trajectory
|KG(s)| = 1 (9)
∠KG(s) = 180° + k 360° (10)
where k is some integer. The root-locus plot lies in the complex plane. The path of roots starts at the open-loop poles and ends at the open-loop zeros, and if part of the path lies on the real axis, then it lies to the left of an odd number of poles [36].
Controllability
Controllability is a concept central in systems theory. It is about the ability of a system to move from any initial state to any final state with final control in finite time. The controllability matrix is defined as H = [B AB A^{2}B⋯] for a linear time invariant system (LTI) of $\dot{x}$ = Ax + Bu where u is an m × 1 vector, x an n × 1 vector, A an n × n matrix, and B an n × m matrix. If the controllability matrix has full rank, then LTI is controllable; otherwise it is uncontrollable. Another way of saying that a matrix is not full rank is that it is singular, and due to numerical inaccuracy of digital computers and model uncertainty, condition number is used to measure how close to singularity a matrix is. The condition number of a matrix is defined to be ||H||·||H^{-1}|| where ||·|| is any matrix norm. We used 2-norm in this report. The condition number of the controllability matrix can be seen as a measure of the degree of controllability. The larger the condition number is, the greater the inputs are needed to reach a target state, even though reaching nearby states requires no great efforts.
Unit-step signal and step-response plots
A unit-step signal is a constant signal of strength one. The step response is the output of a dynamical system in response to a unit-step input. The step-response plot graphically gives much information about the dynamical properties of a system. The most important property the step response manifests is stability. A stable system's plot will converge to a steady state while an unstable system will diverge or oscillate. Step-response plots also show settling time, rise time, and percent overshoot. Settling time measures how fast the system achieves the steady state and rise time how quickly the system responds to perturbation. Rise time is defined to be the time for the output to go from 10% to 90% of the steady state. Settling time is defined to be the time for the output to reach and stay within a 2% neighborhood of the steady-state value. Percent overshoot or undershoot is the percentage of the maximum or minimum minus the steady state and divided by the steady state. Rise time is generally associated with the speed of the dynamics, that is, how fast the system responds to inputs, while overshoot and settling 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 settling time cannot be kept small: decrease in one necessitates increase in the other if nothing else changes. The root-locus technique is one way to use feedbacks to design a closed-loop system with better rise time, better settling time, and better overshoot.
Parameter estimation
We shall also denote all the observations (or outputs) as Y, all the inputs as U, and all the states as X. We will add appropriate subscripts if we mean they are up to a certain time step for a particular time course expressed as a superscript. The E-step needs to estimate the conditional expectation
Q(θ, θ') = E_{ θ' }[log P_{ θ }(X, Y|U)|Y, U] (15)
That constituted the E-step. The M-step, due to the constraints imposed by the network structure, is
[Γ_{new}Σ - Ψ] ∘ M = 0 (22)
Π_{new} = {Φ - ΨΓ_{new}^{ T }- Γ_{new}Ψ^{ T }+ Γ_{new}ΨΓ_{new}T (23)
where Π_{new} and Γ_{new} are the updated parameters, and M is a constraint matrix of the same size as Γ, so that if an entry of Γ is constrained then it is zero and otherwise one. We also assume all noise covariance matrices are diagonal.
Higher order dynamics
If we stick with one gene for one state, then the system will only have first order dynamics, which are either exponential decay or exponential growth, associated with all the genes, but because oscillation is widely observed in biology, at least second order dynamics should be available to models of genetic networks. We will give a simple derivation of how to add second order dynamics for the individual nodes of the genetic networks using the principle of continuous to discrete conversion. This is similar to d'Alché-Buc's method [28]. The third or higher order dynamics can be similarly added but we do not make use of it in this report.
The ones and zeros in equation (28) are fixed except in 1 - λ_{1} where the whole term is variable. An interesting observation is that all interactions and inputs should be on the second order term x_{2}.
Notes
Declarations
Acknowledgements
We wish to thank Sciuto AM, Pilips CS, Orzokek LD, Hege AI, Moran TS, Dillman JF, Yoshizuka N, Yoshizuka-Chadani Y, Krishnan V, Zeichner SL, Rosenfild N, and Alon U for making their data available and Gibson S and Ninness B for making available the Kalman smoother code.
This research was funded in part by Texas A&M University and the Texas Engineering Experiment Station.
Authors’ Affiliations
References
- Wolkenhauer O, Mesarovic M: Feedback dynamics and cell function: Why systems biology is called Systems Biology. Molecular bioSystems. 2005, 1 (1): 14-16. 10.1039/b502088nView ArticlePubMedGoogle Scholar
- Herrero J, Vaquerizas JM, Al-Shahrour F, Conde L, Mateos A, Diaz-Uriarte JS, Dopazo J: New challenges in gene expression data analysis and the extended GEPAS. Nucleic acids research. 2004, 32 (Web Server issue): W485-91. 10.1093/nar/gkh421PubMed CentralView ArticlePubMedGoogle Scholar
- Rosenfeld N, Alon U: Response delays and the structure of transcription networks. Journal of molecular biology. 2003, 329 (4): 645-654. 10.1016/S0022-2836(03)00506-0View ArticlePubMedGoogle Scholar
- Rosenfeld N, Elowitz MB, Alon U: Negative autoregulation speeds the response times of transcription networks. Journal of molecular biology. 2002, 323 (5): 785-793. 10.1016/S0022-2836(02)00994-4View ArticlePubMedGoogle Scholar
- Prill RJ, Iglesias PA, Levchenko A: Dynamic properties of network motifs contribute to biological network organization. PLoS biology. 2005, 3 (11): e343- 10.1371/journal.pbio.0030343PubMed CentralView ArticlePubMedGoogle Scholar
- de Lorenzo V, Serrano L, Valencia A: Synthetic biology: challenges ahead. Bioinformatics (Oxford, England). 2006, 22 (2): 127-128. 10.1093/bioinformatics/btk018View ArticleGoogle Scholar
- Kailath T: Linear Systems. 1980, Englewood Cliffs , Prentice-Hall, Inc.Google Scholar
- Camas FM, Blazquez J, Poyatos JF: Autogenous and nonautogenous control of response in a genetic network. Proceedings of the National Academy of Sciences of the United States of America. 2006, 103 (34): 12718-12723. 10.1073/pnas.0602119103PubMed CentralView ArticlePubMedGoogle Scholar
- Akutsu T, Miyano S, Kuhara S: Identification of genetic networks from a small number of gene expression patterns under the Boolean network model. Pacific Symposium on Biocomputing. 1999, 17-28.Google Scholar
- Martin S, Zhang Z, Martino A, Faulon JL: Boolean Dynamics of Genetic Regulatory Networks Inferred from Microarray Time Series Data. Bioinformatics (Oxford, England). 2007, 23 (7): 86-874. 10.1093/bioinformatics/btm021.View ArticleGoogle Scholar
- Szallasi Z, Liang S: Modeling the normal and neoplastic cell cycle with "realistic Boolean genetic networks": their application for understanding carcinogenesis and assessing therapeutic strategies. Pacific Symposium on Biocomputing. 1998, 66-76.Google Scholar
- Kauffman S, Peterson C, Samuelsson B, Troein C: Random Boolean network models and the yeast transcriptional network. Proceedings of the National Academy of Sciences of the United States of America. 2003, 100 (25): 14796-14799. 10.1073/pnas.2036429100PubMed CentralView ArticlePubMedGoogle Scholar
- Chen T, He HL, Church GM: Modeling gene expression with differential equations. Pacific Symposium on Biocomputing. 1999, 29-40.Google Scholar
- Kikuchi S, Tominaga D, Arita M, Takahashi K, Tomita M: Dynamic modeling of genetic networks using genetic algorithm and S-system. Bioinformatics (Oxford, England). 2003, 19 (5): 643-650. 10.1093/bioinformatics/btg027View ArticleGoogle Scholar
- Kimura S, Ide K, Kashihara A, Kano M, Hatakeyama M, Masui R, Nakagawa N, Yokoyama S, Kuramitsu S, Konagaya A: Inference of S-system models of genetic networks using a cooperative coevolutionary algorithm. Bioinformatics (Oxford, England). 2005, 21 (7): 1154-1163. 10.1093/bioinformatics/bti071View ArticleGoogle Scholar
- Dojer N, Gambin A, Mizera A, Wilczynski B, Tiuryn J: Applying dynamic Bayesian networks to perturbed gene expression data. BMC bioinformatics. 2006, 7: 249- 10.1186/1471-2105-7-249PubMed CentralView ArticlePubMedGoogle Scholar
- Husmeier D: Reverse engineering of genetic networks with Bayesian networks. Biochemical Society transactions. 2003, 31 (Pt 6): 1516-1518.View ArticlePubMedGoogle Scholar
- Wu FX, Zhang WJ, Kusalik AJ: Modeling gene expression from microarray expression data with state-space equations. Pacific Symposium on Biocomputing. 2004, 581-592.Google Scholar
- Rangel C, Angus J, Ghahramani Z, Lioumi M, Sotheran E, Gaiba A, Wild DL, Falciani F: Modeling T-cell activation using gene expression profiling and state-space models. Bioinformatics (Oxford, England). 2004, 20 (9): 1361-1372. 10.1093/bioinformatics/bth093View ArticleGoogle Scholar
- Rangel C, Angus J, Ghahramani Z, Wild DL: Modeling genetic regulatory networks using gene expression profiling and state space models. Applications of Probabilistic Modelling in Medical Informatics and Bioinformatics. Edited by: Husmeier D, Robert S, Dybowski R. 2004, 269-293. Springer-VerlagGoogle Scholar
- Yamaguchi R, Higuchi T: State-space aproach with the maximum likelihood principle to identify the system generating time-course gene expression data of yeast. International Journal of Data Mining and Bioinformatics. 2006, 1 (1): 77-87. 10.1504/IJDMB.2006.009922.View ArticlePubMedGoogle Scholar
- Li Z, Shaw SM, Yedwabnick MJ, Chan C: Using a state-space model with hidden variables to infer transcription factor activities. Bioinformatics (Oxford, England). 2006, 22 (6): 747-754. 10.1093/bioinformatics/btk034View ArticleGoogle Scholar
- Gibson S, Ninness B: Robust Maximum-Likelihood Estimation of Multivariable Dynamic Systems. Automatica. 2005, 41 (10): 1667-1682. 10.1016/j.automatica.2005.05.008.View ArticleGoogle Scholar
- Friedman N, Vardi S, Ronen M, Alon U, Stavans J: Precise temporal modulation in the response of the SOS DNA repair network in individual bacteria. PLoS biology. 2005, 3 (7): e238- 10.1371/journal.pbio.0030238PubMed CentralView ArticlePubMedGoogle Scholar
- B.M. Broom TJMD: Predicting altered pathways using extendable scaffolds. International Journal of Bioinformatics Research and Applications. 2006, 2 (1): 3-18.View ArticleGoogle Scholar
- Sciuto AM, Phillips CS, Orzolek LD, Hege AI, Moran TS, Dillman JF: Genomic analysis of murine pulmonary tissue following carbonyl chloride inhalation. Chemical research in toxicology. 2005, 18 (11): 1654-1660. 10.1021/tx050126fView ArticlePubMedGoogle Scholar
- Yoshizuka N, Yoshizuka-Chadani Y, Krishnan V, Zeichner SL: Human immunodeficiency virus type 1 Vpr-dependent cell cycle arrest through a mitogen-activated protein kinase signal transduction pathway. Journal of virology. 2005, 79 (17): 11366-11381. 10.1128/JVI.79.17.11366-11381.2005PubMed CentralView ArticlePubMedGoogle Scholar
- d'Alché-Buc F, Lahaye PJ, Perrin BE, Ralaviola L, Vujasinovic T, Mazurie A, Bottani S: A Dynamic Model of Gene Regulatory Networks Based on Inertia Principle: Berlin.Edited by: Seiffert U, Jain LC, Schweizer P. 2005, 93-118. SpringerGoogle Scholar
- Kitano H: Opinon - Cancer as a robust system: implications for anticancer therapy. Nature Reviews Cancer. 2004, 4 (3): 227-235. 10.1038/nrc1300View ArticlePubMedGoogle Scholar
- Hornstein E, Shomron N: Canalization of development by microRNAs. Nature genetics. 2006, 38 Suppl: S20-4. 10.1038/ng1803View ArticlePubMedGoogle Scholar
- Kitano H: Biological robustness. Nature Reviews Genetics. 2004, 5 (11): 826-837. 10.1038/nrg1471View ArticlePubMedGoogle Scholar
- Kitano H: Biological robustness in complex host-pathogen systems. Prog Drug Res. 2007, 64: 239, 241-263.Google Scholar
- Ding F, Chen T: Hierarchical least squares identification methods for multivariable systems. IEEE Transactions on Automatic Control. 2005, 50 (3): 397-402. 10.1109/TAC.2005.843856.View ArticleGoogle Scholar
- Wolkenhauer O: Defining systems biology: an engineering perspective. IET systems biology. 2007, 1 (4): 204-206. 10.1049/iet-syb:20079017View ArticlePubMedGoogle Scholar
- 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. 10.1073/pnas.152046799PubMed CentralView ArticlePubMedGoogle Scholar
- Dorf RC, Bishop RH: Modern control systems. 2005, xxv, 881 p.-Upper Saddle River, NJ , Pearson Prentice Hall, 10thGoogle Scholar
- Zhou K, Doyle JC: Essentials of robust control. 1998, xvii, 411 p.-Upper Saddle River, N.J. , Prentice HallGoogle Scholar
- Bhattacharyya SP, Keel LH: Control of uncertain dynamic systems : a collection of papers presented at the International Workshop on Robust Control, San Antonio, Texas, March 1991. 1991, 524 p.-Boca Raton , CRC PressGoogle Scholar
- Kailath T, Sayed AH, Hassibi B: Linear Estimation. Prentice Hall Information and System Sciences Series. 2000, Prentice HallGoogle Scholar
Copyright
This article is published under license to 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.