Network inference via adaptive optimal design
- Johannes D Stigter^{1}Email author and
- Jaap Molenaar^{1}
DOI: 10.1186/1756-0500-5-518
© Stigter and Molenaar; licensee BioMed Central Ltd. 2012
Received: 11 April 2012
Accepted: 19 July 2012
Published: 24 September 2012
Abstract
Background
Current research in network reverse engineering for genetic or metabolic networks very often does not include a proper experimental and/or input design. In this paper we address this issue in more detail and suggest a method that includes an iterative design of experiments based, on the most recent data that become available. The presented approach allows a reliable reconstruction of the network and addresses an important issue, i.e., the analysis and the propagation of uncertainties as they exist in both the data and in our own knowledge. These two types of uncertainties have their immediate ramifications for the uncertainties in the parameter estimates and, hence, are taken into account from the very beginning of our experimental design.
Findings
The method is demonstrated for two small networks that include a genetic network for mRNA synthesis and degradation and an oscillatory network describing a molecular network underlying adenosine 3’-5’ cyclic monophosphate (cAMP) as observed in populations of Dyctyostelium cells. In both cases a substantial reduction in parameter uncertainty was observed. Extension to larger scale networks is possible but needs a more rigorous parameter estimation algorithm that includes sparsity as a constraint in the optimization procedure.
Conclusion
We conclude that a careful experiment design very often (but not always) pays off in terms of reliability in the inferred network topology. For large scale networks a better parameter estimation algorithm is required that includes sparsity as an additional constraint. These algorithms are available in the literature and can also be used in an adaptive optimal design setting as demonstrated in this paper.
Findings
Background
Recent research in Systems Biology shows that the concept of ‘network’ turns out to be highly relevant at all levels of biological complexity [1, 2]. In a modeling approach based on networks, the components of the system under consideration are represented as nodes in a graph and the interactions as edges that connect pairs of nodes. The functioning of a network is determined by its dynamics, i.e., the evolution in time of the nodes. The nodes usually represent the concentrations of some species and their time development is regulated via interactions with other species. At the cellular level, regulatory and metabolic networks are out spoken examples, but also signaling pathways make use of the network concept. Typical examples at the highest complexity levels are predator-prey models and ecological food webs which are in use since long. It is remarkable that not only the network concept forms a unifying element across all biological aggregation levels, but also the mathematics needed to describe the dynamics of the network nodes show great similarities [3]. The central role of network modeling in biology implies that network inference is one of the major challenges of modern biology. The ultimate aim of network inference is to deduce the structure of the network as accurately as possible from data obtained in the experimental practice. To obtain the necessary experimental information, one usually perturbs the network in some way, hoping that the induced response is useful to be explored in some network inference approach. However, if this is not done in a structured way, the harvest might be disappointing. In this paper, we consider the question how we could effectively and fast learn as much as possible of the network structure by designing the experimental setup in an optimal way. It should be realized that the question we put here in a network inference context, is in its generality not new, since similar research is at the core of the mathematical subdiscipline called‘ systems theory’, e.g., [4], and more specifically, in the subdiscipline‘ system identification’ [5]. The fascinating aspect is that insights obtained in the engineering practice may also be of great relevance in a life sciences setting. In an engineering context this kind of research is usually referred to as ‘reverse engineering’. In this paper we also aim at connecting the still too much separated realms of scientists in biology and engineering.
Before getting into details, it is important to remark that network inference may aim at different levels of accuracy [6, 7]. The lowest level approach leads to a relatively rough impression of the network structure. Given the nodes, one uses the data to find out whether a coupling exists between any pair of nodes. In mathematical terms, the network is represented as a graph with undirected edges and the inter action is not made explicit, e.g., in terms of an equation. An undirected edge indicates that the dynamics of the two connected nodes are correlated due to some interaction, but the centeracter of this interaction is not pointed out in detail. For example, the first node could affect the second one, but the interaction could also be the other way around. This case is especially relevant if the data contain steady state levels of the nodes obtained after perturbing the system a considerable number of times, e.g., via a knocking out procedure in gene networks. For this case a number of statistical inference methods have been developed that yield undirected graphs as outcome [8–11]. The next level of accuracy is to strive for a directed graph, in which the edges are arrows. The direction of an arrow then indicates a causal relationship [12]. The most sophisticated level is to deduce more detailed information on the centeracter of the interactions. One then tries to answer questions such as: Is it a promoting or an inhibiting interaction (in regulatory networks), is it a linearly increasing, saturating, or decaying interaction (in metabolic networks), etcetera [13, 14]. In this paper we aim at the third level of accuracy by including the estimation of the strengths of interactions. A common procedure to generate data is to perturb the system under consideration in a more or less random way. The change in the observed data is then analyzed to infer information about its structure. In practice this may lead to very poor results, since it is not assured that the chosen perturbation is really useful for inference purposes. In this paper we propose to follow a much more advanced approach, in which the perturbations are designed such that they contain the optimal amount of information, given the available data. This approach also starts with a more or less randomly chosen perturbation, but we show that the information from such a first step can be effectively used to design a second perturbation that yields enhanced insight in the network structure. So, we maximize in advance the information content of the second perturbation, given the insight deduced from the first perturbation. This maximization is based on improving the condition number of the so-called Fischer Information Matrix of the system. If required, this procedure could be repeated to refine the inference results further.
Our research focuses on networks whose dynamics are described in terms of ordinary differential equations (ODE). We assume that it is possible to observe the network while it develops in time, so that the data consists of time series of some or all of the network components. The general ideas of the present approach are widely applicable and not at all restricted to a special type of network or level of complexity. For illustrational purposes we will make use of a gene is taken from [15]. This system has been used earlier in several studies to explain and illustrate the principles of other inference methods in signaling and gene networks [16–18] and turns out to be very useful for this purpose. Next, we apply our method to the Laub-Loomis model [19], which describes oscillations in excitable cells of Dictyostelium, and show that our new procedure is very effective in reconstructing the underlying network, just by slightly perturbing system and observing its recovery to its natural oscillatory behavior.
The paper is organized as follows. In the next section we will first introduce two motivating examples in which reconstruction of a directed graph representing the interaction nodes of the network is the primary goal. Then we will further elaborate on the details of adaptive optimal input design [20] and explain the method. Finally, the methodology will be applied to the two motivating examples explained earlier and some results will be presented and discussed.
Methods
Two motivating examples
In this section we introduce two biological systems that will be used later on to illustrate how our adaptive optimal design approach works in practice. For both models centeracteristic parameter values are taken from the literature. These values are used to generate artificial data and to mimic real practice these data are spoiled by adding noise. The present aim is to show how the parameters can be efficiently estimated from the data in an adaptive manner.
mRNA synthesis and degradation
Parameter values mRNA synthesis and degradation
V_{1s} = 5 | A_{14} = 4 | K_{14a} = 1.6 | K_{12i} = 0.5 | n_{12} = 1 | V_{1d} = 200 | K_{1d} = 30 |
V_{2s} = 3.5 | A_{24} = 4 | K_{24a} = 1.6 | n_{32} = 2 | n_{14} = 2 | V_{2d} = 500 | K_{2d} = 60 |
V_{3s} = 3 | A_{32} = 5 | K_{32a} = 1.5 | K_{31i} = 0.7 | n_{31} = 1 | V_{3d} = 150 | K_{3d} = 10 |
V_{4s} = 4 | A_{43} = 2 | K_{43a} = 0.15 | K_{12i} = 0.5 | n_{43} = 2 | V_{4d} = 500 | K_{4d} = 50 |
Here, x is the state vector x = (x_{1},x_{2},x_{3},x_{4})^{T} , u is the vector containing the inputs, so u = (V_{1s},V_{2s},V_{3s},V_{4s})^{T} ,θ is the vector containing the remaining parameters in Table 1, and f is the vector valued function defined in equations (5)–(12).
where δ x(t) and δ u are now the deviations from the steady state x_{ ss } and the corresponding (steady) input u_{ ss } that maintains this steady state. The matrix A in equation (14) contains the so-called interaction coefficients of the network and these correspond to a directed graph where there is an interaction from node j to node i if the corresponding matrix element a_{ ij }≠ 0. We note that, of course, this linearisation is a local property that depends on our choice of the linearisation point in state space where we derive the linear system (14). However, it is important to realize that, if nodes i and j are not connected, the matrix A has zeros at positions ij and ji, regardless of the linearisation point that is used! The question of optimal experimental design now comes down to choosing the entries of the vector δ u in such a way that matrix pair [A, B] can be estimated in an optimal manner.
Laub-Loomis model–the no-input case
Parameter values Laub–Loomis
k_{1} = 2.0 | k_{2} = 0.9 | k_{3} = 2.5 | k_{4} = 1.5 | k_{5} = 0.6 | k_{6} = 0.8 | k_{7} = 1.0 |
k_{8} = 1.3 | k_{9} = 0.3 | k_{10} = 0.8 | k_{11} = 0.7 | k_{12} = 4.9 | k_{13} = 23 | k_{14} = 4.5 |
If we assume the perturbation δ u to be a Dirac delta function (or better, a δ-distribution), the question of optimal experimental design now becomes one of finding the values of the entries in the vector b that yield an optimal estimate A, i.e., an estimate Â with lowest possible uncertainty bounds.
Parameter estimation– towards adaptive experimental design
Since there are four parameters{V_{ is }, i = 1,…, 4}, there are four ‘inputs’ at our disposal that we can use as perturbation candidates. In the following we assume that four separate perturbations of all inputs are generated with random values between _{ − }50% and 50% of the steady state input values contained in the vector x_{ ss } (that corresponds to the initial steady state u_{ ss }).
Schmidt et al.[16], inspired by the methods and results of [15], applied linear regression to data obtained by virtual experiments using the original system model as a simulation tool. Based on calculated deviations from the steady state x_{ ss } a linear regression formula was obtained with the coefficients of matrices A and B as unknowns. Schmidt successfully computed the interaction coefficients, including the unknown input perturbations, from step-response data. For clarity it is noted that the step-response data represents the dynamic transient behavior from one steady state to another steady state. This is different from the approach presented by Steinke and co-workers [21], where an interesting statistical methodology was applied with the same idea of iterative experimental design in mind, but for the steady states only, i.e. with no dynamic transient behavior included.
It should further be noted that the regression method applied here only works for relatively small scale problems in which the number of states is not too large. This is mainly because of the limitations in the applied parameter estimation method (see below for details) which easily introduces identifiability problems for large scale networks where the number of parameters to be estimated in the Jacobi matrix is simply too large. The disadvantage of ordinary least squares estimation, however, can be remedied for larger scale networks by more advanced parameter estimation algorithms that include sparsity of the network as an additional constraint (see e.g. [22] for a convex linear-programming approach.)
In summary the parameter estimation method proceeds as follows:
Step 1. Collect N observations of the state vector x in the original non-linear model (13) with a sampling interval Δt. Subtract the steady state x_{ ss } from these observations.
Step 2. Calculate the difference between consecutive time instances of the obtained measurements in step (1) and stack these differences in a (n × (N _{−} 1)) matrix M with n the dimension of the state vector x and N the number of data points.
Step 3. Take M_{ n } as the 2^{ nd } until the last column of M and take M_{ o } as the 1^{ st } until the last but one column of M. Here, the subscripts n and o stand for new and old to denote a cause-effect relation between the matrices M_{ n } and M_{ o }.
Step 4. Augment the matrix M_{ o } with a (1 _{×} (N _{−} 1)) row of 1’s. This is to estimate the input sequence { Δu(k) = Γu(k),k = 0,…,N_{−} 1}. Denote the resulting augmented matrix with M_{ a }.
Step 6. Translate the discrete time result $\widehat{\Phi}\phantom{\rule{1em}{0ex}}\widehat{\Gamma}$ back to continuous-time (see e.g. the standard textbook [23]), yielding the estimates $\widehat{A}$ and $\widehat{B}$.
Having thus progressed to the bottom-right circle in Figure 2, we continue to apply the ‘optimal’ perturbation to the real system, yielding an additional set of observations that can now be used for a second parameter estimation exercise. We have now turned full circle and can start again with an improved estimate of the system A, B, i.e. ${\widehat{A}}^{1},{\widehat{B}}^{1}$. Several iterations of the loop in Figure 2 will now yield parameter estimates that have been inferred from an increasingly rich set of observations. In other words, we have maximized the information content of the model parameters to the best of our (current) knowledge of the network as to converge rapidly to the true parameter values contained in the matrix-pair A, B. Through minimization of the uncertainty bounds on the parameter estimates via optimal input design a better conditioned parameter estimation problem is gained and this is highly preferred above the ad-hoc approach where random perturbations are used as an experimental ‘design’. Although the above methodology (and more specifically the parameter estimates that are obtained using simple linear regression as in [16]) can only be used for relatively small networks, the adaptive input design approach can easily be extended to larger scale problems if other parameter estimation algorithms are utilized that can better handle larger scale networks (see e.g. [22] for a linear-programming approach to the parameter estimation problem for large networks). The message we would like to convey in the current paper is not about a parameter estimation algorithm, but about utilizing such an algorithm in the best possible way as to obtain reliable parameter estimates.
Results
Kholodenko case study
Laub-Loomis case study
The Laub-Loomis case study presents us with a far more challenging problem in comparison to the Kholondenko case study. This is mainly because in this particular case study 49 entries in the Jacobi matrix have to be estimated with linear regression, which is a large number. Furthermore, the non-linearity of the model introduces additional difficulties. It is therefore required to use relatively small perturbations in the initial condition x(0) as to not divert too far from the linearity assumption that underpins the linearized perturbation model. Since small perturbations are difficult to trace back if the noise on the data is too high, we are certainly limited in our possibilities.
Discussion and conclusions
Adaptive optimal experimental design is a natural approach for the recovery of the connections in an interaction network based on a limited number of experiments. Clearly, the methodology as presented in this paper yields promising results as was demonstrated in two case studies. The idea we have pursued is that adaptive or sequential input design closes the model identification loop of experimentation and subsequent calibration of the model, intelligently taking in to account the most recent knowledge that is available. More specifically, this means that the most recent parameter estimate of the Jacobi matrix pair [A, B] are used for subsequent analysis of the uncertainty propagation in the network. An essential feature proposed in our methodology (and also demonstrated in two case studies) is that model based experiments are now prominently included in the loop and here with the best knowledge that is available, i.e. the most recent estimate [$\widehat{A},\widehat{B}$] of the parameter vector θ. This yields a more carefully designed parameter estimation problem that allows a more reliable network reconstruction.
One could argue that the linearized system introduces limitations since, of course, the matrices A and B depend on the linearization point chosen. Especially if the Jacobi matrix structure changes significantly (in terms of zeros and non-zeros) this can introduce a problem. A remedy would then be to repeat the algorithm for different values of the state vector x. However, we emphasize that if there is no connection between node i and j in the network then regardless of the linearization point there will always be a zero in the ij^{ th } entry. If, therefore, a zero is present in more than one linearization point, the chances are increasingly high that there is indeed no connection between node i and j in the network structure. Our approach contributes to the field of systems biology where the recovery of networks based on stimulus-response data is a central topic that has already caught a lot of attention [6, 7]. In systems engineering, model based experimentation has matured substantially but here the design is usually developed for recovery of a set of model parameters in the original non-linear model structure and not for network inference, see e.g. [20] where the input design is formulated in a recursive manner, meaning that the data and subsequent planning of a new experiment are processed/analysed on-line. Of course, the results presented in the current paper do not include large scale networks with several hundreds of differential equations. But, then, this is not our main message. Rather, we have shown that OID based on dynamic time series of simple perturbation experiments leads to a better overall performance of the parameter estimation algorithm. In addition, the idea of adaptive input design may very well be combined with parameter estimation algorithms that are especially geared for such large scale problems and we think this will almost certainly improve the efficiency of subsequent experimentation and calibration as a means to unravel the underlying structure in a network model of a biological system.
Declarations
Authors’ Affiliations
References
- Keurentjes JJB, Angenent GC, Dicke M, Dos Santos VAPM, Molenaar J, vander Putten WH, de Ruiter PC, Struik PC, Thomma BPHJ: Redefining plant systems biology: from cell to ecosystem. Trends Plant Sci. 2011, 16 (4): 183-190. 10.1016/j.tplants.2010.12.002.PubMedView ArticleGoogle Scholar
- Barabasi AL, Oltvai ZN: Network biology: Understanding the cell’s functional organization. Nat Rev Genet. 2004, 5 (2): 101-U15-View ArticleGoogle Scholar
- Lucas M, Laplaze L, Bennett MJ: Plant systems biology: network matters. Plant Cell Environ. 2011, 34 (4): 535-553. 10.1111/j.1365-3040.2010.02273.x.PubMedView ArticleGoogle Scholar
- Dorf RC: Time-Domain Analysis and Design of Control Systems. 1965, Addison-Wesley: Reading, MassGoogle Scholar
- Ljung L: System Identification–Theory for the User. 1987, Prentice HallGoogle Scholar
- Marchal K, De Smet R: Advantages and limitations of current network inference methods. Nat Rev Microbiol. 2010, 8 (10): 717-729.Google Scholar
- Marchal K, Zhao H, Cloots L, Vanden Bulcke T, Wu Y, De Smet R, Storms V, Meysman P, Engelen K: Query-based biclustering of gene expression data using Probabilistic Relational Models. BMC Bioinformatics. 2011, 12 (Suppl. 1): S37-PubMedPubMed CentralGoogle Scholar
- D’Haeseleer P, Liang SD, Somogyi R: Genetic network inference: fromco-expression clustering to reverse engineering. Bioinformatics. 2000, 16 (8): 707-726. 10.1093/bioinformatics/16.8.707.PubMedView ArticleGoogle Scholar
- Margolin A, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Dalla Favera R, Califano A: ARACNE: An algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC Bioinformatics. 2006, 7 (1): S7-10.1186/1471-2105-7-S1-S7.PubMedPubMed CentralView ArticleGoogle Scholar
- Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, Kasif S, Collins JJ, Gardner TS: Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 2007, 5: 54-66. 10.1371/journal.pbio.0050054.View ArticleGoogle Scholar
- Menendez P, Kourmpetis YAI, ter Braak CJF, van Eeuwijk FA: Gene Regulatory Networks from Multifactorial Perturbations Using Graphical Lasso: Application to the DREAM 4 Challenge. PLoS One. 2010, 5 (12): e14147-10.1371/journal.pone.0014147.PubMedPubMed CentralView ArticleGoogle Scholar
- Huynh-Thu VA, Irrthum A, Wehenkel L, Geurts P: Inferring Regulatory Networks from Expression Data Using Tree-Based Methods. PLoS One. 2010, 5 (9): e12776-10.1371/journal.pone.0012776.PubMedPubMed CentralView ArticleGoogle Scholar
- Bonneau R, Reiss DJ, Shannon P, Facciotti M, Hood L, Baliga NS, Thorsson V: The Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets denovo. Genome Biol. 2006, 7 (5): 10.1186/gb-2006-7-5-r36.
- Greenfield A, Madar A, Ostrer H, Bonneau R: DREAM 4: Combining Genetic and Dynamic Information to Identify Biological Networks and Dynamical Models. PLoS One. 2010, 5 (10): e13397-10.1371/journal.pone.0013397. 10.1371/journal.pone.0013397.PubMedPubMed CentralView ArticleGoogle Scholar
- Kholodenko BN, Kiyatkin A, Bruggeman FJ, Sontag E, Westerhoff HV, Hoek JB: Untangling the wires: A strategy to trace functional interactions in signaling and gene networks. Proc Natl Acad Sci USA. 2002, 99 (20): 12841-12846. 10.1073/pnas.192442699.PubMedPubMed CentralView ArticleGoogle Scholar
- Schmid TH, Cho KH, Jacobsen EW: Identification of small scale biochemical networks based on general type system perturbations. FEBS J. 2005, 272 (9): 2141-2151. 10.1111/j.1742-4658.2005.04605.x.View ArticleGoogle Scholar
- Sontag E: Some new directions in control theory inspired by systems biology. Syst Biol. 2004, 1: 9-18. 10.1049/sb:20045006.View ArticleGoogle Scholar
- Gardner TS, Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, Kasif S, Collins JJ: Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 2007, 5: 54-66. 10.1371/journal.pbio.0050054.View ArticleGoogle Scholar
- Laub MT, Loomis WF: A molecular network that produces spontaneous oscillations in excitable cells of Dictyostelium. Mol Biol Cell. 1998, 9 (12): 3521-3532.PubMedPubMed CentralView ArticleGoogle Scholar
- Stigter JD, Vries D, Keesman KJ: On adaptive optimal input design: A bioreactor case study. Aiche J. 2006, 52 (9): 3290-3296. 10.1002/aic.10923.View ArticleGoogle Scholar
- Steinke F, Seeger M, Tsuda K: Experimental Design for Efficient Identification of Gene Regulatory Networks Using Sparse Bayesian Models. BMC Syst Biol. 2007, 1: 51-10.1186/1752-0509-1-51.PubMedPubMed CentralView ArticleGoogle Scholar
- Zavlanos M, Julius AA, Boyd S, Pappas G: Inferring Stable Genetic Networks from Steady-State Data. Automatica. 2011, 47: 1113-1122. 10.1016/j.automatica.2011.02.006.View ArticleGoogle Scholar
- Franklin G, Powell J, Workman M: DigitalControlSystems. 1990, Addison Wesley: Second EditionGoogle Scholar
- Walter E, Pronzato L: Qualitative and Quantitative Experiment Design for Phenomenological Models–ASurvey. Automatica. 1990, 26 (2): 195-213. 10.1016/0005-1098(90)90116-Y.View ArticleGoogle 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.