An Improved Geometric Programming Approach for Optimization of Biochemical Systems

This paper proposes an improved geometric programming approach to address the optimization of biochemical systems. In the proposed method we take advantage of a special and interesting class of nonlinear kinetic models known as generalized mass action (GMA) models. In most situations optimization problems with GMAmodels are nonconvex and difficult problems to solve for global optimality. To deal with this difficulty, in this work, some transformation strategy is first used to convert the optimization problem with GMA models into an equivalent problem. Then a convexification technique is applied to transform this resulting optimization problem into a series of standard geometric programming problems that can be solved to reach a global solution. Two case studies are presented to demonstrate the advantages of the proposed method in terms of computational efficiency.


Introduction
Mathematical optimization of biochemical systems is a key step towards the establishment of rational strategies for yield improvement.In the last decades, much research has been directed toward the development of model-based optimization strategies for biochemical systems .Before optimizing a biochemical system, a key step in these approaches is the selection of the appropriate mathematical model among the different representations available.One type of well-developed representation for biochemical systems is the generalized mass action (GMA) model, which is based on a modeling framework called biochemical systems theory (BST) [22][23][24][25][26][27][28].The hallmark of this representation is that each process is represented separately as power-law formalism so that there are as many contributions as actual fluxes in the real biochemical system [28].The main advantage of GMA models is that they can capture the nonlinear characteristics of real biochemical system while showing some good properties required in the mathematical optimization.In most cases optimization problems of biochemical systems with GMA models belong to a truly nonconvex class of problems, that is, an intrinsically intractable NP-hard problem.Thus, these problems are difficult to solve for global optimality.In recent years, some optimization strategies have been proposed to deal with such difficult problems [6, 8-11, 14, 15].One of these approaches to the optimization of biological systems is the geometric programming (GP) method [6,14].These approaches are very interesting because they take advantage of the special structure features of GMA models; that is, each power-law term of GMA equations is indeed a monomial function required for GP.Several GP techniques including controlled error [6], penalty treatment [6], and successive convex approximation [14] methods have been presented to adapt GP solvers for the treatment of GMA systems.When both controlled error and penalty treatment approaches were used to optimize a biochemical system, a possible outcome is that they cannot find the global optimal solution of optimization problem [14].To deal with this difficulty, Xu presented a successive convex approximation method to solve the optimization problems of biochemical systems under steady-state conditions [14].The main idea of this approach is that the nonconvex optimization task of GMA models is first transformed as a standard GP problem through simple transformation and condense techniques.Then the resulting optimization problem can be solved very efficiently by a series of GPs.In the practical implementation of this approach, however, more iterations are possibly 2 Journal of Applied Mathematics required to find the global optimum of an optimization problem with GMA models (see Section 4).To deal with this issue and enhance the computational efficiency of the successive convex approximation method, it is necessary to make an improvement in its scheme.For this purpose, in the present study, we propose an improved GP method to solve steady-state optimization problems of biochemical systems with GMA models by adopting simple transformation and convexification strategies.Besides, two case studies are conducted to demonstrate the advantages of the proposed method in computational efficiency.Case studies show that the proposed algorithm significantly decreases the iterations to reach a global solution compared with the method in [14].
The rest of the paper is organized as follows.Section 2 describes the steady-state optimization problem of biochemical systems with GMA models.Section 3 presents the global optimization approach for solving steady-state optimization problem of biochemical systems.Section 4 gives two case studies drawn from the literature to demonstrate the advantages of the proposed method in computational efficiency.Finally, brief conclusions are presented in Section 5.

Optimization Problem Statement
Consider the following steady-state optimization problem of biochemical systems: subject to satisfying where  = ( 1 ,  2 , . . .,  + )  ∈  + ; the objective function () is usually a flux or a particular metabolite concentration; the parameters  +  > 0 and  −  > 0 are the stoichiometric coefficients of the metabolite   in the processes  +  and  −  , respectively; the parameters  +  > 0 and  − V > 0 are the positive constants;  +  ,  −  ,  +  , and  − V denote the fluxes of the corresponding reaction processes; constraint (2) ensures that the system will operate under steady-state conditions (i.e.,   / = 0); inequality constraint (3) forces a flux or the ratio of some two fluxes to remain below a certain limit; and constraint (4) imposes both the internal metabolites   ( = 1, 2, . . ., ) and external metabolites   ( =  + 1,  + 2, . . .,  + ) to stay within certain physically and chemically feasible limits (   > 0).

The GMA Formalism. A biochemical process containing
metabolites is usually modeled as a system of differential equations in which the variation in metabolite concentrations   can be represented as the following stoichiometric formalism: In this representation, both of the fluxes  +  and  −  can be expressed in the following power-law functions, respectively [22,23,28]: where the model parameters  +  > 0 and  −  > 0 are the rate constants for fluxes  +  and  −  , respectively, and  +  ∈  and  −  ∈  are the kinetic orders that reflect the direct effects of any given system variable   on fluxes  +  and  −  , respectively.Note that a function with the form of ( 6) or ( 7) is also called a monomial.A sum of one or more monomials is called a posynomial, while a sum of one or more monomials, possibly with negative multiplicative coefficients, is called a signomial.
By introducing ( 6)-( 7) into (5), one obtains the following GMA model of (5): Based on (8) the objective function () and constraint (3) can also be written as the following power-law forms, respectively: where   ∈ , ℎ +  ∈ , and ℎ − V ∈  terms stand for the kinetic orders and  > 0,  +  > 0, and  − V > 0 represent the corresponding rate constants.Now we can rewrite optimization problems (1)-(4) as the following formulation: subject to satisfying: In this representation, both of the constraints ( 11) and ( 12) involve a special structure in the form of signomial function.This kind of optimization problems as shown in ( 10)-( 13) belongs to a truly nonconvex class of problems known as signomial geometric programming (SGP) [29,30], that is, difficult to solve for global optimality.

Improved Geometric Programming Method.
In this subsection, we propose to convert SGP problem ( 10)-( 13) into a sequence of standard GP problems that can be efficiently solved to reach a global solution.
We first denote where  +  (),  −  (),  +  (), and  −  () are posynomial functions.Then optimization problem ( 10)-( 13) can be written as the following formulation: subject to satisfying In this problem, constraint ( 16) has a special structure in the form of a ratio between two posynomials.This type of constraint gives rise to some difficulties in solving problem ( 15)-( 18) for global optimality.To deal with this issue, we rewrite optimization problem ( 15)-( 18) as subject to satisfying by relaxing equality constraints (16) with constraints ( 20)- (26).Here   > 0 are the weighting factors with sufficiently large values, and index sets  1 ,  2 ,  3 ,  4 ,  1 , and  2 are defined, respectively, as with  = {1, 2, . . ., } and  = {1, 2, . . ., }.It can be easily verified that, as the auxiliary variables   decrease from 1 to 0, the feasible region of problem ( 19)-( 30) decreases.When the variables   = 0, the feasible region of relaxed program ( 19)-( 30) is exactly the one of original problem.The variables   in problem ( 19)- (30) do not meet the implicit constraint that the optimization variables are positive in any standard GP [29], but we can use some transformation techniques to transform them into other positive ones.In this study, we propose to transform variable   into a new one in the following representation: where  is a positive constant and  ≤   < 1+.Then we have the following equivalent representation of problem ( 19)-( 30): subject to satisfying It can be observed that optimization problem ( 33)-( 44) is equivalent to original problem (1)-( 4) when   = .In problem ( 33)-( 44), every denominator of constrains ( 36), ( 38)-(40), and (42) is posynomials.These posynomials can be approximated, respectively, with monomials by using the following arithmetic-geometric mean approximation: where the parameters ,  ∈  2 . (46) Applying the method above to each denominator of constrains (36), (38)-(40), and (42), we obtain the following optimization problem: subject to satisfying (48) This is a standard GP that can be transformed into a nonlinear convex problem.
We now present the following successive geometric programming algorithm for steady-state optimization of biochemical systems in this work.
Step 3. Update the weighting coefficients  ()  with where  is a monotonically increasing function of  (−1)
Remark 1.In the similar thought of [14], we can easily prove that the sequent solutions of problem (47)-(48) converge to a point satisfying the KKT conditions of the original problem.
Remark 2. The proposed algorithm requires fewer iterations to obtain the global optimum of a biochemical system than the approach used in [14] does.See Section 4.
Remark 3. We can borrow from the theory and practice of penalty function methods to select an appropriate weighting factor  ()  .

Case Studies
In this section, two case studies are presented to illustrate the calculation algorithm in terms of computational experiments.These systems were optimized using the MATLAB based software GGPLAB [31] on a PC with Intel (R) Core (TM) i5-3230 M Processor 2.60 GHz CPU.The default solver of software GGPLAB was set to be MOSEK [32].

Case Study 1: Tryptophan Biosynthesis in Escherichia coli.
We first apply the proposed method to tryptophan biosynthesis in Escherichia coli.A complete description of this metabolic pathway can be found in [33].The differential equations in dimensionless variables are given as Here,  1 is used for mRNA concentration,  2 is used for enzyme concentration, and  3 is used for tryptophan concentration, respectively.At the basal steady state (see Table 1), the rate equations in (50) are transformed into the following power-law forms, respectively [14]: .

(54)
The following algorithm parameters were chosen in the implementation of the proposed method:  (0) 3 = 1,  () 3 = 1 +  ( ≥ 1),  = 1,  (0) 3 = 1.95, and  = 10 −6 .The proposed algorithm uses about 1.6-second CPU time to perform 5 iterations and reach an optimal solution with the objective value 5.169027 starting from the initial steady state given in Table 1.The optimized results are given in Table 1.The optimal objective value 5.169027 is higher than the results given by [6].Table 2 provides us with the comparison between the proposed approach and the reported method [14] for Case study 1.The calculation results shown in this table were obtained by setting the following algorithm parameters:  (0) 3 = 1,  () 3 = 1 +  ( ≥ 1),  = 1, and  = 10 −6 .As can be seen in Table 2, both optimization strategies obtain a rate of tryptophan production increased more than 3.96 times its initial steady state.But our proposed approach requires fewer iterations to find the global optimum of Case study 1 than the method in [14] does.One can also conclude from Table 2 that, as the value of initial auxiliary variable 3 increases from 1.0 to 1.99, the iterations required for the proposed algorithm decrease.This suggests that we should choose a larger  (0)  value to promote the evolution rate of the proposed approach.

Case Study 2: Maximization of Ethanol Production in
Saccharomyces cerevisiae.In this case study, we will apply the proposed approach to the production of ethanol by Saccharomyces cerevisiae.The dynamics of this system are described as the following formulations [34]: where   represent the following intermediate metabolite concentrations:  1 is the intracellular glucose concentration,  2 represents glucose 6-phosphate,  3 codes for fructose 1, 6-diphosphate,  4 is phosphoenolpyruvate, and  5 represents ATP.The indexed quantities  represent the following fluxes:  in denotes the sugar transport into the cells,  HK summarizes all hexokinases,  PFK is the phosphofructokinase reaction,  GAPD represents glyceraldehyde 3-phosphate dehydrogenase,  PK represents pyruvate kinase,  Pol describes glycogen synthetase,  Gol , the glycerol 3-phosphate dehydrogenase is proportional to  PK , and  ATPase summarizes collectively the use of ATP.At the initial steady state shown in Table 3, these fluxes can be represented as the following power-law expressions, respectively [35]: The performance index describing the rate of ethanol production is given directly by the flux through the pyruvate kinase,  PK .Then we have the following optimization problem [3,36]: (58) This problem has three signomial equality constrains that can be rewritten as follows: Thus, three weighting coefficients  ()  ( ≥ 0 and  = 2, 3, 5) should be chosen in the implementation of the proposed method.To deal with this issue, we first solve two simple GPs as shown in max  = 0.0945 0.
(61) a The method in [14] does not allow  (0) 3 = 0 because of the auxiliary variable  3 > 0 being imposed in the implementation of this algorithm.The following algorithm parameters were firstly used in the implementation of the proposed algorithm:  ()  = 0.0006 ( ≥ 0),  = 1,  (0)  = 1.99, and  = 10 −6 .The proposed method spends about 3.1-second CPU time to perform 5 iterations and reach an optimal solution with the objective value 1577.4227starting from the initial steady state given in Table 3.The detailed results are shown in Table 3. Table 4 presents a comparison of these results and those obtained by other reported methods [14] for Case study 2. The calculation results presented in this table were attained by using the following algorithm parameters:  ()  = 0.0006 ( ≥ 0),  = 1, and  = 10 −6 .From Table 4, it can be seen that both optimization methods yield a rate of ethanol production increased more than 52.38 times its initial steady state.But the proposed approach in this work requires fewer iterations to find the global optimum of Case study 2 than the method in [14] does.
To illustrate the influence of variable  on the performance of our proposed algorithm, ten computational experiments with different  values were performed with the following algorithm parameters:  ()  = 0.0006 ( ≥ 0),  (0)  =  + 0.9, and  = 10 −6 .Table 5 compares the influence of different  values on the performance of the proposed approach.In this table it can be seen that the proposed algorithm with 10 different  values successfully obtains the same rate of ethanol production increased more than 52.38 times its initial steady state, but it requires fewer iterations to find the global optimum of Case study 2, when we select a larger  value.
Next we investigate the influence of weighting coefficients  ()   ( ≥ 0) on the performance of our proposed algorithm.Ten computational experiments with different weighting coefficients  ()  ( ()  ∈ [0.0006, 0.0433] and  ≥ 0) were performed with the following algorithm parameters:  = 1,  (0)  = 1.99, and  = 10 −6 .Table 6 compares the influence of different weighting coefficients  ()   ( ≥ 0) on the performance of the proposed approach.In this table it can be observed that the proposed algorithm with 10 different weighting coefficients attains the global optimum of Case study 2, but it requires fewer iterations to accomplish this optimization task when the weighting coefficient  ()   ( ≥ 0) is chosen as a smaller value.This suggests that we should choose a small  ()   ( ≥ 0) value in the implementation of the proposed algorithm.b The method in [14] does not allow  (0)  = 0 ( = 2, 3, 5) because of the auxiliary variables   > 0 being imposed in the implementation of this algorithm.

Conclusions
In this work, we have presented an improved GP approach for steady-state optimization of biochemical systems.The original nonlinear, nonconvex optimization problem with GMA formalism can be easily transformed into a set of standard GPs that can be solved very efficiently for global optimality.
The proposed algorithm has been applied to two benchmark systems.Compared with the current GP method for steadystate optimization of biochemical systems, the presented approach rapidly and successfully obtains a globally optimal production rate of these systems.This conclusion shows the tractability and effectiveness of the improved GP approach in solving steady-state optimization of nonlinear biochemical systems with GMA models.

Table 2 :
[14]arisons between our proposed algorithm and the method in[14]for Case study 1.

Table 3 :
Optimal solution of Case study 2 obtained by using the proposed approach.

Table 4 :
[14]arisons between the proposed algorithm and the method in[14]for Case study 2.

Table 5 :
Comparisons of effects of 10 different  values on the proposed algorithm.

Table 6 :
Comparisons of effects of 10 different weighting coefficients on the proposed algorithm.