Optimization strategies for metabolic networks

Abstract

Background

The increasing availability of models and data for metabolic networks poses new challenges in what concerns optimization for biological systems. Due to the high level of complexity and uncertainty associated to these networks the suggested models often lack detail and liability, required to determine the proper optimization strategies. A possible approach to overcome this limitation is the combination of both kinetic and stoichiometric models. In this paper three control optimization methods, with different levels of complexity and assuming various degrees of process information, are presented and their results compared using a prototype network.

Results

The results obtained show that Bi-Level optimization lead to a good approximation of the optimum attainable with the full information on the original network. Furthermore, using Pontryagin's Maximum Principle it is shown that the optimal control for the network in question, can only assume values on the extremes of the interval of its possible values.

Conclusions

It is shown that, for a class of networks in which the product that favors cell growth competes with the desired product yield, the optimal control that explores this trade-off assumes only extreme values. The proposed Bi-Level optimization led to a good approximation of the original network, allowing to overcome the limitation on the available information, often present in metabolic network models. Although the prototype network considered, it is stressed that the results obtained concern methods, and provide guidelines that are valid in a wider context.

Background

Current metabolic engineering processes allow to manipulate metabolic networks to improve the desired characteristics of biochemical systems [1]. These manipulations may lead to the maximization of the normal product yield or redirect the production to a flux that was residual or non-significant in the original network. The high level of uncertainty in metabolic network models knowledge makes it extremely difficult to determine what are the required manipulations needed to attain a given objective. Since an heuristic approach to such problems does not allow to explore the maximum potential of metabolic engineering, two approaches are usually considered when modeling metabolic networks. Kinetic models describe the complete dynamics of the network, and have proven useful to implement optimization and control over the network, such as in [2]. The creation of reliable kinetic models involves the estimation of parameters, the complexity of this task increasing with the size of the network considered.

The second approach models the networks on the basis of reaction stoichiometry. Although easier to obtain, these models lack the ability to directly predict the dynamics of the system.

Several techniques have been proposed to optimize and infer network characteristics from these models. In [3] a platform that combines many of these methods is presented. Flux Balance Analysis (FBA) allows the determination of the optimal flux distribution on a network described in terms of the stoichiometry of the reactions and yields reliable results in the study of metabolic systems [4–6]. A review of the method can be seen in [7].

When optimizing a metabolic network for a given objective two distinct problems must be addressed. The first is to find which branch or branches must be manipulated. The second is to determine what type of alterations must be done. Strategies such as OptKnock [8] and the work in [9] address the first problem. In this work a strategy for the second problem is described.

The simulation and engineering of metabolic network models typically involves complex optimization procedures. Geometric Programming (GP), one of the techniques used in this paper, is a powerful mathematical optimization tool that can be applied to problems where the objective and constraint functions have a special form [10]. GP is of particular interest because it can solve large scale problems with extreme efficiency and reliability [11]. Furthermore it has been shown that a problem formulated in S-Systems form can be solved with GP after a minimum adaptation [12].

A common optimization problem is the maximization of the final concentration of a metabolite whose formation competes with the natural objective of the cell (e.g. maximization of biomass). In this work, a prototype network with such behavior is taken as example and the corresponding optimization problem is solved with three alternative methods.

It is stressed that the emphasis of this work is on the methods and not the specific network considered. The key point of the paper consists in establishing properties of a number of optimization methods that may serve as guidelines when considering more complex networks.

This will be further explained in the next section.

Results and Discussion

An overall view of the problem considered and paper contributions is first presented. Details may then be seen in subsequent sections.

The problem to consider consists in finding a control function, defined over a finite interval of operation time, such that the final concentration of a desired product is maximum. This product is yielded by a metabolic network that, depending on the control function, either produces it or a product that favors cell growth. In order to settle ideas, assume that the control variable u is such that it is constrained to be in the interval [0, 1], with u = 0 corresponding to only production of cell growth product and no production of the desired product, and u = 1 corresponds to the inverse situation. Values of u in between 0 and 1 correspond to a mixed production in a way that depends on the network dynamics.

Since the optimization is with respect to a time function, this is an in finite dimensional problem. However we prove in this paper, using Pontryagin's Maximum Principle [13], that the optimal control only assumes values of 0 and 1. This is a priori assumed by other authors [14] and receives now a solid justification. It is a result valid for similar metabolic network problems that aim at optimizing a final yield (e.g. a concentration at the end of the optimization time interval, such as in [15]) and such that the control enters linearly in the network equations.

The significance of this result consists in the fact that, instead of searching the optimal control among piecewise continuous functions assuming values between 0 and 1, one only has to look functions assuming the extreme values of 0 and 1.

Furthermore, in the case study considered, it is shown that the optimum has only one switch between 0 and 1. Therefore the search for the optimum is reduced to find the switching instant, t _reg , that leads to the maximum final yield. Considering the structure of the metabolic network, this is intuitive: the optimum is achieved by first applying all cell resources to population growth and, after t _reg , to redirect them to desired production. If t _reg is too small, the desired production rate is higher during more time, but the cell population to which it applies is small. If t _reg is too big, there are many cells to produce, but they only act during a small time interval. Hence, there is an optimum value for t _reg .

As mentioned in the Background section, a major problem is the high level of uncertainty in the knowledge about metabolic network dynamics. In this respect we consider different optimization algorithms that assume various degrees of information about the system to be optimized.

The first is direct optimization. This assumes complete knowledge about the system and is included to establish a benchmark with which other methods may be compared.

The other two methods are variants of a bi-level algorithm designed in order to accommodate missing information on the network kinetics. Both cases differ from the type of inner-optimization: Geometric Programming in one case and Linear Programming in the other. Both methods lead to good approximations of the optimal control, with a slight advantage of the one relying on Geometric Programming.

Prototype network model

The optimization strategies were tested on a prototype network that is a modified version of a previously one suggested in [16]. The choice of this network was due to its widespread use as a test benchmark for several optimization algorithms. A graphical representation of the network is shown in Figure 1 associated with the following set of ordinary differential equations:

Figure 1

Prototype network. The circles correspond to metabolites and the arrows to fluxes with the reaction rates indicated.

\begin{array}{l}\frac{d{x}_1}{dt}=k-v_1\\ \frac{d{x}_2}{dt}=v_1-v_2(1-u)-v_{3}u\\ \frac{d{x}_3}{dt}=v_2(1-u)\\ \frac{d{x}_4}{dt}=v_{3}u-v_4\\ \frac{d{x}_5}{dt}=v_4\end{array}

(1)

Here the states x _i , i = 1,..., 5 are metabolite concentrations at the network nodes, v _i , i = 1,..., 4 are fluxes associated to the metabolic network branches and k is a constant parameter that represents the uptake of x₁. In the equations, u represents a control function that allows to redirect the flux between the branches x₂ → x₃ and x₂ → x₄. Assuming that x₃ represents a precursor of the cellular objective (such as growth) and x₅ the desired product, if u(t) is biased towards the branch of v₂ this yields the formation of x₃ but little or no production of x₅. If u(t) is biased towards the branch of v₃ the production of x₅ will be affected by the low concentration of x₃ (since there is a forward feedback). Thus, there is an optimal profile for u(t) to maximize the concentration of x₅ at the final time t _final .

In the framework of S-systems [16] the prototype network is described by:

\begin{array}{l}\frac{d{x}_1}{dt}=k-{\beta }_1{x}_1^{h_{11}}\\ \frac{d{x}_2}{dt}={\alpha }_2{x}_1^{g_{21}}-{\beta }_2{x}_3^{h_{23}}{x}_2^{h_{22}}\\ \frac{d{x}_3}{dt}={\alpha }_3{x}_2^{g_{32}}\left(1-u\right)\\ \frac{d{x}_4}{dt}={\alpha }_4{x}_3^{g_{43}}{x}_2^{g_{42}}u-{\beta }_4{x}_4^{h_{44}}\\ \frac{d{x}_5}{dt}={\alpha }_5{x}_4^{g_{54}}\end{array}

(2)

where β _i are the rate constants, g _ij and h _ij are the kinetic orders. Table 1 shows the list of parameters. All the simulations using the prototype network assume x(0) = [0.8 0 1 0 0].

Table 1 Parameters used in the prototype network.

Direct optimization

Direct optimization uses model (2) with the set of parameters from Table 1.

On a first approach, all possible integer values of t _reg in the interval t _reg = [1, 30] were used to compute the final product concentration x₅(t _final ) where t _final = 30s. Figure 2 plots the resulting function J(t _reg ) = x₅(t _final ).

Figure 2

Results of the simulation using Direct optimization. The final product concentration is shown as a function of T _reg . For the value of T _reg corresponding to the dotted line there is a maximum yield.

It is clear from Figure 2 that there is an optimal value for the time of regulation that maximizes the yield of x₅. For the network considered, the optimal time of regulation is t _reg = 9s. If u(t) switches from 0 to 1 before t _reg the formed biomass will not be enough to maximize x₅(t _final ). On the other hand, if u(t) switches from 0 to 1 after t _reg , there will be more biomass but there time will not be enough time to produce the maximum possible amount of x₅.

To illustrate better the behavior of the prototype network, simulations were made for t _reg = 4, t _reg = 9 and t _reg = 14. The obtained optimal t _reg = 9 is compared in Figure 3 with lower and upper values in order to show the different time evolution of the metabolites.

Figure 3

Comparison of three u ( t ) profiles. Three time profiles for the control function u(t) (above) and the corresponding product yield (below). The solid line is the optimal T _reg obtained by Direct Optimization.

Bi-Level Optimization

The Bi-Level optimization was used to test all the possible values of t _reg . Figure 4 plots the normalized curves for J(t _reg ) = x₅(t _final ) for the two optimizations, inner-optimization using Geometric Programming (GP) and inner-optimization using Linear Programming (LP). By comparing Figure 4 with Figure 2 it can be seen that the profiles remain similar. The final product yield, x₅(t _final ), increases with t _reg until the optimal value is reached, then it starts decreasing.

Figure 4

Result of the optimization using the Inner Optimization with Geometric Programming (left) and Linear Programming (right). The profiles of the production of x₅ remain similar to the simulation using Direct optimization.

The optimal time of regulation obtained with both GP and LP on the inner optimization was t _reg = 9. As shown mathematically in the methods section, the optimal control function is either 0 or 1, provided that the dynamics depends linearly on the control and the cost to optimize has only a final term.

In this case the dependency of the Hamiltonian function on u is linear (as given by (8) below). For the prototype considered, \varphi \left(\lambda \left(t\right),x\left(t\right)\right)=\left({\lambda }_4{\alpha }_4{x}_3^{g_{43}}{x}_2^{g_{42}}-{\lambda }_3{\alpha }_3{x}_2^{g_{32}}\right). Figure 5 shows a plot of ϕ(λ(t), x(t)) obtained with a near-optimal control function u(t). As expected, ϕ(λ(t), x(t)) is negative for values smaller than t _reg , leading to an optimal control u(t) = 0 and becomes positive for values larger than t _reg , leading to u(t) = 1. Thus, the optimal control is obtained on the extremes of the allowed interval and furthermore, one single switch (from 0 to 1) is enough to achieve the optimal control. It should be remarked that, since ϕ(λ(t), x(t)) is close to zero around t = t _reg , in practice, when using a numeric method there can be some jittering in the transition of the manipulated variable.

Figure 5

Plot of ϕ ( λ ( t ), x ( t )) (9). This function changes sign at the optimal instant of control switching T _reg .

Conclusions

For a class of networks in which the yield of the product that favors cell population growth (the "natural" product) competes with the desired product yield, with the manipulated variable affecting linearly the fluxes, it has been shown that the optimal control assumes only extreme values. While the implementation of this optimal control poses no challenge on in silico metabolic networks, on real metabolic networks complex bioengineering skills are required. Gene knockout manipulations do not adequate to this kind of control problem due to the long time scale associated with these techniques. The manipulation of specific enzyme levels, controlled by modulating the expression of the corresponding genes using promoter systems and inducers, is a possible solution to this kind of control problem [14].

The use of a bi-level optimization strategy, that maximizes the natural product in the inner level by manipulating the fluxes, leads to a good approximation to the optimal solution, with the advantage of not requiring the full knowledge of the network model. Real networks are extremely complex and exhibit relations between metabolites that are not always expected or fully understood. This gives emphasis to the need of good in silico models and also to the determination of the exact branches to be modified when optimizing a network. Although the example network used is very simple, it has proved to be useful to test the optimization strategies but a more complex network should be used to confirm that the strategy can be scaled to a larger network.

Methods

The optimization problem

The optimization problem consists in relation to a nonlinear state model of a metabolic network like (2), in selecting u(t) for t ∈ [0, t _final ] such that:

J\left(u\right)=x_{\text{5}}\left(t_{final}\right)

(3)

is maximized under the constraint that u(t) ∈ [0, 1], ∇t ≥ 0.

The solution of the optimization problem is obtained using different approaches. Before accomplishing this task, Pontryagin's Maximum Principle is invoked to establish a particular form of the optimal control function for the class of problems at hand.

Optimization

The control function is now optimized in order to obtain a maximum yield of biomass at the end of the run-time (t _final ). Three different methods, assuming various levels of information about the network, are considered in order to attain this goal.

The first method, direct optimization, is used as a benchmark to compare the results of the other methods. The last two methods rely on a Bi-level optimization and illustrate a possible solution to the optimization problem when the information about the network is incomplete.

Direct optimization

The first method, Direct Optimization, is used mainly as a benchmark, to compare the results of the following methods. Since it is assumed that all the information about the network kinetics is known, the system of differential equations, described in (2) is used. Given a function that receives t _reg as input and outputs the final yield of x₅, this optimization tests all the possible values of t_reg and returns the function J(t _reg ) = x₅(t _final ). The value of t _reg that results on a maximum product yield is then determined by solving a simple optimization problem.

The optimization was tested with two MATLAB functions: fmincon, from the standard optimization toolbox, that finds the minimum of a constrained nonlinear multi variable function, and simannealingSB from Systems Biology Toolbox [17] that performs simulated annealing optimization.

Bi-Level Optimization algorithm structure

The Bi-Level optimization algorithm was structured so as to accommodate missing information on the network kinetics. The boxed metabolites and fluxes from Figure 1 are a part of the network that might not be fully described in terms of kinetics. In this approach the missing kinetic information is replaced by stoichiometric data and flux balance analysis is used to obtain the proper flux distribution. Then, an inner optimization determines the fluxes during the batch time. The first step of the inner optimization process is to define the initial conditions of the input x₁ and outputs x₃, x₅. A valid distribution for the fluxes v₁, v₂, v₃ and v₄ is then obtained.

After obtaining the flux distribution, new values for the input/outputs can be calculated by integrating their expressions in the considered time interval. During this time interval the function u(t) and the values of v₁, v₂, v₃ and v₄ are kept constant. This process is repeated along a time grid from t = 0 to t = t _final . The time interval for the integration was defined to be 1 second. The inner optimization process allows us to obtain the product yield, x₅(t _final ), given a certain u(t), taking into account a valid approximation of the network dynamics over the simulation time. The detailed fluxogram of the inner-optimization is shown in Figure 6.

Figure 6

Inner-Optimization algorithm. Block diagram of the Inner-Optimization algorithm.

The bi-level optimization algorithm can be represented schematically as in Figure 7.

Figure 7

Bi-Level optimization formulation. Structure of the Bi-Level optimization.

Inner-optimization using Geometric Programming

On the first implementation of the Bi-Level optimization algorithm the dynamics of the boxed metabolites from Figure 1 are used but, following the algorithm structure, steady-state is assumed. Thus, {\dot{x}}_2 and {\dot{x}}_4 from (2) become:

\begin{array}{l}\frac{d{x}_2}{dt}={\alpha }_2{x}_1^{g_{21}}-{\beta }_2{x}_3^{h_{23}}{x}_2^{h_{22}}=0\\ \frac{d{x}_4}{dt}={\alpha }_4{x}_3^{g_{43}}{x}_2^{g_{42}}u-{\beta }_4{x}_4^{h_{44}}=0\end{array}

(4)

In this algorithm implementation, the inner optimization problem determines the profile of the metabolites, instead of fluxes, due to the nature of the equations. The metabolite concentrations are calculated at the beginning of each time interval, solving a Geometric Programming problem, and used with (2) to integrate the values of x₁, x₃ and x₅ during that interval.

Inner-optimization using Linear Programming

On the second implementation it is assumed that only stoichiometric information is available for the reactions inside the box of Figure 1. Assuming steady state, the equations of {\dot{x}}_2 and {\dot{x}}_4 become:

\begin{array}{l}\frac{d{x}_2}{dt}=v_1-v_2\times \left(1-u\right)-v_3\times u=0\\ \frac{d{x}_4}{dt}=v_3\times u-v_4=0\end{array}

(5)

Figure 1 shows a regulation from x₃ (Biomass) to flux v₃. Since stoichiometric models do not account for feedbacks, the effect of x₃ can not be integrated directly in the equations. Assuming that the forward feedback leads to an over expression of flux v₃, then a valid solution is to model the forward feedback as a variation of the constraints applied to flux v₃. Setting flux v₂ (precursor of Biomass formation) as the objective function, the FBA problem is solved with the previous equations to obtain a valid and unique flux distribution at each time step. In the context of the inner-optimization, these fluxes are then used to calculate the values of the input/outputs.

Pontryagin's Maximum Principle

A general tool to solve dynamic optimization problems such as the one considered here is Pontryagin's Maximum Principle PMP [13].

Let x be the state of a dynamical system with control inputs u such that:

\dot{x}=F\left(x,u\right),x\left(0\right)=x_0,u\left(t\right)\in U,t\in \left[0,T\right]

(6)

where F : ℜⁿ× ℜ → ℜⁿ, U is the set of valid control inputs and T is the final time, assumed here to be constant.

The control function u must be chosen in order to maximize the functional J, defined by:

J(u)=\psi (x_i(T))+{\displaystyle \underset{0}{\overset{T}{\int }}L}(x(t),u(t))dt

(7)

Where ψ is the cost associated with the terminal condition of the system and L the Lagrangian.

According to PMP, a necessary condition for the optimal control is that, along the optimal solution for the state x, co-state λ and control u the Hamiltonian H is maximum with respect to u[13].

Comparing the cost (3) with the generalized case (7) and taking into consideration that, in the case at hand, given by (1), the dynamics vector field depends linearly on the control, it follows that

H(\lambda ,x,u)=\lambda \varphi (\lambda ,x)u

(8)

where ϕ(λ, x) is a function that does not explicitly depend on u. Since, according to (8), the Hamiltonian is linear in u, its maximum is obtained at the boundary of the admissible control set U.

Hence, this shows that, for the metabolic network (1), the control that optimizes (3) only assumes the values u = 0 or u = 1.

In the case at hand, we are interested in maximizing the final value of the state x₅. Since the Lagrangian (L) is zero, (7) becomes J(u) = ψ(x _i (T)). Thus, the functional J to be maximized is:

\psi (x(T))=u_5(T_{final})

(9)

as shown before in (3).

Taking into account that, L = 0 the adjoint equations are reduced to

\dot{\lambda }=-f_x^T\lambda

(10)

The network is described by the system of ordinary differential equations in (2), if we consider the state model in the form of f(x, u), where u is the control function, calculating f _x (x, u) is straightforward.

Thus

\begin{array}{ccc}{\dot{\lambda }}_1& =& {\beta }_1{h}_{11}{x}_1^{h_{11}-1}{\lambda }_1-{\alpha }_2{g}_{21}{x}_1^{g_{21}-1}{\lambda }_2\\ {\dot{\lambda }}_2& =& \begin{array}{l}{\beta }_2{x}_3{h}^{23}{h}_{22}{x}_2^{h_{22}-1}{\lambda }_2-({\alpha }_3{g}_{32}{x}_2^{g_{32}-1})(1-u){\lambda }_3\\ -{\alpha }_4{x}_3^{g43}u{g}_{42}{x}_2^{g_{42}}{}^{-1}{\lambda }_4\end{array}\\ {\dot{\lambda }}_3& =& {\beta }_2{x}_2^{h_{22}}{h}_{23}{x}_3^{h_{23}-1}{\lambda }_2-{\alpha }_{4}u{x}_2^{g}42g_{43}{x}_3^{g_{43}-1}{\lambda }_4\\ {\dot{\lambda }}_4& =& {\beta }_4{h}_{44}{x}_4^{h_{44}-1}{\lambda }_4-{\alpha }_5{g}_{54}{x}_4^{g_{54}-1}{\lambda }_5\\ {\dot{\lambda }}_5& =& 0\end{array}

(11)

The terminal conditions for the co-states λ are

{\lambda }_n(T)=\frac{\partial \psi }{\partial x}{|}_{x=x(T)}=([00001])

(12)

Since L = 0 the Hamiltonian is given by λ^Tf(x).

Substituting in the expression and after some manipulation, becomes:

\begin{array}{l}H(\lambda (t),x(t),u(t),t)=\\ {\lambda }_1\left(k-{\beta }_1-x_1^{h_{11}}\right)+{\lambda }_2\left({\alpha }_2{x}_1^{g_{21}}-{\beta }_2{x}_3^{h_{23}}{x}_2^{h_{22}}\right)+{\lambda }_5({\alpha }_5{x}_4^{g_{54}})+\\ {\lambda }_3{\alpha }_3{x}_2^{g_{32}}-{\alpha }_4{\beta }_4{x}_4^{h_{44}}+({\lambda }_4{\alpha }_4{x}_2^{g_{43}}{x}_2^{g_{42}}-{\lambda }_3{\alpha }_3{x}_2^{g_{32}})u\end{array}

(13)

that depends linearly on the control function u, as expected.

The derivative of the Hamiltonian in order to the control function is:

\begin{array}{ccc}{H}_u& =& {\lambda }^T{f}_u\\ =& -{\lambda }_3{\alpha }_3{x}_2^{g_{32}}+{\lambda }_4{\alpha }_4{x}_3^{g_{43}}{x}_2^{g_{42}}\end{array}

(14)

This expression is the one that determines the number of switches between 0 and 1 of the control variable.