Modeling metabolic networks in C. glutamicum: a comparison of rate laws in combination with various parameter optimization strategies
- Andreas Dräger^{1}Email author,
- Marcel Kronfeld^{1}Email author,
- Michael J Ziller^{1}Email author,
- Jochen Supper^{1},
- Hannes Planatscher^{1},
- Jørgen B Magnus^{2, 3},
- Marco Oldiges^{3},
- Oliver Kohlbacher^{1} and
- Andreas Zell^{1}
DOI: 10.1186/1752-0509-3-5
© Dräger et al; licensee BioMed Central Ltd. 2009
Received: 01 June 2008
Accepted: 14 January 2009
Published: 14 January 2009
Abstract
Background
To understand the dynamic behavior of cellular systems, mathematical modeling is often necessary and comprises three steps: (1) experimental measurement of participating molecules, (2) assignment of rate laws to each reaction, and (3) parameter calibration with respect to the measurements. In each of these steps the modeler is confronted with a plethora of alternative approaches, e. g., the selection of approximative rate laws in step two as specific equations are often unknown, or the choice of an estimation procedure with its specific settings in step three. This overall process with its numerous choices and the mutual influence between them makes it hard to single out the best modeling approach for a given problem.
Results
We investigate the modeling process using multiple kinetic equations together with various parameter optimization methods for a well-characterized example network, the biosynthesis of valine and leucine in C. glutamicum. For this purpose, we derive seven dynamic models based on generalized mass action, Michaelis-Menten and convenience kinetics as well as the stochastic Langevin equation. In addition, we introduce two modeling approaches for feedback inhibition to the mass action kinetics. The parameters of each model are estimated using eight optimization strategies. To determine the most promising modeling approaches together with the best optimization algorithms, we carry out a two-step benchmark: (1) coarse-grained comparison of the algorithms on all models and (2) fine-grained tuning of the best optimization algorithms and models. To analyze the space of the best parameters found for each model, we apply clustering, variance, and correlation analysis.
Conclusion
A mixed model based on the convenience rate law and the Michaelis-Menten equation, in which all reactions are assumed to be reversible, is the most suitable deterministic modeling approach followed by a reversible generalized mass action kinetics model. A Langevin model is advisable to take stochastic effects into account. To estimate the model parameters, three algorithms are particularly useful: For first attempts the settings-free Tribes algorithm yields valuable results. Particle swarm optimization and differential evolution provide significantly better results with appropriate settings.
Background
The metabolism of whole cells can be described as a network of metabolites and reactions interconverting these metabolites. To understand cellular systems, dynamic modeling of cellular processes has become an important task in systems biology [1–4]. Dynamic models describe the whole system and the state of each reacting species therein in a time-dependent manner [2]. Once such a model is constructed, several network properties can be derived, for instance stability, robustness, or the long-term behavior [5]. Furthermore, a well developed model provides a basis for predictions under different perturbations or varied environmental circumstances and can be applied to enhance the yield of desired metabolic products like certain amino acids [3].
To set up such models, appropriate rate laws have to be assigned to each reaction within the network. From these, a differential equation system that characterizes the rates of change of each reactant can be derived. However, setting up model equations is a difficult task. For many larger networks available in databases like KEGG[6, 7] or META CYC[8] the reaction mechanism remains unknown. In many cases, reliable rate equations for the reactions are not known because these actually have to be derived for each catalyzing enzyme individually [9]. It is therefore a common approach to apply approximative rate laws, which characterize the most important features of the reaction rate. Many rate laws, which are either continuous or discrete, and either deterministic or stochastic [2], have been proposed for this purpose. Several examples of each group exist such as probabilistic [10, 11], phenomenological [5, 12–15], or semi-mechanistic approaches [5, 16].
A second problem arises whenever a dynamic model of biochemical systems is created, because any such model contains a certain number of parameters like the reaction rates, Michaelis constants or the limiting rate as well as constants describing the influence of certain inhibitors [17–19] or, in stochastic systems, the reaction propensity [20]. Except for phenomenological models like power law approximations [21], linlog [12, 22], or loglin kinetics [13, 23], the parameter values can often be measured. However, this procedure is time-consuming, expensive, and often impractical. Online databases like the Brunswick Enzyme Database (BRENDA) [24–26] provide measured parameter values for many enzymes, but variations in the experimental settings, under which these and the time series measurements for the system under study were obtained, limit the applicability of these values for modeling purposes. In addition, it was observed that there are differences between parameters measured in vivo and in vitro [[27], p. 461]. The application of computational methods to optimize model parameters regarding the fit error has therefore become an important task in the model identification process [28–31]. In this connection, the optimizer tries to minimize the distance between measured values or values created in silico and the simulated time course for each reacting species by varying the model parameters. The smaller this distance is, the higher is the quality of a possible solution for one parameter set. This quality measure is often called the "fitness" of the parameter set. As an exhaustive search for the best solution is computationally not possible, heuristic optimization methods try to find the global optimum of the system. Usually, metabolic systems are analytically hard or infeasible to solve. Often those systems are non-convex or multimodal, i. e., contain numerous local optima, and the gradient cannot be computed easily. Biologically inspired optimization procedures like Evolutionary Algorithms (EAs) are known to handle even highly nonlinear optimization problems [32–34]. Many such optimization algorithms are freely available in several software packages [35–38] or included in commercial toolboxes [39].
During the last few decades, manifold derivatives of EAs have been proposed. Each of them has certain advantages and is therefore more or less appropriate for a special problem. Their development was driven by analogies to natural phenomena such as Darwinian evolution (genetic algorithm [40], evolution strategy [41], differential evolution [42]), hill climbing [43], the formation of crystal structures in metallurgy (simulated annealing [44]), or the swarm intelligence idea (particle swarm optimization [45, 46]). Each one of these optimization procedures provides several settings that influence its performance; for instance, the temperature in simulated annealing, the crossover probability in genetic algorithms, or the population size in particle swarm optimization. For a detailled introduction to all heuristic optimization procedures used in this work [see Additional file 1].
In several studies, heuristic optimization procedures like EAs have been applied successfully to biochemical systems after these have been translated into sets of differential equations [32–34, 47–50]. Most of the time standard settings for the optimization procedures are used. An analysis of these settings does, in fact, enable enhancement of the performance of the optimization process [33, 47, 48]. In many cases, a lack of time prevents researchers from systematically benchmarking these settings. In order to improve the model quality, however, this would be necessary. The resulting model systems, including the identified parameters, are often used to derive network properties like the long-term behavior or to perform steady-state analyses [49–53] of the system, but only in few cases detailed and specially derived rate laws could be applied to deduce a model system [49, 50]. Many studies are available, in which a single type of approximative rate equation is applied to set up a model – in many cases without a comparison with alternative approaches. Guthke et al. modeled the amino acid metabolism of primary human liver cells using a phenomenological approach [51, 53] whereas Liu and Wang used S-systems [21, 54, 55] for their biochemical models [56]. Magnus et al. applied linlog kinetics [12, 22] to model the valine and leucine metabolism in C. glutamicum [52].
Very few studies compare alternative modeling approaches to investigate their applicability for the specific problem [9, 34]. While Spieth et al. studied whether in silico time series data generated with certain model systems can be reproduced crosswise with other ones [34], Bulik et al. analyzed properties like stability when detailed kinetic equations within a system are replaced by approximative ones [9].
An investigation of both, alternative modeling approaches together with a systematical benchmark of the settings of optimization procedures, was rarely done. The disposability of software tools, which assign rate equations more or less automatically [57–61], requires the user to be especially aware of the properties of different modeling approaches and the possible quality that can be achieved with a certain type of kinetic models. Since automatically created models have already been published [62], and attempts have been made to scale up this modeling process to derive even genome-scale kinetic models automatically [63], these properties are of paramount importance. Summarizing, to construct a mathematical description of a biochemical reaction system the modeler has to consider at least two central questions:
1. Which rate laws are the most appropriate ones for the specific purpose?
2. Which optimization procedure performs best on the problem class of parameter inference?
Once a model has been created, the choice of initial conditions represents a further important question for an appropriate simulation of the model. Fundamentally, the system can be written as an initial value problem, sometimes referred to as a single-shoot approach [30], or as a multiple-shooting problem [31, 64]. Here, a single-shoot strategy is employed using the steady-state concentrations of the participating metabolites as initial values.
The reaction system in more detail
No | Reaction | Enzyme | Inhibitor |
---|---|---|---|
R _{1} | 2 Pyr → AcLac + CO_{2} | AHAS | Val |
R _{2} | AcLac + NADPH_{2} ⇌ DHIV + NADP^{+} | AHAIR | Val |
R _{3} | DHIV → KIV + H_{2}O | DHAD | Val |
R _{4} | KIV + Glut → Val + α KG | BCAAT_{ValB} | |
R _{5} | KIV + Ala → Val + Pyr | BCAAT_{ValC} | |
R _{6} | Val → Val_{ext} | Trans_{Val} | Leu |
R _{7} | KIV + AcCoA → IPM + CoA | IPMS | Leu |
R _{8} | IPM + NAD^{+} → KIC + NADH_{2} + CO_{2} | IPMDH | |
R _{9} | KIC + Glut ⇌ Leu + α KG | BCAAT_{LeuB} | |
R _{10} | Leu → Leu_{ext} | Trans_{Leu} | Val |
Results and discussion
Mathematical models
Figure 1 shows the biosynthesis of Val and Leu in C. glutamicum, starting with pyruvate (Pyr). The reactions of this pathway are summarized in Table 1. In R_{3}, one substrate (DHIV) turns into one product (KIV) when neglecting the influence of the second product, water, which is plentiful in the cell. Therefore, we can apply a Michaelis-Menten equation to model the kinetics of this reaction. The transport of Val and Leu through the cell wall (R_{6} and R_{10}) has exactly one substrate and one product and can, therefore, also be modeled using the Michaelis-Menten approach. All other reactions in the network cannot be modeled using this classical approach because multiple substrate or product molecules are involved. Therefore, approximative rate laws can be applied to these reactions. Here we employ one stochastic and three deterministic rate equations. Approximative rate laws can be used for the three reactions with a Michaelis-Menten mechanism as well.
The reversibility of the reactions constitutes another important question to be solved before modeling. As transport through the cell wall removes both products from the cellular system, we assume that there is no reverse reaction (an uptake of Val or Leu) and, hence, consider both reactions irreversible. The two reactions R_{2} and R_{9} are reversible [6–8, 52]. We let the optimization procedure "decide" if the remaining reactions should be modeled in a reversible or irreversible way. To this end, we construct one "reversible" and one "irreversible" alternative model for each rate law, keeping only the two known reactions R_{2} and R_{9} reversible.
In this way we derive the following seven models based on four approaches for the reaction velocity on this pathway. Details and the formulas can be found in Methods and [see Additional file 2].
GMAKr Pure generalized mass action kinetics, in which all reactions apart from R_{6} and R_{10} are modeled reversibly, with 24 parameters.
GMAKi Pure generalized mass action kinetics, in which only the two reactions R_{2} and R_{9} are considered reversible, with 18 parameters.
GMMr Like GMAKr but with three Michaelis-Menten equations for reactions R_{3}, R_{6} and R_{10}. This model contains 31 parameters. In the GMAKr model the influences of all enzymes are neglected and hidden in the rate constants, which is an oversimplification of the biochemical process. The model comes closer to the biochemical process when inserting Michaelis-Menten equations for the three reactions R_{3}, R_{6}, and R_{10}.
GMMi Like GMMr but with only two reversible reactions R_{2} and R_{9}, leading to 24 parameters.
CKMMr Convenience kinetics with three Michaelis-Menten equations as in GMMr. All reactions apart from R_{6} and R_{10} are considered reversible leading to 59 parameters. The convenience rate law is also an approximation of the biochemical process. Therefore, we do not construct a pure convenience kinetics model of the whole system but apply Michaelis-Menten kinetics whenever this is possible (R_{3}, R_{6} and R_{10}).
CKMMi Convenience kinetics with three Michaelis-Menten equations as in model GMMi, in which only the two reactions R_{2} and R_{9} are considered reversible, with 41 parameters.
LANG To demonstrate the possibility of large-scale parameter optimization, even for stochastic models, and to model the effects of random fluctuations in the metabolite concentrations, we consider a stochastic description as well. Based on the Langevin equation, this system contains 24 parameters.
In a glucose stimulus-response experiment 47 measurements are taken for 13 metabolites on this pathway (for details see Methods). The parameters of all models are calibrated with regard to these data.
Fine-tuned optimization algorithms and models
Settings for the standard algorithms in detail
Algorithm | Population | Mutation | Crossover | Selection |
---|---|---|---|---|
Monte Carlo | 50 | no | no | Best |
Hill Climber | 1, 10, 25, 50, 100, 250 | Fixed Step σ = 0.2, p_{ m }= 1 | no | Best |
binGA | 250 | one-point, p_{ m }= 0.1 | one-point, p_{ c }= 0.7 | Tournament, group of 8 |
realGA | 250 | global, p_{ m }= 0.1 | UNDX, p_{ c }= 0.8 | Tournament, group of 8 |
stdES | (5, 25) | global, p_{ m }= 0.8 | discrete one-point, p_{ c }= 0.2 | Best |
cmaES | (5,^{+}25) | CMA, p_{ m }= 1 | no | Best |
SA | 250 | linear annealing schedule, α = 0.1, initial T = 5 | Best | |
DE | 100 | current-to-best/1, λ = F = 0.8, CR = 0.5 | ||
PSO | 100 | star topology, ϕ_{1} = ϕ_{2} = 2.05, χ = 0.73 |
Preliminary test results
GMAK, reversible | GMAK, irreversible | ||||||||
---|---|---|---|---|---|---|---|---|---|
Min. | Algorithm | Average | Std. Dev. | Algorithm | Min. | Algorithm | Average | Std. Dev. | Algorithm |
20.334 | PSO | 21.190 | 0.576 | PSO | 24.587 | PSO | 25.967 | 1.171 | Tribes |
20.335 | DE | 21.228 | 0.756 | DE | 25.006 | Tribes | 29.502 | 9.610 | DE |
21.401 | Tribes | 21.725 | 0.275 | Tribes | 25.683 | DE | 33.169 | 10.143 | PSO |
23.097 | binGA | 26.106 | 2.204 | binGA | 25.981 | binGA | 35.670 | 3.125 | cmaESplus |
24.321 | cmaESplus | 27.598 | 2.091 | cmaESplus | 30.704 | cmaESplus | 50.663 | 4.138 | HC MS 10 |
GMM, reversible | GMM, irreversible | ||||||||
Min. | Algorithm | Average | Std. Dev. | Algorithm | Min. | Algorithm | Average | Std. Dev. | Algorithm |
20.312 | PSO | 21.272 | 0.461 | PSO | 24.477 | Tribes | 24.654 | 0.282 | Tribes |
20.407 | DE | 21.711 | 1.153 | DE | 24.499 | DE | 37.696 | 19.374 | DE |
21.590 | Tribes | 21.887 | 0.243 | Tribes | 24.553 | PSO | 41.529 | 21.768 | binGA |
22.913 | binGA | 26.742 | 2.711 | binGA | 25.266 | binGA | 31.136 | 9.052 | cmaESplus |
23.890 | cmaESplus | 26.624 | 1.851 | cmaESplus | 25.812 | stdES | 31.338 | 0.546 | HC MS 1 |
CKMM, reversible | CKMM, irreversible | ||||||||
Min. | Algorithm | Average | Std. Dev. | Algorithm | Min. | Algorithm | Average | Std. Dev. | Algorithm |
20.882 | PSO | 21.773 | 0.352 | PSO | 21.632 | PSO | 23.968 | 0.931 | DE |
21.821 | DE | 22.633 | 0.562 | DE | 22.651 | DE | 24.624 | 0.315 | Tribes |
22.258 | Tribes | 23.079 | 0.464 | Tribes | 24.191 | Tribes | 25.761 | 0.331 | binGA |
22.829 | cmaES | 24.341 | 1.026 | cmaES | 25.152 | binGA | 26.434 | 0.339 | cmaESplus |
23.687 | binGA | 24.736 | 0.557 | cmaESplus | 25.738 | cmaESplus | 26.539 | 0.210 | HC MS 10 |
Due to the generally better performance of the three reversible models, we examine alternative settings for the optimization procedures only for these models and subsequently apply the best setting found to each alternative irreversible model. The most promising settings are used to optimize the stochastic model as well.
Comparison of the performance of the optimization algorithms
Statistics on the most successful runs of each main optimizer
Model | reversible | irreversible | Algorithm | Population | ||||
---|---|---|---|---|---|---|---|---|
Min | Avg | Std. Dev. | Min | Avg | Std. Dev. | |||
GMAK | 20.326 | 20.742 | 0.501 | 25.745 | 34.472 | 15.285 | PSO linear 3, ϕ_{1} = ϕ_{2} = 2.05 | 25 |
20.403 | 21.787 | 1.297 | 25.183 | 31.694 | 13.255 | DE, f = 0.8, λ = 0.5, CR = 0.3 | 100 | |
21.975 | 23.812 | 1.604 | 24.741 | 49.045 | 21.285 | binGA, adaptive MUT, no CROSS | 250 | |
24.321 | 27.598 | 2.091 | 30.704 | 35.670 | 3.125 | cmaES | (5+25) | |
GMM | 20.280 | 22.818 | 2.186 | 24.857 | 27.978 | 9.146 | DE, f = 0.8, λ = 0.5, CR = 0.5 | 100 |
20.312 | 21.272 | 0.461 | 24.553 | 58.957 | 17.253 | PSO star, ϕ_{1} = ϕ_{2} = 2.05 | 25 | |
21.649 | 24.628 | 1.801 | 24.616 | 40.896 | 25.881 | binGA, adaptive MUT, one-point CROSS | 250 | |
23.890 | 26.624 | 1.851 | 26.414 | 31.136 | 9.052 | cmaES | (5+25) | |
CKMM | 20.100 | 21.434 | 0.563 | 21.511 | 26.077 | 7.729 | PSO grid3, ϕ_{1} = ϕ_{2} = 2.05 | 25 |
20.862 | 22.499 | 1.119 | 21.763 | 23.603 | 1.268 | DE, f = λ = 0.8, CR = 0.3 | 100 | |
21.431 | 23.222 | 1.066 | 25.632 | 28.055 | 2.697 | cmaES | (10, 50) | |
22.092 | 23.353 | 0.666 | 25.040 | 25.346 | 0.176 | binGA, adaptive MUT, no CROS | 100 |
Comparison of the modeling approaches
To evaluate which model is the most promising one we consider the following three criteria:
1. Fit of the model to the data
2. Number of model parameters
3. Computational time for simulation (computational complexity, weak criterion)
Each reversible deterministic model can be fitted to the data with a similar deviation from the measurements. The irreversible alternatives show a significantly higher deviation. Only the irreversible CKMM model is able to fit the data almost as well as the reversible models. However, most curves resulting from all of the irreversible models tend to become straight lines through the measurements, and thus behave in a biologically implausible manner [see Additional file 2]. This suggests that the irreversible models are unable to follow the dynamics of the system due to their fewer degrees of freedom. The rather abstract reversible models are able to tackle possible side effects of reactions not included in this reaction system and simulate them in terms of the reverse reaction. These models also consider the fact that, in biological systems, reaction products are normally not completely absent. Their concentration may be low, but they still take part in the reaction, in some cases giving a kind of feedback to the reactants [[19], pp. 312–313]. Therefore, the irreversible models are generally not competitive with respect to the data fit. From the three remaining reversible models, the CKMMr model achieves the best fit to the data, with 20.100, followed by GMMr model at 20.280 (worse by 0.180), and the GMAKr model at 20.326 (slightly worse by a further 0.046). The difference in the best model fit between the CKMMr and the GMAKr models is only 0.303. Hence, all three reversible models can be fitted to the data with a similarly small relative squared error (RSE, see Methods). A comparison of the best model fit 20.100 (CKMMr) to the fitness of the independently computed splines (19.670) evinces a difference of only 0.430. When considering the number of parameters (criterion 2), the GMAKr model shows a clear advantage with its 24 parameters compared to the 31 of the GMMr model or even 59 of the CKMMr model.
When choosing the parameter values for the GMAKr model completely randomly, this system can hardly be integrated without step size adaptation. Therefore, it is necessary to identify meaningful ranges of kinetic parameters within BRENDA [24–26]. A low-value initialization is necessary to assure numerically stable initial populations for the optimization procedures. For the other models this happens only if the parameters are chosen from implausibly large ranges.
The last criterion, the computation time, also depends on the complexity of the model but not necessarily on the number of parameters. The GMAKr model requires the smallest number of mathematical operations, followed by the GMMr model. The most complicated model is the CKMMr model. The average evaluation time over 100,000 repeats with randomly chosen parameters is 49 ms for the GMAKr model, 101 ms for the GMMr model and 151 ms for the CKMMr model. For hardware details see the Hardware Configuration section.
In order to take the effects of random fluctuations into account, one has to use stochastic models of the chemical reaction system. However, the most general approach, the chemical master equation [67], can hardly be solved numerically for larger systems [20]. Taking the number of parameters (criterion 2) and the computational costs (criterion 3) into account, the Langevin model is the most suitable formalism to consider the effects of random fluctuations while still providing an acceptable performance. Since the model can be stated in a way that allows it to be integrated with standard solvers for ordinary differential equations, the computational costs are of the same order of magnitude as the solution of the GMAKr model. However, care must be taken with respect to justification of the underlying simplifying assumptions [10]. The stochastic model equations of the biological system under consideration show no qualitatively different behavior in comparison to the deterministic model. The main reasons for this observation are found in the rather large molecule populations and the absence of points of instability in the allowed phase-space region.
Parameter space analysis
Once a model has been optimized several times and locally optimal parameter sets for the model are available, an analysis of the space of potential solutions becomes possible. This allows deducing characteristics of the solution space aiming to reduce the model complexity and enhance the optimization performance. If, for instance, two parameters show a linear dependency or correlation to each other, one of these can be removed from the model. Another interesting experiment would be to determine new ranges for each parameter. This can be done if a certain parameter varies only in a very small range compared to its maximal possible range.
For each of the three models, GMAKr, GMMr, and CKMMr, we gather all parameter values from all optimization runs that lead to a fitness less than 25. In this way, we obtain one parameter matrix for each model, in which each column corresponds to one optimization run and each row stands for one parameter. We conduct three analyses on the best parameters on each reversible model (see Methods):
1. Clustering to identify groups of similar ranges or almost constant values and to visualize the values of each parameter.
2. Variance analysis to visualize the scattering of each parameter.
3. Multiple correlation analysis aiming at finding highly correlated parameters, some of which can be eliminated.
For Histograms showing the parameter distribution of the GMAKr, GMMr, and CKMMr models [see Additional file 2].
Reduced reaction system and optimization results
Model | Replaced Parameter | Linear Regression Model | Mean Fitness ± Std. Dev. |
---|---|---|---|
GMAKr | k _{+1} | a_{1} · ${K}_{3}^{\text{I}}$ | 21.4782 ± 0.3254 |
k _{+3} | a_{2} · ${K}_{3}^{\text{I}}$ | ||
k _{-1} | a_{3} · ${K}_{3}^{\text{I}}$ | ||
k _{-3} | a_{4} · ${K}_{3}^{\text{I}}$ | ||
k _{+4} | a_{5} · k_{-4} | ||
k _{+7} | a_{6} · k_{-7} | ||
k _{+9} | a_{7} · K_{-9} | ||
GMMr | k _{-4} | b_{1} · k_{+4} + b_{2} · k_{+5} | 21.9538 ± 0.0832 |
k _{+1} | b_{3} · k_{-1} + b_{4} · ${K}_{1}^{\text{I}}$ | ||
k _{+7} | b_{5} · k_{-7} | ||
k _{+8} | b_{6} · k_{-8} | ||
k _{+9} | b_{7} · k_{-9} | ||
CKMMr | ${k}_{+1}^{\text{cat}}$ | c_{1} · ${K}_{{[\text{Pyr}]}_{1}}^{\text{M}}$ | 22.7219 ± 0.6626 |
${K}_{{[\text{DHIV}]}_{1}}^{\text{M}}$ | c_{2} · ${K}_{{[{\text{NADP}}^{\text{+}}]}_{1}}^{\text{M}}$ | ||
${k}_{\text{+3}}^{\text{cat}}$ | c_{3} · ${K}_{{\text{[KIV]}}_{\text{1}}}^{\text{M}}$ | ||
${K}_{{\text{[KIC]}}_{\text{2}}}^{\text{M}}$ | c_{4} · ${K}_{\text{[KIV]}}^{\text{M}}$ + c_{5} · ${K}_{{\text{[Glut]}}_{\text{2}}}^{\text{M}}$ |
Conclusion
The purpose of this study is to identify both the most suitable modeling approach and the best-performing optimization algorithm to calibrate the parameters contained in metabolic network models. To this end, we constructed one probabilistic and six deterministic mathematical descriptions of valine and leucine biosynthesis in C. glutamicum. The parameters of each model were optimized with respect to in vivo measurements for the reacting species within the system. In this way we compared eight optimization procedures. We systematically benchmarked both the algorithms and the alternative models to highlight their advantages and drawbacks. In the following paragraphs, we draw several conclusions from the comparison of these seven variants of a realistic reaction system, and we assume them to hold for similar systems. Thus, if no prior knowledge about a comparable metabolic system is available, our results can serve as a starting point for model construction and calibration.
Let us consider the capabilities of the modeling approaches in more detail, when taking into account the ability to approximate measured data, the hybrid model for the reversible system based on convenience rate laws and Michaelis-Menten equations (CKMMr) has the best performance. At the same time, this is the most complex model with respect to the number of parameters and computational costs for each simulation. The acceptable parameter values for this model, found by multiple optimization runs, varied over several orders of magnitude. This corresponds to the fact that the optimization problem shows a large number of local optima, which is often referred to as multimodal behavior. Furthermore, this model is integrable with parameter values selected by chance from an almost arbitrarily wide range. This means that no preliminary data analysis in enzyme databases is required to obtain an integrable start population for the CKMMr model.
On the other hand, a simplified deterministic description of the reaction system based on the generalized mass action rate law (GMAKr) yields good performance as well. Its advantage over the CKMMr model is its small number of parameters. Its major disadvantage is the strong tendency to become non-integrable when selecting parameter values by chance from a larger range. A restriction of the parameter space is required to ensure numerical stability when integrating metabolic network models. Some of the parameters showed almost no variance among the results of multiple optimization runs.
For models of biochemical systems with low metabolic concentrations or systems operating close to the point of instability, stochastic effects should be considered. Generally, simulating large stochastic models is computationally not feasible. The Langevin approach is a simplified stochastic description which facilitates taking these effects into account at acceptable computational costs. For the specific biochemical example network considered in this study, the stochastic effects are negligible as the behavior of the Langevin model is similar to the GMAKr model due to the rather large molecular populations. However, the benchmark showed that this approach is suitable for large-scale parameter optimization and model inference.
Modeling certain reactions in a non-reversible way as was done in all remaining models leads to a significantly worse ability to fit the measured data. We conclude that possible side effects are compensated by means of the reverse reaction. When modeling multi-enzyme systems all reactions should be treated reversibly, unless there is significant biological evidence to introduce irreversible rate laws. If neither the kinetic equations nor meaningful ranges of the parameter space are known, the model should be constructed using convenience kinetics. If the parameter space can be restricted using prior knowledge and the number of parameters matters, a model based on generalized mass-action rate laws constitutes an appropriate choice.
The second aim of this study is to identify the best-performing optimization algorithm for parameter estimation. The ability to find good local optima for the parameter values is the first quality measure for the algorithms. All five evolutionary algorithms tested yielded reasonable performance. From a user perspective, these algorithms differ in the number of settings which influence their behavior and are therefore more or less easy to apply. Hence, the effort to find a good configuration for an algorithm constitutes the second criterion of quality. The Tribes algorithm was among the best-performing algorithms in our benchmarks. As a settings-free optimization procedure, it is the most user-friendly method. However, other algorithms are able to yield even better results after fine-tuning. Particularly, DE and PSO provided the best performance while keeping the effort necessary for their fine-tuning within reasonable limits. ES and binGA are also able to identify valuable local optima for all systems but require examining a large number of well-established alternative operators for their crossover and mutation steps.
For first optimization attempts, the easy to use Tribes algorithm is a good choice. With slightly more effort, the user can adjust the algorithms PSO and DE to yield even better results. Combined with the convenience kinetics modeling approach, these algorithms provide a suitable choice to model unknown systems of metabolic reactions.
Methods
The biochemical example network
Figure 1 illustrates the biosynthesis of Val and Leu in a process diagram [68, 69] according to the META CYC[8] and KEGG[6, 7] databases. The pathway starts with pyruvate (Pyr), from which Val and Leu are produced. Both products are used for biomass production or can be transported out of the cell if not needed. It is important to note that this pathway is regulated by both products in six feedback inhibition mechanisms. The transport of Leu and Val across the cell wall is actually performed by the same enzyme, so that both substrates compete with each other. However, for modeling purposes two distinct reactions are necessary in which the competition is included as inhibition. Some reactions are lumped together (Table 1) as suggested by Magnus et al. [52]. Since the reaction 2 IPM) ⇌ 3 IPM is fast, it is assumed to be in equilibrium. This and the two following reactions 3 IPM + NAD^{+} → 2 I_{3}OS + NADH_{2} as well as (2S)-2-isopropyl-3-oxosuccinate (2 I_{3}OS) → 2-ketoisocaproate (KIC) +CO_{2} that only depend on the concentration of 2 IPM are lumped together, introducing the symbol IPM for both derivatives.
Glucose stimulus-response experiment
Steady-state concentrations of the reacting species
Metabolite | Concentration (mM) |
---|---|
AcLac | 0.236 |
DHIV | 0.132 |
IPM | 0.0227 |
KIC | 0.0741 |
Leu | 0.209 |
Val | 29.4 |
Mathematical models
For the resulting seven-dimensional differential equation system [see Additional file 2]. The databases KEGG[6, 7] and META CYC[8] do not indicate if the reaction network pictured in Figure 1 contains irreversible reactions besides the draining of both products as listed in Table 1. One way to study the influence of the existence or non-existence of reverse reactions on the dynamics of the whole system is to derive alternative models and investigate their ability to approximate the data. The simpler irreversible reactions are favored if they are able to fit the data. The following paragraphs present the general equations for all rate laws, which are inserted for each v_{ i }and investigated in this study. For the resulting differential equation systems [see Additional file 2]. An implementation of the fourth-order Runge-Kutta method [70] solves the ordinary differential equation systems. The stochastic Langevin system is adapted to be integrated using the MAT LAB™ integrator ode15s [71, 72].
Generalized Mass Action Kinetics (GMAKr)
with ${K}_{j}^{\text{I}}\u2a7e0$. While Equation (4) is derived intuitively, driven by the assumption that the exponential function constitutes an important growth and shrinkage function in biology, Equation (3) can be derived from the competing reactions of the enzyme with its substrate or inhibitor. The first equation can be derived using the equilibrium constant for the inhibition reaction and the conservation law of the enzyme as well as the enzyme-inhibitor complex concentrations. Applying Equation (2) combined with Equation (3) to reaction system R_{1} through R_{10} leads to an Ordinary Differential Equation (ODE) system with 24 parameters k_{± j}, ${K}_{j}^{\text{I}}$.
Irreversible GMAK with exp inhibition (GMAKi)
By setting all product concentrations apart from R_{2} and R_{9} to zero and applying Equation (4) to Equation (2) we obtain the irreversible version of this equation system with 18 parameters k_{± j}, ${K}_{j}^{\text{I}}$.
Michaelis-Menten equations (GMMr)
In the case of R_{3} there might be a reverse reaction. R_{6} and R_{10} are assumed to be irreversible because they describe the transport of Val and Leu out of the cell. We further assume that the constants v^{m} in R_{6} and R_{10} are allowed to be zero so that there is no need to export Val or Leu if it is needed for biomass formation. All other reactions are modeled using the GMAKr approach including Equation (3) for inhibition. The complete GMMr model contains 31 parameters to be estimated. To avoid numerical problems, the inhibition constants in Michaelis-Menten kinetics are transformed into their reciprocals ${K}^{\text{I}a|b}\text{'}=\frac{1}{{K}^{\text{I}a|b}}$. This modification allows the optimization procedure to "decide" which kind of inhibition occurs [5] [see Additional file 2].
Irreversible Michaelis-Menten Model (GMMi)
From the GMMr model an irreversible alternative is established by setting all product concentrations to zero. The resulting system contains 24 parameters ${K}_{j}^{\text{Ia}|\text{b}}$, k_{± j}, ${K}_{ij}^{\text{M}}$ to be estimated.
Reversible Convenience Kinetics (CKMMr)
has been suggested and herein applied. Besides the reciprocal constant this approach equals Equation (3). The product [E_{ j }]${k}_{\pm j}^{\text{cat}}$ is lumped into one parameter ${V}_{\pm j}^{\text{m}}$ for all j assuming that all enzyme concentrations remain constant during the 25 s. No enzyme concentrations have been measured, so that an optimizer cannot distinguish between the product of two parameters or one parameter. The three reactions that follow the Michaelis-Menten mechanism are modeled using Equation (5). The reactions R_{6} and R_{10} are considered irreversible as described before. Applying Equation (6) to all remaining reactions in the system R_{1} through R_{10} yields an equation system with 59 parameters. The stoichiometric matrix has full column rank. Hence, the parameters ${k}_{\pm j}^{\text{cat}}$ can be estimated directly without violating thermodynamic constraints [16].
Irreversible Convenience Kinetics (CKMMi)
By setting all product concentrations apart from R_{2} and R_{9} to zero, we obtain an irreversible version of this model containing 41 parameters.
Stochastic modeling based on the Langevin equation (LANG)
allowing adaptive step-size control of the ODE solver. The reaction propensities are calculated according to Gillespie [11]. This leads to an equation system with 24 parameters.
Representing external metabolites with splines
As suggested by Magnus et al. [52], metabolites whose concentrations cannot be explained in terms of the model are considered external, i. e., they are an input to the model but involved in numerous other reactions that are not reflected by this system (Figure 1). We approximate these using cubic splines which smooth the measurements. To weight all measurements equally, all ω_{ i }are set to 1. Due to the different ranges of the concentrations of the six metabolites, it is not possible to find one appropriate degree of smoothness λ that leads to equally smooth curves. Hence, we transform all concentrations into the range [0, 1], set λ = 1, compute the spline coefficients and re-transform them back into the original range of the specific metabolite. For details and a figure showing the resulting splines [see Additional file 2].
Fitness function and search-space restrictions
This fitness function was used in several publications for similar problems [33, 47, 48].
To restrict the search space for the optimizers, we limit parameter values to the range [0, 2000], covering 98.748% of all known kinetic parameters in BRENDA [24–26], as suggested in [47, 48]. All known parameters in BRENDA are greater than or equal to zero. To avoid division by zero in some parameters, the range is set to [ε, 2000] with ε = 10^{-8}. For parameters transformed in Michaelis-Menten equations (as described above), the range is limited to [0, 10^{8}], resulting in a search space from 10^{-8} through ∞. Only 0.962% of all K^{I} and 0.004% of all K^{M}values in BRENDA are reported to be lower than 10^{-8}. This ε is chosen to guarantee numerical stability. These restrictions are applied to all parameters. In cases where no division by the parameter value is necessary, ε = 0 is allowed to be a lower bound of the parameter range.
All parameters are initialized with low values to avoid obtaining unstable initial populations and it is assumed that large parameter values are rather infrequent in nature [47, 48]. A Gaussian distribution with μ = σ = 1 guarantees low initial values and ensures stable initial populations. Each parameter is set to the boundary values if it breaks any of the search-space restrictions. We limit the initialization procedure for all models to a low value initialization to ensure equal conditions in all comparisons.
Systematically improving the performance of the optimization procedures
Standard settings for the optimization algorithms
For a comprehensive introduction to all optimization algorithms used in this study [see Additional file 1], in which we also explain the specific settings and operators for each method in detail.
Using the open-source framework EVA2 for nature-inspired optimization procedures [37, 38], we test the following standard settings of the algorithms (Table 2) on the inference problem, of which the following are evolutionary optimization procedures:
• Binary Genetic Algorithm (binGA) with one-point mutation, p_{ m }= 0.1, and one-point crossover, p_{ c }= 0.7.
• Real-valued Genetic Algorithm (realGA) with global mutation, p_{ m }= 0.1, and UNDX crossover, p_{ c }= 0.8. Both GAs use tournament selection with a group size of 8 in a population of 250 individuals.
• Standard Evolution Strategy (stdES) as (5,25)-ES with global mutation, p_{ m }= 0.8, and discrete one-point crossover, p_{ c }= 0.2.
• Evolution Strategy with covariance matrix adaption with and without elitism, i. e., "plus" strategy, (cmaES/cmaESplus) with μ = 5 and λ = 25, p_{ m }= 1.0, and no crossover. All ESs use deterministic best-first selection to choose the next generation.
• Differential Evolution (DE) with the scheme DE/current-to-best/1 [42] setting λ = F = 0.8, CR = 0.5 and a population size of 100.
Two optimization strategies are swarm intelligence optimization procedures:
• Constricted Particle Swarm Optimization (PSO) [45, 76]: setting ϕ_{1} = 2.05, ϕ_{2} = 2.05, χ = 0.73 using star topology and a population size of 100 as well as its derivative
• Tribes [46]
We also study the following classical non-evolutionary methods:
• Monte Carlo Optimization (MCO) with 50 multi-runs
• (Multi-start) Hill Climber (HC), the number of starts varying from 1, 10, 25, 50, 100 to 250 using Gaussian mutation with a fixed standard deviation of σ = 0.2 and a mutation probability p_{ m }= 1.0.
• Simulated Annealing (SA) with α = 0.1 and an initial temperature of T = 5 using a linear annealing schedule and a population size of 250.
For all algorithms with population sizes lower than 250 individuals, a pre-population with 250 parameter vectors is generated and the best are selected to create the initial population. This step is crucial to obtain comparable results for algorithms with different population sizes [33]. Every setting described in this and all following paragraphs is repeated 20 times with 100,000 fitness evaluations per run on the deterministic models. The results of this preliminary trial are shown in Table 3 and Figure 3. The methods found to be successful in these analyses are utilized to optimize the Langevin model as well.
Alternative settings for the best optimization algorithms
We first study the influence of various mutation and crossover operators on the three deterministic reversible models for binGA and ES. In a grid search, all combinations of the following operators are evaluated, and we exclude the combinations of no crossover with no mutation. For binGA, no mutation, one-point, and adaptive mutation, an operator which modifies individual mutation probabilities similar to ES step-size adaptation, is tested, paired with one- and n-point (n = 3), uniform, and bit-simulated crossover, each with p_{ m }= 0.1, p_{ c }= 0.7. On ES we use the mutation operators covariance matrix adaptation (CMA), the 1/5^{th} success rule as well as correlated, global, local, and no mutation paired with the crossover operators one- and n-point (n = 3), UNDX, and no crossover, each with p_{ m }= 0.8, p_{ c }= 0.2. To study the influence of the mutation or the crossover operator alone, the values for p_{ m }and p_{ c }are set to zero or one, depending on which influence is investigated. The population size for binGA is set to 100 and for ES (5, 25) (Figure 4). For the most successful operator combination (adaptive mutation with bit-simulated crossover) on the GMAKr model, we study the influence of the probabilities p_{ m }and p_{ c }, with which the respective operator is invoked by binGA. All pairs of pm and pc are evaluated from 0.0 through 1.0 in 0.1 steps and a population size of 100, excluding p_{ m }= p_{ c }= 0.0 (Figure 5(a)). Subsequently, the population sizes 50, 100, 250, 500, 1,000, and 2,000 are tested for binGA with p_{ m }= 0.2 and p_{ c }= 1.0 (Figure 5(b)). For the cmaESplus, we evaluate all combinations of μ ∈ {5, 10, 25, 50, 75}, and λ ∈ {10, 25, 50, 75, 100, 125, 150}, excluding cases where μ > λ and keeping p_{ m }at 1.0, and p_{ c }at 0.0 (Figure 6).
For the DE approach, another grid search is performed on the three reversible deterministic models, altogether testing values for F, λ ∈ {0.5, 0.8}, and CR ∈ {0.3, 0.5, 0.9} (Figure 7). For the most promising parameter set (F = 0.8, λ = 0.5, and CR = 0.3 for GMAKr and GMMr, and CR = 0.9 for the CKMMr), the population size is additionally varied over {50, 250, 500, 1000} for each model (Figure 8).
The PSO with star topology and standard settings for ϕ_{1} = ϕ_{2} = 2.05 is compared with a star topology and settings ϕ_{1} = 2.8, ϕ_{2} = 1.3, as suggested in [77] as well as a linear 3 and grid 3 topology with standard parameters on all three reversible deterministic models (Figure 9). The population size is set to 25. Additional population sizes, namely {50, 250, 500}, are evaluated using a grid 3 topology and standard values for ϕ_{1} and ϕ_{2} on these models (Figure 10).
Parameter-space analysis
The completion of the large-scale parameter optimization study in the first part yields a large set of different optimal parameters. We select all parameter sets with a fitness less than 25 for each reversible deterministic model (GMAKr, GMMr, and CKMMr).
Cluster analysis
To cluster these best parameter sets (Figure 12), we apply the agglomerative nesting algorithm (AGNES) [[78], pp. 199–251], using a Euclidean metric, which is implemented in the cluster package [79] of the R-project [80].
Variance analysis to visualize the scattering of each parameter
The variances of each parameter among the best optimization results of the three reversible deterministic models are calculated and plotted on a logarithmic scale due to the large differences in their orders of magnitude (Figure 13).
Multiple correlation analysis
measures how well the paremeter Y can be explained in a linear sense by the other model parameters X_{1},...,X_{ p }. After identification of highly correlated parameters, a subset of replacement candidates is selected. A parameter X is considered a replacement candidate for parameter Y if r_{X,Y}≥ 0.7 and the p-value of the correlation, computed using a t-test, is close to zero. The degrees of freedom selected for replacement are subsequently substituted by a linear regression model of their correlated parameters.
Hardware configuration
All experiments are run on a cluster with 16 AMD dual Opteron CPUs with 2.4 GHz, 1 MB level 2 cache, and 2 GB RAM per node under the Sun Grid Engine, and JVM 1.5.0 with Scientific Linux 5 as operating system. An optimization of one model in 20 parallel multi-runs takes approximately 1.5 h.
Availability of models and optimization procedures
All models investigated in this study are included as optimization problems in EvA2, a Java™-based workbench for heuristic optimization [37, 38], which can be downloaded at http://www.ra.cs.uni-tuebingen.de/software/EvA2.
Declarations
Acknowledgements
We are grateful to Ralf Takors and Klaus Beyreuther. This work was funded by the German Federal Ministry of Education and Research (BMBF) in the National Genome Research Network (NGFN-II EP) under Project Number 0313323 and the German federal state Baden-Württemberg in the project "Identifikation und Analyse metabolischer Netze aus experimentellen Daten" under Project Number 7532.22-26-18. Conflict of Interest: none declared.
Authors’ Affiliations
References
- Kitano H: Computational systems biology. Nature. 2002, 420 (6912): 206-210.View ArticlePubMedGoogle Scholar
- Albert R: Network Inference, Analysis, and Modeling in Systems Biology. Plant Cell. 2007, 19 (11): 3327-3338.PubMed CentralView ArticlePubMedGoogle Scholar
- Gombert AK, Nielsen J: Mathematical modelling of metabolism. Current Opinion in Biotechnology. 2000, 11 (2): 180-186.View ArticlePubMedGoogle Scholar
- Covert MW, Schilling CH, Famili I, Edwards JS, Goryanin II, Selkov E, Palsson BO: Metabolic modeling of microbial strains in silico. Trends in Biochemical Sciences. 2001, 26 (3): 179-186.View ArticlePubMedGoogle Scholar
- Heinrich R, Schuster S: The Regulation of Cellular Systems. 115 Fifth Avenue New York, NY 10003. 1996, Chapman and HallGoogle Scholar
- Kanehisa M, Goto S, Hattori M, Aoki-Kinoshita KF, Itoh M, Kawashima S, Katayama T, Araki M, Hirakawa M: From genomics to chemical genomics: new developments in KEGG. Nucl Acids Res. 2006, 34: D354-357.PubMed CentralView ArticlePubMedGoogle Scholar
- Ogata H, Goto S, Sato K, Fujibuchi W, Bono H, Kanehisa M: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Research. 2000, 27: 29-34. http://nar.oxfordjournals.org/cgi/content/abstract/27/1/29View ArticleGoogle Scholar
- Caspi R, Foerster H, Fulcher CA, Hopkinson R, Ingraham J, Kaipa P, Krummenacker M, Paley S, Pick J, Rhee SY, Tissier C, Zhang P, Karp PD: MetaCyc: a multiorganism database of metabolic pathways and enzymes. Nucleic Acids Research. 2006, D511-D516. 34 DatabasePubMed CentralView ArticlePubMedGoogle Scholar
- Bulik S, Grimbs S, Selbig J, Holzhütter HG: Combining mechanistic and simplified enzymatic rate equations: A promising approach for speeding up the kinetic modeling of complex metabolic networks. FEBS Journal. 2009, 276: 410-524.View ArticlePubMedGoogle Scholar
- Gillespie DT: The chemical Langevin equation. Journal of Chemical Physics. 2000, 113: 297-306.http://link.aip.org/link/?JCPSA6/113/297/1View ArticleGoogle Scholar
- Gillespie DT: A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. Journal of Computational Physics. 1976, 22 (4): 403-434. http://www.sciencedirect.com/science/article/B6WHY-4DD1NC9-CP/2/43ade5f11fb949602b3a2abdbbb29f0eView ArticleGoogle Scholar
- Visser D, Heijnen J: The Mathematics of Metabolic Control Analysis Revisited. Metabolic Engineering. 2002, 4: 114-123. http://www.sciencedirect.com/science/article/B6WN3-45V802C-3/2/d624a20d0e70ca2a1058359d7fd00cb0View ArticlePubMedGoogle Scholar
- Hatzimanikatis V, Floudas CA, Bailey JE: Analysis and design of metabolic reaction networks via mixed-integer linear optimization. AIChE. 1996, 42 (5): 1277-1292.http://cat.inist.fr/?aModele=afficheN&cpsidt=3105838View ArticleGoogle Scholar
- Hatzimanikatis V, Bailey JE: Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log)linear kinetic models. Biotechnology and Bioengineering. 1997, 54 (2): 91-104. http://www3.interscience.wiley.com/journal/71003853/abstractView ArticlePubMedGoogle Scholar
- Hatzimanikatis V, Emmerling M, Sauer U, Bailey JE: Application of mathematical tools for metabolic design of microbial ethanol production. Biotechnology and Bioengineering. 1998, 58 (2): 154-161. http://www3.interscience.wiley.com/journal/71002326/abstractView ArticlePubMedGoogle Scholar
- Liebermeister W, Klipp E: Bringing metabolic networks to life: convenience rate law and thermodynamic constraints. Theor Biol Med Model. 2006, 3: 41-PubMed CentralView ArticlePubMedGoogle Scholar
- Segel IH: Enzyme Kinetics – Behavior and Analysis of Rapid Equilibrium and Steady-State Enzyme Systems. 1993, Wiley Classics Library EditionGoogle Scholar
- Bisswanger H: Enzymkinetik – Theorie und Methoden. 2000, Weinheim, Germany: Wiley-VCH, 3View ArticleGoogle Scholar
- Cornish-Bowden A: Fundamentals of Enzyme Kinetics. 2004, 59 Portland Place, London: Portland Press Ltd, 3Google Scholar
- Gillespie DT: Stochastic simulation of chemical kinetics. Annu Rev Phys Chem. 2007, 58: 35-55.View ArticlePubMedGoogle Scholar
- Savageau MA: Biochemical systems analysis. 3. Dynamic solutions using a power-law approximation. J Theor Biol. 1970, 26 (2): 215-226.View ArticlePubMedGoogle Scholar
- Visser D, Heijnen JJ: Dynamic simulation and metabolic re-design of a branched pathway using linlog kinetics. Metab Eng. 2003, 5 (3): 164-176.View ArticlePubMedGoogle Scholar
- Hatzimanikatis V, Bailey J: MCA Has More to Say. Journal of theoretical Biology. 1996, 233-342.Google Scholar
- Barthelmes J, Ebeling C, Chang A, Schomburg I, Schomburg D: BRENDA, AMENDA and FRENDA: the enzyme information system in 2007. Nucl Acids Res. 2007, 35 (suppl_1): D511-514. http://nar.oxfordjournals.org/cgi/content/abstract/35/suppl_1/D511PubMed CentralView ArticlePubMedGoogle Scholar
- Schomburg I, Chang A, Schomburg D: BRENDA, enzyme data and metabolic information. Nucl Acids Res. 2002, 30: 47-49.PubMed CentralView ArticlePubMedGoogle Scholar
- Schomburg I, Chang A, Ebeling C, Gremse M, Heldt C, Huhn G, Schomburg D: BRENDA, the enzyme database: updates and major new developments. Nucleic Acids Research. 2004, 32 (Database Issue): D431-433.PubMed CentralView ArticlePubMedGoogle Scholar
- Metzler DE: Biochemistry. 2001, Harcourt/Academic PressGoogle Scholar
- Banga JR: Optimization in computational systems biology. BMC Systems Biology. 2008, 2: 47-PubMed CentralView ArticlePubMedGoogle Scholar
- Rodriguez-Fernandez MR, Mendes P, Banga JR: A hybrid approach for efficient and robust parameter estimation in biochemical pathways. Biosystems. 2006, 83: 248-265.View ArticlePubMedGoogle Scholar
- Rodriguez-Fernandez M, Egea JA, Banga JR: Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems. BMC Bioinformatics. 2006, 7: 483-PubMed CentralView ArticlePubMedGoogle Scholar
- Balsa-Canto E, Peifer M, Banga JR, Timmer J, Fleck C: Hybrid optimization method with general switching strategy for parameter estimation. BMC Systems Biology. 2008, 2: 26-PubMed CentralView ArticlePubMedGoogle Scholar
- Spieth C, Streichert F, Speer N, Zell A: Optimizing Topology and Parameters of Gene Regulatory Network Models from Time-Series Experiments. Proceedings of the Genetic and Evolutionary Computation Conference (GECCO 2004), (Part I) of LNCS. 2004, 3102: 461-470.http://www.springerlink.com/content/cx9mmtkl2ca0fcx4/View ArticleGoogle Scholar
- Spieth C, Worzischek R, Streichert F, Supper J, Speer N, Zell A: Comparing Evolutionary Algorithms on the Problem of Network Inference. Proceedings of the Genetic and Evolutionary Computation Conference (GECCO 2006). 2006,http://portal.acm.org/citation.cfm?id=1143997.1144052Google Scholar
- Spieth C, Hassis N, Streichert F, Supper J, Beyreuther K, Zell A: Comparing Mathematical Models on the Problem of Network Inference. Proceedings of the Genetic and Evolutionary Computation Conference (GECCO 2006). 2006,http://portal.acm.org/citation.cfm?id=1143997.1144045Google Scholar
- Weise T: Global Optimization Algorithms – Theory and Application. 2008, Thomas Weise, 2 http://www.it-weise.de/Google Scholar
- Charbonneau P, Knapp B: A User's Guide to PIKAIA 1.0. 1995, http://download.hao.ucar.edu/archive/pikaia/userguide.psGoogle Scholar
- Streichert F, Ulmer H: JavaEvA – A Java Framework for Evolutionary Algorithms. Technical Report WSI-2005–06, Center for Bioinformatics Tübingen. 2005, University of Tübingen,http://w210.ub.uni-tuebingen.de/dbt/volltexte/2005/1702/Google Scholar
- Kronfeld M: EvA2 Short Documentation. 2008, University of Tübingen, Dept. of Computer Architecture, Sand 1, 72076 Tübingen,http://www.ra.cs.uni-tuebingen.de/software/EvA2Google Scholar
- Mathtools.net : MATLAB/Optimization.http://www.mathtools.net/MATLAB/Optimization
- Holland JH: Adaptation in Natural and Artificial Systems. 1975, The University of Michigan PressGoogle Scholar
- Rechenberg I: Evolutionsstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution. 1973, Fromman-Holzboog, StuttgartGoogle Scholar
- Storn R: On the Usage of Differential Evolution for Function Optimization. 1996 Biennial Conference of the North American Fuzzy Information Processing Society. 1996, 519-523. Berkeley: IEEE, New York, USA,http://www.icsi.berkeley.edu/~storn/bisc1.ps.gzGoogle Scholar
- Tovey CA: Hill climbing with multiple local optima. Alg Disc Meth. 1985, 6 (3): 384-393.http://link.aip.org/link/?SML/6/384/1View ArticleGoogle Scholar
- Kirkpatrick S, Gelatt CD, Vecchi MP: Optimization by Simulated Annealing. Science. 1983, 220 (4598): 671-680. http://www.sciencemag.org/cgi/content/abstract/220/4598/671View ArticlePubMedGoogle Scholar
- Clerc M, Kennedy J: The Particle Swarm – Explosion, Stability, and Convergence in a Multidimensional Complex Space. IEEE Transactions on Evolutionary Computation. 2002, 6: 58-73.http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=985692View ArticleGoogle Scholar
- Clerc M: Particle Swarm Optimization. 2005, ISTE LtdGoogle Scholar
- Dräger A, Supper J, Planatscher H, Magnus JB, Oldiges M, Zell A: Comparing Various Evolutionary Algorithms on the Parameter Optimization of the Valine and Leucine Biosynthesis in Corynebacterium glutamicum. 2007 IEEE Congress on Evolutionary Computation. Edited by: Srinivasan D, Wang L. 2007, 620-627. IEEE Computational Intelligence Society, Singapore: IEEE Press,http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4424528View ArticleGoogle Scholar
- Dräger A, Kronfeld M, Supper J, Planatscher H, Magnus JB, Oldiges M, Zell A: Benchmarking Evolutionary Algorithms on Convenience Kinetics Models of the Valine and Leucine Biosynthesis in C. glutamicum. 2007 IEEE Congress on Evolutionary Computation. Edited by: Srinivasan D, Wang L. 2007, 896-903. IEEE Computational Intelligence Society, Singapore: IEEE Press,http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4424565View ArticleGoogle Scholar
- Chassagnole C, Noisommit-Rizzi N, Schmid JW, Mauch K, Reuss M: Dynamic Modeling of the Central Carbon Metabolism of Escherichia coli. 2002, 54-73. Wiley Periodicals, Inc,http://www3.interscience.wiley.com/journal/93519745/abstractGoogle Scholar
- Klipp E, Nordlander B, Kruger R, Gennemark P, Hohmann S: Integrative model of the response of yeast to osmotic shock. Nature Biotechnology. 2005, 23 (8): 975-982.View ArticlePubMedGoogle Scholar
- Guthke R, Schmidt-Heck W, Pless G, Gebhardt R, Pfaff M, Gerlach JC, Zeilinger K: Dynamic Model of Amino Acid and Carbohydrate Metabolism in Primary Human Liver Cells. VII International Symposium on Biological and Medical Data Analysis. 2006,http://www.springerlink.com/content/f78167p04n1426w6/Google Scholar
- Magnus JB, Hollwedel D, Oldiges M, Takors R: Monitoring and Modeling of the Reaction Dynamics in the Valine/Leucine Synthesis Pathway in Corynebacterium glutamicum. Biotechnology Progress. 2006, 22 (4): 1071-1083.View ArticlePubMedGoogle Scholar
- Guthke R, Zeilinger K, Sickinger S, Schmidt-Heck W, Buentemeyer H, Iding K, Lehmann J, Pfaff M, Pless G, Gerlach JC: Dynamics of amino acid metabolism of primary human liver cells in 3D bioreactors. Bioprocess Biosystem Engineering. 2006, 28 (5): 331-340.View ArticleGoogle Scholar
- Savageau MA: Biochemical systems analysis. I. Some mathematical properties of the rate law for the component enzymatic reactions. J Theor Biol. 1969, 25 (3): 365-369.View ArticlePubMedGoogle Scholar
- Savageau MA: Biochemical systems analysis. II. The steady-state solutions for an n-pool system using a power-law approximation. J Theor Biol. 1969, 25 (3): 370-379.View ArticlePubMedGoogle Scholar
- Liu PK, Wang FS: Inference of Biochemical Network Models in S-System Using Multi-Objective Optimization Approach. Bioinformatics. 2008, 24 (8): 1085-1093.View ArticlePubMedGoogle Scholar
- Vera J, Sun C, Oertel Y, Wolkenhauer O: PLMaddon: a power-law module for the Matlab™ SBToolbox. Bioinformatics. 2007, 23 (19): 2638-2640.View ArticlePubMedGoogle Scholar
- Shapiro BE, Levchenko A, Meyerowitz EM, Wold BJ, Mjolsness ED: Cellerator: extending a computer algebra system to include biochemical arrows for signal transduction simulations. Bioinformatics. 2002, 19 (5): 677-678.View ArticleGoogle Scholar
- Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI-a COmplex PAthway SImulator. Bioinformatics. 2006, 22 (24): 3067-3074.View ArticlePubMedGoogle Scholar
- Dräger A, Hassis N, Supper J, Schröder A, Zell A: SBMLsqueezer: a CellDesigner plug-in to generate kinetic rate equations for biochemical networks. BMC Systems Biology. 2008, 2: 39-PubMed CentralView ArticlePubMedGoogle Scholar
- Yang CR, Shapiro BE, Mjolsness ED, Hatfield GW: An enzyme mechanism language for the mathematical modeling of metabolic pathways. Bioinformatics. 2004, 21 (6): 774-780.View ArticlePubMedGoogle Scholar
- Borger S, Liebermeister W, Uhlendorf J, Klipp E: Automatically generated model of a metabolic network. International Conference on Genome Informatics. 2007, 18: 215-224. http://eproceedings.worldscinet.com/9781860949920/9781860949920_0021.htmlGoogle Scholar
- Jamshidi N, Palsson BO: Formulating genome-scale kinetic models in the post-genome era. Molecular Systems Biology. 2008, 4:Google Scholar
- Voss HU, Timmer J, Kurths J: Nonlinear Dynamical System Identification from Uncertain and Indirect Measurements. International Journal of Bifurcation and Chaos. 2004, 14 (6): 1905-1933. http://www.worldscinet.com/cgi-bin/details.cgi?id=pii:S0218127404010345&type=htmlView ArticleGoogle Scholar
- Eggeling L, Bott M: Handbook of Corynebacterium glutamicum. 2005, Boca Raton: Taylor & FrancisView ArticleGoogle Scholar
- Hansen N, Ostermeier A: Completely Derandomized Self-Adaptation in Evolution Strategies. Evolutionary Computation. 2001, 9 (2): 159-195.View ArticlePubMedGoogle Scholar
- Gillespie D: A rigorous derivation of the chemical master equation. Physica A. 1992, 188: 404-425. http://www.comp.nus.edu.sg/~cs6280/Materials/06-gillespie92.pdfView ArticleGoogle Scholar
- Kitano H, Funahashi A, Matsuoka Y, Oda K: Using process diagrams for the graphical representation of biological networks. Nature Biotechnology. 2005, 23 (8): 961-966.View ArticlePubMedGoogle Scholar
- Funahashi A, Tanimura N, Morohashi M, Kitano H: CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BioSilico. 2003, 1 (5): 159-162.http://www.sciencedirect.com/science/article/B75GS-4BS08JD-5/2/5531c80ca62a425f55d224b8a0d3f702View ArticleGoogle Scholar
- Spieth C, Supper J, Streichert F, Speer N, Zell A: JCell – a Java-based framework for inferring regulatory networks from time series data. Bioinformatics. 2006, 22 (16): 2051-2052. http://bioinformatics.oxfordjournals.org/cgi/content/abstract/22/16/2051View ArticlePubMedGoogle Scholar
- Shampine LF, Reichelt MW: The MATLAB ODE Suite. Tech rep. 2007, http://www.mathworks.com/access/helpdesk/help/pdf_doc/otherdocs/ode_suite.pdfGoogle Scholar
- Shampine LF, Reichelt MW, Kierzenka JA: Solving Index-1 DAEs in MATLAB and Simulink. SIAM Rev. 1999, 41 (3): 538-552.http://www.mathworks.com/support/solutions/files/s8314/dae.psView ArticleGoogle Scholar
- Schauer M, Heinrich R: Quasi-steady-state approximation in the mathematical modeling of biochemical reaction networks. Mathematical Bioscience. 1983, 65: 155-171.http://cat.inist.fr/?aModele=afficheN&cpsidt=9308909View ArticleGoogle Scholar
- Kloeden PE, Platen E: Numerical Solution of Stochastic Differential Equations. 1992, Applications of Mathematics, Berlin: Springer-VerlagView ArticleGoogle Scholar
- Bentele M: Stochastic Simulation and System Identification of large Signal Transduction Networks in Cells. PhD thesis. 2004, Combined Faculties for the Natural Sciences and for Mathematics of the Ruperto-Carola University Heidelberg, Germany, http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=488968Google Scholar
- Kennedy J, Eberhart R: Particle Swarm Optimization. IEEE Int Conf on Neural Networks, Perth, Australia. 1995Google Scholar
- Carlisle A, Dozier G: An off-the-shelf PSO. Proceedings of the Workshop on Particle Swarm Optimization. 2001, Indianapolis: Purdue School of Engineering and Technology, IndianapolisGoogle Scholar
- Kaufman L, Rousseeuv PJ: Finding Groups in Data: An Introduction to Cluster Analysis. 1990, Probability and Mathematical Statistics, New York: John Wiley and Sons, IncView ArticleGoogle Scholar
- Maechler M, Rousseeuw P, Struyf A, Hubert M: Cluster Analysis Basics and Extensions. 2005Google Scholar
- R Development Core Team : R: A Language and Environment for Statistical Computing. 2007, , R Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.orgGoogle Scholar
- Hartung J: Statistik. 2002, München: Oldenburg Wissenschaftsverlag GmbHView 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.