| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Estimating parameters and hidden variables in non-linear state-space models based on ODEs for biological networks inference
IBISC FRE CNRS 2873, University of Evry and Genopole® 523, place des terrasses 91025 Evry, France
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
Motivation: Statistical inference of biological networks such as gene regulatory networks, signaling pathways and metabolic networks can contribute to build a picture of complex interactions that take place in the cell. However, biological systems considered as dynamical, non-linear and generally partially observed processes may be difficult to estimate even if the structure of interactions is given.
Results: Using the same approach as Sitz et al. proposed in another context, we derive non-linear state-space models from ODEs describing biological networks. In this framework, we apply Unscented Kalman Filtering (UKF) to the estimation of both parameters and hidden variables of non-linear state-space models. We instantiate the method on a transcriptional regulatory model based on Hill kinetics and a signaling pathway model based on mass action kinetics. We successfully use synthetic data and experimental data to test our approach.
Conclusion: This approach covers a large set of biological networks models and gives rise to simple and fast estimation algorithms. Moreover, the Bayesian tool used here directly provides uncertainty estimates on parameters and hidden states. Let us also emphasize that it can be coupled with structure inference methods used in Graphical Probabilistic Models.
Availability: Matlab code available on demand.
Contact: florence.dalche{at}ibisc.univ-evry.fr
Supplementary information: Supplementary data are available from http://amisbio.ibisc.fr/dm
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
Cellular networks (Elowitz and Leibler, 2000; Klipp and Liebermeister, 2006) implement complex mechanisms that enable the cell to respond through time to input signals. Identifying the structure and the parameters of these networks from experimental data is undoubtedly one of the most important challenges in systems biology. Recently, several directions have been simultaneously but independently explored in reverse engineering of metabolic networks, signaling pathways and transcriptional regulatory. This diversity of approaches is essentially due to the kind of available data.
In the case of metabolic pathways, various modeling frameworks from Flux Balance Analysis (FBA) to Ordinary Differential Equations (ODEs) have been developed. As modelers in this domain often benefit from an important background knowledge, much of the current work is focused on model refinement (Herrgard et al., 2006) using perturbation data. Only few works concern parameter estimation of ODEs with a strong assumption about the network structure.
Regarding signaling pathways, most of the related work in the modeling literature (Klipp and Liebermeister, 2006) also considers that the structure of signaling pathways is known. In this case, the most important task becomes the estimation of (generally) non-linear models based on Hill or mass action kinetics.
On the contrary, in modeling transcriptional regulatory networks, the availability of mRNA concentrations has led researchers to develop algorithms for the estimation of both parameters and structure from prior knowledge and experimental data (Nachman et al., 2004; Perrin et al., 2003; Rangel et al., 2004).
While differences can be stressed between gene regulatory networks, signaling pathways and metabolic networks (Klipp and Liebermeister, 2006), we adopt here a transversal point of view and propose to solve in a unique framework the parameter estimation task when the structure of the network is known. We notice that whichever biological network is under study (metabolic, signaling or regulatory), some of its variables may not be observed, increasing the difficulty of the estimation task. We thus suggest the use of a general framework based on state-space models that accounts not only for non-linear dynamics but also for partially observed systems. Similarly to Sitz et al.'s (2002) work developed in a non-biological context, we derive such models from well-grounded Ordinary Differential Equations used in systems biology. This broadens our ability to cover a large variety of biological systems and establishes a bridge between dynamical graphical models and ODEs used in Systems Biology. In non-linear systems, the statistical learning problem is no longer solved in closed form as opposed to linear systems and this raises computational difficulties. In this work, we have chosen to use Unscented Kalman Filtering to tackle non-linearities. Among other Bayesian approaches, this one is fast and relatively easy to implement. It can be considered as a first step towards more sophisticated approaches such as particle filtering. We illustrate the efficiency of this framework by estimating the parameters and hidden variables of two different systems. The first system is the Repressilator, a synthetic transcriptional regulatory network proposed by Elowitz and Leibler (2000) in order to exhibit how sustained oscillations can be obtained through a simple system of three repressors. We use Michaelis–Menten kinetics with Hill curves to describe it.
The second system under consideration is the JAK-STAT signaling pathway (Swameye et al., 2003; Zi and Klipp, 2006) which takes part in the regulation of cellular responses to cytokines and growth factors. Following the setting introduced by Zi and Klipp (2006), we study the parameter and hidden variables estimation problem when using mass action dynamics.
We show by our derivations that the same state-space model can encompass these two kinds of modeling. Both on artificial and experimental data, the estimation method performs successfully.
|
2 MODELING BIOLOGICAL NETWORKS WITH NON-LINEAR STATE-SPACE MODELS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
Let us consider a biological network composed of px variables evolving with time, denoted
. The vector x(t) is supposed to represent the state of the network at time t, which is observed at N + 1 times t0 = 0 < t1 < ... < tN, The graph structure of the network is known. Thus the functional nature of interactions is contained in a parameter 
and the state evolves in the following way for n 
0
|
| (1) |
through the observation functions On, n = 0, ... , N
|
| (2) |
is a measurement noise chosen as a centered Gaussian noise with covariance Rn. The model defined by Equations (1) and (2) is a state-space model, frequently encountered in engineering science. In general state-space models, x can have a stochastic evolution, so equation (1) may be replaced by the more general one
, with
being a white noise. This assumption also has a biological motivation; for instance, McAdams and Arkin (1999) have shown the intrinsic randomness of gene regulatory networks, where x represents gene expression levels and concentrations of transcription factors in the cell. Our goal is to provide a general learning framework in which parameters and hidden variables can be estimated from a time series y0:N = (y0, ... ,yN).
2.1 Deriving non-linear state-space models from ODEs
Quantitative models of biological networks are usually based on Ordinary Differential Equations (ODE), which means that the state of the networks is supposed to satisfy the following ODE
|
| (3) |
n N), its evolution can be cast into (1) with functions Fn and x(tn) = xn for n 0
|
| (4) |
) = o(·) being a function independent of n and . In the following, we show that this framework encompasses both models of transcriptional regulatory networks and models of signaling pathways.
2.2 Modeling transcriptional regulatory networks with Hill equations
In transcriptional regulatory networks, variables of interest are mRNA and protein concentrations, denoted respectively by ri and pi, i = 1, ... ,d. Let us make the assumption here that one gene can only produce one protein. We consider transcription and translation as dynamical processes, in which the production of mRNAs depends on the concentrations of protein transcription factors (TFs) and the production of proteins depends on the concentrations of mRNAs. Hence, we have x(t) = (r(t)
, p(t)) with r(t) = (r1(t), ... , rd(t)) and p(t) = (p1(t), ... ,pd(t)). Equation (3) can be split into the following equations:
|
| (5) |
|
| (6) |
and
are respectively the degradation rates of mRNA i and protein i. The function gi describes how TFs regulate the transcription of gene i and Equation (6) describes the production and the degradation of protein i as linear functions where ki is the translational constant for gene i (Chen et al., 1999).
Various forms have been proposed to model gi(p), ranging from linear (Chen et al., 1999) to non-linear approaches (de Jong, 2002; Elowitz and Leibler, 2000; Nachman et al., 2004; Smolen et al., 2000). Experimental evidence has suggested that the response of mRNA to TFs concentrations has a Hill curve form (de Jong, 2002; Elowitz and Leibler, 2000). The regulation function of transcription factor pj on its target gene i can be described by
for the activation case and
for the inhibition case.
is the maximum rate of transcription of gene i, kij is the concentration of protein pj at which gene i reaches half of its maximum transcription rate and n is a steepness parameter describing the shape of sigmoid responses. The parameter
for i, j = 1, ... , d is the set of kinetic constants to be estimated. Note that if a gene has several regulators, the regulatory part of the Equation (5) can be extended into a product of functions g+ and g– that expresses the combined effect of regulators. However, we consider here examples where the genes have only one regulator.
Elowitz and Leibler (2000) introduced the Repressilator, a synthetic network based on three transcriptional repressors in order to implement the desired dynamical behavior (sustained oscillations) as illustrated in Figure 1. The system was also built experimentally by genetic engineering with mutated Escherichia coli strains. Despite the simplicity of the transcriptional regulation model, the negative feedback loop leads to oscillating concentrations confirmed by experiments. The kinetics of the system can be described by six coupled ODEs, which exactly fit the framework previously described and can be translated into a discrete-time state-space model of form (1). The hidden variables are the protein concentrations with evolution as in Equation (6), and the observations are the (noisily) observed mRNA concentrations ri, i = 1, ... , 3 which satisfies:
|
|
|
where [i + 1] equals i + 1 modulo 3.
2.3 Modeling signaling pathways: the JAK-STAT example
Signaling pathways usually involve numerous and various intermediate products in a complex sequence of transformations hence it is quite difficult to describe them in a general way. Depending on the types of signals and intermediary components, and the localization of the pathways, there exist several relevant types of ODEs. Consequently, it seems rather difficult to give the same wide picture as for transcriptional regulatory networks, but most of the time we can say that the ODEs system involves non-linear reaction rates derived from mass action law and Michaelis–Menten (or Hill) kinetics. We focus here on the JAK-STAT signaling pathway involved in the cellular response to cytokines and growth factors, which involves Janus kinases (JAKs) and Signal Transducers and Activators of Transcription (STATs), see the graph on the right of Figure 1. This pathway transduces the signal carried by these extracellular polypeptides to the cell nucleus, where activated STAT proteins modify gene expression. In both cases, there may be some difficulties to observe the variables of the pathway, and this gives rise to different observation functions o from gene regulatory networks. This is particularly emphasized in the JAK-STAT pathway, for which it is difficult to discriminate between several intermediates in the pathway. Swameye et al. (2003) have suggested an ODE linking the Erythropoietin receptor (EpoR) to the various forms of the STAT5 protein: dephosphorylated STAT5 monomer (x1) and phosphorylated STAT5 dimer (x2) in the cytoplasm, phosphorylated STAT5 dimer (x3) and STAT5 monomer (x4) in the nucleus. In Swameye et al. (2003), the concentration of EpoR is considered as an exogenous variable of the system. The evolution of this network can be described by the following system of coupled differential equations with an input variable u(·) (EpoR), which is an adaptation of the system proposed in Swameye et al. (2003) by Zi and Klipp (2006):
|
| (7) |
} denotes the indicator function, equal to 0 for t and equal to 1 otherwise. The concentrations and constants ai, i = 1, 3,4 in (7) stand for normalized quantities (Zi, 2006). = (a1, a3, a4) is the parameter to be estimated. As pointed out by Swameye et al., the individual STAT5 population is difficult to access experimentally, and only the following variables could be measured: y1 = (x2 + 2x3), the concentration of phosphorylated STAT5 in the cytoplasm and y2 = (x1 + x2 + 2x3), the total amount of STAT5 in the cytoplasm. Thus, the model and data obtained fit the framework of the state space model described in Section 2.1. |
3 ESTIMATION WITH UNSCENTED KALMAN FILTERING |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
3.1 Bayesian estimation
In a Bayesian framework, parameters as well as hidden states are random variables. The goal of inference is to compute the posterior distribution of parameters and initial state, p(
, x0|y0:N), given a prior distribution 
(,x0). It is then possible to estimate variances for parameters. Moreover, we benefit from the large family of algorithms developed for this framework: the state-space form given by equations (1) and (2) is exploited for deriving a sequential estimation procedure based on filtering. We then use an extension of Kalman filtering that computes an approximation of the posterior probability and which also gives an approximation of the Minimum Mean Squared Error estimator (MMSE). This method makes it possible to deal with hidden variables and to estimate them quite simply. We first recall the general principle of filtering and describe the adaptation of the Kalman filter to the case of non-linear evolution equations, using the Unscented Transformation (UT). We then derive the estimation method for unknown parameters and partially observed states and describe it for the biological models considered in the previous section.
3.2 Filtering
In this section, we remove for sake of clarity parameter in Equations (1) and (2), and we describe only the estimation of hidden states. Filtering is the sequential computation of the posterior (or filtering) probability 
n(x) = p(xn|y0:n) for n = 0, ... , N (Cappé, et al., 2005). Without loss of generality, the complete process x = (xn)0 n N may be a Markov (non-deterministic) chain, with values in
(here
). The computation of the filtering probability consists of the alternate and sequential computation of the prediction probability p(xn|y0:n–1), n0, the so-called prediction step:
|
| (8) |
n(·) (the so-called correction step) by:
|
| (9) |
. Note that at n = 0, the prediction step is replaced by setting p(x0|y0:–1)
(x0), where (x0) is our prior distribution on the initial state.
3.3 Approximate filtering
When X is a Gaussian and linear Markov process (F and o are linear), the prediction-correction algorithm is the well-known Kalman filter that consists of a recursive computation of the mean and covariance of the (Gaussian) distribution n(·). This algorithm is no longer valid when the process X is not Gaussian, nor when the function F is non-linear, but several extensions have been proposed to tackle non-linearity or non-gaussianity: Extended Kalman Filtering (EKF, Wan and van der Merwe, 2001), Kernel Kalman Filtering (KKF, Ralaivola and d'Alché-Buc, 2004), Unscented Kalman Filtering (UKF, Julier et al., 2000; Wan and van der Merwe, 2001) and Particle Filtering (PF, Doucet et al., 2001). Among the previously published methods, we have focused on UKF that is an approximation of the posterior distributions n(·) by Gaussian distributions with mean mn(x) and covariance 
n(x). The UKF has the same computational complexity as EKF but offers a better approximation of the true covariance n(x) and does not require the derivatives of F to be computed. Compared to the other filtering methods, KKF requires us to define a kernel function that may be difficult to choose if one wants to stick to classical equations of kinetics. PF involves several hyperparameters and, unlike UKF, is based on clouds of randomly generated points of important size that induces numerous ODE integrations in prediction-correction steps and leads to a higher computational cost. UKF relies on small deterministic sets of appropriately chosen points used in order to mimic the non-linear evolution of the state variable: the so-called sigma points 
0,n, ... , 2px, n (px is the dimension of x). The key idea in UKF lies in the prediction step, where the ‘unscented transformation’ allows one to compute an approximation of the mean mn + 1|n(x) and covariance n+1|n(x) of the prediction probability. The mean and covariance of the transformed variable F(xn) (when xn has the posterior distribution n) can indeed be approximated simply by using the first empirical moments of transformed sigma points chosen as 0,n = mn(x), i,n = mn(x) + Qi,n and
, where Qn is a square root matrix of (2px + 1/2) n(x). Other interesting choices of sigma points are given in Julier et al. (2000).
Then, the correction step is carried out in a way similar to Kalman filtering using the classical Kalman gain matrix Kn + 1 and the (approximate) covariance n+1|n(y) of the pdf p(yn+1| y0:n) (see Section 1 in the Supplementary Material for a complete description of the algorithm). The sequence of estimates (filtered process)
of the hidden variables is the sequence of means mn(x).
3.4 Bayesian estimation of parameters
We adapt here the previous general setting to the joint estimation of hidden states and parameters. This can be accomplished using the augmented state vector approach, which consists in rewriting the dynamical system (1), (2):
|
| (10) |
|
| (11) |
|
| (12) |
n, xn| y0:n) and we can derive a sequence of (improving) estimates
.
The minimizer of the squared error is approximated by mn(). Nevertheless, in practice mn() can be a spurious minimizer, and it is recommended to perform several sweeps of the algorithm on the data, i.e. the parameters estimated in the previous sweep are chosen as the initial value for the next sweep.
As described for hidden states, at n = 0, the prediction step is replaced by setting p(0, x0 | y0:–1) = (,x0). We propose using the rather non-informative hierarchical prior (,x0) = ()(x0), where
, with
, i.e. all of the components of the vector are independent and Gaussian, with a mean drawn according to a uniform distribution whose support is determined by a hyperparameter 
i computed from the data, and the variance
is a fixed value (depending on the data). However, if a certain constraint concerning the initial value is made available, more informative prior could be used.
The estimation procedure for the Repressilator and the JAK-STAT model is done by considering the stacked state variable
. The state components of the points i,n are computed using by Fn(·, n), i.e. by integration of the ODE described in Section 2.1. For the repressilator, the 
i,n are vectors of dimension 3 corresponding to the gene expression levels. For the JAK-STAT model, the i,n are vectors of dimension 2 computed from i,n by the linear combinations x2 + 2x3 and x1 + x2 + 2x3.
|
4 RESULTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
We first illustrate our approach on artificial data generated from the Repressilator model and second, on experimental data of the JAK-STAT pathway. Simulation studies are useful in providing insights into the strengths and weaknesses of learning algorithms, such as robustness against numerous choices of settings, including the quantity of observed data, the sampling interval for observing the data, the number of time points in the observed time series.
The size of the systems used in the experiments are quite representative of the size of the models that the proposed method (based on UKF) can efficiently handle, i.e. around 10 variables and parameters. In higher dimensions, the recursive optimization using the UKF approximation can lead to spurious minimizers, but such a limitation can be partly overcome by using a better approximation, such as the Particle Filter. The main limitation of the approach depends on the respective sizes of the observed and hidden parts of the system.
4.1 Parameter estimation of the Repressilator
4.1.1 Simulated data
We start from the equations given in (2.2) and fix the following values of the parameters according to the stability study presented in Elowitz and Leibler (2000):
,
,
,
and n = 3. The components of the initial state are drawn independently from a uniform distribution on [0, 100] (arbitrary units). Simulations are performed using the MATLAB numerical integrator ode45 over the time interval [0, T], with T = 20. The observation noises
are added to three observed variables to mimic gene expression data and the SD of
shown in the experiments is chosen to be equal to 20% of the SD of the states. The robustness of the method has been tested with respect to a higher noise level (30%, 40%), and similar results for the estimation for the states and parameters have been obtained. The estimated predicted variance and the variance of the estimators increase, although no systematic divergence of the method has been detected.
During the simulation, measurements are sampled at a fixed interval 
t, so that for each experiment a time series containing T/t time points is collected. We assume that the learning problem consists in identifying the following six parameters:
while the degradation rates for proteins and mRNAs are known. In order to learn the true parameters, we use a multi-start approach by sampling I = 50 different initial states and parameters from our prior (, x0), so that we compute 50 filters in parallel. Our final state and parameter estimates are simply the mean of the prediction of the 50 different filters (an alternative way to combine the different filters would be to select the filter with the lowest prediction error). The Gaussian priors for the parameter are such that
and
, and for the unobserved variables
and
is set to 20% of the SD of the state xi. For the observed variables, the prior is also Gaussian with mean
and the same formula as for the unobserved variables is used for the SD.
4.1.2 Evaluation of parameter estimation
The filtered protein concentrations and parameters using UKF are shown in Figure 2 for t = 0.2, (see also Fig. 1 in Supplementary Material for the case of one experiment: one time series corresponding to one initial condition). The filtered concentrations of proteins quickly adjust to the true protein trajectories. The SDs of the estimators are estimated using the square-root of the diagonal of the matrix n(x). Parameter estimates start far from the true value with high SDs, but they gradually converge to the true parameters. The small values of the final SDs of the estimates point out the convergence of the learning algorithm. Finally, since a single sweep of UKF on a time series containing 100 observations takes less 5 s, our multi – start approach gives an estimation within 
250 s.
|
4.1.3 Dependency between the prediction error of the hidden states and the sampling interval
In order to analyze the dependency between the prediction error and the sampling interval, we used different sampling intervals (respectively)
t = 0.2, 0.4, 0.8, 1 and 2, corresponding to (respectively) 100, 50, 25, 20 and 10 time points. The estimates for all parameters are reported in Table 1 of the Supplementary Material. We plot the estimation results for parameter
in Figure 3. Obviously, the estimates are closer to the true values with smaller SD when there are more time points available. In order to compare errors for the different cases, we introduce the normalized Mean Squared Error between the true and estimated trajectories: |
|
n is the covariance matrix of the UKF estimate at time n and diag(A) is the diagonal matrix equal to the diagonal of a square matrix A. The errors for each component of the state variable are rescaled in order for them to be comparable (though the covariance between the components is not taken into account). The MSE is plotted in Figure 2 of the Supplementary Material and increases when there are less time points. However, a relevant result is still obtained for only 10 observed time points. As we can see in the final column of Table 1 of the Supplementary Material, the true parameters stands within ±
confidence interval of the estimated parameters.
|
versus the length of observed time series. The dash line shows the true parameter. The bars show means and SDs of estimated parameter.4.1.4 Dependency between the prediction error of the parameters and the number of repeated experiments
We show here that the influence of the number of different experiments (i.e. time series corresponding to the observation of the same system but with different initial conditions). The learning algorithm can be adapted to this setting in a straightforward manner, and we show that it is possible to deal with more difficult situations. As an illustration, we assume that the parameters k1, k2, k3 are also unknown, so that we have to estimate nine parameters from several time series (with 50 observations each). The MSE's obtained for a number of experiments varying from 1 to 6 are reported in Figure 4 of the Supplementary material. The MSE decreases when more experiments are available. One notes that three different experiments provides the same mean level of MSE as six experiments but with higher variance. This may help in designing experiments and gives some hints for obtaining accurate estimates of the parameters from only a few experiments. We also plot the estimation of parameter
versus the number of repeated experiments in Figure 3 of the Supplementary Material. The estimated parameter tends to the true value with smaller SDs when the number of experiment increases.
|
4.2 Parameter estimation for the JAK-STAT pathway model using experimental data
Experimental data of JAK-STAT pathways from Swameye et al. (2003) was used. Time series of two observed variables y1 (the total concentration of phosphorylated STAT5) and y2 (the total concentration of STAT5 in the cytoplasm) are measurable. Each time series contains 16 time points sampled at t = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 60] min. Data for the input EPoR phosphorylation is also available. Here, we use a linear interpolation in order to obtain a continuous time input. We initialize the parameters a1, a3, a4 and the initial condition x1 independently with a uniform distribution on [0,5]. We plot the normalized MSE between the predicted time series and the data in Figure 4 of the Supplementary Material. The convergence of this curve shows the stability of the learning algorithm, and ensures that we have reached a local minima. Eventually, the parameter estimates (with SDs) are
,
and
, and the prediction for the observed variables y1 and y2 are shown in Figure 5, which shows a good fit of the learned model. We also check the coherence of the estimation by simulating the JAK-STAT pathway with these estimates. A new time series x
is simulated from (7) with initial conditions x1 = 0.2, x2 = x3 = x4 = 0 and the estimated parameters. The result in Figure 5 showed that the learned model is able to predict well the four unobserved components of x, so we may have a higher confidence in the prediction of the unobserved variables.
|
|
5 DISCUSSION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
The approach we derived for non-linear biological systems has already been proposed in another context by Sitz et al. (2002). However, the present work is a novel application of this approach in Systems Biology, which opens new perspectives in estimating non-linear biological systems. So far, linear state-space models have been mainly used for transcriptional regulatory networks. Perrin et al. (2003) proposed IDBN (Inertial Dynamic Bayesian Network), a second-order linear ODE model that accounts for inertia and allows us to represent damped oscillations: in this model hidden variables are the true mRNA concentrations and their derivative. Parameters and hidden variables have been estimated in this context by linear Kalman filtering and smoothing (d'Alché-Buc et al., 2005). Rangel et al. (2004) introduced also a linear dynamical probabilistic model for which the biological interpretation is less explicit but this work has been interestingly followed by a study on hidden variable estimation in Beal et al. (2005). However, a higher level of detail is often required in order to improve the biological relevance of the models. This is why Nachman et al. (2004) suggested a non-linear state-space model based on Michaelis–Menten equations but only to cope with transcription factors. In their model as well as in Rogers et al. (2007), the equation taking into account protein production and degradation is not used. Our model in the case of regulatory networks can thus be seen as an extension to a full modeling of protein and mRNA concentrations through time.
Another interest in the state-space interpretation of ODE model lies in the access to new estimation methods, which could be faster than classical ones and could also be able to incorporate priors on the system. Indeed, most of the methods proposed so far for signaling pathway consist of the minimization of a least squares criterion (Li et al., 2005; van Riel and Sontag, 2006) without any regularization term. Even when the model is slightly complex (around 10 variables), the minimization step requires special care in order to reach the true maxima, because classical local optimization methods (gradient-based or Newton-like) are doomed to reach local minima. In the end, global optimization techniques have been proposed in order to solve these problems such as simulated annealing, evolutionary algorithms, (Koh et al., 2006; Moles et al., 2003; Polisetty et al., 2006). These techniques are batch methods that do not use the recursive structure induced by the ODE model. In contrast to the former, the alternate regression (Chou et al., 2006) or the system perturbation method (Schmidt et al., 2005) exploit the particular structure of the learning problem in order to derive relatively simple algorithms. Our method takes advantage of the dynamical nature of the model by implementing a recursive optimization. Though the UKF is only able to reach a local minimum of the posterior probability, it delivers a sequence of estimated states
and parameters
which are likely for all the intermediate estimation problems with observations y0:n. This sequence of estimates remains plausible, at least when the initial conditions are correctly chosen (which can be done in some case with the literature). We have used in our experiments simple priors (flat or Gaussian priors), but Bayesian estimation may benefit from more elaborated prior distributions in order to favor meaningful regions of the parameter and state spaces. Moreover, the most striking property of this estimator is its ease of implementation and above all its speed, which is an advantage over global optimization methods. One should also note that the variance of the estimator
is simultaneously computed, whereas it is not straightforward to compute it in batch methods since the model can be too complex to be done analytically or approximately. Moreover, our estimation method can still be enhanced by the use of smoothing probabilities and the promising results of UKF calls for more sophisticated filtering approaches such as particle filtering.
|
6 CONCLUSION AND PERSPECTIVES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
We have presented non-linear state-space models for describing biological networks and non-linear filtering approaches to estimate both parameters and hidden variables. As the models are built up from ODEs, they benefit from all the existing background in biological modeling with ODEs and thus are ensured to exhibit high biological relevance. This point was illustrated on two different kinds of networks models: a transcriptional regulatory model based on Hill kinetics and a signaling pathway model based on mass action kinetics. Let us notice that, given the type of equations we already dealt with, there is no reason not to apply this approach to model metabolic networks and to estimate their parameters as soon as experimental data are available. Moreover, this work raises several issues that encourage further works. First, our estimation algorithm like others in literature requires time series of sufficient length to be efficient. In this case, we need to reduce the complexity of the parameter space by introducing relevant biological priors. The Bayesian framework we use is appropriate for this. Second, large networks with large number of parameters may not be identifiable. In order to overcome this limitation, we suggest applying the decompositional scheme developed in Koh et al. (2006) in order to work only on subnetworks. Third, it should be emphasized that our parameter estimation method can be coupled with any of the classical structure learning schemes used in graphical probabilistic models (MCMC, evolutionary approaches) in order to fully reverse-engineer biological networks. Finally, it should also be stressed that this framework could account for joint modeling of metabolic, signaling and regulatory networks if one can deal with the various time scales and has access to appropriately observed time series.
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
This research has been mainly funded by Genopole® through an ATIGE grant and a postdoc grant. ANR (National Agency for Research) has also contributed to fund the project through the ARA 2005 call for projects.
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
Associate Editor: Chris Stoecker
Received on June 15, 2007; revised on October 6, 2007; accepted on October 7, 2007
|
REFERENCES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 MODELING BIOLOGICAL NETWORKS...3 ESTIMATION WITH UNSCENTED...4 RESULTS5 DISCUSSION6 CONCLUSION AND PERSPECTIVESACKNOWLEDGEMENTSREFERENCES |
|---|
Beal M, et al. A Bayesian approach to reconstructing genetic regulatory networks with hidden factors. Bioinformatics (2005) 21:349–356.
Cappé O, et al. Inference in Hidden Markov Models (2005) Springer-Verlag.
Chen T, et al. Modeling gene expression with differential equations. In: Pacific Symposium of Biocomputing (1999).
Chou I-C, et al. Parameter estimation in biochemical systems models with alternating regression. Theor. Biol. Med. Model (2006) 3:1–15.[CrossRef][Medline]
d'Alché-Buc F, et al. Dynamic model of gene regulatory networks based on inertia principle. Bioinformatics using Computational Intelligence Paradigms, volume 176 of Studies in Fuzziness and Soft Computing. Springer (2005) 93–117.
de Jong H. Modeling and simulation of genetic regulatory systems: a literature review. J. Comput. Biol (2002) 9:67–103.[CrossRef][Web of Science][Medline]
Doucet A, et al. Sequential Monte Carlo Methods in Practice (2001) Springer-Verlag.
Elowitz M, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature (2000) 403:335–338.[CrossRef][Medline]
Herrgard MJ, et al. Identification of genome-scale metabolic network models using experimentally measured flux profiles. PLoS Comput. Biol (2006) 2:676–686.[Web of Science]
Julier S, et al. A new method for the nonlinear transformation of means and covariances in filters and estimators. IEEE Trans. Automat. Contr (2000) 45:477–482.[CrossRef]
Klipp E, Liebermeister W. Mathematical modeling of intracellular signaling pathways. BMC Neurosci (2006) 7(Suppl. 1).
Koh G, et al. A decompositional approach to parameter estimation in pathway modeling: a case study of the Akt and MAPK pathways and their crosstalk. Bioinformatics (2006) 22:271–280.[CrossRef]
Li Z, et al. Parameter estimation of ordinary differential equations. IMA J. Numer. Anal (2005) 25:264–285.
McAdams H, Arkin A. It's a noisy business! Genetic regulation at the nanomolar scale. Trends Genet (1999) 15:65–69.[CrossRef][Web of Science][Medline]
Moles C, et al. Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Res (2003) 13:2467–2474.
Nachman I, et al. Inferring quantitative models of regulatory networks from expression data. Bioinformatics (2004) 20(Suppl. 1):248–256.[CrossRef]
Perrin B-E, et al. Inference of gene regulatory network with Dynamic Bayesian Network. Bioinformatics (2003) 19:i386–i49.
Polisetty PK, et al. Identification of metabolic system parameters using global optimization methods. Theor. Biol. Med. Model (2006) 3:1–15.[CrossRef][Medline]
Ralaivola L, d'Alché-Buc F. Dynamical modeling with kernels for nonlinear time series prediction. In: Advances in Neural Information Processing Systems (NIPS'04) (2004).
Rangel C, et al. Modeling T-cell activation using gene expression profiling and state-space models. Bioinformatics (2004) 20:1361–1372.
Rogers S, et al. Bayesian model-based inference of transcription factor activity. BMC Bioinformatics (special issue from PMSB'06) (2007).
Schmidt H, et al. Identification of small scale biochemical networks based on general type system perturbations. FEBS J (2005) 272:2141–2151.[CrossRef][Medline]
Sitz A, et al. Estimation of parameters and unobserved components for nonlinear systems from noisy time series. Phys. Rev. E (2002) 66:016210.[CrossRef]
Smolen P, et al. Modeling transcriptional control in gene networks - methods, recent results, and future directions. Bull. Math. Biol (2000) 62:247–292.[CrossRef][Web of Science][Medline]
Swameye I, et al. Identification of nucleocytoplasmic cycling as remote sensor in signaling by databased modeling. Proc. Natl. Acad. Sci. USA (2003) 100:1028–1033.
van Riel N, Sontag E. Parameter estimation in models combining signal transduction and metabolic pathways: the dependent input approach. IEE Proc. Syst. Biol (2006) 153:263–274.[Medline]
Wan EA, van der Merwe R. The Unscented Kalman Filter. Series in adaptive and Learning Systems for Signal Processing, Communications, and Control. Haykin S, ed. (2001) Kalman Filtering and Neural Networks. Wiley, pp. 221–283.
Zi Z. User Guide: SBML-PET, A Sytems Biology Markup Language Based Estimation Parameter Tool. Max Planck Institute for Molecular Genetics (2006) Available from www.sysbio.molgen.mpg.de/SBML-PET/.
Zi Z, Klipp E. SBML-PET: a systems biology markup language-based parameter estimation tool. Bioinformatics (2006) 22:2704–2705.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  B. Anchang, M. J. Sadeh, J. Jacob, A. Tresch, M. O. Vlad, P. J. Oefner, and R. Spang Special Feature: Complex Systems: From Chemistry to Systems Biology Special Feature: Modeling the temporal interplay of molecular signaling and gene expression by using dynamic nested effects models PNAS, April 21, 2009; 106(16): 6447 - 6452. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 R. Yoshida, M. Nagasaki, R. Yamaguchi, S. Imoto, S. Miyano, and T. Higuchi Bayesian learning of biological pathways on genomic data assimilation Bioinformatics, November 15, 2008; 24(22): 2592 - 2601. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||