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.

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:

$\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].

Param.	Value	Param.	Value
α 2	8	h 11	0.5
α 3	4.0556	h 22	1.4224
α 4	1.8397	h 23	0.6109
α 5	4.0556	h 44	0.5829
β 1	1	g 21	0.5
β 2	5.1179	g 32	0.4171
β 4	4.0556	g 42	2.8274
k	0.8	g 43	1.4646
		g 54	0.5

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 ).

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.

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.

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.

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.