Systematic calibration of a cell signaling network model
- Kyoung Ae Kim^{1},
- Sabrina L Spencer^{2, 3},
- John G Albeck^{2, 3},
- John M Burke^{2, 3},
- Peter K Sorger^{2, 3},
- Suzanne Gaudet^{2, 3} and
- Do Hyun Kim^{1}Email author
DOI: 10.1186/1471-2105-11-202
© Kim et al; licensee BioMed Central Ltd. 2010
Received: 12 June 2009
Accepted: 23 April 2010
Published: 23 April 2010
Abstract
Background
Mathematical modeling is being applied to increasingly complex biological systems and datasets; however, the process of analyzing and calibrating against experimental data is often challenging and a rate limiting step in model development. To address this problem, we developed a systematic methodology for calibrating quantitative models of dynamic biological processes and illustrate its utility by validating a model of TRAIL (Tumor necrosis factor Related Apoptosis-Inducing Ligand)-induced cell death.
Results
We propose a serial framework integrating analysis and calibration modules and we compare various methods for global sensitivity analysis and global parameter estimation. First, adequacy of the network structure is checked by global sensitivity analysis to changes in concentrations of molecular species, validating that the model can reproduce qualitative features of the system behavior derived from experiments or literature surveys. Second, rate parameters are ranked by importance using gradient-based and variance-based sensitivity indices, and we systematically determine the optimal number of parameters to include in model calibration. Third, deterministic, stochastic and hybrid algorithms for global optimization are applied to estimate the values of the most important parameters by fitting to time series data. We compare the performance of these three optimization algorithms.
Conclusions
Our proposed framework covers the entire process from validating a proto-model to establishing a realistic model for in silico experiments and thereby provides a generalized workflow for the construction of predictive models of complex network systems.
Background
Comprehensive and predictive models of biological systems are expected to improve our ability to analyze complex systems, from molecular pathways to populations of organisms. Thus, there is much interest in sophisticated computational modeling techniques and high-throughput data generation [1]. One of the major difficulties in modeling cell signaling networks is the identification of the directionality and strength of relationship between molecular species in specific pathways. However, once this has been done, the knowledge can be formalized in mathematical models based on various computational methods. In particular, differential equations are widely used in biological modeling to describe dynamic processes in terms of rates of change [2–4]. The variables in these models represent the concentrations of molecular species and the directionality and strength of their relationships are encoded in the rate parameters governing their interactions. Following the construction of a mathematical representation, cycles of experimental validation and model improvement are essential for generating a predictive model, by ensuring that all required molecular species are adequately represented and that the parameter values are accurate. However, calibration of the mathematical model is not trivial because non-linearity and feedback/feedforward connections commonly found in cell signaling pathways make the analysis difficult [5, 6]. Here, we develop a systematic methodology for validating quantitative models of biological processes and apply our methodology to an existing model of TRAIL-induced apoptosis [7].
Systematic procedure of model calibration
Model calibration or regression by data fitting is necessary for computational modeling in any field of science or engineering. Systems biology faces the same challenge to construct experimentally validated models. However, formal tools for quantitative biological models have not been established yet and manual analysis is common in practice. In fact, manual fitting has the advantage that researchers may apply their experimental intuition or prior knowledge to the model relatively easily with minimal aid of mathematical or computational skills. However, the structural complexity of signaling pathways makes it difficult to fit the model heuristically based on intuition or simple analyses only. There are three dominant differences between manual fitting and systematic calibration: (1) As in Yang's work [8], manual fitting is attempted to estimate uncertain parameter values which cannot be decided directly by experimental measurement or literature. On the other hand, the systematic calibration in our study aims principally to estimate, among uncertain parameters, only the most important. We investigated the individual effect of parameters and focused on the dominant parameters to calibrate the model. (2) Manual fitting is carried out mainly by a trial-and-error process that does not guarantee optimal fit of the model. On the other hand, our systematic calibration method approaches the problem globally over the multi-dimensional domain of important uncertain parameters. Thus, it has higher probability of finding the optimal solution. (3) Manual fitting ends with what are, at the time, the best parameter values, while systematic calibration provides additional information, such as important subsets of pathways in a network or possible local optimum solutions.
The computational methods used for analyzing the network model
Computational methods | Searching principle | Reference | |
---|---|---|---|
Sensitivity analysis | Local sensitivity | Local | |
Average of local sensitivities | Global | [17] | |
Sobol's method | [11] | ||
First-order sensitivity | Global | ||
Total-effect sensitivity | Global | ||
Parameter estimation | Local estimation | Local, deterministic | |
Multi-start of local estimations | Global, deterministic | ||
Evolutionary strategy using stochastic ranking (SRES) | Global, stochastic | [22] |
In the third step, we perform a quantitative fit, or calibration, of the model to experimental data, to determine parameter values that minimize the deviations between experimental results and model simulations (Figure 1; Estimation of globally optimal fit). Parameter estimation by global optimization has been developed for engineering optimization problems [15, 16]. Below we investigate the advantages and disadvantages of three methods for biological applications, including computational efficiency, and compare the results (Table 1).
Lastly, as the model evolves in light of newly available data, the overall procedure should be iterated. We believe that by implementing the intermediary steps where sensitivity analyses are used both to assess the qualitative behavior of the model and determine which parameters to optimize, our systematic method will significantly facilitate model calibration.
Results and discussion
Qualitative analysis of a proto-model
Analysis of sensitivity with respect to initial species concentrations provides a criterion for the qualitative correctness of a cell signaling model. Sensitivity analysis assesses how changes in model inputs contribute to model output variability, and its ability to deduce model input-output relationships makes sensitivity analysis one of the critical parts of model development, verification, and evaluation. Changes in initial species concentrations can mimic the effects of mutations or changes in the expression level of the molecular players involved, and the sensitivity of the model output to changes in initial species concentrations should match the expected change in system behavior.
where model outputs and inputs are represented as y_{ i }and p_{ j }respectively. This method is often called local sensitivity because it reflects output variability accurately near a given nominal input value, p*. However, most kinetic parameters are quite uncertain and a range rather than a single parameter value is available, either from the literature or from biophysical constraints on the reactions. Thus, "model-independent", or more precisely, parameter-independent, global sensitivity analysis techniques have generated great interest[10]. Averaging of local sensitivities over a range of plausible values for uncertain parameters is one possible method for global sensitivity. Local sensitivities are calculated with multiple parameter choices that are selected randomly or regularly within parameter ranges. The sensitivities for those parameter choices, integrated over the time interval of interest in the monitoring of the output, ∫|S_{ ij }(t)|dt, are then averaged to determine global sensitivity [17]. Importantly, because integration over time and averaging are necessary, a compromise must be made between accurately calculating the magnitude of the effects by using absolute sensitivity values, and assessing the directionality of the effects, by maintaining the sign of the values.
The model on which we applied our methodology simulates the response of a single cell to TRAIL. TRAIL is a protein ligand which triggers the process of programmed cell death, or apoptosis. This model of TRAIL-induced cell death signaling network encompasses the activation of initiator (caspase-8 or C8) and effector (caspase-3 or C3) caspases, the onset of mitochondrial outer membrane permeabilization and the death of the cell, as marked by cleavage of the caspase-3 substrate, PARP. According to a recent study [7], extrinsic apoptosis shows a specific behavior of all-or-none effector caspase activation at the single-cell level. As the authors termed it, the process shows "variable-delay, snap-action": a long, variable delay between TRAIL stimulation and effector caspase activation is followed by rapid and sudden progression to completion. The original model is composed of 58 ordinary differential equations based on mass action kinetics. Eighteen out of 58 protein species have non-zero initial concentrations, and 70 rate constants regulate the reactions in the model network. The original parameters were determined from the literature and manual fitting. In this study, we applied our methods to analyze qualitative properties of the model and fit the model to dynamic quantitative experimental data in a systematical and computationally effective way. Hereafter, the original model in [7] will be referred to as the manually calibrated model, to distinguish it from our improved model.
Second, the absolute value of sensitivity provides a measure of how strongly the perturbation of a single species' concentration affects the model output. The sensitivity with respect to perturbation of XIAP was found to be relatively high on average, implying that the model output can be changed dramatically by small changes in XIAP (Figure 2). This prediction is supported by biological evidence that XIAP directly inhibits the enzymatic activity of caspases and the degree of inhibition is highly dependent on the concentration of XIAP[18]. Cleavage of PARP is insensitive to the concentration of caspase-6 (C6), in agreement with experiments in which reducing the expression of caspase-6 by ~90% did not affect TRAIL-induced cell death (Figure 2; [7]). Overall, our sensitivity analysis agreed with the known effects of varying protein concentration. If, however, the signs or strengths of the sensitivities in our analysis had not agreed with experimental results, the model construction would have to be re-examined. Modification of the model and this qualitative analysis would be done iteratively until a satisfactory result could be reached. Although the TRAIL model study does not provide us with an example of failure of qualitative agreement at this step of the procedure, it is still worth noting that qualitative agreement with known experimental system behavior can be a strong preliminary criteria for adequacy of the model structure. In effect, it sets a minimal qualification that must be met before more computationally intensive methods are applied to improve the proto-model by quantitative fitting.
In a third type of assessment of the results our sensitivity analysis of PARP cleavage, we analyzed the influence of the uncertainty of rate constants on the sensitivity with respect to initial species concentrations. Sensitivities that are not affected by parameter values will have narrow distributions, and by consequence, their sensitivity value is very reliable. The sensitivities related to the perturbation of some species like XIAP and caspase-8 were found to be broadly distributed and thus to be relatively uncertain (Figure 2). Particularly interesting is the fact that the sensitivity of PARP cleavage to caspase-8 is negative in some cases, even if it is known to have a pro-apoptotic function, invalidating certain parameter sets.
Dominant parameter selection
Comparison of results from different global sensitivity indices
Average of local sensitivities | Sobol's first order sensitivity | Sobol's total effect sensitivity | |
---|---|---|---|
8 dominant parameters* | k_{8}, k_{9}, k_{ 5 }, kc_{1}, k_{3}, k_{ 1 }, k_{-_1}, k_{4} | kc_{1}, k_{ 12 }, k_{3}, k_{ 13 }, k_{9}, k_{24}, k_{ 10 }, k_{-1} | kc_{1}, k_{3}, k_{ 12 }, k_{9}, k_{ 13 }, k_{ 1 }, k_{ 10 }, k_{ 5 } |
CPU time | 24 hours | 160 hours | 160 hours |
Objective function^{†} | 3.546 | 8.364 | 7.745 |
Importantly, there are often biologically meaningful quantities of interest for which partial derivatives cannot be defined, and these may be the outputs for which the dominant parameters need to be identified. For example, for TRAIL-induced apoptosis we can define biologically meaningful features of the dynamic behavior of cell death. One example is the delay time (t_{delay}) that measures how long it takes from the time of TRAIL addition to the time at which 50% of PARP is cleaved. Another is the switching time (t_{switch}), which measures the rapidity of PARP cleavage after caspase-3 (C3) activation. These features are variables that are discontinuous with respect to input parameter variation, and to determine the dominant parameters in controlling t_{delay,} we therefore explored other sensitivity analysis methods to replace gradient-based sensitivity analysis.
To determine which parameters dominate the control of PARP cleavage dynamics and t_{delay}, the model parameters were ranked by highest to lowest amplitude in sensitivities (Figure 4) and the eight most dominant parameters from each of the three sensitivity indices are listed in Table 2. The parameters that are commonly selected by all three methods are bolded, and those selected by two are underlined; the nomenclature of the parameters follows that of Albeck et al[7]. For example, k_{9}, which is the forward reaction rate constant of PARP cleavage by caspase-3, is ranked within the eight dominant parameters by all three sensitivity indices. k_{3} and kc_{1,} relevant to caspase-8 activation and death ligand binding to the receptor respectively, are also dominant by all three methods. Even though all the reactions in the network play a role in cell death signaling, the sets of reactions rate constants listed in Table 2 were identified as the most critical in regulating the dynamic of PARP cleavage. This prediction, that reactions relevant to caspase-8 activation are critical in regulating the delay time to death was arrived at by our computational sensitivity analysis, but, importantly, it is supported by experimental evidence: the reactions involved in caspase-8 substrate cleavage strongly influence t_{delay}[19].
Once the ranking of parameters has been determined, the next question is how many parameters to target during a calibration to accurately capture network behavior. While there are no general and definitive criteria, it should be noted that estimation of too many parameters increases the number of degrees of freedom and the probability that inadequate local optima are detected. On the other hand, choosing too few parameters decreases fitting performance as well as the reliability of the optimal solution. To address this trade-off, we used the ranked parameters to determine the optimal cut-off for the calibration of the model of TRAIL-induced cell death. In Figure 4b, the 70 parameters on the x-axis were ordered by their ranking number (determined from Figure 4a). We observed that the sensitivities dropped off sharply after a few steps - ranked sensitivities generated L-shaped curves. The three different sensitivity algorithms have a common property that the parameter of 8^{th} highest sensitivity was approximately at the border between horizontal and vertical lines. This analysis suggested that for this particular model, the eight most sensitive parameters can cover much of the variation in PARP cleavage, and should be sufficient to include in model calibration.
Parameter estimation by global optimization
Among the various approaches for global parameter estimation, the simplest one is the deterministic multi-start method where a large number of local estimations start from different initial parameter combinations (Figure 5(a); red circles). The logarithmic space of parameters is divided uniformly in a grid and deterministic local estimation starts from every grid point, comparing the fit of nearby points. Because the entire parameter space is explored, the guarantee for finding the global optimum is high, as long as the grid samples the space sufficiently well. In Figure 5(a), parameter sets starting from initial grid points converge to the points aligned along the valley after local estimations have terminated. However the computational load increases exponentially with the number of parameters, as dimensions are added to the sampling grid. To overcome this difficulty, random sampling in a Latin hypercube of parameter space[20] or parallel computing with cluster processors could be utilized.
Stochastic methods on the other hand, can find the global solution with relatively less computational effort. These methods start with parameter values that are randomly sampled in parameter space, and, according to a set of rules, explore new solutions in the neighborhood of the initial point looking for a better solution and repeat until no further improvement of fit is found. Genetic algorithms and simulated annealing are well known examples of stochastic methods[21]. In a comparative study of various optimization methods, Stochastic Ranking Evolutionary Strategy (SRES) showed the best performance [15]. In SRES, a "population" composed of randomly selected "elements", or sets of parameter values, is generated. The elements are ranked by their fit to the data using a bubble-sort procedure[22]. Only highly ranked elements are retained as ancestors for the next generation, which are used to probabilistically produce a new population of random elements with a better fit, on average. The source code of SRES is available in the public domain[23].
It is not surprising that many local minima were detected using a multi-start method because nonlinear and complex models like cell signaling networks usually exhibits objective function surfaces with multiple local optima. The plateaus near the global optimum and around the objective function values of 40 and 100 in Figure 6(a) could be due to: 1) a wide well on the hypothetical surface of parameter space so that estimations from many nearby starting points converge to a single minimal solution, allowing us to easily arrive at the optimal solution or 2) a valley-shaped local optima on the surface of objective function. Because a wide well on the surface of parameter space is rare in network models, the most likely causes of the plateaus were valley-shaped optima. Along a valley, solutions may be distinct if they are located far apart from one another, but nevertheless fit the model in similarly well. Ideally, when constructing predictive models, this situation should be avoided by reducing valleys to more focused wells on the surface of parameter space by adding constraints to the optimization problem.
Comparison of estimation performance by different optimization algorithms
Local optimization | Multi-start of local optimization | SRES^{†}(G^{‡} = 33) | SRES (G = 300) | SRES (G = 33) & local optimization | |
---|---|---|---|---|---|
Number of surveyed parameter combinations | 1 | 6561 | 6600 | 60000 | 6600 |
CPU time | ~1 minute | 2600 hours | 2 hours | 18 hours | 2 hours |
Objective function | 45.04 | 0.2447 | 3.546 | 1.739 | 0.2679 |
Influence of dominant parameter choice on optimization
Conclusion
In this report, we proposed a framework for efficiently calibrating computational models of biological systems, and applied it to a model of TRAIL-induced apoptosis while comparing several sensitivity analysis methods and model optimization algorithms. Importantly, we showed how sensitivity analysis can be used to rapidly test whether the model structure adequately allows qualitative matching to the behavior of the biological system. This step implements a minimal qualification, focusing the initial search on the qualitative performance of the proto-model. Within our framework, this validation step is required before proceeding to quantitative optimization of the model, ensuring that computationally costly optimization algorithms are used effectively. Furthermore, we showed how global sensitivity analysis methods can be used to identify the parameters that dominantly regulate the dynamics of the output of interest. With the application of Sobol's algorithms, we were also able to identify parameters that control the TRAIL-induced delay time to cell death (t_{delay}), a biologically relevant quantity that is not a state variable of the model. Undoubtedly, this type of sensitivity analysis will prove useful within our outlined framework for other models as well, for example in models of oscillatory systems where, in certain cases, the period of the oscillations is more meaningful than their amplitude. Finally, while comparing different model calibration algorithms, we showed that global sensitivity analysis could successfully identify the parameters to include in quantitative optimization, allowing great computational savings by constraining the search to the important model dimensions. In the future, we foresee that the predictive quality of models would be further improved by repeating this cycle of model validation, identification of dominant parameters and optimization with different model outputs that are controlled by other parameters, allowing the determination of more and more parameter values.
Methods
Mathematical model and experimental data
The serial methodology was applied to fit a recently described mathematical model of TRAIL-induced cell death signaling [7]. This model is composed of ordinary differential equations based on mass action kinetics. Although most ODE models assume simply that the inside of cell is a mixed soup and do not include spatial information, the original model describes reactions and transport of molecular species in two compartments; cytoplasm and mitochondria. The authors verified that sudden activation of effector caspase after a long delay is related to permeabilization of the mitochondrial membrane and relocalization of certain proteins. To evaluate our methodology, we carried out model calibration using the original ODE model and the experimental data. The live-cell imaging data were obtained by microscopy monitoring of a population of HeLa cells treated with 10 ng/ml, 50 ng/ml, or 250 ng/ml of TRAIL and 2.5 μg/ml cycloheximide [7]. Although the cells were isogenic, the delay period until sudden death (t_{delay}) varied from cell to cell due to inherent fluctuations in cell state [19]. Figure 2B in [7] actually shows 5 examples of single-cell dynamic data for each condition; these examples were chosen after monitoring well over 100 cells in each condition. However we focus here on deterministic modeling at the single-cell level. Thus, for each condition we chose a single representative cell whose t_{delay} is the median in the population of more than 100 cells. The cleavage of effector caspase reporter protein (EC-RP) was quantified at 3-min intervals, and used for fitting the model output corresponding to cleavage of PARP, the effector caspase substrate. In addition, once the rapid cleavage of EC-RP is complete, then the output signal cannot be accurately measured by microscopy - it becomes extremely noisy as the cells detaches from the surface and moves out of the focal plane, and thus poorly reports cellular activity. Therefore we neglected the fluctuations in the experimental data after cleavage of PARP reaches a value of 1. Instead, we fit the model to a plateau with a value of one.
where p is the set of parameters, N_{exp} is the number of experiments, w_{ i }is the weight associated with the measurement of the i^{th} experiment, , and y_{ i }(p) is the corresponding value computed from the model. The weight may be given differently depending on reliability of specific experimental measurement. If there are uncertain or less confident data points, those should take less part in evaluating the objective function by using smaller weight than other data points. Since we took all the data with equal importance, the weights were all set equal to 1 in this case. The estimation calculation was stopped if the normalized difference of objective function values between two successive iterations was less than 1E-6.
Parameter space
The plausible range for uncertain rate constant parameter values was set around a nominal value of the corresponding parameter. For most parameters, the upper and lower bounds were set as 100 and 1/100 times the nominal value, respectively. For several parameters whose nominal values are considered to be relatively certain by previous experiences we set narrower ranges to minimize the effect of uncertain parameter values in the global sensitivity analysis. For instance, forward, backward and catalytic rate constants relevant to caspase-3 cleavage by caspage-8 (k_{5}, k_{–5}, kc_{5}), caspase-6 cleavage by caspage-3 (k_{6}, k_{–}_{6}, kc_{6}), PARP cleavage by caspase-3 (k_{9}, k_{–9}, kc_{9}), Bid activation by cleaved caspase-8 (k_{10}, k_{–}_{10}, kc_{10}), Bcl2 reacting with activated Bax in the mitochontrial compartment (k_{14}, k_{–14}), Bax_{2} reacting with Bcl2 (k_{16}, k_{–16}), Bax4 reacting with Bcl2 (k_{18}, k_{–}_{18}), Cytochrome c and Smac release from mitochondria (k_{20}, k_{–}_{20}, kc_{20}, k_{21}, k_{–21}, kc_{21}) had a range set to between 10 and 1/10 times of their nominal value. The rate constants regarding ubiquitination of cleaved caspase-3 by XIAP (k_{8}, _{k–}_{8}, kc_{8}) and XIAP reacting with the apoptosome and Smac (k_{27}, k_{–}_{27}, k_{28}, k_{–28}) have an even narrower range between 2 and 1/2 times the nominal value. The nomimal values were either obtained from the literature or set by trial and error to allow the model to reproduced experimental data, as previously described [7].
Computations
All computations were performed using JACOBIAN^{®} 4.0, a dynamic modeling software provided by Numerica Technology, LLC. The local estimation was executed by the built-in JACOBIAN^{®} function of Limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) estimation solver and a weighted least square objective function. The High-Performance Computing facility at Harvard Medical School was utilized for intensive computations. The repetitive jobs of the multi-start estimation as well as Sobol's sensitivity analysis were parallelized and distributed to over 200 computing nodes (AMD dual core processors). For comparisons between different algorithms, the CPU usage time of each node was summed as if a single computing machine was utilized.
Declarations
Acknowledgements
The authors thank Glen Ko and Taeshin Park in Numerica Technology for their support concerning JACOBIAN^{®} software.
Authors’ Affiliations
References
- Palsson B: The challenges of in silico biology. Nat Biotech 2000, 18: 1147–1150. 10.1038/81125View ArticleGoogle Scholar
- Wolkenhauer O, Ullah M, Wellstead P, Cho K-H: The dynamic systems approach to control and regulation of intracellular networks. FEBS Lett 2005, 579: 1846–1853. 10.1016/j.febslet.2005.02.008View ArticlePubMedGoogle Scholar
- Aldridge BB, Burke JM, Lauffenburger DA, Sorger PK: Physicochemical modelling of cell signalling pathways. Nat Cell Biol 2006, 8: 1195–1203. 10.1038/ncb1497View ArticlePubMedGoogle Scholar
- Polynikis A, Hogan SJ, di Bernardo M: Comparing different ode modelling approaches for gene regulatory networks. J Theor Biol 2009, 261: 511–530. 10.1016/j.jtbi.2009.07.040View ArticlePubMedGoogle Scholar
- Asthagiri AR, Lauffenburger DA: Bioengineering models of cell signaling. Annu Rev Biomed Eng 2000, 2: 31–53. 10.1146/annurev.bioeng.2.1.31View ArticlePubMedGoogle Scholar
- Tyson JJ, Chen KC, Novak B: Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol 2003, 15: 221–231. 10.1016/S0955-0674(03)00017-6View ArticlePubMedGoogle Scholar
- Albeck JG, Burke JM, Spencer SL, Lauffenburger DA, Sorger PK: Modeling a snap-action, variable-delay switch controlling extrinsic cell death. PLoS Biol 2008, 6: 2831–2852. 10.1371/journal.pbio.0060299View ArticlePubMedGoogle Scholar
- Yang K, Ma W, Liang H, Ouyang Q, Tang C, Lai L: Dynamic simulations on the arachidonic acid metabolic network. PLoS Comput Biol 2007, 3: 523–530.Google Scholar
- Frey HC, Patil SR: Identification and review of sensitivity analysis methods. Risk Anal 2002, 22: 553–578. 10.1111/0272-4332.00039View ArticlePubMedGoogle Scholar
- Saltelli A, Tarantola S, Chan KP-S: A quantitative model-independent method for global sensitivity analysis of model output. Technometrics 1999, 41: 39–56. 10.2307/1270993View ArticleGoogle Scholar
- Sobol IM: Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates. Math Comput Simulat 2001, 55: 271–280. 10.1016/S0378-4754(00)00270-6View ArticleGoogle Scholar
- Zheng Y, Rundell A: Comparative study of parameter sensitivity analyses of the tcr-activated erk-mapk signalling pathway. IEE P Syst Biol 2006, 153: 201–211. 10.1049/ip-syb:20050088View ArticleGoogle Scholar
- Cho K-H, Shin S-Y, Kolch W, Wolkenhauer O: Experimental design in systems biology, based on parameter sensitivity analysis using a monte carlo method: A case study for the tnf{alpha}-mediated nf-{kappa} b signal transduction pathway. Simulation 2003, 79: 726–739. 10.1177/0037549703040943View ArticleGoogle Scholar
- Zi Z, Cho K-H, Sung M-H, Xia X, Zheng J, Sun Z: In silico identification of the key components and steps in ifn-[gamma] induced jak-stat signaling pathway. FEBS Lett 2005, 579: 1101–1108. 10.1016/j.febslet.2005.01.009View ArticlePubMedGoogle Scholar
- Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: A comparison of global optimization methods. Genome Res 2003, 13: 2467–2474. 10.1101/gr.1262503View ArticlePubMedPubMed CentralGoogle Scholar
- Mendes P, Kell D: Non-linear optimization of biochemical pathways: Applications to metabolic engineering and parameter estimation. Bioinformatics 1998, 14: 869–883. 10.1093/bioinformatics/14.10.869View ArticlePubMedGoogle Scholar
- Bentele M, Lavrik I, Ulrich M, Stosser S, Heermann DW, Kalthoff H, Krammer PH, Eils R: Mathematical modeling reveals threshold mechanism in cd95-induced apoptosis. J Cell Biol 2004, 166: 839–851. 10.1083/jcb.200404158View ArticlePubMedPubMed CentralGoogle Scholar
- Deveraux QL, Takahashi R, Salvesen GS, Reed JC: X-linked iap is a direct inhibitor of cell-death proteases. Nature 1997, 388: 300–304. 10.1038/40901View ArticlePubMedGoogle Scholar
- Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK: Non-genetic origins of cell-to-cell variability in trail-induced apoptosis. Nature 2009, 459: 428–432. 10.1038/nature08012View ArticlePubMedPubMed CentralGoogle Scholar
- Stein M: Large sample properties of simulations using latin hypercube sampling. Technometrics 1987, 29: 143–151. 10.2307/1269769View ArticleGoogle Scholar
- Edgar TF, Himmelblau DM, Lasdon LS: Optimization of chemical processes. second edition. Singapore: McGraw-Hill; 2001.Google Scholar
- Runarsson TP, Xin Y: Stochastic ranking for constrained evolutionary optimization. IEEE T Evolut Comput 2000, 4: 284–294. 10.1109/4235.873238View ArticleGoogle Scholar
- Stochastic ranking with evolution strategy for matlab[http://www3.hi.is/~tpr/index.php?page=software/sres/sres]
- Rodriguez-Fernandez M, Mendes P, Banga JR: A hybrid approach for efficient and robust parameter estimation in biochemical pathways. Biosystems 2006, 83: 248–265. 10.1016/j.biosystems.2005.06.016View ArticlePubMedGoogle 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.