 Methodology article
 Open Access
 Published:
Steadystate global optimization of metabolic nonlinear dynamic models through recasting into powerlaw canonical models
BMC Systems Biology volume 5, Article number: 137 (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.
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
For a proof of concept of the difficulties of global optimizing nonlinear models and of the use of recasting for attaining practical solutions, we shall start by defining a reference biochemical network that corresponds to the reaction scheme in Figure 1. This hypothetical system has a source metabolite X_{5} and four internal metabolites. The network includes six reactions and a branch point. X_{3} acts as a feedback inhibitor of the synthesis of X_{2}, while X_{1} is an activator of the synthesis of X_{4}.
The generic model for this system is:
Each of the velocities is a nonlinear function of the involved metabolites. The SC representation, provides a systematic way for defining a functional model of this pathway. As a demonstrative example, let us suppose that the numerical model is:
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
We first address the optimization of the aforementioned model in their original SC form using state of the art global optimization techniques. The model was coded in the algebraic modeling system GAMS 23.0.2 and solved with the commercial global optimization package BARON v.8.1.5. on an Intel 1.2 GHz machine. An optimality gap (i.e., tolerance) of 0.2% was set in all the instances. As can be seen in Table 1, BARON produce results with an optimality gap (OG) below the specified tolerance.
Table 1 only shows one solution for each particular instance. However, BARON identified in each case a set of equivalent optima (i.e, solutions with the same objective function value) involving different changes in enzyme activities, which indicates that the optimization problem is somehow degenerated. This redundancy is a consequence of the system's structure and has practical implications. As an example, we have calculated some of these equivalent points for case O1v_{4} using the NumSol option of BARON (see Figure 2). In particular, a well defined triangular region containing the changes in k_{2} and k_{5}, and k_{1} and k_{2} that lead to the same objective function value is identified. Within these regions, one can decide which combination of changes should be selected based on additional cost arguments, as they all show the same performance in terms of the predefined objective function. This region could be further reduced by imposing additional constraints to the optimization.
2.4 Recasting SC models into GMA models
Any SC model can be recast into a GMA canonical model by introducing the auxiliary variables ${z}_{rj}={K}_{rj}+{X}_{j}^{{n}_{rj}}$. Substitution and differentiation generates the following recast GMA (rGMA) model:
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}}$.
For simulation purposes, model (3) is equivalent to the original SC model. As discussed in [36], a model recast into a GMA model has the same steadystate that the original nonlinear model. The steadystate equations of the rGMA model can be expressed as:
2.5 Steadystate optimization of SC models through recasting
The steadystate solutions of Eqn. (4b) satisfy also Eqn. (4a). Thus, for optimization purposes, the steadystate constraints of interest are:
According to these results, the optimization problem can be stated as:
In our reference model, we shall consider the following constraints:
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.
As can be seen in Table 2, the optimization of the rGMA formulation using our customized OA yields similar results to those obtained when BARON is applied to the original SC model. In some cases, significant reductions in computational time are attained with our OA algorithm. While BARON took a total time of 407.68 CPU seconds for solving the 8 instances, the customized OA algorithm solved the same problems in 8.5 CPU seconds.
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.
The region illustrated in Figure 2 should not be misunderstood as a feasibility region. In fact, solutions do exist outside this region, but they lead to worse objective function values. To further clarify this issue, we consider a grid of values for k_{2} and k_{5} in the region defined by constraints 4 ≤ k_{2} ≤ 5 and 0.2 ≤ k_{5} ≤ 0.8, and solve the optimization problem within each cell applying BARON to the SC model, and our OA to the rGMA model. Recall that these linear constraints define a region that contains that in Figure 2. The results obtained in this optimization are illustrated in Tables 3 and 4, and are exactly equal for both methods. However, the CPU time is much lower when using our OA algorithm applied to rGMA (11,811 CPU seconds for generating all the points with BARON applied to the SC model vs 17 CPU seconds with the customized OA applied to the rGMA model; as shown in Tables 5 and 6).
2.6 Difficult optimization tasks can be solved via recasting
The reference model can be optimized either by general purpose techniques or by rGMA specific methods such as the customized OA. However, even with this simple example, we may encounter instances that are hard to solve using standard techniques. Consider, for instance, the same reaction scheme as before but this time with the alternative parameters indicated in the following model:
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)
When BARON is employed to solve this case using the native SC form, it cannot reduce the optimality gap below the specified tolerance after 1 hour of CPU time. In contrast, when the model is recast into its rGMA form and our OA method is applied, the global optimum can be determined with an optimality gap of 2% in 10.95 seconds (see Table 7). This illustrates both, the utility of using the rGMA as a canonical form for dealing with the optimization of SC models, and the computational efficiency of our global optimization methods specifically designed to take advantage of the mathematical structure of the GMA.
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.
The three strategies have as a starting point a set of ordinary differential equations, in which the dependent variables or nodes are the chemical species whose dynamical behavior one is interested in studying. For a system with n dependent variables, p processes and m independent variables, the node equations are written as follows:
μ_{ 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.
At this stage, the various strategies begin to differ in the way that they implement and analyze the equations. Typically, Flux balance analysis (FBA) and related techniques consider only the steady state behavior of the system, and treat v_{ r } as a variable whose value can be changed in order to optimize specific steady state constraints. To accomplish this, FBAlike methods attempt to find solutions for the following system of linear equations:
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.
As an alternative, theoretically well supported canonical representations can be derived using approximation theory. One type of such representations are powerlaw models. In a powerlaw model, each v_{ r } in Eqn. (10) is approximated as [19, 21]:
This approximation is derived at a given operating point $\left({X}_{{1}_{0}},{X}_{{2}_{0}},..,{X}_{{\left(n+m\right)}_{0}}\right)$ as a firstorder Taylor series representation of the target function in loglog space. This approximation can generate models with different representations. The two that are most commonly used are the Ssystem representation and the GMA representation. The Ssystem representation is obtained by lumping the various processes that contribute to the synthesis of a given metabolite into a global process of synthesis ${V}_{i}^{+}$ and those that contribute to the utilization of a given metabolite into a global degradation process ${V}_{i}^{}$:
Then, the aggregated processes are represented by powerlaw functions:
Alternatively, the GMA form is obtained representing each individual v_{ r } as a powerlaw:
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.
To complement the powerlaw approach, the Saturable and Cooperative (SC) formalism was introduced by Sorribas et al. [23] as an extension of the ideas that led to the powerlaw formalism. The SC representation of a given velocity is:
This representation can be obtained from a powerlaw model defined at a given operating point X_{0} = (X_{10},.., X_{(n + m)0}) through the following relationships:
Thus SC uses the same information as the powerlaw except for the new parameters p_{ rj } (saturation fractions), which are defined as:
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.
As a very simple introductory example, consider a linear pathway with two internal metabolites X_{1} and X_{2} and a source metabolite X_{3} (Figure 3). In this pathway, X_{2} is a competitive inhibitor of the synthesis of X_{1} from the source metabolite. A generic model using MichaelisMenten kinetic functions, assuming a competitive inhibition of the first reaction by X_{2}, can be written as:
in which X_{3} is an externally fixed variable.
Recasting this model as a rGMA can be done as follows. First, let us define three new variables:
We can now write the model in 20 as:
with initial conditions ${X}_{1}\left(0\right)={X}_{{1}_{0}}$ and ${X}_{2}\left(0\right)={X}_{{2}_{0}}$.
To complete the recasting we must now provide the equations that follow the change in the new variables over time. These are given by the following equations:
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.
References
 1.
Voit EO: Optimization in integrated biochemical systems. Biotechnol Bioeng. 1992, 40 (5): 57282. 10.1002/bit.260400504
 2.
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;2V
 3.
Banga JR: Optimization in computational systems biology. BMC Syst Biol. 2008, 2: 47 10.1186/17520509247
 4.
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.22271
 5.
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.
 6.
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.
 7.
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.
 8.
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.
 9.
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.
 10.
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.
 11.
Fonseca LL, Sanchez C, Santos H, Voit EO: Complex coordination of multiscale cellular responses to environmental stress. Molecular bioSystems. 2010,
 12.
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.
 13.
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.002
 14.
Grossmann I, Biegler LT: Part II. Future perspective on optimization. Computers and Chemical Engineering. 2004, 28: 11931218.
 15.
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.37
 16.
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.60
 17.
Voit EO: Design principles and operating principles: the yin and yang of optimal functioning. Math Biosci. 2003, 182: 8192. 10.1016/S00255564(02)001621
 18.
Jamshidi N, Palsson BO: Formulating genomescale kinetic models in the postgenome era. Molecular systems biology. 2008, 4: 171
 19.
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)800263
 20.
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)800275
 21.
Savageau MA: Biochemical systems analysis. 3. Dynamic solutions using a powerlaw approximation. Journal of theoretical biology. 1970, 26 (2): 215226. 10.1016/S00225193(70)800133
 22.
Savageau MA: Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology. 1976, Reading, Mass.: AddisonWesley,
 23.
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.21316
 24.
Liebermeister W, Klipp E: Bringing metabolic networks to life: convenience rate law and thermodynamic constraints. Theoretical biology & medical modelling. 2006, 3: 41 10.1186/17424682341
 25.
Goel G, Chou IC, Voit EO: System Estimation from Metabolic Time Series Data. Bioinformatics (Oxford, England). 2008,
 26.
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.0072
 27.
Arkin A, Ross J, McAdams H: Stochastic kinetic analysis of developmental pathway bifurcation in phage lambdainfected Escherichia coli cells. Genetics. 1998, 149 (4): 16331648.
 28.
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.21679
 29.
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:
 30.
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.
 31.
Voit EO: Computational Analysis of Biochemical Systems. A Practical Guide for Biochemists and Molecular Biologists. 2000, Cambridge, U.K.: Cambridge University Press,
 32.
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,
 33.
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.026
 34.
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.10288
 35.
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.010
 36.
Savageau MA, Voit EO: Recasting nonlinear differential equations as Ssystems: a canonical nonlinear form. Mathematical Biosciences. 1987, 87: 83115. 10.1016/00255564(87)900356.
 37.
Voit EO: Recasting nonlinear models as Ssystems. Mathematical and Computer Modelling. 1988, 11 (C): 140145.
 38.
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.
 39.
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.10676
 40.
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/147121057184
 41.
Ehrgott M: Multicriteria Optimization. 2005, Springer,
 42.
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.2235812100
 43.
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.050385
 44.
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)000301
 45.
Savageau MA: Influence of fractal kinetics on molecular recognition. Journal of Molecular Recognition: JMR. 1993, 6 (4): 149157. 10.1002/jmr.300060403
 46.
Savageau MA: MichaelisMenten mechanism reconsidered: implications of fractal kinetics. Journal of theoretical biology. 1995, 176: 115124. 10.1006/jtbi.1995.0181
 47.
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)000203
 48.
HernandezBermejo B, Fairen V, Sorribas A: Powerlaw modeling based on leastsquares minimization criteria. Mathematical biosciences. 1999, 161 (12): 8394. 10.1016/S00255564(99)000358
 49.
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)000390
 50.
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/bger251
 51.
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)000464
 52.
MarinSanguino A, Voit EO, GonzalezAlcon C, Torres NV: Optimization of biotechnological systems through geometric programming. Theor Biol Med Model. 2007, 4: 38 10.1186/17424682438
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).
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
AMS suggested the potential utility of recasting for optimizing nonlinear kinetic models. AS and AMS elaborate on the recasting of SC models and planned the work. CP, GGG and LJ implemented the OA algorithm and worked out the technical solution for applying it to a rGMA model. CP and GGG performed the optimization tasks. AS and RA defined the reference model and obtained the numerical parameters used in the paper. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Pozo, C., MarínSanguino, A., Alves, R. et al. Steadystate global optimization of metabolic nonlinear dynamic models through recasting into powerlaw canonical models. BMC Syst Biol 5, 137 (2011). https://doi.org/10.1186/175205095137
Received:
Accepted:
Published:
Keywords
 Global Optimization
 Flux Balance Analysis
 Optimization Task
 Global Optimization Method
 Efficient Global Optimization