 Research article
 Open Access
 Published:
Global dynamic optimization approach to predict activation in metabolic pathways
BMC Systems Biology volume 8, Article number: 1 (2014)
Abstract
Background
During the last decade, a number of authors have shown that the genetic regulation of metabolic networks may follow optimality principles. Optimal control theory has been succesfully used to compute optimal enzyme profiles considering simple metabolic pathways. However, applying this optimal control framework to more general networks (e.g. branched networks, or networks incorporating enzyme production dynamics) yields problems that are analytically intractable and/or numerically very challenging. Further, these previous studies have only considered a singleobjective framework.
Results
In this work we consider a more general multiobjective formulation and we present solutions based on recent developments in global dynamic optimization techniques. We illustrate the performance and capabilities of these techniques considering two sets of problems. First, we consider a set of singleobjective examples of increasing complexity taken from the recent literature. We analyze the multimodal character of the associated non linear optimization problems, and we also evaluate different global optimization approaches in terms of numerical robustness, efficiency and scalability. Second, we consider generalized multiobjective formulations for several examples, and we show how this framework results in more biologically meaningful results.
Conclusions
The proposed strategy was used to solve a set of singleobjective case studies related to unbranched and branched metabolic networks of different levels of complexity. All problems were successfully solved in reasonable computation times with our global dynamic optimization approach, reaching solutions which were comparable or better than those reported in previous literature. Further, we considered, for the first time, multiobjective formulations, illustrating how activation in metabolic pathways can be explained in terms of the best tradeoffs between conflicting objectives. This new methodology can be applied to metabolic networks with arbitrary topologies, nonlinear dynamics and constraints.
Background
Optimality principles have been successfully used to describe the design, organization and behavior of biological systems at different levels. Sutherland [1] argues that optimization can in fact play an even major role, allowing biology to move from just describing mechanisms to being able to predict, from first principles, how organisms should be designed. In the context of cell biology, mathematical optimization has been the underlying hypothesis in applications such as flux balance analysis [2], the analysis of activation of metabolic pathways [3–7], model identification (including parameter estimation and optimal experimental design) [8] and optimal stimulation procedures to achieve a desired biological behavior [9–11].
Here we focus on the problem of enzyme activation in metabolic networks, which has received substantial attention in the recent literature. Several authors have shown that the genetic regulation of metabolic networks may follow an optimality principle, such as the minimization of the transition time or the maximization of the production of a given metabolite. In their seminal work, Klipp et al. [4] showed how sequential gene expression appears in unbranched metabolic networks under the hypothesis of minimum transition time. The sequence matched the enzyme order in the pathway. This “justintime” activation pattern was experimentally confirmed later for the case of the aminoacid biosynthesis systems of E. coli[5]. These authors also showed that in the arginine, methionine and serine systems, the earlier the enzyme is involved in a pathway, the shorter is the response time and the higher the maximal promoter activity of the corresponding gene. Evidences of temporal distribution of enzymes were also found in the lysine biosynthesis pathway by E. coli[12]. Chechik et al. [13] proposed the notion of activity motif to systematically study the dynamic behavior of metabolic networks. As a case study they considered the transcriptional response of metabolic genes after a sudden change in environmental or nutritional condition in S. cerevisiae. These authors also found that enzymes in a metabolic chain are induced in the same order they are used in the pathway in both directions forward and backward. All these results support the idea that “justintime” activation in metabolic pathways is a widespread phenomena.
These works indicate that, starting with a mathematical model of a given metabolic network, it is in principle possible to anticipate activation profiles for specific cellular objectives solving the corresponding dynamic optimization problem. Needless to say, the results obtained must be compared with existing or new experimental measurements. In this dynamic optimization framework, the objective is to compute time varying control profiles, usually enzyme concentrations or the corresponding expression rates, so as to minimize a given cost function, such as the time needed to reach a given amount of product, subject to the system dynamics (the model) and algebraic constraints, for example in the maximum amount of enzyme available.
However solving this class of optimization problems for arbitrary metabolic networks of certain complexity is still a challenge. So far, the literature reports applications in very simple scenarios. Oyarzún et al. [6, 14] suggested the application of classical optimal control theory, giving an analytical proof of sequential activation for unbranched networks for the case the transition time is to be minimized subject to specific constraints. Almost in parallel, Bart et al. [7, 15] also considered the maximum principle of Pontryagin to solve related examples. Although successful for simple unbranched pathways, these works suggest that the use of indirects approaches from optimal control theory, such as the maximum principle of Pontryagin or the use of the HamiltonJacobiBellman equation, may not be appropriate to deal with general problems (such as those related to branched and/or large networks) due to the complexity of the resulting numerical problems.
It is also important to highlight that all these authors have studied enzyme activation in metabolic networks considering that the metabolism is optimal with respect to a singleobjective. We believe that a more realistic approach is to take more than one objective into account. Multiobjective optimization has been successfully applied in several classes of biological problems [16], including metabolic engineering [17, 18], the heatshock response of bacteria [19], and the allosteric regulation of enzymes [20]. These and other recent works [21–23] illustrate the importance of studying optimality as the tradeoff between conflicting objectives, such as economy (cost) and effectiveness (benefit), in order to obtain results with more biological meaning. It should be noted that multiobjective optimization problems do not have a unique solution, but a set of optimal solutions known as the Pareto front. All Pareto solutions are optimal but represent different tradeoffs between objectives.
Here we consider, for the first time, a multiobjective formulation to explain activation in metabolic networks. It should be noted that several authors [6, 14, 24] have considered composed objective functions, such us combinations of transition time and enzyme consumption, thus obtaining only one of the possible Pareto solutions instead of a Pareto front. In this work we adopt a general multicriteria framework and we propose the use of advanced numerical dynamic optimization techniques to study/predict enzyme activation in general pathways. The underlying idea is to combine the control vector parameterization (CVP) approach with adequate global optimization techniques. This new methodology can be applied to metabolic networks with arbitrary topologies, non linear dynamics and constraints.
We then illustrate the approach with two sets of problems. First, a set of singleobjective examples of increasing complexity taken from the literature. The multimodality of these problems is evaluated by means of multistart local deterministic methods. The need of global optimization methods arises from the nonconvexity of the general problem, due to the frequent bilinear terms for the controls, the nonlinear character of the systems dynamics, and the presence of (nonlinear) constraints [25]. As a consequence, the use of standard local optimization methods results in local solutions. To surmount these difficulties, we also present a metaheuristic approach which is compared with several other stochastic global optimization methods. Second, we consider the multiobjective formulation of several problems, illustrating how activation in metabolic pathways can be explained in terms of the best tradeoffs between conflicting objectives.
Methods
Problem formulation
The aim of the general multiobjective dynamic optimization problem is to compute the time varying control profiles and the final time (u(t), t _{ f }) so as to minimize (or maximize) a given set of cost functions (J) subject to the system dynamics and possibly algebraic constraints. Mathematically this reads:
Where:
Subject to:
where:

The costs functional (1) corresponds, for example, to the time needed to reach a given state of the system or the enzyme cost. In a singleobjective problem J[x,u]=J _{1}[x,u];

x is the vector of state variables, typically metabolite concentrations;

u is the vector of control variables: enzyme concentrations (e) or expression rates (r);

Equation (2) represents the system dynamics (dynamic mathematical model of the network);

Equations (3)–(4) represent equality and inequality path constraints, such as the total amount of enzyme available during the process;

Equations (5)–(6) represent equality and inequality point constraints, i.e. those that must be verified at a given time, for example the amount of metabolite at the end of the process;

Equation (7) corresponds to the upper and lower bounds for the control variables, for instance the minimum and maximum enzyme available through out the process.
Numerical methods
Roughly speaking, numerical methods for the solution of dynamic optimization problems can be classified into two groups: direct and indirect methods. Indirect methods solve the optimization problem using the Pontryagin’s maximum principle. This method is based on the transformation of the original problem using the necessary optimality conditions of Pontryagin. This results in a two or multipoint boundary value problem, in the presence the constraints, to be solved for state and costate variables.
Direct approaches, such us complete parameterization (CP, [26]), multiple shooting (MS, [27]) or control vector parameterization (CVP, [28]) transform the original problem into a nonlinear programming problem (NLP) by means of the discretization and approximation, either of the control variables or both the control and state variables. The main differences among these methods are the number of decision variables, the presence or absence of parameterization constraints and the necessity of using an initial value problem solver. The use of CP or MS results in larger NLP problems with junction constraints to handle the system dynamics, which requires the use of specific optimization methods and may be computationally intensive for large scale dynamic systems. Therefore the CVP method is selected here.
Control vector parametrization
In the CVP approach the time horizon is divided into a number of ρ time intervals, with variable or fixed length. The control variables are then approximated within each interval by means of low order polynomials. With this parametrization the original dynamic optimization problem is transformed into a nonlinear programming problem with dynamic and algebraic constraints. The nonlinear programming problem obtained must be solved by employing a suitable NLP method and an initial value problem solver.
Nonlinear programming methods
Nonlinear programming methods are basically classified in local and global methods. Local methods are designed to converge to the closest solution to the initial guess. The iterates are computed by means of the cost function value or the gradient and/or the Hessian of the cost function with respect to the decision variables. In addition there are implementations that are able to automatically handle constraints. In the context of dynamic optimization, Sequential Quadratic Programming (SQP) methods, which guarantee feasibility of the solution at convergence, or feasible SQP methods, which guarantee feasibility throughout the optimization, have shown to be efficient in the solution of largescale constrained NLPs.
Unfortunately, NLPs arising from the application of direct approaches (such as CVP) are frequently multimodal. This is particularly the case in the presence of nonlinear or bilinear dynamic constraints (as in the problems considered in this work) or when variable length elements are used in the CVP approach. Therefore, local optimization techniques may converge to local optima. A multistart of local methods may help to analyze multimodality of the optimization problem. Interested readers will find illustrative examples in the Additional file 1.
Stochastic global optimization methods are robust alternatives for the search of the global solution. However these methods are well known to be particularly inefficient in the refining of the solutions. In this concern, BalsaCanto et al. [29], presented the combination of a stochastic global solver with a deterministic local method in a so called sequential hybrid approach for the solution of dynamic optimization problems. That work highlights the role of the tuning of the hybrid, i.e. the selection of the methods and the switching point, in the overall performance of the method. However several examples illustrated the benefits of the combination with respect to the individual solvers in the sense of robustness and efficiency.
Recent works [30, 31] show how the scatter search metaheuristic in combination with local methods outperforms previous sequential hybrid methods in the solution of complex optimization problems. The key property of this approach is that the switching points are automatically selected by the algorithm based on a balance between quality and diversity among the solutions generated in every iteration.
Based on the above, the following methods were selected attending to their previously reported performance:

SRES, the stochastic ranking evolution strategy method [32]. An evolutionary approach that is able to handle constraints by means of the stochastic ranking approach.

DE, differential evolution [33]. It is based on the generation of new solutions by adding the weighted difference between two population vectors to a third vector. The method was complemented with a death penalty function approach to handle constraints.

Sequential hybrids of SRES or DE with different implementations of successive quadratic programming methods.

eSS, scatter search [34] in combination with different implementations of successive quadratic programming methods.
Multiobjective optimization methods
In contrast to the singleobjective case, the aim in a multiobjective case is to find the optimal tradeoffs between conflicting objectives. The notion of a single optimal solution is replaced by the notion of a Pareto front, i.e. a set of optimal tradeoff solutions. All solutions in this set are optimal in the sense that it is not possible to improve one of the objectives without degrading one or more of the others. The weighted sum or the ε−constraint approach transforms the multiobjective problem into a set of singleobjectives problems whose solutions result in the Pareto front. However the use of the weighted sum presents a number of drawbacks such as the impossibility to recover nonconvex area in the Pareto front or the selection of weigths to obtain a uniform distribution of Pareto solutions. In view of this, the ε−constraint method [18] was selected here. The underlying idea is to solve a set of singleobjective problems where one objective (J _{ p }) is to be minimized while the others (J _{ i }) are incorporated as inequality constraints. Mathematically this reads as:
Subject to:
The Pareto front is obtained modifying the upperbounds of the constraints. Each of the singleobjective problems was solved using control vector parametrization combined with the eSS scatter search solver.
Results and discussion
This work considers the solution of two sets of examples. The first set consists of singleobjective problems taken from the literature. All these problems were initially solved with a multistart of local methods to analyze their possible multimodal nature. Since all problems appeared to be multimodal, a selection of global and hybrid methods were then used to find the global solution for each case study. The details are listed in Additional file 1, showing that scatter search (eSS) offered the best compromise between efficiency and robustness. A second set of problems considered a more general multiobjective formulation. The ε−constraint approach and scatter search (eSS) were subsequently used to solve this second set of problems, enabling us to obtain the optimal tradeoffs between different objectives. A detailed analysis of the resulting enzyme activation profiles revealed the “justintime” activation is still present in multicriteria optimality. The examples considered are summarized in Table 1. Results are presented and discussed below.
Threestep linear pathway with mass action kinetics (LPN3B)
This example was originally formulated by Bartl et al. [7] and solved by means of the maximum principle of Pontryagin. The pathway (Figure 1A) consists of three enzymatic reactions with mass action kinetics, each reaction catalyzed by a specific enzyme (e _{ i }). S _{1} corresponds to the substrate, S _{2} and S _{3} to the intermediate metabolites and S _{4} to the product. Metabolites and enzymes are expressed in concentration units and time in seconds.
In this work we considered the objective of the minimization of the time needed to reach a certain amount of product (Equation (10)), when the substrate remains unalterable (buffered substrate concentration) as represented in matrix N. This objective reflects a biological situation where there is no need to convert all the substrate into the product but a certain amount is essential.
The mathematical statement of the singleobjective problem is as follows.
Find e over t∈[ t _{ o },t _{ f }] to minimize:
Subject to the system dynamics:
Where N:
And:
With the following endpoint constraint:
and path constraint:
with: E _{ T }=1M, k _{ i }=1 s ^{−1} i=2,3,4, S _{1}(t _{0})=1M, S _{ i }(t _{0})=0 for i=2,3,4 and P(t _{ f })=0.9 M.
The above problem was solved with several methods, with the eSS scatter search solver as the best performer. Detailed results, including a discussion on multimodality, are given in Additional file 1: Tables S.6 and S.7.
Figure 1B presents the corresponding metabolite dynamics and optimal profiles of the enzymes for the singleobjective optimization problem presented in Equation (10). The solution was obtained under the assumption of no substrate consumption (for example, in a constant culture medium). The first phase of the process is devoted to produce S _{2} as fast as possible, the enzyme is fully activated and since substrate is not consumed a high amount of S _{2} is produced and accumulated. Similarly the second and third phases correspond to a full activation of enzymes e _{2} and e _{3} respectively. The process time is t _{ f }=4.2s. It should be remarked that the optimal solution for this case was achieved by eSS in less than 30 s of computation in a standard PC. Sequential hybrids of SRES and DE achieved similar values but requiring longer computation times.
The optimal solution of the above singleobjective problem (see Equation (10)) results in a significant accumulation of intermediates at the end of the process which may be harmful for the cell. Thus, we formulated a generalized problem so as to find the best compromise between the time required to achieve a certain amount of product and the intermediates accumulation. Mathematically the formulation of this multiobjective problem is:
Find e over t∈ [ t _{ o },t _{ f }] to minimize:
subject to the system dynamics presented in Equation (11), and with the constraints presented in Equations (13)–(14). The obtained Pareto front for the multiobjective problem (Equation (15)) is presented in Figure 1C. It can be observed that for process durations longer than 8 seconds the effect on the intermediate concentration is negligible, which means that no further improvements in the reduction of intermediates concentration can be achieved to satisfy the required amount of product.
To illustrate the differences in enzyme activation profiles for different tradeoffs, Figure 1D presents the optimal profiles for three solutions from the Pareto front (P _{1} and P _{3} as extreme points, and P _{2} as a balanced tradeoff between objectives). It can be noticed that the first enzyme is always fully activated, and its activation time is reduced when both objectives are considered (P _{2} and P _{3}). It can also be noted that an intermediate zone is generated where both enzymes are activated. In this region S _{2} and S _{3} are produced and consumed in order to avoid their accumulation in the pathway. Additionally it is observed that in the final stage the third enzyme is fully activated to obtain the desired amount of product (S _{4}).
Fourstep linear pathway with MichaelisMenten kinetics (LPDN4)
This example considers a four step linear pathway with MichaelisMenten kinetics and was originally considered by Oyarzún et al. [6]. In this case, in order to obtain more realistic results, enzymes are assumed not to become activated instantaneously, but to follow first order kinetics. The pathway (Figure 2A) consists in four enzymatic reactions catalyzed by a specific enzyme (e _{ i }) where S _{1} corresponds to the substrate, S _{2}−S _{3} to the intermediate metabolites and S _{4} to the product. The objective of the system is to minimize the enzyme consumption and the time needed to reach a given steady state.
This singleobjective problem reads mathematically as follows:
Find r over t∈[ t _{ o },t _{ f }] to minimize:
Subject to the system dynamics:
where the enzyme dynamics is considered to be linear with the expression rate (r) and λ=0.5.
with:
and:
with the following end point constraints, which describe the given steady state.
In addition the following constraints are imposed to limit the amount of enzymes and their rates:
With: E _{ T }=1 m M, k _{ c a t 1}=1 s ^{−1}, k _{ c a t 2}=2 s ^{−1}, k _{ c a t 3}=4 s ^{−1}, k _{ c a t 4}=3 s ^{−1}, K _{ M }=1 s ^{−1}, V=0.2 m M/s, S _{1}(t _{0})=5mM, S _{ i }(t _{0})=0 for i=2,3,4, e _{ SS }= [0.24 0.26 0.21 0.29] mM for i=1,2,3,4.
Note that in this case the expression rates are computed via dynamic optimization and the optimal enzyme activation profiles are obtained from the Equation (18).
The optimal value J ^{∗} = 6.154 was obtained using eSS. Note that this result is in good agreement with the one previously reported (J ^{∗} = 6.298) [35]. Figure 2B presents the optimal metabolite and enzyme dynamics for the singleobjective case (Equation (16)). As expected, the incorporation of the enzyme dynamics slows down the entire process. The optimal profiles for the expression rates follow a switching pattern that matches the pathway topology, leading to enzyme profiles that follow a sequential activation in agreement with previous results. The enzyme profiles show that for the synthesis of one enzyme the degradation of the previous enzyme is needed. The constraints imposed in the final amount of metabolites require a high accumulation of metabolites in the system and this could be lethal for the system [36].
It is interesting to note that in two previous works [6, 35] a combined objective function including the transition time and the enzyme cost was proposed, obtaining a single tradeoff solution. In other words, they obtained one of all the possible Pareto solutions. Here we go a step a further and formulate and solve the multiobjective optimization problem in order to obtain the full Pareto front. Mathematically the formulation of the multiobjective problem is:
Find e over t∈[ t _{ o },t _{ f }] to minimize:
Subject to Equations (17)–(18) and with the constraints presented in Equations (21)–(22).
Due to the final time constraints used to define the steady state the system keeps several enzymes activated during the process, which implies an extra effort for the system as discussed in literature [7]. In this case end point constraints (Equation (20)) are replaced by a constraint on the final amount of product:
The Pareto front for the multiobjective case (Equation (23)) is presented in Figure 2C. It can be noticed that for large times (around 15 seconds) the effect on the reduction of enzyme concentration is negligible, which means that the enzyme consumption can not be further reduced to achieve the desired amount of product. In Figure 2D optimal profiles for expression rates and enzyme dynamics are shown in order to check how enzyme activation is affected by the different tradeoff. P _{1} corresponds to the minimization of the time and P _{2} and P _{3} to different tradeoff between objectives.
We also notice that the optimal profiles for the expression rates follow a switching pattern that matches with the pathway topology, i.e. is reflected in a sequential activation of the enzymes. Before activating one enzyme, the degradation process of the previous one has started (i.e. the cell has only a certain amount of protein available). This situation is more relevant when the enzyme consumption is reduced. Note that when the process time increases the activation time of the enzymes is slightly reduced (this decrease is more relevant for short process times).
Glycolysis inspired network (GBD)
The original problem was considered by Bartl et al. [15]. Results achieved for that formulation are discussed in Additional file 1. Here we propose a more realistic formulation that incorporates the enzyme dynamics. The enzyme dynamics (Equation (27)) are considered to be linear with the expression rate (r _{ i }). The pathway (Figure 3A) consists of four enzymatic reactions with one branch where S _{1} corresponds to the substrate, S _{2}−S _{4} to the intermediate metabolites and S _{5} to the product. The hypothesis in this problem is that the pathway activation minimizes the time needed to transform the substrate (S _{1}) into a fixed amount of product (S _{5}). Note that the substrate S _{1} is not consumed during the process.
The activation profiles for the singleobjective case may be acquired by computing r over t∈ [ t _{ o },t _{ f }] to minimize the cost function.
Subject to the system dynamics:
Where:
With:
and the following end point constraints:
and the following path constraint:
With E _{ T }=1mM, P(t _{ f })=0.75mM, k _{ cati }=1 s ^{−1}, K _{ M }=1 s ^{−1}, S _{1}(t _{0})=1mM, S _{ i }(t _{0})=e _{ i }(t _{0})=0 for i=2,3,4,5, λ=0.5 s ^{−1}.
The optimal value (J ^{∗}=9.5) was obtained using scatter search (eSS) in less than 50 s of CPU time. Detailed optimization results are presented in Additional file 1 (Table S.13). The corresponding optimal enzyme activation profiles for the singleobjective case (Equation (25)) are shown in Figure 3B. Again the optimal profiles for the expression rate follow a switching pattern that matches with the pathway topology leading to enzyme profiles that follow a sequential activation profile.
The enzyme profiles show that for the synthesis of one enzyme the degradation of the previous is needed, meeting the constraint on the total amount of enzyme (Equation (30)). As in previous cases, there is a high accumulation of metabolites, which could be harmful for the cell. Note that the problem formulation could be modified so to predict the scenario with no accumulation of intermediates. In addition, these calculations were preformed for a fixed λ value, but of course the problem can be easily solved to consider other cases where genes might have different expression rates and proteins might have different degradation rates.
We now consider a multiobjective formulation, extending the objective function and keeping unchanged the rest of the problem. Mathematically, the multiobjective problem is formulated as follows:
Find r over t∈ [ t _{ o },t _{ f }] to minimize:
The Pareto front is presented in Figure 3C for the objectives in Equation (31). For long process times, no significant improvement in the enzyme cost is achieved. The optimal enzyme profiles corresponding to 3 points in the Pareto front are depicted in Figure 3D. P _{1} corresponds to the singleobjective case (time minimization), and P _{2} and P _{3} represent different tradeoffs between process time and enzyme cost minimization. The optimal profiles for the expression rate follow a switching pattern that matches with the pathway topology. In P _{1} we observed that for the synthesis of one enzyme the degradation of the previous is needed. In fact, from P _{2} and P _{3} it can be noted that such situation gains relevance as we move in the Pareto front.
One interesting and common situation in branched pathways is that the system could have two different outputs (p.e. produce several amino acids in sensible ratios for protein synthesis), which in practice means that the cell resources are distributed accordingly. The introduction of a second output in the singleoutput pathway was considered to study the behavior of the system if several products are generated. This new pathway is presented in Figure 4A.
The new multiobjective problem mathematically reads as: find activation profiles by computing r over t∈[ t _{ o },t _{ f }] to maximize:
Subject to the system dynamics presented in Equations (26)–(27), with:
And with the path constraint presented in (30). In this case the final time was considered constant (t _{ f }=15 s). The Pareto front for the problem formulated in Equation (32) is presented in Figure 4B, where point P _{1} and P _{3} stand for singleobjective situations were only one product is produced, and P _{2} is a balanced tradeoff. In Figure 4C the optimal profiles of the selected Pareto points (P _{1}, P _{2} and P _{3}) are depicted. In all cases we can see a sequential enzymatic activation. In P _{1} four enzymes are activated to convert the substrate (S _{1}) in the product S _{5}. In P _{3} only 3 enzymes are activated to obtain the product S _{6}. Regarding P _{2}, it can be noted that e _{1} and e _{2} are sequentially activated during the early process. Later on, in the intermediate stage, e _{4} and e _{5} are simultaneously activated. And finally, at the end of the process, only e _{4} is active.
Central metabolism of Saccharomyces cerevisiae during diauxic shift (SC)
The diauxic shift is characterized by decreased growth rate and by switching metabolism from glycolysis to aerobic utilization of ethanol under conditions of glucose depletion. This allows the maintenance of the cellular redox potential, NADH/NAD and ATP levels, enabling the cell to survive over longer periods without nutrients, in a so called quiescent state. The idea is to explain this particular behavior by computing the enzymatic activation profile that maximizes the survival time.
This case was considered by Klipp et al. [3, 4] using orthogonal collocation in combination with a genetic algorithm. The pathway (Figure 5A) consists in a simplified model of the central metabolism of Saccharomyces cerevisiae. The system complies six reactions: upper glycolysis, lower glycolysis, ethanol formation, ethanol consumption TCA cycle and respiratory chain. The model has arbitrary units.
Mathematically the singleobjective problem is formulated as follows:
Find e over t∈ [ t _{ o },t _{ f }] to maximize:
Subject to the system dynamics:
Where S and N correspond to:
And:
With a constraint in the total amount of enzyme:
And:
With: X _{1}(t _{0})=1, X _{2}(t _{0})=1, X _{3}(t _{0})=1, X _{4}(t _{0})=10, N A D H(t _{0})=0.7, N A D(t _{0})=0.3, A T P(t _{0})=0.8, A D P(t _{0})=0.2, NADH _{ c }=0.5, ATP _{ c }=0.7, k _{7}=3 and k _{8}=0.1.
It should be noted that constraints (38)–(39) indicate that the cells will survive provided the concentrations of NADH and ATP are above some given critical values (NADH _{ c } and ATP _{ c } respectively).
In order to compare our results with the ones reported in the work by Klipp et al. [4] and the available experimental data, two types of enzyme CVP approximations were considered: i) a stepwise approximation with 120 constant elements for each enzyme (i.e. 721 decision variables) and ii) a piecewise linear approximation with 8 variable length linear elements (i.e. 55 decision variables).
The optimal value (J* = 95.8) was again obtained using eSS scatter search. Note that the optimal solution is highly sensitive to the switching times (details are given in Additional file 1). The optimal profiles and the metabolite dynamics for the singleobjective problem (Equation (33)) are depicted in Figure 5B, showing a higher activity of the lower part of the glycolysis (e _{2}) and in the ethanol formation (e _{3}), as expected. During the last part of the process the activity of (e _{4}) is increased allowing the formation of NADH. Also there is an increase in the tricarboxylic acid cycle (e _{5}) and in the respiratory chain (e _{6}) which compensate the decline in supply of NADH and ATP.
In the above singleobjective formulation (Equation (33)), we can interpret that the system tries to survive as long as possible without any limitation on the investment of resources used. In the next case, we try to find the best compromise between the survival of the system and the resources used to do that. The aim in this case is to maximize the survival time with a low investment on the enzyme cost. We formulate the optimization problem as a multiobjective problem (with objectives time and the enzyme cost). Mathematically the multiobjective problem is formulated as follows:
Find e over t∈ [ t _{ o },t _{ f }] to maximize:
The set of solutions for the maximization problem (Equation (40)) is presented in Figure 5C. P _{1} corresponds to the singleobjective situation above, while P _{2}, P _{3} and O are different tradeoffs among the objectives. There, we can observe that the survival time of the system is not significantly reduced by a decrease on the enzyme consumption until the reduction on the enzyme consumption is over 60%. For low enzyme consumptions (under 60%), the survival time is clearly affected and diminishes rapidly.
The reduction on enzyme consumption is reflected as a reduction on the initial value of enzymes and their dynamics during the process. In Figure 5E optimal profiles for e _{5} obtained for four different points of the Pareto front are compared. Besides the reduction on the initial concentration of the enzyme, it can be observed that in P _{1} and O there is a sharp increase on the activity of the enzyme towards then end of the process. This dynamic is significantly reduced in P _{2}, and disappears in P _{3}, where the enzyme has a constant activation value. Complete optimal profiles are presented in Additional file 1.
We computed, for all the Pareto points, the distance between experimental data (taken from [4]) and model predictions, finding that point “O” in the Pareto corresponds to the best fit. Figure 5D presents the comparison between experimental data and model predictions for this optimal tradeoff solution. The numerical profiles are in good agreement with the general tendency of the experimental data. e _{1} and e _{3} slowly decrease the activity whereas e _{4} and e _{5} slowly increase their activity arriving to high activation values. e _{2} keeps an almost constant activation value along the time horizon.
Conclusions
This work proposes the use of multiobjective optimization, combined with an advanced dynamic optimization method, as a systematic way to predict or explain genetic metabolic regulation. The presented methodology can handle very general formulations, including several objectives, arbitrary network topologies and nonlinear kinetics, large numbers of control variables (enzymes or activation rates), and general path and point constraints on controls and states (metabolites concentrations).
The methodology is based on the transformation of the original multiobjective dynamic optimization problem into a set of nonlinear programming (NLP) problems by means of the control vector parametrization. The NLP problems can be solved using the combination of an initial value problem solver and an efficient global optimizer. The need of global optimization arises from the nonconvexity of these problems.
To illustrate this methodology, two sets of problems were considered. A first set of singleobjective examples were taken from the literature. The obtained results are comparable or better than those previously published and were obtained with a minimum computational cost, between a few seconds and a few minutes in a standard PC. Our results show clearly that a hybrid metaheuristic based on scatter search was the most robust and efficient solver when dealing with this class of problems.
A second set of multiobjective formulations were considered, for the first time, to achieve more biologically meaningful results. The Pareto front which represents the different tradeoffs between objectives was obtained by transforming the multicriteria problem into a set of singleobjective problems, each of them solved by means of the previously mentioned methodology. Interestingly, the optimal activation profiles computed exhibited the “justintime” behavior in most cases and for both single and multicriteria formulations. By definition, all the solutions in a Pareto front are equally optimal. However, if one introduces an additional requirement or constraint (e.g. fit to experimental data), it is then possible to select a single solution from the front which meets such additional criteria (as illustrated with the S. cerevisiae case study).
The methodology presented here opens up new possibilities, especially with respect to the handling of medium and large scale complex networks. It will also enable further research, including experimental validation, in order to asses biological optimality in terms of Paretooptimal tradeoffs.
References
 1.
Sutherland WJ: The best solution. Nature. 2005, 435: 56910.1038/435569a.
 2.
Orth JD, Thiele I, Palsson BO: What is flux balance analysis?. Nat Biotechnol. 2010, 28 (3): 245248. 10.1038/nbt.1614.
 3.
Klipp E, Holzhütter HG, Heinrich R: Reprogramming of the metabolic system by altering gene expression. Proceedings of the 9th Int. Meeting on BioThermoKinetics: Animating the Cellular Map. 2000, South Africa: Stellenbosch University Press, 6570.
 4.
Klipp E, Heinrich R, Holzhütter HG: Prediction of temporal gene expression. Metabolic optimization by redistribution of enzyme activities. Eur J Biochem. 2002, 269: 54065413. 10.1046/j.14321033.2002.03223.x.
 5.
Zaslaver A, Mayo AE, Rosenberg R, Bashkin P, Sberro H, Tsalyuk M, Surette M, Alon U: Justintime transcription program in metabolic pathways. Nat Genet. 2004, 36: 486491. 10.1038/ng1348.
 6.
Oyarzún DA, Ingalls BP, Middleton RH, Kalamatianos D: Sequential activation of metabolic pathways: a dynamic optimization approach. Bull Math Biol. 2009, 71: 18511872. 10.1007/s1153800994275.
 7.
Bartl M, Li P, Schuster S: Modelling the optimal timing in metabolic pathway activationuse of Pontryagin’s Maximum Principle and role of the golden section. Biosystems. 2010, 101: 6777. 10.1016/j.biosystems.2010.04.007.
 8.
Banga JR, BalsaCanto E: Parameter estimation and optimal experimental design. Essays Biochem. 2008, 45: 195210. 10.1042/BSE0450195.
 9.
Lebiedz D, Maurer H: External optimal control of selforganisation dynamics in a chemotaxis reaction diffusion system. IET Syst Biol. 2004, 2: 222229.
 10.
Salby O, Sager S, Shaik OS, Kummer U, Lebiedz D: Optimal control of selforganized dynamics in cellular signal transduction. Math Comput Model Dynamical Syst. 2007, 13: 487502. 10.1080/13873950701243969.
 11.
Vilas C, BalsaCanto E, Garcia MSG, Banga JR, Alonso AA: Dynamic optimization of distributed biological systems using robust and efficient numerical techniques. BMC Syst Biol. 2012, 6: 7910.1186/17520509679.
 12.
Ou J, Yamada T, Nagahis K, Hirasawa T, Furusawa C, Yomo T, Shimizu H: Dynamic change in promoter activation during lysine biosynthesis in Escherichia coli cells. Mol Biosyst. 2008, 4: 128134. 10.1039/b711035a.
 13.
Chechik G, Oh E, Rando O, Weissman J, Regev A, Koller D: Activity motifs reveal principles of timing in transcriptional control of the yeast metabolic network. Nat Biotechnol. 2008, 26: 12511259. 10.1038/nbt.1499.
 14.
Oyarzún DA: Optimal control of metabolic networks with saturable enzyme kinetics. IET Syst Biol. 2011, 5 (2): 110119. 10.1049/ietsyb.2010.0044.
 15.
Bartl M, Kötzing M, Kaleta C, Schuster S, Li P: Justintime activation of a glycolysis inspired metabolic network  solution with a dynamic optimization approach. Proceedings 55nd International Scientific Colloquium. 2010, Ilmenau (Germany): Universitätsbibliothek Ilmenau
 16.
Handl J, Kell DB, Knowles J: Multiobjective optimization in bioinformatics and computational biology. Computat Biol Bioinformatics, IEEE/ACM Trans. 2007, 4 (2): 279292.
 17.
Sendin JOH, Vera J, Torres NV, Banga JR: Model based optimization of biochemical systems using multiple objectives: a comparison of several solution strategies. Math Comput Model Dynamical Syst. 2006, 12 (5): 469487. 10.1080/13873950600723442.
 18.
Sendin JOH, Exler O, Banga JR: Multiobjective mixed integer strategy for the optimisation of biological networks. IET Syst Biol. 2010, 4 (3): 236248. 10.1049/ietsyb.2009.0045.
 19.
ElSamad H, Khammash M, Homescu C, Petzold L: Optimal performance of the heatshock gene regulatory network. Proceedings 16th IFAC World Congress, Volume 16. 2005, Prague, Czech Republic: Elsevier, 22062206.
 20.
Higuera C, Villaverde AF, Banga JR, Ross J, Moran F: Multicriteria optimization of regulation in metabolic networks. PLoS ONE. 2012, 7 (7): e4112210.1371/journal.pone.0041122.
 21.
Shoval O, Sheftel H, Shinar G, Hart Y, Ramote O, Mayo A, Dekel E, Kavanagh K, Alon U: Evolutionary tradeoffs, pareto optimality, and the geometry of phenotype space. Science. 2012, 336 (6085): 11571160. 10.1126/science.1217405.
 22.
Noor E, Milo R: Efficiency in evolutionary tradeoffs. Science. 2012, 336 (6085): 11141115. 10.1126/science.1223193.
 23.
Szekely P, Sheftel H, Mayo A, Alon U: Evolutionary tradeoffs between economy and effectiveness in biological homeostasis systems. PLOS Comput Biol. 2013, 9: e100316310.1371/journal.pcbi.1003163.
 24.
Wessely F, Bartl M, Guthke R, Li P, Schuster S, Kaleta C: Optimal regulatory strategies for metabolic pathways in Escherichia coli depending on protein costs. Mol Syst Biol. 2011, 7: 515528.
 25.
Banga JR, BalsaCanto E, Moles C, Alonso AA: Dynamic optimization of bioprocesses: efficient and robust numerical strategies. J Biotechnol. 2005, 117: 407419. 10.1016/j.jbiotec.2005.02.013.
 26.
Biegler LT, Cervantes AM, Wätcher A: Advances in simultaneous strategies for dynamic process optimization. Chem Eng Sci. 2002, 57 (4): 575593. 10.1016/S00092509(01)003761.
 27.
Bock HG, Plitt KJ: A multiple shooting algorithm for direct solution of optimal control problems. Proceedings 9th IFAC World Congress. 1984, New York: Pergamon Press, 242247.
 28.
Vassiliadis VS, Sargent RWH, Pantelides CC: Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints. Ind Eng Chem Res. 1994, 33 (9): 21112122. 10.1021/ie00033a014.
 29.
BalsaCanto E, Vassiliadis VS, Banga JR: Dynamic optimization of single and multistage systems using a hybrid stochasticdeterministic method. Ind Eng Chem Res. 2005, 44 (5): 15141523. 10.1021/ie0493659.
 30.
Egea JA, RodriguezFernandez M, Banga JR, Marti R: Scatter search for chemical and bioprocess optimization. J Glob Optimization. 2007, 37 (3): 481503. 10.1007/s1089800690753.
 31.
Egea JA, BalsaCanto E, Garcia MG, Banga JR: Dynamic optimization of nonlinear processes with an enhanced scatter search method. Ind Eng Chem Res. 2009, 48 (9): 43884401. 10.1021/ie801717t.
 32.
Runarsson TP, Yao X: Stochastic ranking for constrained evolutionary optimization. IEEE Trans Evol Comput. 2000, 564: 284294.
 33.
Storn R, Price K: Differential evolution  a simple and efficient heuristic for global optimization over continuous spaces. J Glob Optimization. 1997, 11: 341359. 10.1023/A:1008202821328.
 34.
Egea JA, Martí R, Banga JR: An evolutionary method for complexprocess optimization. Comput Oper Res. 2010, 37 (2): 315324. 10.1016/j.cor.2009.05.003.
 35.
Oyarzún DA: A controltheoretic approach to dynamic optimization of metabolic networks. PhD thesis. PhD thesis. National University of Ireland Maynooth, Ireland 2010
 36.
Cascante M, Lloréns M, MeléndezHevia E, Puigjaner J, Montero F, Martí E: The metabolic productivity of the cell factory. J Theor Biol. 1996, 153: 317325.
Acknowledgements
This research received financial support from the Spanish Ministerio de Economía y Competitividad (and the FEDER) through the project “MultiScales” (DPI201128112C0403), and from the CSIC intramural project “BioREDES” (PIE201170E018). Gundián M. de Hijas Liste acknowledges financial support from the MICINNFPI programme.
Author information
Affiliations
Corresponding authors
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
GMHL performed the computations. EBC implemented the dynamic optimization code. EBC and JRB contributed to the conception and design of the work. EBC, JRB and EK supervised the project and helped to draft the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
Further details on the use of global dynamic optimization to predict the activation in metabolic pathways.
Additional file 1: The Additional file 1 presents a more detailed description of the numerical approaches used in this work as well as a comparative study of the results achieved. (PDF 488 KB)
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
de HijasListe, G.M., Klipp, E., BalsaCanto, E. et al. Global dynamic optimization approach to predict activation in metabolic pathways. BMC Syst Biol 8, 1 (2014). https://doi.org/10.1186/1752050981
Received:
Accepted:
Published:
Keywords
 Dynamic optimization
 Global optimization
 Multiobjective optimization
 Pareto optimality
 Metabolic pathways
 Gene expression