Abstract
Background
An efficient and reliable parameter estimation method is essential for the creation of biological models using ordinary differential equation (ODE). Most of the existing estimation methods involve finding the global minimum of data fitting residuals over the entire parameter space simultaneously. Unfortunately, the associated computational requirement often becomes prohibitively high due to the large number of parameters and the lack of complete parameter identifiability (i.e. not all parameters can be uniquely identified).
Results
In this work, an incremental approach was applied to the parameter estimation of ODE models from concentration time profiles. Particularly, the method was developed to address a commonly encountered circumstance in the modeling of metabolic networks, where the number of metabolic fluxes (reaction rates) exceeds that of metabolites (chemical species). Here, the minimization of model residuals was performed over a subset of the parameter space that is associated with the degrees of freedom in the dynamic flux estimation from the concentration timeslopes. The efficacy of this method was demonstrated using two generalized mass action (GMA) models, where the method significantly outperformed singlestep estimations. In addition, an extension of the estimation method to handle missing data is also presented.
Conclusions
The proposed incremental estimation method is able to tackle the issue on the lack of complete parameter identifiability and to significantly reduce the computational efforts in estimating model parameters, which will facilitate kinetic modeling of genomescale cellular metabolism in the future.
Keywords:
Incremental parameter estimation; Kinetic modeling; Metabolic network; GMA modelBackground
The estimation of unknown kinetic parameters from timeseries measurements of biological molecules is a major bottleneck in the ODE model building process in systems biology and metabolic engineering [1]. The majority of current estimation methods involve simultaneous (singlestep) parameter identification, where model prediction errors are minimized over the entire parameter space. These methods often rely on global optimization methods, such as simulated annealing, genetic algorithms and other evolutionary approaches [13]. The problem of obtaining the bestfit parameter estimates however, is typically illposed due to issues related with data informativeness, problem formulation and parameter correlation, all of which contribute to the lack of complete parameter identifiability. Not to mention, finding the global minimum of model residuals over highly multidimensional parameter space is challenging and can become prohibitively expensive to perform on a computer workstation, even for tens of parameters.
Here, we consider the modeling of cellular metabolism using the canonical powerlaw formalism, specifically the generalized mass action (GMA) systems [4,5]. The powerlaw formalism has many advantages, which have been detailed elsewhere [1,6]. Notably, power laws have a relatively simple structure that permits algebraic manipulation in the logarithmic scale, but nonetheless is capable of describing essentially any nonlinearity. Regulatory interactions among metabolites can also be described straightforwardly through the kinetic order parameters, establishing an equivalence between structural identification and parametric estimation. However, the number of parameters increases proportionally with the number of metabolites and fluxes, leading to a largescale parameter identification problem, one where singlestep estimation methods often struggle to converge.
The integration of ODE often constitutes a major part of the computational cost in the parameter estimation, especially when the ODE model is stiff [7]. While stiffness can genuinely arise due to a large time scale separation of the reaction kinetics in the real system, stiff ODEs could also result from unrealistic combinations of parameter values during the parameter optimization procedure, especially when a global optimizer is used. The parameter estimation of ODE models using powerlaw kinetics is particularly prone to stiffness problem since many of the unknown parameters are the exponents of the concentrations. For this reason, alternative formulations have been proposed that avoid these ODE integrations either completely [7,8] or partially [911]. Particularly, computational cost could be significantly reduced by decomposing the estimation problem into two phases, starting with the calculation of dynamic reaction rates or fluxes from the slopes of concentration data, followed by the least square regressions of kinetic parameters [1214]. In this case, the final parameter estimation is done one flux at a time, each involving only a handful of parameters and thus, the global minimum solution can be either computed analytically (for example, when using loglinear powerlaw flux functions) or determined efficiently. Moreover, as the first estimation phase (flux estimation) depends only on the assumption of the topology of the metabolic network, the flux estimates can subsequently be used to guide the selection of the most appropriate flux functions for the second phase or to detect inconsistencies in the assumed topology of the network separately from the flux equations [14]. However, the application of this method requires the number of metabolites to be equal to or larger than that of fluxes, so that the flux estimation can result in a unique solution. Since the reverse situation is more commonly encountered in the typical metabolic networks, a generalization of this incremental estimation approach becomes the main focus in this study.
As noted above, the new parameter estimation method in this work is built on the concept of incremental identification [12,13] or dynamical flux estimation (DFE) method [14,15]. The proposed method provides two new contributions: (1) an ability to handle the more general scenario, where the number of reactions exceeds that of the metabolites and (2) high numerical efficiency through the reduction of the parameter search space. Specifically, two parameter estimation formulations are proposed with objective functions that depend on model prediction errors of metabolite concentrations and of concentration timeslopes. An extension of this strategy to circumstances where concentration data of some metabolites are missing is also presented. The proposed method is applied to two previously published GMA models and compared with singlestep estimation methods, in order to demonstrate its efficacy.
Methods
The generalized mass action model of cellular metabolism describes the mass balance of metabolites, taking into account all metabolic influxes and effluxes and their stoichiometric ratios, as follows:
where X(t,p) is the vector of metabolic concentration time profiles, S ∈ R^{m × n} is the stoichiometric matrix for m metabolites that participate in n reactions, and v(X,p) denotes the vector of metabolic fluxes (i.e. reaction rates). Here, each flux is described by a powerlaw equation:
where γ_{j} is the rate constant of the jth flux and f_{ji} is the kinetic order parameter, representing the influence of metabolite X_{i} on the jth flux (positive: X_{i} is an activating factor or a substrate, negative: X_{i} is an inhibiting factor). In incremental parameter identification, a data preprocessing
step (e.g. smoothing or filtering) is usually applied to the noisy timecourse concentration
data X_{m}(t_{k}), in order to improve the timeslope estimates
A unique set of dynamic flux values v(t_{k}) can only be computed from
Specifically, we start by decomposing the fluxes into two groups: v(t_{k}) = [ v_{I}(t_{k})^{T}v_{D}(t_{k})^{T} ]^{T} , where the subscripts I and D denote the independent and dependent subset, respectively. Then, the parameter vector
p and the stoichiometric matrix S can be structured correspondingly as p = [ p_{I}p_{D} ] and S = [ S_{I}S_{D} ]. The relationship between the independent and dependent fluxes can be formulated
by rearranging
In this case, given p_{I}, one can compute the independent fluxes v_{I}(X_{m}(t_{k}),p_{I}) using the concentration data X_{m}(t_{k}), and subsequently obtain v_{D}(t_{k}) from Equation (3). Finally, p_{D} can be estimated by a simple least square fitting of v_{D}(X_{m}(t_{k}),p_{D}) to the computed v_{D}(t_{k}) one flux at a time, when there are more time points than the number of parameters in each flux.
In this study, two formulations of the parameter estimation of ODE models in Equation (1) are investigated, involving the minimization of concentration and slope errors. The objective function for the concentration error is given by
and that for the slope error is given by
where K denotes the total number of measurement time points and X(t_{k},p) is the concentration prediction (i.e. the solution to the ODE model in Equation (1)). Figure 1 describes the formulation of the incremental parameter estimation and the procedure for computing the objective functions. Note that the computation of Φ_{C} requires an integration of the ODE model and thus, the estimation using this objective function is expected to be computationally costlier than that using Φ_{S}. On the other hand, metabolic mass balance is only approximately satisfied at discrete time points t_{k} during the parameter estimation using Φ_{S}, as the ODE model is not integrated.
Figure 1. Flowchart of the incremental parameter estimation.
There are several important practical considerations in the implementation of the proposed method. The first consideration is on the selection of the independent fluxes. Here, the set of these fluxes is selected such that (i) the m × m submatrix S_{D} is invertible, (ii) the total number of the independent parameters p_{I} is small, and (iii) the prior knowledge of the corresponding p_{I} is maximized. The last two aspects should lead to a reduction in the parameter search space and the cost of finding the global optimal solution of the minimization problem in Figure 1. The second consideration is regarding constraints in the parameter estimation. Biologically relevant values of parameters are often available, providing lower and/or upper bounds for the parameter estimates. In addition, enzymatic reactions in the ODE model are often assumed to be irreversible and thus, dynamic flux estimates are constrained to be positive. Hence, the parameter estimation involves a constrained minimization problem, for which many global optimization algorithms exist.
So far, we have assumed that the timecourse concentration data are available for all metabolites. However, the method above can be modified to accommodate more general circumstances, in which data for one or several metabolites are missing. In this case, the ODE model is first rewritten to separate the mass balances associated with measured and unmeasured metabolites, such that
where the subscripts M and U refer to components that correspond to measured and unmeasured metabolites, respectively. Again, if the fluxes are split into two categories v_{I} and v_{D} as above, the following relationship still applies for the measured metabolites:
Naturally, the degree of freedom associated with the dynamic flux estimation is higher
by the number of component in X_{U} than before. Figure
2 presents a modification of the parameter estimation procedure in Figure
1 to handle the case of missing data, in which an additional step involving the simulation
of unmeasured metabolites
Figure 2. Flowchart of the incremental parameter estimation when metabolites are not completely measured.
Results
Two case studies: a generic branched pathway [7] and the glycolytic pathway of L. lactis[16], were used to evaluate the performance of the proposed estimation method. In addition, simultaneous estimation methods employing the same objective functions in Equations (4) and (5) were applied to these case studies, to gauge the reduction in the computational cost from using the proposed strategy. In order to alleviate the ODE stiffness issue, parameter combinations that lead to a violation in the MATLAB (ode15s) integration time step criterion is assigned a large error value (Φ_{C} = 10^{3} for the branched pathway and 10^{5} for the glycolytic pathway). Alternatively, one could also set a maximum allowable integration time and penalize the associated parameter values upon violation, as described above. In this study, the optimization problems were solved in MATLAB using publicly available eSSM GO (Enhanced Scatter Search Method for Global Optimization) toolbox, a populationbased metaheuristic global optimization method incorporating probabilistic and deterministic strategies [17,18]. The MATLAB codes of the case studies below are available in Additional file 1. Each parameter estimation was repeated five times to ensure the reliability of the global optimal solution. Unless noted differently, the iterations in the optimization algorithm were terminated when the values of objective functions improve by less than 0.01% or the runtime has exceeded the maximum duration (5 days).
Additional file 1. Incremental Estimation Code. Additional file 1 contains MATLAB codes for the parameter estimations in the two case studies: branched pathway model and L. lactis pathway model.
Format: ZIP Size: 271KB Download file
A generic branched pathway
The generic branched pathway in this example consists of four metabolites and six fluxes, describing the transformations among the metabolites (doubleline arrows), with feedback activation and inhibition (dashed arrows with plus or minus signs, respectively), as shown in Figure 3A. The GMA model of this pathway is given in Figure 3B, containing a total of thirteen rate constants and kinetic orders. This model with the parameter values and initial conditions reported previously [7] were used to generate noisefree and noisy timecourse concentration data (i.i.d additive noise from a Gaussian distribution with 10% coefficient of variation). The noisy data were smoothened using a 6th order polynomial, which provided the best relative goodness of fit among polynomials according to Akaike Information Criterion (AIC) [19] and adjusted R^{2}[20]. Subsequently, timeslopes of noisefree and smoothened noisy data were computed using the central finite difference approximation.
Here, v_{1} and v_{6} were chosen as the independent fluxes as they comprise the least number of kinetic parameters and lead to an invertible S_{D}. The two rate constants and two kinetic orders were constrained to within [0,25] and [0,2], respectively. In addition, all the reactions are assumed to be irreversible.
Table 1 compares simultaneous and incremental parameter estimation runs using noisefree data, employing the two objective functions above. Regardless of the objective function, the proposed incremental approach significantly outperformed the simultaneous estimation. When using the concentrationerror minimization, simultaneous optimization met great difficulty to converge due to stiff ODE integrations. Only one out of five repeated runs could complete after relaxing the convergence criteria of the objective function to 1%, while the others were prematurely terminated after the prescribed maximum runtime of 5 days. In contrast, the proposed incremental estimation was able to find a minima of Φ_{C} in less than 96 seconds on average with good concentration fit and parameter accuracy (see Figure 4A and Table 1). By avoiding ODE integrations using Φ_{S}, the simultaneous estimation of parameters could be completed in roughly 10 minutes duration, but this was much slower than the incremental estimation using Φ_{C}. In this case, the incremental method was able to converge in below 2 seconds or over 250 times faster. The goodness of fit to concentration data and the accuracy of parameter estimates were relatively equal for all three completed estimations (see Figure 4B and Table 1). The parameter inaccuracy in this case was mainly due to the polynomial smoothing of the concentration data, since the same estimations using the analytical values of the slopes (by evaluating the right hand side of the ODE model in Equation (1)) could give accurate parameter estimates (see Additional file 2: Table S1).
Additional file 2. Supplementary Tables. Additional file 2 contains the parameter estimation results of the branched pathway model using noisefree data and analytical slopes, the parameter estimates of the two case studies, and the parameter estimation results of five repeated runs.
Format: PDF Size: 233KB Download file
This file can be viewed with: Adobe Acrobat Reader
Table 1. Parameter estimations of the branched pathway model using noisefree data
Figure 4. Simultaneous and incremental estimation of the branched pathway using in silico noisefree data (×). (A) concentration predictions using parameter estimates from incremental method by Φ_{C} minimization (–––); (B) concentration predictions using parameter estimates from simultaneous method (○) and proposed method () by Φ_{S} minimization.
Table 2 provides the results of the same estimation procedures as above using noisy data. Data noise led to a loss of information and an expected decline in the parameter accuracy. Like before, the simultaneous estimation using Φ_{C} met stiffness problem and three out of five runs did not finish within the fiveday time limit. The incremental approach using either one of the objective functions offered a significant reduction in the computational time over the simultaneous estimation using Φ_{S}, while providing comparable parameter accuracy and concentration and slope fit (see Figure 5 and Table 2). In this example, data noise did not affect the computational cost in obtaining the (global) minimum of the objective functions.
Table 2. Parameter estimations of the branched pathway model using noisy data
Figure 5. Simultaneous and incremental estimation of the branched pathway using in silico noisy data (×). (A) concentration predictions using parameter estimates from incremental method by Φ_{C} minimization (–––); (B) concentration predictions using parameter estimates from simultaneous method (○) and proposed method () by Φ_{S} minimization.
Finally, the estimation strategy described in Figure
2 was applied to this example using noisefree data and assuming X_{3} data were missing. Fluxes v_{3} and v_{4} that appear in
Table 3 summarizes the parameter estimation results. Four out of five repeated runs of Φ_{C} simultaneous optimization were again prematurely terminated after 5 days. Meanwhile, the rest of the estimations could provide reasonably good data fitting with the exception of fitting to X_{3} data as expected (see Figure 6). Like data noise, missing data led to increased inaccuracy of the parameter estimates, regardless of the estimation methods. Finally, the computational speedup by using the incremental over the simultaneous estimation was significant, but was lower than in the previous runs due to the additional integration of X_{U} and the larger number of independent parameters. The detailed values of the parameter estimates in this case study can be found in the Additional file 2: Tables S2 and S3.
Table 3. Parameter estimations of the branched pathway model using noisefree data with X_{3 }missing
Figure 6. Simultaneous and incremental estimation of the branched pathway with missing X_{3}: in silico noisyfree data (×). (A) concentration predictions using parameter estimates from incremental method by Φ_{C} minimization (); (B) concentration predictions using parameter estimates from simultaneous method (○) and proposed method (–––) by Φ_{S} minimization.
The glycolytic pathway in Lactococcus. lactis
The second case study was taken from the GMA modeling of the glycolytic pathway in L. lactis[16], involving six internal metabolites: glucose 6phosphate (G6P) – X_{1}, fructose 1, 6biphosphate (FBP) – X_{2}, 3phosphoglycerate (3PGA) – X_{3}, phosphoenolpyruvate (PEP)  X_{4}, Pyruvate – X_{5}, Lactate – X_{6}, and nine metabolic fluxes. In addition, external glucose (Glu), ATP and Pi are treated as offline variables, whose values were interpolated from measurement data. The pathway connectivity is given in Figure 7A, while the model equations are provided in Figure 7B.
The timecourse concentration dataset of all metabolites were measured using in vivo NMR [21,22], and smoothened data used for the parameter estimations below were shown in Figure 8. The raw data has been filtered previously [16], and these smoothened data for all metabolites but X_{6}, were directly used for the concentration slope calculation in this case study. In the case of X_{6}, a saturating Hilltype equation: k_{1}t^{n} / (k_{2} + t^{n}) where t is time and the constants k_{1}, k_{2}, n are smoothing parameters, was fitted to the filtered data to remove unrealistic fluctuations. The central difference approximation was also adopted to obtain the timeslope data.
Figure 8. Incremental estimation of the L. lactis model: Experimental data (×) compared with model predictions using parameters from concentration error minimization (–––) and slope error minimization ().
Fluxes v_{4}, v_{7} and v_{9} were selected as the DOF, again to give the least number of p_{I} and to ensure that S_{D} is invertible. All rate constants were constrained to within [0, 50], while the independent and dependent kinetic orders were allowed within [0, 5] and [5, 5], respectively. The difference between the bounds for the independent and dependent kinetic orders was done on purpose to simulate a scenario where the signs of the independent kinetic orders were known a priori.
Table 4 reports the outcome of the singlestep and incremental parameter estimation runs using Φ_{C} and Φ_{S}. The values of the parameter estimates are given in the Additional file 2: Table S4. Like in the previous case study, there was a significant reduction in the estimation runtime by using the proposed method over the simultaneous estimation, with comparable goodness of fit in concentration and slope. None of the five repeats of Φ_{C} simultaneous minimization converged within the fiveday time limit, even after relaxing the convergence criteria of the objective function to 1%. On the other hand, the incremental estimation using Φ_{C} was not only able to converge, but was also faster than the simultaneous estimation of Φ_{S} that did not require any ODE integration. The incremental estimation using Φ_{C} was able to provide parameters with the best overall concentration fit (see Figure 8), despite having a large slope error. Finally, minimizing Φ_{S} does not guarantee that the resulting ODE is numerically solvable, as was the case of simultaneous estimation, due to numerical stiffness. But the incremental parameter estimation from minimizing Φ_{S} can produce solvable ODEs with good concentration and slope fits.
Table 4. Parameter estimations of the L. lactis model
Discussion
In this study, an incremental strategy is used to develop a computationally efficient method for the parameter estimation of ODE models. Unlike most commonly used methods, where the parameter estimation is performed to minimize model residuals over the entire parameter space simultaneously, here the estimation is done in two incremental steps, involving the estimation of dynamic reaction rates or fluxes and fluxbased parameter regressions. Importantly, the proposed strategy is designed to handle systems in which there exist extra degrees of freedom in the dynamic flux estimation, when the number of metabolic fluxes exceeds that of metabolites. The positive DOF means that there exist infinitely many solutions to the dynamic flux estimation, which is one of the factors underlying the parameter identifiability issues plaguing many estimation problems in systems biology [23,24].
The main premise of the new method is in recognizing that while many equivalent solutions
exist for the dynamic flux estimation, the subsequent fluxbased regression will give
parameter values with different goodnessoffit, as measured by Φ_{C} or Φ_{S}. In other words, given any two dynamic flux vectors v(t_{k}) satisfying
There are many factors, including datarelated, modelrelated, computational and mathematical issues, which contribute to the difficulty in estimating kinetic parameters of ODE models from timecourse concentration data [1]. Each of these factors has been addressed to a certain degree by using the incremental identification strategy presented in this work. For example, in datarelated issues, the proposed method can be modified to handle the absence of concentration data of some metabolites, as shown in Figure 2. Nevertheless, the method is neither able nor expected to resolve the lack of complete parameter identifiability due to insufficient (dynamical) information contained in the data [23,24]. As illustrated in the first case study, singlestep and incremental approaches provided parameter estimates with similar accuracies, which expectedly deteriorated with noise contamination and loss of data.
The appropriateness of using a particular mathematical formulation, like power law, is an example of modelrelated issues. As discussed above, this issue can be addressed after the dynamic fluxes are estimated, where the chosen functional dependence of the fluxes on a specific set of metabolite concentrations can be tested prior to the parameter regression [14]. Next, the computational issues associated with performing a global optimization over a large number of variables and the need to integrate ODEs have been mitigated in the proposed method by performing optimization only over the independent parameter subset and using a minimization of slope error, respectively. Finally, in this work, we have also addressed a mathematical issue related to the degrees of freedom that exist during the inference of dynamic fluxes from slopes of concentration data. However, extra degrees of freedom (mathematical redundancies) are also expected to influence the second step of the method, i.e. onefluxatatime parameter estimation. For (log)linear regression of parameters in GMA models, such redundancy will lead to a lack of full column rank of the matrix containing the logarithms of concentration data X_{m}(t_{k}) and thus, can be straightforwardly detected.
The proposed estimation method has several weaknesses that are common among incremental estimation methods. As demonstrated in the first case study, the accuracy of the identified parameter relies on the ability to obtain good estimates of the concentration slopes. Direct slope estimation from the raw data, for example using central finite difference approximation, is usually not advisable due to high degree of noise in the typical biological data. Hence, presmoothing of the timecourse data is often required, as done in this study. Many algorithms are available for such purpose, from simplistic polynomial regression and splines to more advanced artificial neural network [7,25] and WhittakerEilers smoother [26,27]. If reliable concentration slope estimates are not available, but bounds for the slope values can be obtained, then one can use interval arithmetic to derive upper and lower limits for the dependent fluxes and parameters using Equation (3) (or Equation (7) [28]. When the objective function involves integrating the model, validated solution to ODE with interval parameters can be used to produce the corresponding upper and lower bounds of concentration predictions [29]. Finally, the estimation can be reformulated, for example by minimizing the upper bound of the objective.
In addition to the drawback discussed above, the proposed strategy requires a priori knowledge about the topology of the network. For cellular metabolism, such information has become more readily available as genomescale metabolic network of many important organisms, including human, E. coli and S. cereviseae, have been and are continuously being reconstructed [30]. For other networks, many algorithms also exist for the estimation of network topology based on timeseries concentration data, including Bayesian network inference, transfer entropy, and Granger causality [3133].
Conclusions
The estimation of kinetic parameters of ODE models from timecourse concentration
data remains a key bottleneck in model building in systems biology. The lack of complete
parameter identifiability has been blamed as the root cause of the difficulty in such
estimation. In this study, a new incremental estimation method is proposed that is
able to overcome the existence of extra degrees of freedom in the dynamic flux estimation
from concentration slopes and to significantly reduce the computational requirements
in finding parameter estimates. The method can also be applied, after minor modifications,
to circumstances where concentration data for a few molecules are missing. While the
present work concerns with the GMA modeling of metabolic networks, the estimation
strategies discussed in this work have general applicability to any kinetic models
that can be written as
Competing interest
The authors declare that they have no competing interests.
Authors’ contributions
GJ conceived of the study, carried out the parameter estimation and wrote the manuscript. GS participated in the design of the study. RG conceived and guided the study and wrote the manuscript. All authors have read and approved the final manuscript.
Funding
SingaporeMIT Alliance and ETH Zurich.
References

Chou IC, Voit EO: Recent developments in parameter estimation and structure identification of biochemical and genomic systems.

Mendes P, Kell D: Nonlinear optimization of biochemical pathways: applications to metabolic engineering and parameter estimation.

Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: a comparison of global optimization methods.

Savageau MA: Biochemical systems analysis. I. Some mathematical properties of the rate law for the component enzymatic reactions.

Savageau MA: Biochemical systems analysis. II. The steadystate solutions for an npool system using a powerlaw approximation.

Voit EO: Computational analysis of biochemical systems: a practical guide for biochemists and molecular biologists. New York: Cambridge University Press; 2000.

Voit EO, Almeida J: Decoupling dynamical systems for pathway identification from metabolic profiles.

Tsai KY, Wang FS: Evolutionary optimization with data collocation for reverse engineering of biological networks.

Kimura S, Ide K, Kashihara A, Kano M, Hatakeyama M, Masui R, Nakagawa N, Yokoyama S, Kuramitsu S, Konagaya A: Inference of Ssystem models of genetic networks using a cooperative coevolutionary algorithm.

Maki Y, Ueda T, Masahiro O, Naoya U, Kentaro I, Uchida K: Inference of genetic network using the expression profile time course data of mouse P19 cells.

Jia G, Stephanopoulos G, Gunawan R: Parameter estimation of kinetic models from metabolic profiles: twophase dynamic decoupling method.

Bardow A, Marquardt W: Incremental and simultaneous identification of reaction kinetics: methods and comparison.

Marquardt W, Brendel M, Bonvin D: Incremental identification of kinetic models for homogeneous reaction systems.

Goel G, Chou IC, Voit EO: System estimation from metabolic timeseries data.

Voit EO, Goel G, Chou IC, Fonseca LL: Estimation of metabolic pathway systems from different data sources.

Voit EO, Almeida J, Marino S, Lall R, Goel G, Neves AR, Santos H: Regulation of glycolysis in Lactococcus lactis: an unfinished systems biological case study.

Egea JA, RodriguezFernandez M, Banga JR, Marti R: Scatter search for chemical and bioprocess optimization.

RodriguezFernandez M, Egea JA, Banga JR: Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems.

Montgomery DC, Runger GC: Applied statistics and probability for engineers. 4th edition. Hoboken, NJ: Wiley; 2007.

Neves AR, Ramos A, Costa H, van Swam II, Hugenholtz J, Kleerebezem M, de Vos W, Santos H: Effect of different NADH oxidase levels on glucose metabolism by Lactococcus lactis: kinetics of intracellular metabolite pools determined by in vivo nuclear magnetic resonance.

Neves AR, Ramos A, Nunes MC, Kleerebezem M, Hugenholtz J, de Vos WM, Almeida J, Santos H: In vivo nuclear magnetic resonance studies of glycolytic kinetics in Lactococcus lactis.

Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmuller U, Timmer J: Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood.

Srinath S, Gunawan R: Parameter identifiability of powerlaw biochemical system models.

Almeida JS: Predictive nonlinear modeling of complex data by artificial neural networks.

Vilela M, Borges CC, Vinga S, Vasconcelos AT, Santos H, Voit EO, Almeida JS: Automated smoother for the numerical decoupling of dynamics models.

Jaulin L, Kieffer M, Didrit O, Walter E: Applied interval analysis: with examples in parameter and state estimation, robust control and robotics. London: Springer; 2001.

Lin YD, Stadtherr MA: Validated solution of ODEs with parametric uncertainties.
16th European Symposium on Computer Aided Process Engineering and 9th International Symposium on Process Systems Engineering 2006, 21:167172.

Latendresse M, Paley S, Karp PD: Browsing metabolic and regulatory networks with BioCyc.

Imoto S, Kim S, Goto T, Miyano S, Aburatani S, Tashiro K, Kuhara S: Bayesian network and nonparametric heteroscedastic regression for nonlinear modeling of genetic network.

Nagarajan R, Upreti M: Comment on causality and pathway search in microarray time series experiment.

Tung TQ, Ryu T, Lee KH, Lee D: Inferring gene regulatory networks from microarray time series data using transfer entropy. In Proceedings of the Twentieth IEEE International Symposium on ComputerBased Medical Systems:2022 June 2007; Maribor, Slovenia. Edited by Kokol P, Los A. Los Alamitos: IEEE Computer Society; 2007:383388.