Steadystate global optimization of metabolic nonlinear dynamic models through recasting into powerlaw canonical models
 Carlos Pozo^{1},
 Alberto MarínSanguino^{2},
 Rui Alves^{3},
 Gonzalo GuillénGosálbez^{1},
 Laureano Jiménez^{1} and
 Albert Sorribas^{3}Email author
https://doi.org/10.1186/175205095137
© Pozo et al; licensee BioMed Central Ltd. 2011
Received: 5 April 2011
Accepted: 25 August 2011
Published: 25 August 2011
Abstract
Background
Design of newly engineered microbial strains for biotechnological purposes would greatly benefit from the development of realistic mathematical models for the processes to be optimized. Such models can then be analyzed and, with the development and application of appropriate optimization techniques, one could identify the modifications that need to be made to the organism in order to achieve the desired biotechnological goal. As appropriate models to perform such an analysis are necessarily nonlinear and typically nonconvex, finding their global optimum is a challenging task. Canonical modeling techniques, such as Generalized Mass Action (GMA) models based on the powerlaw formalism, offer a possible solution to this problem because they have a mathematical structure that enables the development of specific algorithms for global optimization.
Results
Based on the GMA canonical representation, we have developed in previous works a highly efficient optimization algorithm and a set of related strategies for understanding the evolution of adaptive responses in cellular metabolism. Here, we explore the possibility of recasting kinetic nonlinear models into an equivalent GMA model, so that global optimization on the recast GMA model can be performed. With this technique, optimization is greatly facilitated and the results are transposable to the original nonlinear problem. This procedure is straightforward for a particular class of nonlinear models known as Saturable and Cooperative (SC) models that extend the powerlaw formalism to deal with saturation and cooperativity.
Conclusions
Our results show that recasting nonlinear kinetic models into GMA models is indeed an appropriate strategy that helps overcoming some of the numerical difficulties that arise during the global optimization task.
Keywords
1 Background
Identifying optimization strategies for increasing strain productivity should be possible by applying optimization methods to detailed kinetic models of the target metabolism. Thus, a rational approach would pinpoint the changes to be done  e.g. by modulating gene expression  in order to achieve the desired biotechnological goals [1–4]. To build such models we can either start from a detailed description of the underlying processes (bottomup strategy) or we can fit kinetic models to experimental data obtained in vivo (topdown strategy).
The bottomup approach was the original strategy for model building in the biological sciences. Bottomup kinetic models require information that is seldom available, despite the increasing amount of kinetic data contained in a growing set of databases (for example see [5, 6] and the webpage http://kinetics.nist.gov/kinetics/index.jsp). Even in the best case scenarios where kinetic data are available, the data have often been obtained in different labs and under in vitro conditions that are hardly ever comparable or representative of the situation in vivo. In addition, models built using this strategy often fail to adequately reproduce the known behavior of the target system [7–10]. With the accumulation of timeseries data, which were originally generated from the accurate measurement of transient responses, topdown modeling became viable as an alternative to the bottomup strategy [11]. However, topdown modeling also faces important difficulties. For example, regulatory interactions between metabolites and enzymes are very poorly characterized and most metabolic maps lack such crucial information. Therefore, for a given network structure (i.e. a stoichiometric description) obtained from databases, a large number of alternative regulatory patterns may exist that account for the observed experimental data [12]. Model discrimination among the alternative regulatory patterns requires appropriate experimental design. However, this is seldom considered when performing the time series measurements. Last, but not the least, parameter identifiability in highly nonlinear models can be problematic (for a review see [13]).
An additional issue that is common to models built using both strategies is that such detailed kinetic models include nonconvexities that lead to the existence of multiple local optima in which standard nonlinear optimization algorithms may get trapped during the search. Several stochastic and deterministic global optimization methods have been proposed to overcome this limitation [14]. Deterministic methods, which are the only ones that can rigorously guarantee global optimality, rely on the use of convex envelopes or underestimators to formulate lowerbounding convex problems that are typically combined with spatial branch and bound strategies. Most of these methods are general purpose and assume special structures in the continuous terms of the mathematical model. Because of this, they can encounter numerical difficulties in specific metabolic engineering systems that require the optimization of kinetic models with a large number of nonconvexities of different nature.
Given all these issues, it is hardly surprising that linear stoichiometric models have emerged as the most popular tool to analyze genomewide metabolic networks using optimization techniques. Linear optimization problems can be solved using very fast and efficient algorithms [15, 16] that are implemented in almost every kind of computer, ranging from laptops to cloud computing centers. In addition, such models require a relatively small amount of information.
The possibility of condensing information about a very large network in a compact form enabled stoichiometric models to provide interesting insights in many different cases. However, the apparent simplicity in building and analyzing stoichiometric models comes at the cost of neglecting regulatory signals, metabolite levels and dynamic constraints. Accounting for these features in a dynamic way requires using more detailed, nonlinear, mathematical models [17, 18].
These models go a step further than stoichiometric models by incorporating regulatory influences through a set of ordinary differential equations that can account for the system's dynamics. Building such models is often impossible because the appropriate functional form that needs to be used to describe the dynamical behavior of specific processes is in general unknown. Modeling strategies based on systematic approximated kinetic representations, such as powerlaws [19–22], Saturating and Cooperative [23], or convenience kinetics [24], sidestep this issue by providing uniform forms that are guaranteed to be accurate over a range of conditions and reduce the amount of information required to build the models. Because of the regularity in the form of the mathematical function, models based on approximate formalisms can be automatically built from the reaction scheme of the target system. The model parameters can subsequently be estimated from experimental data using different procedures [13, 25].
Although building and analyzing of comprehensive genomewide detailed models is still not viable in most cases (see however [26, 27]), developing ways to extend large scale optimization analysis to larger and more realistic nonlinear kinetic models is an important part of the future of systems biology [18]. In fact, the optimization of certain types of nonlinear problems can already be solved very efficiently and geometric programming problems with up to 1,000 variables and 10,000 constraints can be solved in minutes on a personal computer.
Efficient global optimization techniques are available for powerlaw models [1, 28–30], either in Ssystem form or in Generalized Mass Action (GMA) form (for a review see [31]). In the case of Ssystem models, a simple logarithmic transformation brings the model to a linear form [1]. In the case of GMA models, the problem can be efficiently solved using branchandbound [28, 32] and outerapproximation techniques [29, 30].
The usefulness of the global optimization techniques developed for GMA models has been shown in the analysis of the adaptive response of yeast to heat shock [29, 33]. In essence, starting with a GMA model and considering a set of constraints on flux and metabolite values, we can obtain: (i) The pattern of enzyme activities that maximizes a given objective, (ii) The region of feasible changes in enzyme activities so that the model fulfills a set of constraints on fluxes, metabolites, maximum allowable change in activity, etc., and (iii) A heat map of how the objective function changes within the feasible region. These results share some similarities with those produced with stoichiometric models, but incorporate many additional features.
Based on ideas similar to those that led to the development of the powerlaw formalism, Sorribas et al. [23] proposed a new Saturable and Cooperative (SC) formalism, that extends the powerlaw representation to include cooperativity and saturation. Although models built using this new formalism loses some of the simplicity inherent to the analysis of Ssystems and GMA models, they tend to be accurate over a wider range of conditions than both the SSystem and GMA representations [23]. Thus, it is important to enlarge the scope of global optimization methods developed for powerlaw models in order to deal with the SC formalism and analyze under which situations the later models behave better than the former.
Optimization of SC models faces a number of practical problems common to kinetic nonlinear models [34, 35]. To sidestep these problems, and in order to be able to use global optimization methods developed for powerlaw models, we will use a technique called recasting. Recasting permits the exact transformation of a continuous nonlinear model with an arbitrary form into a canonical GMA model [36, 37]. This transformation is typically performed by increasing the number of variables of the original model. Through this technique, arbitrary nonlinear models can be represented using a canonical form such as GMA or Ssystem that can be used for simulation and optimization purposes, which opens the door for effectively extending the optimization and feasibility analysis originally devised for GMA models to other detailed kinetic models.
In this paper, and as a first step to define a framework for optimization of nonlinear models with arbitrary form and extend FBA and related approaches to detailed kinetic models, we shall show the practical utility of recasting SC models into GMA models for optimization purposes. This technique is similar to the symbolic reformulation algorithm proposed by Smith and Pantelides [38]. Our method, however, focuses on obtaining a powerlaw representation that greatly facilitates global optimization, instead of continuing with the recasting until converting the model to a standard form containing linear constraints and a set of nonlinearities corresponding to bilinear product, linear fractional, simple exponentiation and univariate function terms. After recasting the model to the canonical form, we can apply any of the optimization strategies we have presented for GMA models [29, 32] to obtain the global optimum of the original SC problem.
2 Results
2.1 Global optimization of nonlinear models through recasting
In these equations, k_{ r }, r = 1,.., 6 is an auxiliary variable used to model changes in the enzyme activity. At the basal level, k_{ r } = 1 for all the reactions. During the optimization tasks, it is possible to limit the maximum change in gene expression by imposing a maximum allowable change in k_{ r }.
We shall now address the following questions:
(i) To what extent can general purpose global optimization methods be applied to SC models?, (ii) Given that a SC model can be recast as a GMA (rGMA), is this useful for optimization of the original SC model?, (iii) Are the results obtained with the rGMA equivalent to the results of the original SC model?, and (iv) What are the practical advantages of optimizing a rGMA model?.
2.2 Optimization goals
In order to address the questions posed at the end of the previous section we shall define the following optimizations tasks (note that changes in enzyme activities and metabolite concentrations are constrained between 0.2 ≤ k_{ r } ≤ 5.0 and 0.1 ≤ X_{ i } ≤ 10.0 respectively in all the instances unless otherwise specified):

O1: What is the optimal pattern of changes in enzyme activities that maximizes the objective function in the new steadystate for a fixed value of X_{5}?

O2: What is the optimal pattern of changes in enzyme activities that maximizes the objective function in the new steadystate for a fixed value of X_{5} considering a maximum allowable variation of 10% in the steadystate values of the intermediaries?

O3: What is the optimal pattern of changes in enzyme activities that maximizes the objective function in the new steadystate for a fixed value of X_{5} considering changes in the output flux from X_{4} of less than 10% with respect to its reference value?

O4: What is the best set of changes, assuming that we can only manipulate three enzymes, that maximizes the objective function in the new steadystate for a fixed value of X_{5} considering a maximum variation of 10% in the steadystate values of the intermediaries?
Two different objective functions (OF), steadystate concentration of X_{3} and flux v_{4}, have been considered for each optimization case, except for O3. This latter case has been optimized in terms only of the first objective (i.e., steadystate concentration of X_{3}), because limits on v_{4} are already included in the formulation of the optimization problem.
2.3 Global optimization of SC models using BARON
Results for the maximization of X_{3} and v_{4} and optimization goals O1O4 using BARON v.8.1.5. for a tolerance of 0.2%.
O  k _{1}  k _{2}  k _{3}  k _{4}  k _{5}  k _{6}  X _{3}  OG (%)  CPU (s) 

1  0.26  5.00  4.97  0.20  0.20  0.54  8.30  0.20  136.17 
2  0.20  0.24  0.22  0.20  0.21  0.20  1.10  0.00  0.06 
3  0.60  5.00  5.00  0.53  0.20  0.27  5.39  0.20  96.39 
4  0.99  1.15  1.00  0.96  1.00  1.00  1.10  0.00  1.42 
O  k _{ 1 }  k _{ 2 }  k _{ 3 }  k _{ 4 }  k _{ 5 }  k _{ 6 }  v _{ 4 }  OG (%)  CPU (s) 
1  4.61  5.00  5.00  5.00  0.72  1.20  37.40  0.20  157.83 
2  3.22  3.73  5.00  4.99  0.21  0.22  31.33  0.00  1.67 
3  0.88  0.94  0.88  0.96  0.23  3.00  6.60  0.00  10.53 
4  1.16  1.00  1.34  1.34  1.00  1.00  7.61  0.00  3.61 
2.4 Recasting SC models into GMA models
with appropriate initial conditions ${X}_{j}\left(0\right)={X}_{{j}_{0}}$ and ${z}_{rj}\left(0\right)={K}_{rj}+{X}_{{j}_{0}}^{{n}_{rj}}$.
2.5 Steadystate optimization of SC models through recasting
Once the problem has been recast into a rGMA, its mathematical structure can be exploited in order to improve the efficiency of the solution procedure, as demonstrated by the authors in previous works. This problem has a GMA form except for the auxiliary constraint 5b, which is required to recast the SC into the rGMA. This constraint can be easily handled by means of relaxation techniques and exponential transformations similar to those used by the authors in their global optimization algorithms for pure GMA models [32, 33]. In particular, two algorithms were developed for the global optimization of GMA models: a customized outerapproximation (OA, [30]) and a tailored spatial branchandbound (sBB, [32]). The authors showed that the numerical performance of these methods depends on the specific problem being solved, and that none of them is clearly better than the other one. Here, we use the OA algorithm to solve 6, as this method proved to be faster than sBB for problems of smaller size ([32]). Again, the main body of the algorithm was coded in GAMS 23.0.2, using CPLEX 11.2.1 as MILP solver for the master subproblems and CONOPT 3.14 s as NLP solver for the slave subproblems of the algorithm. For a fair comparison, we also set a tolerance of 0.2%, the same as when using BARON.
Results for the maximization of X_{3} and v_{4} using the rGMA model and optimization goals O1O4 using the customized OA for a tolerance of 0.2%.
O  k _{1}  k _{2}  k _{3}  k _{4}  k _{5}  k _{6}  X _{3}  OG (%)  CPU (s) 

1  0.26  5.00  5.00  0.20  0.20  0.20  8.30  0.20  2.94 
2  0.21  0.22  0.21  0.20  0.20  0.20  1.10  0.00  0.06 
3  0.60  5.00  5.00  0.53  0.20  0.24  5.40  0.13  2.35 
4  1.00  1.05  0.97  0.92  1.00  1.00  1.10  0.00  0.23 
O  k _{ 1 }  k _{ 2 }  k _{ 3 }  k _{ 4 }  k _{ 5 }  k _{ 6 }  v _{ 4 }  OG (%)  CPU (s) 
1  3.96  5.00  5.00  5.00  0.20  2.99  37.47  0.00  0.16 
2  3.22  3.55  5.00  4.99  0.20  0.21  31.33  0.17  0.66 
3  0.68  1.79  1.12  1.27  0.20  0.21  6.60  0.00  0.12 
4  1.16  1.00  1.34  1.34  1.00  1.00  7.61  0.11  1.98 
Note that the objective function values obtained with the SC and rGMA models only differ within the tolerance imposed. In some cases, discrepancies regarding the enzymatic profiles calculated are observed mainly due to the system's structure, that is, to the fact that the problem contains multiple solutions attaining the same performance in terms of objective function value but involving different enzymatic configurations, as discussed in section 2.3.
To further investigate this issue, we apply the multisolution capability of BARON to the rGMA model (Figure 2). Again, different equivalent optima are obtained, but this time the dispersion of the equivalent solutions generated for a given case tend to concentrate either in the center or in the extremes of the region containing the solutions with the same objective function value calculated with the SC model.
Results (objective function) of the optimization of case O1 v_{4} for specific regions of k_{2} and k_{5} obtained with BARON for the SC model.
k_{5}/k_{2}  1  2  3  4  5  6  7  8 

8  36.50  36.71  36.90  37.08  37.24  37.37  37.47  37.47 
7  36.62  36.83  37.02  37.19  37.34  37.46  37.47  37.47 
6  36.75  36.95  37.14  37.31  37.44  37.47  37.47  37.47 
5  36.88  37.08  37.26  37.41  37.47  37.47  37.47  37.47 
4  37.02  37.21  37.38  37.47  37.47  37.47  37.47  37.47 
3  37.15  37.34  37.47  37.47  37.47  37.47  37.47  37.47 
2  37.29  37.46  37.47  37.47  37.47  37.47  37.47  37.47 
1  37.43  37.47  37.47  37.47  37.47  37.47  37.47  37.47 
Results (objective function) of the optimization of case O1v_{4} for specific regions of k_{2} and k_{5} obtained with the customized OA for the rGMA model.
k_{5}k_{2}  1  2  3  4  5  6  7  8 

8  36.50  36.71  36.90  37.08  37.24  37.37  37.47  37.47 
7  36.62  36.83  37.02  37.19  37.34  37.46  37.47  37.47 
6  36.75  36.95  37.14  37.31  37.44  37.47  37.47  37.47 
5  36.88  37.08  37.26  37.41  37.47  37.47  37.47  37.47 
4  37.02  37.21  37.38  37.47  37.47  37.47  37.47  37.47 
3  37.15  37.34  37.47  37.47  37.47  37.47  37.47  37.47 
2  37.29  37.46  37.47  37.47  37.47  37.47  37.47  37.47 
1  37.43  37.47  37.47  37.47  37.47  37.47  37.47  37.47 
Results (CPU time in seconds) of the optimization of case O1 v_{4} for specific regions of k_{2} and k_{5} obtained with BARON for the SC model.
k_{5}/k_{2}  1  2  3  4  5  6  7  8 

8  212.53  308.53  185.64  201.80  222.30  201.53  139.16  178.31 
7  194.81  161.16  215.80  196.81  344.73  243.02  0.03  174.81 
6  234.30  203.75  147.08  180.69  328.34  254.42  304.11  280.53 
5  212.08  282.41  329.33  237.34  208.02  292.27  200.00  154.62 
4  288.00  160.14  92.94  235.80  172.69  147.14  56.11  150.28 
3  125.56  111.17  150.27  187.52  337.97  158.16  112.66  264.12 
2  239.70  190.59  100.03  138.47  106.38  205.14  119.39  246.34 
1  140.42  102.12  80.45  21.69  73.12  96.61  89.94  80.03 
Results (CPU time in seconds) of the optimization of case O1v_{4} for specific regions of k_{2} and k_{5} obtained with the customized OA for the rGMA model.
k_{5}/k_{2}  1  2  3  4  5  6  7  8 

8  0.13  0.27  0.23  0.18  0.17  0.19  0.28  0.28 
7  0.26  0.28  0.28  0.26  0.28  0.23  0.32  0.25 
6  0.32  0.30  0.28  0.28  0.27  0.23  0.19  0.25 
5  0.31  0.21  0.25  0.25  0.26  0.28  0.27  0.29 
4  0.25  0.27  0.32  0.30  0.25  0.27  0.26  0.28 
3  0.20  0.22  0.28  0.28  0.29  0.30  0.19  0.53 
2  0.28  0.25  0.19  0.19  0.22  0.17  0.30  0.25 
1  0.23  0.24  0.26  0.27  0.23  0.21  0.24  0.31 
2.6 Difficult optimization tasks can be solved via recasting
The optimization task of interest being:

O5: Which is the optimal pattern of changes in enzyme activities that maximize v_{6} in the new steadystate for a fixed value of X_{5} and considering the following constraints?$\begin{array}{c}\hfill 0.3\le {X}_{1}\le 30\hfill \\ \hfill 0.1\le {X}_{2}\le 10\hfill \\ \hfill 0.1\le {X}_{3}\le 10\hfill \\ \hfill 0.6\le {X}_{4}\le 50\hfill \\ \hfill 0.1\le {k}_{r}\le 20\hfill & \hfill r=1,...,p\hfill \end{array}$(9)
Results of the optimization of model 8 with BARON (SC model) and the customized OA (rGMA model).
Solver  k _{1}  k _{2}  k _{3}  k _{4}  k _{5}  k _{6}  OF  OG (%)  CPU (s) 

BARON (SC)  6.24  5.16  0.46  0.6  8.46  9.09  60.36  45.18  3600 
OA (rGMA)  6.25  5.17  0.45  0.6  8.44  9.1  60.46  2.18  10.95 
3 Discussion
While experimental tools to manipulate gene expression are already available, there is no established set of guidelines on how these tools can be used to achieve a certain goal. So far, two main difficulties have prevented model driven optimization from becoming a standard in providing such guidelines: (i) the lack of information to build detailed kinetic models and (ii) the computational difficulties that arise upon the optimization of such models. The latter can be exemplified by the application of mixed integer nonlinear optimization techniques (MINLP) in the context of kinetic models presented in [34, 35]. In such cases, the optimization task showed to be computationally very demanding and global optimality could not be guaranteed in many cases. We propose that using models with a standardized structure may offer a solution to both problems. On one hand, approximate kinetics, such as the SC formalism, can provide very accurate approximations and retain key features of the system like saturation and cooperativity. On the other hand, these formalisms can be automatically recast into GMA form and using efficient global optimization methods developed specifically for this canonical representation. Although this technique will certainly have limitations, our previous results indicate that it can be applied to models of moderate complexity without major problems [32]. Optimization of GMA models comprising up to 60 reactions and 40 metabolites offer no limitation to our technique. We have shown how these methods can be easily used to optimize SC via recasting into rGMA models while still being quite efficient.
Our results can be of particular interest for dealing with multicriteria optimization on realistic models. This kind of problems are relevant when exploring the adaptive response to changing conditions, were conflictive goals may be at play [39, 40]. Particularly, we should notice that several multiobjective optimization techniques, such as the weighted sum or epsilon constraint methods [41] are based on solving a set of auxiliary singleobjective problems. These approaches could directly benefit from the numerical advances presented in this work. This kind of problems are relevant when exploring the adaptive response to changing conditions, were conflictive goals may be on play [39, 40]. The highly efficient OA algorithm applied to rGMA models provide a practical way for extending multicriteria optimization methods, for instance as used in [39], to nonlinear kinetic models. It is in principle possible to make use of methods such as ours to analyze the optimality of large scale dynamic systems much in the same way that Flux Balance Analysis can be applied to analyze the stoichiometry of an organism on a genomic scale. To make this possible, however, extensive experimental and modeling efforts would be required to characterize the most important properties of the involved processes. In fact, we anticipate that practical limitations to apply the techniques presented here in solving larger problems will be dominated by the lack of information about the component processes and metabolites rather than by the technical capacity of the optimization technique presented here. Although a complete kinetic characterization of the processes in a complete metabolic network may yet be far, information on kinetic orders and saturation fractions is easier to obtain. In this context, the SC formalism provides a sound approximation that results in a mathematical representation useful for simulation and optimization through recasting.
4 Conclusions
We expect that the possibility of building models using nonlinear approximate formalisms and of subsequently optimizing these models will trigger interest in the experimental characterization of the components of cellular metabolism. After the genomic explosion, we need to step back and begin to measure enzyme activities, metabolite levels, and regulatory signals on a larger scale than we used to do before, if we want to understand the emergence of the dynamic properties of biological systems and to be able to develop successful biotechnological applications.
5 Methods
5.1 Modelling strategies
The process of model building and optimization can be used to understand how a system should be changed in order to achieve specific biotechnological goals or how the same system has evolved in order to more efficiently execute a given biological function. Different tradeoffs are considered during the modeling process. On the one hand, one wants to use models that are as simple as possible to guarantee numerical tractability. Unfortunately simplifications may lead to models whose accuracy is only ensured for a limited range of physiological conditions. On the other hand, models that are very detailed and accurate over a wide range of physiological conditions are typically more difficult to analyze and optimize. Needless to say, the type of modeling strategy and the model one chooses to implement have a large impact on the results of the analysis. The most widely used strategies in the context of optimization are: (1) Stoichiometric models, (2) Kinetic models, and (3) Approximated models.
μ_{ ir } stands for the stoichiometry of each metabolite X_{ i } in each reaction r with respect to metabolite i and can be derived from the reaction scheme.
This system of equations is solved under different assumptions. A typical problem is that of understanding the effect of knocking out different genes from the system. This analysis can be performed by setting v_{ r } = 0 for the process(es) that depend on the product of the genes that are knocked out. Once these constraints are set, linear optimization techniques can be used to identify the region of the variable space that satisfies the steady state and optimizes at the same time a set of specific measurable aspects of the systems [42–44]. It must be noted that FBA analysis of Eqn. (11) does not account for the regulatory effects that can result from gene knockout and it cannot be used to predict changes in metabolic concentrations and temporal responses. Thus, optimization constraints are limited to steadystate fluxes [15].
To overcome these limitations, we must use more complex kinetic models where the effect of changing the values of the variables on the fluxes is taken into account. This requires defining a functional form for each v_{ r } in Eqn. (10). Often, this functional form is drawn from a number of classical enzyme kinetic ratelaws. As a result, we use an approximate expression for the kinetic behavior of each elementary process whose form depends on the underlying mechanism of the process. The reason for this is that the classical rate laws are rational functions of the variables and they are built upon different types of simplifying assumptions on the detailed mechanism of the reactions. Such assumptions range from considering that the elementary chemical steps of the catalytic process occur at very different timescales to assuming that the concentration of the catalyst and of the reactants differ in orders of magnitude. Thus, rate laws such as the popular MichaelisMenten are approximations to the actual mechanism in specific conditions. However, more often than not, one does not have enough information to judge if such conditions meet those one is trying to model. Thus, using rational enzyme kinetics in models lacks a sound theoretical ground. In fact, within the complex architecture of the intracellular milieu, many of the assumptions that justify these classical ratelaws may not hold [45–47]. Even in the best case scenario where a detailed kinetic model using classical enzyme kinetics can be derived and numerically identified, it may be hard to globally optimize that model using general purpose algorithms. As we will show here, available optimization techniques may fail to solve fairly trivial optimization tasks even in simple models. These numerical difficulties can be overcome by defining reformulated models based on canonical representations that are easier to handle using customized global optimization algorithms devised for specific canonical functional forms.
The parameters in these representations have a clear physical interpretation. Kinetic orders, the exponents in the powerlaws, are local sensitivities of the fluxes, either individual (f_{ rj } for v_{ r }) or aggregated (g_{ ij } for ${V}_{i}^{+}$ and h_{ ij } for ${V}_{i}^{}$), with respect to X_{ j }. Rateconstants (α_{ i }, β_{ i } and γ_{ r }) are parameters that are computed so that the flux in the model at steady state is equal to the operating flux at the operating point for the metabolites. Parameter estimation techniques have been developed so that powerlaw parameters can be calculated from experimental measurements [13]. It should also be noted that the use of estimation procedures (i.e., leastsquares), alternate regression or similar procedures to estimate powerlaw parameters from dynamic curves lead to a powerlaw representation that is no longer local according to the classical definition [48–50]. Those models may, by definition, slightly improve their accuracy over strictly local models.
where v_{ r0 } = v_{ r }(X_{10},.., X_{ n0 },... X_{(n + m)0}) and V_{ rj } is either the limit velocity (saturation) when X_{ j } → ∞ if n_{ rj } > 0, or the limit velocity when X_{ j } → 0 if n_{ rj } < 0.
Using SC models for global optimization can raise some numerical issues. These difficulties can be avoided to a large extent by recasting SC models into a canonical GMA model, through the introduction of auxiliary variables, as will be shown in the next section.
5.2 Recasting nonlinear models into powerlaw canonical models by increasing the number of variables
Nonlinear models can be exactly recast into GMA or Ssystem models through the use of auxiliary variables [36]. As a result, the final model is an exact representation of the original model, written in a canonical form. In other words, the resulting GMA model is not an approximation to the original model: it is an exact replica of it. To avoid confusion, hereafter, we refer to a GMA model that exactly recasts another as an rGMA model.
in which X_{3} is an externally fixed variable.
with initial conditions ${X}_{1}\left(0\right)={X}_{{1}_{0}}$ and ${X}_{2}\left(0\right)={X}_{{2}_{0}}$.
with initial conditions ${X}_{4}\left(0\right)={K}_{1}\left(1+Ki\mathsf{\text{/}}{X}_{{2}_{0}}\right)+{X}_{{3}_{0}}$, ${X}_{5}\left(0\right)={K}_{2}+{X}_{{1}_{0}}$, and ${X}_{6}\left(0\right)={K}_{3}^{2}+{X}_{{2}_{0}}^{2}$.
The resulting rGMA model (2223) is an exact representation of model in (20). Hence, for a set of appropriate initial conditions, the simulation of the dynamic response using either the model recast as a rGMA or the original model will produce the same trajectory. In principle, any nonlinear model can be recast into a rGMA following a similar procedure [36]. This can be extremely useful, because it allows for the application of tailored global optimization procedures originally devised for GMA models [28–30, 32, 51, 52] to generic nonlinear models.
Declarations
Acknowledgements
AS is funded by MICINN (Spain) (BFU20080196). RA is partially supported by MICINN (Spain) through Grants BFU200762772/BMC and BFU201017704). AS and RA are members of the 2009SGR809 research group of the Generalitat de Catalunya. GGG and CP acknowledges support from the Spanish Ministry of Science and Innovation (Projects DPI200804099 and CTQ200914420C0201) and the Spanish Ministry of External Affairs and Cooperation (Projects A/023551/09 and A/031707/10).
Authors’ Affiliations
References
 Voit EO: Optimization in integrated biochemical systems. Biotechnol Bioeng. 1992, 40 (5): 57282. 10.1002/bit.260400504View ArticlePubMedGoogle Scholar
 AlvarezVasquez F, GonzalezAlcon C, Torres NV: Metabolism of citric acid production by Aspergillus niger: model definition, steadystate analysis and constrained optimization of citric acid production rate. Biotechnol Bioeng. 2000, 70: 82108. 10.1002/10970290(20001005)70:1<82::AIDBIT10>3.0.CO;2VView ArticlePubMedGoogle Scholar
 Banga JR: Optimization in computational systems biology. BMC Syst Biol. 2008, 2: 47 10.1186/17520509247PubMed CentralView ArticlePubMedGoogle Scholar
 RodriguezPrados JC, de Atauri P, Maury J, Ortega F, Portais JC, Chassagnole C, Acerenza L, Lindley ND, Cascante M: In silico strategy to rationally engineer metabolite production: A case study for threonine in Escherichia coli. Biotechnology and bioengineering. 2009, 103 (3): 609620. 10.1002/bit.22271View ArticlePubMedGoogle Scholar
 Scheer M, Grote A, Chang A, Schomburg I, Munaretto C, Rother M, Stelzer M, Thiele J, Schomburg D: BRENDA, the enzyme information system in 2011. Nucleic Acids Res. 2011, 39: 670676. 10.1093/nar/gkq1089.View ArticleGoogle Scholar
 Rojas I, Golebiewski M, Kania R, Krebs O, Mir S, Weidemann A, Wittig U: SABIORK: a database for biochemical reactions and their kinetics. BMC Systems Biology. 2007, 1: S610.1186/175205091S1S6.View ArticleGoogle Scholar
 Shiraishi F, Savageau MA: The tricarboxylic acid cycle in Dictyostelium discoideum. I. Formulation of alternative kinetic representations. The Journal of biological chemistry. 1992, 267 (32): 2291222918.PubMedGoogle Scholar
 Shiraishi F, Savageau MA: The tricarboxylic acid cycle in Dictyostelium discoideum. II. Evaluation of model consistency and robustness. The Journal of biological chemistry. 1992, 267 (32): 2291922925.PubMedGoogle Scholar
 Shiraishi F, Savageau MA: The tricarboxylic acid cycle in Dictyostelium discoideum. III. Analysis of steady state and dynamic behavior. The Journal of biological chemistry. 1992, 267 (32): 2292622933.PubMedGoogle Scholar
 Shiraishi F, Savageau MA: The tricarboxylic acid cycle in Dictyostelium discoideum. IV. Resolution of discrepancies between alternative methods of analysis. The Journal of biological chemistry. 1992, 267 (32): 2293422943.PubMedGoogle Scholar
 Fonseca LL, Sanchez C, Santos H, Voit EO: Complex coordination of multiscale cellular responses to environmental stress. Molecular bioSystems. 2010,Google Scholar
 Sorribas A, Cascante M: Structure identifiability in metabolic pathways: parameter estimation in models based on the powerlaw formalism. The Biochemical journal. 1994, 298 (Pt 2): 303311.PubMed CentralView ArticlePubMedGoogle Scholar
 Chou IC, Voit EO: Recent developments in parameter estimation and structure identification of biochemical and genomic systems. Mathematical Biosciences. 2009, 219 (2): 5783. 10.1016/j.mbs.2009.03.002PubMed CentralView ArticlePubMedGoogle Scholar
 Grossmann I, Biegler LT: Part II. Future perspective on optimization. Computers and Chemical Engineering. 2004, 28: 11931218.View ArticleGoogle Scholar
 Terzer M, Maynard ND, Covert MW, Stelling J: Genomescale metabolic networks. Wiley interdisciplinary reviews. Systems biology and medicine. 2009, 1 (3): 285297. 10.1002/wsbm.37View ArticlePubMedGoogle Scholar
 Gianchandani EP, Chavali AK, Papin JA: The application of flux balance analysis in systems biology. Wiley interdisciplinary reviews. Systems biology and medicine. 2010, 2 (3): 372382. [JID: 101516550; ppublish], 10.1002/wsbm.60View ArticlePubMedGoogle Scholar
 Voit EO: Design principles and operating principles: the yin and yang of optimal functioning. Math Biosci. 2003, 182: 8192. 10.1016/S00255564(02)001621View ArticlePubMedGoogle Scholar
 Jamshidi N, Palsson BO: Formulating genomescale kinetic models in the postgenome era. Molecular systems biology. 2008, 4: 171PubMed CentralView ArticlePubMedGoogle Scholar
 Savageau MA: Biochemical systems analysis. I. Some mathematical properties of the rate law for the component enzymatic reactions. Journal of theoretical biology. 1969, 25 (3): 365369. 10.1016/S00225193(69)800263View ArticlePubMedGoogle Scholar
 Savageau MA: Biochemical systems analysis. II. The steadystate solutions for an npool system using a powerlaw approximation. Journal of theoretical biology. 1969, 25 (3): 370379. 10.1016/S00225193(69)800275View ArticlePubMedGoogle Scholar
 Savageau MA: Biochemical systems analysis. 3. Dynamic solutions using a powerlaw approximation. Journal of theoretical biology. 1970, 26 (2): 215226. 10.1016/S00225193(70)800133View ArticlePubMedGoogle Scholar
 Savageau MA: Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology. 1976, Reading, Mass.: AddisonWesley,Google Scholar
 Sorribas A, HernandezBermejo B, Vilaprinyo E, Alves R: Cooperativity and saturation in biochemical networks: a saturable formalism using Taylor series approximations. Biotechnology and bioengineering. 2007, 97 (5): 12591277. 10.1002/bit.21316View ArticlePubMedGoogle Scholar
 Liebermeister W, Klipp E: Bringing metabolic networks to life: convenience rate law and thermodynamic constraints. Theoretical biology & medical modelling. 2006, 3: 41 10.1186/17424682341View ArticleGoogle Scholar
 Goel G, Chou IC, Voit EO: System Estimation from Metabolic Time Series Data. Bioinformatics (Oxford, England). 2008,Google Scholar
 Ni T, Savageau M: Model assessment and refinement using strategies from biochemical systems theory: Application to metabolism in human red blood cells. Journal of Theoretical Biology. 1996, 179 (4): 329368. 10.1006/jtbi.1996.0072View ArticlePubMedGoogle Scholar
 Arkin A, Ross J, McAdams H: Stochastic kinetic analysis of developmental pathway bifurcation in phage lambdainfected Escherichia coli cells. Genetics. 1998, 149 (4): 16331648.PubMed CentralPubMedGoogle Scholar
 Polisetty PK, Gatzke EP, Voit EO: Yield optimization of regulated metabolic systems using deterministic branchandreduce methods. Biotechnol Bioeng. 2008, 99 (5): 115469. 10.1002/bit.21679View ArticlePubMedGoogle Scholar
 GuillénGosálbez G, Sorribas A: Identifying quantitative operation principles in metabolic pathways: A systematic method for searching feasible enzyme activity patterns leading to cellular adaptive responses. BMC Bioinformatics. 2009, 10:Google Scholar
 GuillénGosálbez G, Pozo C, Jiménez L, Sorribas A: A global optimization strategy to identify quantitative design principles for gene expression in yeast adaptation to heat shock. Computer Aided Chemical Engineering. 2009, 26: 10451050.View ArticleGoogle Scholar
 Voit EO: Computational Analysis of Biochemical Systems. A Practical Guide for Biochemists and Molecular Biologists. 2000, Cambridge, U.K.: Cambridge University Press,Google Scholar
 Pozo C, GuillénGosálbez G, Sorribas A, Jiménez L: A Spatial BranchandBound Framework for the Global Optimization of Kinetic Models of Metabolic Networks. Industrial and Engineering Chemistry Research. 2010,Google Scholar
 Sorribas A, Pozo C, Vilaprinyo E, GuillénGosálbez G, Jiménez L, Alves R: Optimization and evolution in metabolic pathways: global optimization techniques in Generalized Mass Action models. Journal of Biotechnology. 2010, 149 (3): 141153. 10.1016/j.jbiotec.2010.01.026View ArticlePubMedGoogle Scholar
 Chassagnole C, NoisommitRizzi N, Schmid J, Mauch K, Reuss M: Dynamic modeling of the central carbon metabolism of Escherichia coli. Biotechnol Bioeng. 2002, 79: 5373. 10.1002/bit.10288View ArticlePubMedGoogle Scholar
 Nikolaev E: The elucidation of metabolic pathways and their improvements using stable optimization of largescale kinetic models of cellular systems. Metab Eng. 2010, 12: 2638. 10.1016/j.ymben.2009.08.010View ArticlePubMedGoogle Scholar
 Savageau MA, Voit EO: Recasting nonlinear differential equations as Ssystems: a canonical nonlinear form. Mathematical Biosciences. 1987, 87: 83115. 10.1016/00255564(87)900356.View ArticleGoogle Scholar
 Voit EO: Recasting nonlinear models as Ssystems. Mathematical and Computer Modelling. 1988, 11 (C): 140145.View ArticleGoogle Scholar
 Smith EMB, Pantelides CC: A symbolic reformulation/spatial branchandbound algorithm for the global optimisation of nonconvex MINLPs. Computers and Chemical Engineering. 1999, 23 (45): 457478. 10.1016/S00981354(98)002865.View ArticleGoogle Scholar
 Vera J, de Atauri P, Cascante M, Torres NV: Multicriteria optimization of biochemical systems by linear programming: application to production of ethanol by Saccharomyces cerevisiae. Biotechnol Bioeng. 2003, 83 (3): 33543. 10.1002/bit.10676View ArticlePubMedGoogle Scholar
 Vilaprinyo E, Alves R, Sorribas A: Use of physiological constraints to identify quantitative design principles for gene expression in yeast adaptation to heat shock. BMC Bioinformatics. 2006, 7: 184 10.1186/147121057184PubMed CentralView ArticlePubMedGoogle Scholar
 Ehrgott M: Multicriteria Optimization. 2005, Springer,Google Scholar
 Famili I, Forster J, Nielsen J, Palsson BO: Saccharomyces cerevisiae phenotypes can be predicted by using constraintbased analysis of a genomescale reconstructed metabolic network. Proc Natl Acad Sci USA. 2003, 100 (23): 131349. 10.1073/pnas.2235812100PubMed CentralView ArticlePubMedGoogle Scholar
 Famili I, Mahadevan R, Palsson BO: kCone analysis: determining all candidate values for kinetic parameters on a network scale. Biophys J. 2005, 88 (3): 161625. 10.1529/biophysj.104.050385PubMed CentralView ArticlePubMedGoogle Scholar
 Price ND, Papin JA, Schilling CH, Palsson BO: Genomescale microbial in silico models: the constraintsbased approach. Trends Biotechnol. 2003, 21 (4): 1629. 10.1016/S01677799(03)000301View ArticlePubMedGoogle Scholar
 Savageau MA: Influence of fractal kinetics on molecular recognition. Journal of Molecular Recognition: JMR. 1993, 6 (4): 149157. 10.1002/jmr.300060403View ArticlePubMedGoogle Scholar
 Savageau MA: MichaelisMenten mechanism reconsidered: implications of fractal kinetics. Journal of theoretical biology. 1995, 176: 115124. 10.1006/jtbi.1995.0181View ArticlePubMedGoogle Scholar
 Savageau MA: Development of fractal kinetic theory for enzymecatalysed reactions and implications for the design of biochemical pathways. Bio Systems. 1998, 47 (12): 936. 10.1016/S03032647(98)000203View ArticlePubMedGoogle Scholar
 HernandezBermejo B, Fairen V, Sorribas A: Powerlaw modeling based on leastsquares minimization criteria. Mathematical biosciences. 1999, 161 (12): 8394. 10.1016/S00255564(99)000358View ArticlePubMedGoogle Scholar
 HernandezBermejo B, Fairen V, Sorribas A: Powerlaw modeling based on leastsquares criteria: consequences for system analysis and simulation. Mathematical biosciences. 2000, 167 (2): 87107. 10.1016/S00255564(00)000390View ArticlePubMedGoogle Scholar
 Alves R, Vilaprinyo E, HernandezBermejo B, Sorribas A: Mathematical formalisms based on approximated kinetic representations for modeling genetic and metabolic pathways. Biotechnology and Genetic Engineering Reviews. 2008, 25: 140. 10.5661/bger251View ArticlePubMedGoogle Scholar
 MarinSanguino A, Torres NV: Optimization of biochemical systems by linear programming and general mass action model representations. Math Biosci. 2003, 184 (2): 187200. 10.1016/S00255564(03)000464View ArticlePubMedGoogle Scholar
 MarinSanguino A, Voit EO, GonzalezAlcon C, Torres NV: Optimization of biotechnological systems through geometric programming. Theor Biol Med Model. 2007, 4: 38 10.1186/17424682438PubMed CentralView ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.