Prediction of activation of metabolic pathways via dynamic optimization

21st European Symposium on Computer Aided Process Engineering – ESCAPE 21 E.N. Pistikopoulos, M.C. Georgiadis and A. Kokossis (Editors) © 2011 Elsevier B.V. All rights reserved.

Prediction of activation of metabolic pathways via dynamic optimization Gundian M. De Hijas-Listea, Eva Balsa-Cantoa, Julio R. Bangaa. a

(Bio)Process Engineering Group, IIM-CSIC, Eduardo Cabello 6, Vigo 36208, Spain.

Abstract Genetic regulation of metabolic networks may be driven by optimality principles. Dynamic optimization can be used to either predict (a priori) or explain (a posteriori) the activation profile of genes and/or enzymes in metabolic pathways. In the proposed approach, the use of a suitable dynamic optimization direct method is combined with adequate global optimization solvers, implemented in the DOTcvpSB toolbox. The main advantage of this approach is that arbitrarily complex pathways can be considered, including branched ones, while avoiding convergence to local solutions. In order to illustrate the capabilities of the proposed methodology, several examples are presented and solved. Keywords: Metabolic pathways, dynamic optimization, global optimization.

1. Introduction The optimal control of metabolic systems is a key question in bioprocess systems engineering. 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. Klipp et al. [1] considered the case of linear pathways following mass action kinetics to illustrate the enzymatic switching profile that would lead to minimum transition time. Similarly Zaslaver et al. [2] theoretically confirmed the experimental observation of what they called “just-in-time” activation pattern in enzyme expression for the case of an unbranched pathway responsible of amino acid synthesis. More recently, optimal control theory has been used to deal with these problems. In particular, Oyarzún et al. [3] demonstrated by means of the maximum principle of Pontryagin the optimality of the bang-bang switching enzymatic activation profiles for the case of instantaneous enzyme activation. Bart et al. [4] generalized the results for the case of linear dynamics for the enzyme activation. However both works state the difficulty (or impossibility) of analytically solving other more realistic cases, such as those considering non-linear dynamics for the enzyme expression, or other arbitrarily complex networks, like e.g. branched pathways. For those cases, one must rely on advanced numerical techniques, and this is the objective of this contribution. This work formulates the problem as a general nonlinear dynamic optimization problem. The numerical solution is then approached through the combination of the control vector parameterization scheme with adequate non-linear optimization techniques, as provided in the recently developed toolbox DOTcvpSB [5]. It should be noted that the non-linear character of the models at hand, together with the usual presence of path and point constraints related to the metabolic and the enzymes concentrations or fluxes, usually induce multimodality, thus requiring the use of global optimization methods, as illustrated below.

Prediction of activation of metabolic pathways via dynamic optimization

1387

2. Problem formulation The objective of the general dynamic optimization problem is to compute the time varying control profiles and the final time (u(t), tf) so as to minimize (or maximize) a given cost function (J) subject to the system dynamics and possibly algebraic constraints. Mathematically this reads: tf

~ ~ min J [ x , u] min T [ x{tf }]  ³ I [ x{t}, u{t}, t ]dt

u ( t ),t f

u ( t ),tf

(1)

t0

Subject to:

dx ~ \ [ x{t}, u{t}, t ] , xt0 x 0 dt

>

@

h xt , ut

0,

>

@ d ut d uU

g xt , ut d 0

(2) (3), (4)

(5), (6) x d xt d x , u where x is the vector of state variables, u the vector of control variables, eqns. (2) represent the system dynamics, eqns. (3) and (4) regard the path or point equality and inequality constraints and eqns. (5) and (6) correspond to the upper and lower bounds for the state and control variables, respectively. L

U

L

3. Numerical methods. The methods employed for solving general dynamic optimization problems may be classified into two main groups: 1) Indirect approaches, which are based in the maximum principle of Pontryagin; 2) Direct approaches, which transform the original problem into a non-linear programming (NLP) by means of some type of discretization and interpolation scheme for the control variables or the controls and the states. 3.1. Control Vector Parametrization In the CVP approach, the time horizon is divided into a number of U time intervals. The control variables are then approximated within each interval by means of usually low order polynomials, so that the original (infinite dimensional) dynamic optimization problem is transformed into a NLP problem with dynamic and algebraic constraints. The non-linear programming problem (NLP) obtained must be then solved by employing a suitable solver. 3.2. Need of Global Optimization The non-linear character of the models at hand, together with the presence of path and point constraints related to the metabolic and the enzymes concentrations or fluxes usually induce multimodality. To illustrate this point, we considered the computation of the optimal enzyme activation profile for a given metabolic pathway so as to minimize a combined cost of the time taken to reach a given steady state, and a measure of the enzyme usage (Example 1 below) using a multi-start of a local method. From 500 starts, a local SQP method was able to converge to feasible solutions in only 48 runs (Fig. 1). Figure 1: Histogram of solutions obtained Most solutions are completely different from . with multi-start SQP method (500 runs)

G. M. De Hijas-Liste et al

1388

each other, a clear consequence of the multimodal nature of the problem.Therefore, robust and efficient global optimization (GO) methods are required in order to find the best solution. There are a large number of possibilities, from deterministic global methods to stochastic and hybrid global-local methods. Here we make use of certain stochastic or hybrid global methods which have proven to be successful in the solution of challenging dynamic optimization methods [6,7]. 3.3. DOTcvpSB This work makes use of DOTcvpSB [5] a MATLAB® toolbox devoted to solve dynamic optimization problems in systems biology (http://www.iim.csic.es/~dotcvpsb/). It combines in an efficient way the control vector parameterization approach with several state of the art local and global optimization techniques. In particular, popular global optimization solvers such as Differential Evolution (DE, [8]) or hybrid metaheuristics such as an ant colony optimization based method (ACOmi, [9]) are available.

4. Case studies. 4.1. Example 1: Metabolic pathway with Michaelis-Menten kinetics. In this first example, a simple linear pathway with Michaelis-Menten kinetics is considered. The objective is to minimize a combined cost of the time taken to reach a given steady state, and a measure of the enzyme usage. Such given steady state is characterized pre-defined values for end point constraints and total enzyme abundance [3]. The mathematical statement is as follows: E3 E2 E1 E0 S3 S0 S1 S2 Find E(t) over t[t0,tf] to minimize:

J

tf

3

0

i 0

³ (1  ¦ Ei (t )) ˜ dt

(7)

Subject to the system dynamics:

dSi dt

Xi 1 ( Si 1 (t ), Ei 1 (t ))  Xi ( Si (t ), Ei (t ))

Xi ( Si (t ), Ei (t ))

kcat i ˜ Si (t ) K mi  Si (t )

Ei (t )

With the following end point constraints: Xi (t ) V for t t t f , i 0, 1, ..., 3

i 1, ..., 3

(8) (9)

(10)

And the following path constraint: 3

¦ E (t ) d E i

i 0

T

(t )

(11)

With: kcat0=1s-1, kcat1=2s-1, kcat2=4s-1, kcat3=3s-1, Kmi=1mM, S0=5mM, Si=0 for i=1,…3, ET=1mM, V=0.2 The optimal solution corresponds to a bang-bang profile (Jopt=4.757) for each of the enzymes following a temporal sequence that matches with the pathway topology (see Figure 2). It should be noted that the best solution found had not been achieved with the

Prediction of activation of metabolic pathways via dynamic optimization

1389

multistart in Figure 1. Note that at final time certain enzyme levels are required so as to reach the steady state given by the flux endpoint constraints. Due to the minimization of the enzyme usage together with the path constraint on the total amount of enzyme, the initial substrate is not fully converted into final product. Figure 2: Dynamic profiles and optimal control values obtained with DOTcvpSB (ACOmi solver)

4.2. Example 2: Metabolic pathway with enzyme dynamics In nature, the synthesis of enzymes is not instantaneous as assumed in example 1. In order to obtain more realistic results, the model can be extended by considering the dynamics of the enzyme synthesis as a linear expression. This example modifies previous one by adding the following expressions for the enzyme dynamics:

dEi dt

ri (t )  O ˜ Ei (t ) with O

0.5 s 1

(12)

In this example the aim is once again to minimize a combined cost of the time taken to reach a given steady state, and a measure of the enzyme usage. In order to compare results with the previous example, flux end point constraints are replaced by concentration end point constraints corresponding to the steady state found in example 1. In addition the following constraints are imposed so as to limit the amount of enzymes ant their rates throughout the process.

0 d Ei (t ) d ET , 0 d ri (t ) d 1mM ˜ s 1

(13)

Figure 3: Dynamic profiles for the metabolite concentrations and for the enzymes obtained with DOTcvpSB (DE solver).

1390

G. M. De Hijas-Liste et al

As expected, the incorporation of the enzyme dynamics slows down the entire process. The optimal profiles (Jopt=6.164 ) for the expression rates follow a switching pattern that matches the pathway topology, leading to enzyme profiles that follow a sequential activation, similarly to the previous case (another example of just-in-time scheduling), in agreement with previous results from other authors.

5. Conclusions. Here we have shown how nonlinear dynamic optimization can be used as an efficient way to predict/explain regulation patterns in metabolic networks. We have also illustrated the need of global optimization methods in order to properly solve these problems. Computations for several case studies of increasing complexity were successfully carried out with the DOTcvpSB toolbox, recently developed by our group.

Acknowledgments The authors acknowledge financial support from Spanish MICINN project MultiSysBio (DPI2008-06880-C03-02) and CSIC project PIE 200730I002. Gundian M. De HijasListe acknowledges financial support from the MICINN FPI programme.

References [1] E. Klipp, R. Heinch,H.G. Holzhütter, 2002, Prediction of temporal gene expression. Metabolic opimization by re-distribution of enzyme activities, Eur. J.Biochem. 269:5406 [2] A. Zaslaver, A. Mayo, R. Rosenberg, P. Bashkin, H. Sberro, M. Tsalyuk, M. Surette, U. Alon, 2004, Just-in-time transcription program in metabolic pathways, Nat. Genet. 36:486. [3] D. A. Oyarzún,B.P. Ingalls, R.H. Middleton, D. Kalamatianos, 2009, Sequential activation of metabolic pathways: a dynamic optimization approach, Bulletin of Mathematical Biology 71:1851 [4] M. Bartl, P. Li, S. Schuster, 2010, Modelling the optimal timing in metabolic pathway activation-use of Pontryagin's Maximum Principle and role of the Golden section, Biosystems,101(1):67-77 [5] T. Hirmajer, E. Balsa-Canto, J. R., Banga, 2009, DOTcvpSB, a software toolbox for dynamic optimization on systems biology, BMC Bioinformatics. 10:199 [6] J.R. Banga and W.D. Seider. Global optimization of chemical processes using stochastic algorithms. In State of the Art in Global Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996; pp 563-583 [7] E. Balsa-Canto, V. S. Vassiliadis, J. R. Banga (2005) Dynamic Optimization of Single and Multi-Stage Systems Using a Hybrid Stochastic-Deterministic Method. Ind. Eng. Chem. Res., 44(5): 1514–1523. [8] R. Storn and K. Price (1997) Differential Evolution- a Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. J. Glob. Opt. 11:341-359 [9] M. Schlueter and J.R. Banga (2007) Ant Colony Optimization for Mixed Integer Nonlinear Programming with Real-World applications. 7th International Conference on Optimization: Techniques and applications ICOTA7, Kobe (Japan), 12-15 Dec. 2007.