BACKGROUND
the estimation of unknown kinetic parameters from time-series measurements of biological molecules is a major bottleneck in the ode model building process in systems biology and metabolic engineering
<cit> .
the majority of current estimation methods involve simultaneous  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
<cit> .
the problem of obtaining the best-fit parameter estimates however, is typically ill-posed 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 power-law formalism, specifically the generalized mass action  systems
<cit> .
the power-law formalism has many advantages, which have been detailed elsewhere
<cit> .
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 large-scale parameter identification problem, one where single-step 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
<cit> .
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 power-law 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
<cit>  or partially
<cit> .
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
<cit> .
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  or determined efficiently.
moreover, as the first estimation phase  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
<cit> .
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
<cit>  or dynamical flux estimation  method
<cit> .
the proposed method provides two new contributions:  an ability to handle the more general scenario, where the number of reactions exceeds that of the metabolites and  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 time-slopes.
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 single-step 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:
dxt,p/dt=x˙t,p=svx,p,
where x is the vector of metabolic concentration time profiles, s ∈ rm × n is the stoichiometric matrix for m metabolites that participate in n reactions, and v denotes the vector of metabolic fluxes .
here, each flux is described by a power-law equation:
vjx,p=γj∏ixifji,
where γj is the rate constant of the j-th flux and fji is the kinetic order parameter, representing the influence of metabolite xi on the j-th flux .
in incremental parameter identification, a data pre-processing step  is usually applied to the noisy time-course concentration data xm, in order to improve the time-slope estimates
x˙mtk.
subsequently, the dynamic metabolic fluxes v are estimated from equation  by substituting
x˙t with
x˙mtk.
finally, the kinetic parameters associated with the j-th flux  can be calculated using a least square regression of the power law flux function in equation  against the estimated vj.
note that for gma models, the least square parameter regressions in the last step are linear in the logarithmic scale and thus, can be performed very efficiently.
a unique set of dynamic flux values v can only be computed from
x˙mtk=svtk, when the number of metabolites exceeds that of fluxes.
however, a metabolite in general can participate in more than one metabolic flux .
in such a situation, there exist an infinite number of dynamic flux combinations v that satisfy
x˙mtk=svtk.
the dimensionality of the set of flux solutions is equal to the degree of freedom , given by the difference between the number of fluxes and the number of metabolites: ndof = n-m > <dig> ).
the positive dof means that the values of ndof selected fluxes can be independently set, from which the remaining fluxes can be computed.
this relationship forms the basis of the proposed estimation method, in which the model goodness of fit to data is optimized by adjusting only a subset of parameters associated with the independent fluxes above.
specifically, we start by decomposing the fluxes into two groups: v = 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 =  and s = .
the relationship between the independent and dependent fluxes can be formulated by rearranging
x˙mtk=svtk into:
vdtk=sd−1x˙mtk−sivixmtk,pi.
in this case, given pi, one can compute the independent fluxes vi,pi) using the concentration data xm, and subsequently obtain vd from equation .
finally, pd can be estimated by a simple least square fitting of vd,pd) to the computed vd 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  are investigated, involving the minimization of concentration and slope errors.
the objective function for the concentration error is given by
Φcp,x=1mk∑k=1kxmtk−xtk,ptxmtk−xtk,p
and that for the slope error is given by
Φsp,x=1mk∑k=1kx˙mtk−svxmtk,ptx˙mtk−svxmtk,p,
where k denotes the total number of measurement time points and x is the concentration prediction ).
figure
<dig> 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 tk during the parameter estimation using Φs, as the ode model is not integrated.
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  the m × m submatrix sd is invertible,  the total number of the independent parameters pi is small, and  the prior knowledge of the corresponding pi 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
<dig>  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 time-course 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
x˙t,p=x˙mx˙ut,p=smsuvxm,xu,p
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 vi and vd as above, the following relationship still applies for the measured metabolites:
vdtk=sd,m−1x˙mtk−si,mvitk
naturally, the degree of freedom associated with the dynamic flux estimation is higher by the number of component in xu than before.
figure
<dig> presents a modification of the parameter estimation procedure in figure
<dig> to handle the case of missing data, in which an additional step involving the simulation of unmeasured metabolites
x˙u=suvxm,xu,p will be performed.
in this integration, xm is set as an external variable, whose time-profiles are interpolated from the measured concentrations.
the set of independent fluxes vi are now selected to include all fluxes that appear in
x˙u and those that lead to a full column ranked sd,m.
if sd,m is a non-square matrix, then a pseudo-inverse will be done in equation .
of course, the same considerations mentioned above are equally relevant in this case.
note that the initial conditions of xu will also need to be estimated.
RESULTS
two case studies: a generic branched pathway
<cit>  and the glycolytic pathway of l. lactis <cit> , were used to evaluate the performance of the proposed estimation method.
in addition, simultaneous estimation methods employing the same objective functions in equations  and  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  integration time step criterion is assigned a large error value .
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  toolbox, a population-based metaheuristic global optimization method incorporating probabilistic and deterministic strategies
<cit> .
the matlab codes of the case studies below are available in additional file
<dig>  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  <dig> % or the runtime has exceeded the maximum duration .
a generic branched pathway
the generic branched pathway in this example consists of four metabolites and six fluxes, describing the transformations among the metabolites , with feedback activation and inhibition , 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
<cit>  were used to generate noise-free and noisy time-course concentration data .
the noisy data were smoothened using a 6-th order polynomial, which provided the best relative goodness of fit among polynomials according to akaike information criterion
<cit>  and adjusted r <dig> <cit> .
subsequently, time-slopes of noise-free and smoothened noisy data were computed using the central finite difference approximation.
<cit> .
here, v <dig> and v <dig> were chosen as the independent fluxes as they comprise the least number of kinetic parameters and lead to an invertible sd.
the two rate constants and two kinetic orders were constrained to within  <cit>  and  <cit> , respectively.
in addition, all the reactions are assumed to be irreversible.
table
<dig> compares simultaneous and incremental parameter estimation runs using noise-free data, employing the two objective functions above.
regardless of the objective function, the proposed incremental approach significantly outperformed the simultaneous estimation.
when using the concentration-error 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  <dig> days.
in contrast, the proposed incremental estimation was able to find a minima of Φc in less than  <dig> seconds on average with good concentration fit and parameter accuracy .
by avoiding ode integrations using Φs, the simultaneous estimation of parameters could be completed in roughly  <dig> 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  <dig> seconds or over  <dig> times faster.
the goodness of fit to concentration data and the accuracy of parameter estimates were relatively equal for all three completed estimations .
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 ) could give accurate parameter estimates .
a. cpu time was based on a workstation with dual intel quad-core  <dig>  ghz processors.
b. only one out of five runs completed with a relative improvement of the objective function below 1% between iterations.
the rest did not converge within the 5-day time limit after iterating for  <dig>   <dig>   <dig>  and  <dig> times.
the corresponding Φc at termination were  <dig> × 10- <dig>   <dig>  × 10- <dig>   <dig>  × 10- <dig> and  <dig>  × 10- <dig>  respectively.
c. mean value ± standard deviation out of five repeats.
d. root mean square error of model predictions, where the underlined value refers to the objective function of the minimization.
table
<dig> 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 five-day 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 .
in this example, data noise did not affect the computational cost in obtaining the  minimum of the objective functions.
a. two out of five runs completed with a relative improvement of the objective function below 1% between iterations.
the rest did not converge within the 5-day time limit after iterating for  <dig>   <dig>  and  <dig> times.
the corresponding Φc at termination were  <dig>  × 10- <dig>   <dig>  × 10- <dig> and  <dig>  × 10- <dig>  respectively.
finally, the estimation strategy described in figure
<dig> was applied to this example using noise-free data and assuming x <dig> data were missing.
fluxes v <dig> and v <dig> that appear in
x˙ <dig> were chosen to be among the independent fluxes and flux v <dig> was also added to the set such that the dependent fluxes can be uniquely determined from equation .
in addition to the parameters associated with the aforementioned fluxes, the initial condition x <dig> was also estimated.
the bounds for the rate constants and kinetic orders were kept the same as above, while the initial concentration was bounded within  <cit> .
table
<dig> summarizes the parameter estimation results.
four out of five repeated runs of Φc simultaneous optimization were again prematurely terminated after  <dig> days.
meanwhile, the rest of the estimations could provide reasonably good data fitting with the exception of fitting to x <dig> data as expected .
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 xu 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 s <dig> and s <dig>
a. only one out of five runs completed with a relative improvement of the objective function below 1% between iterations.
the rest did not converge within the 5-day time limit after iterating for  <dig>   <dig>   <dig> and  <dig> times.
the corresponding Φc at termination were  <dig> × 10- <dig>   <dig>  × 10- <dig>   <dig>  × 10- <dig> and  <dig>  × 10- <dig>  respectively.
the glycolytic pathway in lactococcus.
lactis
the second case study was taken from the gma modeling of the glycolytic pathway in l. lactis <cit> , involving six internal metabolites: glucose 6-phosphate  – x <dig>  fructose  <dig>  6-biphosphate  – x <dig>  3-phosphoglycerate  – x <dig>  phosphoenolpyruvate  - x <dig>  pyruvate – x <dig>  lactate – x <dig>  and nine metabolic fluxes.
in addition, external glucose , atp and pi are treated as off-line 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.
<cit> .
the time-course concentration dataset of all metabolites were measured using in vivo nmr
<cit> , and smoothened data used for the parameter estimations below were shown in figure
<dig>  the raw data has been filtered previously
<cit> , and these smoothened data for all metabolites but x <dig>  were directly used for the concentration slope calculation in this case study.
in the case of x <dig>  a saturating hill-type equation: k1tn /  where t is time and the constants k <dig>  k <dig>  n are smoothing parameters, was fitted to the filtered data to remove unrealistic fluctuations.
the central difference approximation was also adopted to obtain the time-slope data.
fluxes v <dig>  v <dig> and v <dig> were selected as the dof, again to give the least number of pi and to ensure that sd is invertible.
all rate constants were constrained to within  <cit> , while the independent and dependent kinetic orders were allowed within  <cit>  and , 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
<dig> reports the outcome of the single-step and incremental parameter estimation runs using Φc and Φs.
the values of the parameter estimates are given in the additional file
2: table s <dig>  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 five-day 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 , 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.
a. none of five runs finished with a relative improvement of the objective function below 1% within the 5-day time limit, after iterating for  <dig>   <dig>   <dig>   <dig> and  <dig> times.
the corresponding Φc at termination were  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig> , respectively.
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 flux-based 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
<cit> .
the main premise of the new method is in recognizing that while many equivalent solutions exist for the dynamic flux estimation, the subsequent flux-based regression will give parameter values with different goodness-of-fit, as measured by Φc or Φs.
in other words, given any two dynamic flux vectors v satisfying
x˙mtk=svtk, the associated parameter pairs  may not predict the slope or concentration data equally well, due to differences in the quality of parameter regression for each v. also, because of the dof, the minimization of model residuals needs to be done only over a subset of parameters that are associated with the flux degrees of freedom, resulting in much reduced parameter search space and correspondingly much faster convergence to the  optimal solution.
the superior performance of the proposed method over simultaneous estimation was convincingly demonstrated in the two gma modeling case studies in the previous section.
the minimization of slope error, also known as slope-estimation-decoupling strategy method
<cit> , is arguably one of the most computationally efficient simultaneous methods.
in this strategy, the parameter fitting essentially constitutes a zero-finding problem and the estimation can be done without having to integrate the odes.
yet, the incremental estimation could offer more than two orders of magnitude reduction in the computational time over this strategy.
there are many factors, including data-related, model-related, computational and mathematical issues, which contribute to the difficulty in estimating kinetic parameters of ode models from time-course concentration data
<cit> .
each of these factors has been addressed to a certain degree by using the incremental identification strategy presented in this work.
for example, in data-related issues, the proposed method can be modified to handle the absence of concentration data of some metabolites, as shown in figure
<dig>  nevertheless, the method is neither able nor expected to resolve the lack of complete parameter identifiability due to insufficient  information contained in the data
<cit> .
as illustrated in the first case study, single-step 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 model-related 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
<cit> .
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  are also expected to influence the second step of the method, i.e.
one-flux-at-a-time parameter estimation.
for 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 xm 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, pre-smoothing of the time-course 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
<cit>  and whittaker-eilers smoother
<cit> .
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
<cit> .
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
<cit> .
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 genome-scale metabolic network of many important organisms, including human, e. coli and s. cereviseae, have been and are continuously being reconstructed
<cit> .
for other networks, many algorithms also exist for the estimation of network topology based on time-series concentration data, including bayesian network inference, transfer entropy, and granger causality
<cit> .
CONCLUSIONS
the estimation of kinetic parameters of ode models from time-course 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
x˙tk=svtk.
the creation of computationally efficient parameter estimation methods, such as the one presented here, represents an important step toward genome-scale kinetic modeling of cellular metabolism.
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
singapore-mit alliance and eth zurich.
supplementary material
additional file 1
incremental estimation code.
additional file
<dig> contains matlab codes for the parameter estimations in the two case studies: branched pathway model and l. lactis pathway model.
click here for file
additional file 2
supplementary tables.
additional file  <dig> contains the parameter estimation results of the branched pathway model using noise-free data and analytical slopes, the parameter estimates of the two case studies, and the parameter estimation results of five repeated runs.
click here for file

 the estimation of unknown kinetic parameters from time-series measurements of biological molecules is a major bottleneck in the ode model building process in systems biology and metabolic engineering. the majority of current estimation methods involve simultaneous 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. the problem of obtaining the best-fit parameter estimates however, is typically ill-posed 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 power-law formalism, specifically the generalized mass action systems. the power-law formalism has many advantages, which have been detailed elsewhere. 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 large-scale parameter identification problem, one where single-step 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. 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 power-law 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 or partially. 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. 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 or determined efficiently. moreover, as the first estimation phase 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. 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 or dynamical flux estimation method. the proposed method provides two new contributions: an ability to handle the more general scenario, where the number of reactions exceeds that of the metabolites and 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 time-slopes. 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 single-step 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: dxt,p/dt=x˙t,p=svx,p, where x is the vector of metabolic concentration time profiles, s ∈ rm × n is the stoichiometric matrix for m metabolites that participate in n reactions, and v denotes the vector of metabolic fluxes. here, each flux is described by a power-law equation: vjx,p=γj∏ixifji, where γj is the rate constant of the j-th flux and fji is the kinetic order parameter, representing the influence of metabolite xi on the j-th flux. in incremental parameter identification, a data pre-processing step is usually applied to the noisy time-course concentration data xm, in order to improve the time-slope estimates x˙mtk. subsequently, the dynamic metabolic fluxes v are estimated from equation by substituting x˙t with x˙mtk. finally, the kinetic parameters associated with the j-th flux can be calculated using a least square regression of the power law flux function in equation against the estimated vj. note that for gma models, the least square parameter regressions in the last step are linear in the logarithmic scale and thus, can be performed very efficiently. a unique set of dynamic flux values v can only be computed from x˙mtk=svtk, when the number of metabolites exceeds that of fluxes. however, a metabolite in general can participate in more than one metabolic flux. in such a situation, there exist an infinite number of dynamic flux combinations v that satisfy x˙mtk=svtk. the dimensionality of the set of flux solutions is equal to the degree of freedom, given by the difference between the number of fluxes and the number of metabolites: ndof = n-m > ). the positive dof means that the values of ndof selected fluxes can be independently set, from which the remaining fluxes can be computed. this relationship forms the basis of the proposed estimation method, in which the model goodness of fit to data is optimized by adjusting only a subset of parameters associated with the independent fluxes above. specifically, we start by decomposing the fluxes into two groups: v = 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 =  and s =. the relationship between the independent and dependent fluxes can be formulated by rearranging x˙mtk=svtk into: vdtk=sd−1x˙mtk−sivixmtk,pi. in this case, given pi, one can compute the independent fluxes vi,pi) using the concentration data xm, and subsequently obtain vd from equation. finally, pd can be estimated by a simple least square fitting of vd,pd) to the computed vd 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 are investigated, involving the minimization of concentration and slope errors. the objective function for the concentration error is given by Φcp,x=1mk∑k=1kxmtk−xtk,ptxmtk−xtk,p and that for the slope error is given by Φsp,x=1mk∑k=1kx˙mtk−svxmtk,ptx˙mtk−svxmtk,p, where k denotes the total number of measurement time points and x is the concentration prediction ). figure 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 tk during the parameter estimation using Φs, as the ode model is not integrated. 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 the m × m submatrix sd is invertible, the total number of the independent parameters pi is small, and the prior knowledge of the corresponding pi 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 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 time-course 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 x˙t,p=x˙mx˙ut,p=smsuvxm,xu,p 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 vi and vd as above, the following relationship still applies for the measured metabolites: vdtk=sd,m−1x˙mtk−si,mvitk naturally, the degree of freedom associated with the dynamic flux estimation is higher by the number of component in xu than before. figure presents a modification of the parameter estimation procedure in figure to handle the case of missing data, in which an additional step involving the simulation of unmeasured metabolites x˙u=suvxm,xu,p will be performed. in this integration, xm is set as an external variable, whose time-profiles are interpolated from the measured concentrations. the set of independent fluxes vi are now selected to include all fluxes that appear in x˙u and those that lead to a full column ranked sd,m. if sd,m is a non-square matrix, then a pseudo-inverse will be done in equation. of course, the same considerations mentioned above are equally relevant in this case. note that the initial conditions of xu will also need to be estimated.  two case studies: a generic branched pathway and the glycolytic pathway of l. lactis, were used to evaluate the performance of the proposed estimation method. in addition, simultaneous estimation methods employing the same objective functions in equations and 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 integration time step criterion is assigned a large error value. 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 toolbox, a population-based metaheuristic global optimization method incorporating probabilistic and deterministic strategies. the matlab codes of the case studies below are available in additional file 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 % or the runtime has exceeded the maximum duration. a generic branched pathway the generic branched pathway in this example consists of four metabolites and six fluxes, describing the transformations among the metabolites, with feedback activation and inhibition, 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 were used to generate noise-free and noisy time-course concentration data. the noisy data were smoothened using a 6-th order polynomial, which provided the best relative goodness of fit among polynomials according to akaike information criterion and adjusted r. subsequently, time-slopes of noise-free and smoothened noisy data were computed using the central finite difference approximation.. here, v and v were chosen as the independent fluxes as they comprise the least number of kinetic parameters and lead to an invertible sd. the two rate constants and two kinetic orders were constrained to within and, respectively. in addition, all the reactions are assumed to be irreversible. table compares simultaneous and incremental parameter estimation runs using noise-free data, employing the two objective functions above. regardless of the objective function, the proposed incremental approach significantly outperformed the simultaneous estimation. when using the concentration-error 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 days. in contrast, the proposed incremental estimation was able to find a minima of Φc in less than seconds on average with good concentration fit and parameter accuracy. by avoiding ode integrations using Φs, the simultaneous estimation of parameters could be completed in roughly 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 seconds or over times faster. the goodness of fit to concentration data and the accuracy of parameter estimates were relatively equal for all three completed estimations. 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 ) could give accurate parameter estimates. a. cpu time was based on a workstation with dual intel quad-core ghz processors. b. only one out of five runs completed with a relative improvement of the objective function below 1% between iterations. the rest did not converge within the 5-day time limit after iterating for and times. the corresponding Φc at termination were × 10-  × 10-  × 10- and  × 10- respectively. c. mean value ± standard deviation out of five repeats. d. root mean square error of model predictions, where the underlined value refers to the objective function of the minimization. table 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 five-day 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. in this example, data noise did not affect the computational cost in obtaining the minimum of the objective functions. a. two out of five runs completed with a relative improvement of the objective function below 1% between iterations. the rest did not converge within the 5-day time limit after iterating for and times. the corresponding Φc at termination were  × 10-  × 10- and  × 10- respectively. finally, the estimation strategy described in figure was applied to this example using noise-free data and assuming x data were missing. fluxes v and v that appear in x˙ were chosen to be among the independent fluxes and flux v was also added to the set such that the dependent fluxes can be uniquely determined from equation. in addition to the parameters associated with the aforementioned fluxes, the initial condition x was also estimated. the bounds for the rate constants and kinetic orders were kept the same as above, while the initial concentration was bounded within. table summarizes the parameter estimation results. four out of five repeated runs of Φc simultaneous optimization were again prematurely terminated after days. meanwhile, the rest of the estimations could provide reasonably good data fitting with the exception of fitting to x data as expected. 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 xu 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 s and s a. only one out of five runs completed with a relative improvement of the objective function below 1% between iterations. the rest did not converge within the 5-day time limit after iterating for and times. the corresponding Φc at termination were × 10-  × 10-  × 10- and  × 10- respectively. the glycolytic pathway in lactococcus. lactis the second case study was taken from the gma modeling of the glycolytic pathway in l. lactis, involving six internal metabolites: glucose 6-phosphate – x fructose 6-biphosphate – x 3-phosphoglycerate – x phosphoenolpyruvate - x pyruvate – x lactate – x and nine metabolic fluxes. in addition, external glucose, atp and pi are treated as off-line 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 time-course concentration dataset of all metabolites were measured using in vivo nmr, and smoothened data used for the parameter estimations below were shown in figure the raw data has been filtered previously, and these smoothened data for all metabolites but x were directly used for the concentration slope calculation in this case study. in the case of x a saturating hill-type equation: k1tn / where t is time and the constants k k n are smoothing parameters, was fitted to the filtered data to remove unrealistic fluctuations. the central difference approximation was also adopted to obtain the time-slope data. fluxes v v and v were selected as the dof, again to give the least number of pi and to ensure that sd is invertible. all rate constants were constrained to within, while the independent and dependent kinetic orders were allowed within and, 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 reports the outcome of the single-step and incremental parameter estimation runs using Φc and Φs. the values of the parameter estimates are given in the additional file 2: table s 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 five-day 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, 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. a. none of five runs finished with a relative improvement of the objective function below 1% within the 5-day time limit, after iterating for and times. the corresponding Φc at termination were,,, and, respectively. 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 flux-based 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. the main premise of the new method is in recognizing that while many equivalent solutions exist for the dynamic flux estimation, the subsequent flux-based regression will give parameter values with different goodness-of-fit, as measured by Φc or Φs. in other words, given any two dynamic flux vectors v satisfying x˙mtk=svtk, the associated parameter pairs may not predict the slope or concentration data equally well, due to differences in the quality of parameter regression for each v. also, because of the dof, the minimization of model residuals needs to be done only over a subset of parameters that are associated with the flux degrees of freedom, resulting in much reduced parameter search space and correspondingly much faster convergence to the optimal solution. the superior performance of the proposed method over simultaneous estimation was convincingly demonstrated in the two gma modeling case studies in the previous section. the minimization of slope error, also known as slope-estimation-decoupling strategy method, is arguably one of the most computationally efficient simultaneous methods. in this strategy, the parameter fitting essentially constitutes a zero-finding problem and the estimation can be done without having to integrate the odes. yet, the incremental estimation could offer more than two orders of magnitude reduction in the computational time over this strategy. there are many factors, including data-related, model-related, computational and mathematical issues, which contribute to the difficulty in estimating kinetic parameters of ode models from time-course concentration data. each of these factors has been addressed to a certain degree by using the incremental identification strategy presented in this work. for example, in data-related issues, the proposed method can be modified to handle the absence of concentration data of some metabolites, as shown in figure nevertheless, the method is neither able nor expected to resolve the lack of complete parameter identifiability due to insufficient information contained in the data. as illustrated in the first case study, single-step 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 model-related 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. 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 are also expected to influence the second step of the method, i.e. one-flux-at-a-time parameter estimation. for 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 xm 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, pre-smoothing of the time-course 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 and whittaker-eilers smoother. 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. 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. 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 genome-scale metabolic network of many important organisms, including human, e. coli and s. cereviseae, have been and are continuously being reconstructed. for other networks, many algorithms also exist for the estimation of network topology based on time-series concentration data, including bayesian network inference, transfer entropy, and granger causality.  the estimation of kinetic parameters of ode models from time-course 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 x˙tk=svtk. the creation of computationally efficient parameter estimation methods, such as the one presented here, represents an important step toward genome-scale kinetic modeling of cellular metabolism. 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 singapore-mit alliance and eth zurich. supplementary material additional file 1 incremental estimation code. additional file contains matlab codes for the parameter estimations in the two case studies: branched pathway model and l. lactis pathway model. click here for file additional file 2 supplementary tables. additional file contains the parameter estimation results of the branched pathway model using noise-free data and analytical slopes, the parameter estimates of the two case studies, and the parameter estimation results of five repeated runs. click here for file
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 time-slopes. most of the existing estimation methods involve finding the global minimum of data fitting residuals over the entire parameter space simultaneously. an efficient and reliable parameter estimation method is essential for the creation of biological models using ordinary differential equation. in this work, an incremental approach was applied to the parameter estimation of ode models from concentration time profiles. incremental parameter estimationkinetic modelingmetabolic networkgma model unfortunately, the associated computational requirement often becomes prohibitively high due to the large number of parameters and the lack of complete parameter identifiability. 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 genome-scale cellular metabolism in the future. particularly, the method was developed to address a commonly encountered circumstance in the modeling of metabolic networks, where the number of metabolic fluxes exceeds that of metabolites. in addition, an extension of the estimation method to handle missing data is also presented.
 an efficient and reliable parameter estimation method is essential for the creation of biological models using ordinary differential equation. 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.  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 exceeds that of metabolites. 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 time-slopes. the efficacy of this method was demonstrated using two generalized mass action models, where the method significantly outperformed single-step estimations. in addition, an extension of the estimation method to handle missing data is also presented.  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 genome-scale cellular metabolism in the future. incremental parameter estimationkinetic modelingmetabolic networkgma model 