Abstract

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 GMA models 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.

1. 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 [121]. 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) [2228]. 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, 811, 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 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.

2. Optimization Problem Statement

Consider the following steady-state optimization problem of biochemical systems:

subject to satisfying where ; the objective function is usually a flux or a particular metabolite concentration; the parameters and are the stoichiometric coefficients of the metabolite in the processes and , respectively; the parameters and are the positive constants; , , , and denote the fluxes of the corresponding reaction processes; constraint (2) ensures that the system will operate under steady-state conditions (i.e., ); 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 () and external metabolites () to stay within certain physically and chemically feasible limits ().

3. Optimization Methods

3.1. 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 and 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 terms stand for the kinetic orders and , , and 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.

3.2. 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 are the weighting factors with sufficiently large values, and index sets , , , , , and are defined, respectively, as with and . 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 , 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 . 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 (),  (), , (), (), (),   (), and () can be computed as follows for some given and () values:

Applying the method above to each denominator of constrains (36), (38)–(40), and (42), we obtain the following optimization problem:

subject to satisfying 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  0. Given initial steady-state , positive constant , initial values of auxiliary variables (, ), initial weights , and solution accuracy , set iteration counter .

Step  1. For and , solve the standard GP (47)-(48) to obtain and with weighting factors .

Step  2. If and ( is a vector whose elements are ), then stop.

Step  3. Update the weighting coefficients with where is a monotonically increasing function of . Set and continue from Step  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 .

4. 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].

4.1. 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, is used for mRNA concentration, is used for enzyme concentration, and 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]:

The reaction can be selected as the optimized objective function [5]. Then we have the following steady-state optimization problem:

subject to satisfying

The following algorithm parameters were chosen in the implementation of the proposed method: , , , , and . 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: , , , and . 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 increases from 1.0 to 1.99, the iterations required for the proposed algorithm decrease. This suggests that we should choose a larger value to promote the evolution rate of the proposed approach.

4.2. 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: is the intracellular glucose concentration, represents glucose 6-phosphate, codes for fructose 1, 6-diphosphate, is phosphoenolpyruvate, and represents ATP. The indexed quantities represent the following fluxes: denotes the sugar transport into the cells, summarizes all hexokinases, is the phosphofructokinase reaction, represents glyceraldehyde 3-phosphate dehydrogenase, represents pyruvate kinase, describes glycogen synthetase, , the glycerol 3-phosphate dehydrogenase is proportional to , and 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, . Then we have the following optimization problem [3, 36]:

subject to satisfying This problem has three signomial equality constrains that can be rewritten as follows: Thus, three weighting coefficients ( and ) should be chosen in the implementation of the proposed method. To deal with this issue, we first solve two simple GPs as shown in The optimal objective values of optimization problems (60) can be easily obtained as follows, respectively: Then we get an approximate interval for objective function of original problem (58). So we choose the weighting coefficients () as where .

The following algorithm parameters were firstly used in the implementation of the proposed algorithm: , , , and . The proposed method spends about 3.1-second CPU time to perform 5 iterations and reach an optimal solution with the objective value 1577.4227 starting 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: , , and . 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: , , and . 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 on the performance of our proposed algorithm. Ten computational experiments with different weighting coefficients ( and ) were performed with the following algorithm parameters: , , and . Table 6 compares the influence of different weighting coefficients 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 is chosen as a smaller value. This suggests that we should choose a small value in the implementation of the proposed algorithm.

5. 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 steady-state 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (nos. 11101051 and 11371071) and the Program for Liaoning Excellent Talents in University (no. LJQ2013115).