# An Improved Geometric Programming Approach for Optimization of Biochemical Systems

Authors: Gongxian Xu; Lei Wang
Journal: Journal of Applied Mathematics
 (2014)

---

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


---

## Body

## 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 [1–21]. 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–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 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. Section2 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:(1)
max
⁡
⁡
F
(
X
)
,subject to satisfying(2)
∑
j
=
1
p
i
+
μ
i
j
+
V
i
j
+
(
X
)
-
∑
k
=
1
p
i
-
μ
i
k
-
V
i
k
-
(
X
)
=
0
,
i
=
1,2
,
…
,
n
,
(3)
∑
u
=
1
q
l
+
λ
l
u
+
V
l
u
+
(
X
)
-
∑
v
=
1
q
l
-
λ
l
v
-
V
l
v
-
(
X
)
≤
0
,
l
=
1,2
,
…
,
d
,
(4)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
n
+
m
,
where X
=
(
X
1
,
X
2
,
…
,
X
n
+
m
)
T
∈
R
n
+
m; the objective function F
(
X
) is usually a flux or a particular metabolite concentration; the parameters μ
i
j
+
>
0 and μ
i
k
-
>
0 are the stoichiometric coefficients of the metabolite X
i in the processes V
i
j
+ and V
i
k
-, respectively; the parameters λ
l
u
+
>
0 and λ
l
v
-
>
0 are the positive constants; V
i
j
+, V
i
k
-, V
l
u
+, and V
l
v
- denote the fluxes of the corresponding reaction processes; constraint (2) ensures that the system will operate under steady-state conditions (i.e., d
X
i
/
d
t
=
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 X
e (e
=
1,2
,
…
,
n) and external metabolites X
e (e
=
n
+
1
,
n
+
2
,
…
,
n
+
m) to stay within certain physically and chemically feasible limits (X
e
L
>
0).


## 3. Optimization Methods

### 3.1. The GMA Formalism

A biochemical process containingn metabolites is usually modeled as a system of differential equations in which the variation in metabolite concentrations X
i can be represented as the following stoichiometric formalism:
(5)
d
X
i
d
t
=
∑
j
=
1
p
i
+
μ
i
j
+
V
i
j
+
(
X
)
-
∑
k
=
1
p
i
-
μ
i
k
-
V
i
k
-
(
X
)
,
i
=
1,2
,
…
,
n
.
In this representation, both of the fluxes V
i
j
+ and V
i
k
- can be expressed in the following power-law functions, respectively [22, 23, 28]:
(6)
V
i
j
+
(
X
)
=
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
,
j
=
1,2
,
…
,
p
i
+
,
(7)
V
i
k
-
(
X
)
=
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
,
k
=
1,2
,
…
,
p
i
-
,
where the model parameters α
i
j
+
>
0 and α
i
k
-
>
0 are the rate constants for fluxes V
i
j
+ and V
i
k
-, respectively, and g
i
j
e
+
∈
R and g
i
k
e
-
∈
R are the kinetic orders that reflect the direct effects of any given system variable X
e on fluxes V
i
j
+ and V
i
k
-, 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):
(8)
d
X
i
d
t
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
-
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
,
i
=
1,2
,
…
,
n
.Based on (8) the objective function F
(
X
) and constraint (3) can also be written as the following power-law forms, respectively:
(9)
F
(
X
)
=
γ
∏
e
=
1
n
+
m
X
e
f
e
,
∑
u
=
1
q
l
+
(
λ
l
u
+
β
l
u
+
∏
e
=
1
n
+
m
X
e
h
l
u
e
+
)
-
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
≤
0
,
l
=
1,2
,
…
,
d
,
where f
e
∈
R, h
l
u
e
+
∈
R, and h
l
v
e
-
∈
R terms stand for the kinetic orders and γ
>
0, β
l
u
+
>
0, and β
l
v
-
>
0 represent the corresponding rate constants.Now we can rewrite optimization problems (1)–(4) as the following formulation:
(10)
max
⁡
⁡
F
(
X
)
=
γ
∏
e
=
1
n
+
m
X
e
f
e
,subject to satisfying:(11)
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
-
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
=
0
,
i
=
1,2
,
…
,
n
,
(12)
∑
u
=
1
q
l
+
(
λ
l
u
+
β
l
u
+
∏
e
=
1
n
+
m
X
e
h
l
u
e
+
)
-
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
≤
0
,
l
=
1,2
,
…
,
d
,
(13)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
n
+
m
.
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(14)
H
i
+
(
X
)
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
,
i
=
1,2
,
…
,
n
,
H
i
-
(
X
)
=
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
,
i
=
1,2
,
…
,
n
,
G
l
+
(
X
)
=
∑
u
=
1
q
l
+
(
λ
l
u
+
β
l
u
+
∏
e
=
1
n
+
m
X
e
h
l
u
e
+
)
,
l
=
1,2
,
…
,
d
,
G
l
-
(
X
)
=
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
,
l
=
1,2
,
…
,
d
,
where H
i
+
(
X
), H
i
-
(
X
), G
l
+
(
X
), and G
l
-
(
X
) are posynomial functions. Then optimization problem (10)–(13) can be written as the following formulation:
(15)
max
⁡
⁡
F
(
X
)
=
γ
∏
e
=
1
n
+
m
X
e
f
e
,subject to satisfying(16)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
=
1,2
,
…
,
n
,
(17)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
=
1,2
,
…
,
d
,
(18)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
n
+
m
.
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
(19)
min
⁡
1
γ
∏
e
=
1
n
+
m
X
e
-
f
e
+
∑
i
∈
I
2
∪
I
3
∪
I
4
w
i
t
i
,subject to satisfying(20)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
∈
I
1
,
(21)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
2
,
(22)
H
i
+
(
X
)
H
i
-
(
X
)
≥
1
-
t
i
,
i
∈
I
2
,
(23)
H
i
-
(
X
)
H
i
+
(
X
)
≤
1
,
i
∈
I
3
,
(24)
H
i
-
(
X
)
H
i
+
(
X
)
≥
1
-
t
i
,
i
∈
I
3
,
(25)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
4
,
(26)
H
i
+
(
X
)
H
i
-
(
X
)
≥
1
-
t
i
,
i
∈
I
4
,
(27)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
1
,
(28)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
2
,
(29)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
m
+
n
,
(30)
0
≤
t
i
<
1
,
i
∈
I
2
∪
I
3
∪
I
4
,
by relaxing equality constraints (16) with constraints (20)–(26). Here w
i
>
0 are the weighting factors with sufficiently large values, and index sets I
1, I
2, I
3, I
4, L
1, and L
2 are defined, respectively, as
(31)
I
1
=
{
i
∣
i
∈
I
,
p
i
+
=
1
,
p
i
-
=
1
}
,
I
2
=
{
i
∣
i
∈
I
,
p
i
+
≥
2
,
p
i
-
=
1
}
,
I
3
=
{
i
∣
i
∈
I
,
p
i
+
=
1
,
p
i
-
≥
2
}
,
I
4
=
{
i
∣
i
∈
I
,
p
i
+
≥
2
,
p
i
-
≥
2
}
,
L
1
=
{
l
∣
l
∈
L
,
q
l
-
=
1
}
,
L
2
=
{
l
∣
l
∈
L
,
q
l
-
≥
2
}
,
with I
=
{
1,2
,
…
,
n
} and L
=
{
1,2
,
…
,
d
}. It can be easily verified that, as the auxiliary variables t
i decrease from 1 to 0, the feasible region of problem (19)–(30) decreases. When the variables t
i
=
0, the feasible region of relaxed program (19)–(30) is exactly the one of original problem.The variablest
i 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 t
i into a new one in the following representation:
(32)
s
i
=
t
i
+
θ
,
i
∈
I
2
∪
I
3
∪
I
4
,
where θ is a positive constant and θ
≤
s
i
<
1
+
θ. Then we have the following equivalent representation of problem (19)–(30):
(33)
min
⁡
⁡
1
γ
∏
e
=
1
n
+
m
X
e
-
f
e
+
∑
i
∈
I
2
∪
I
3
∪
I
4
w
i
s
i
,subject to satisfying(34)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
∈
I
1
(35)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
2
(36)
(
1
+
θ
)
H
i
-
(
X
)
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
≤
1
,
i
∈
I
2
(37)
H
i
-
(
X
)
H
i
+
(
X
)
≤
1
,
i
∈
I
3
(38)
(
1
+
θ
)
H
i
+
(
X
)
H
i
-
(
X
)
+
s
i
H
i
+
(
X
)
≤
1
,
i
∈
I
3
(39)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
4
(40)
(
1
+
θ
)
H
i
-
(
X
)
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
≤
1
,
i
∈
I
4
(41)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
1
(42)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
2
(43)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
m
+
n
(44)
θ
≤
s
i
<
1
+
θ
,
i
∈
I
2
∪
I
3
∪
I
4
.
It can be observed that optimization problem (33)–(44) is equivalent to original problem (1)–(4) when s
i
=
θ. 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:
(45)
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
+
μ
i
1
-
α
i
1
-
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
-
≥
∏
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
ρ
i
j
+
)
ρ
i
j
+
·
(
μ
i
1
-
α
i
1
-
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
-
ρ
i
1
-
)
ρ
i
1
-
=
H
^
2
i
+
(
X
,
s
i
)
,
i
∈
I
2
,
H
i
-
(
X
)
+
s
i
H
i
+
(
X
)
=
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
+
μ
i
1
+
α
i
1
+
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
+
≥
∏
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
ρ
i
k
-
)
ρ
i
k
-
·
(
μ
i
1
+
α
i
1
+
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
+
ρ
i
1
+
)
ρ
i
1
+
=
H
^
3
i
-
(
X
,
s
i
)
,
i
∈
I
3
,
H
i
-
(
X
)
=
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
≥
∏
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
ρ
i
k
-
)
ρ
i
k
-
=
H
^
4
i
-
(
X
,
s
i
)
,
i
∈
I
4
,
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
+
s
i
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
≥
∏
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
ρ
i
j
+
)
ρ
i
j
+
·
∏
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
s
i
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
ρ
i
k
s
-
)
ρ
i
k
s
-
=
H
^
4
i
+
(
X
,
s
i
)
,
i
∈
I
4
,
G
l
-
(
X
)
=
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
≥
∏
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
ρ
l
v
-
)
ρ
l
v
-
=
G
^
l
-
(
X
)
,
l
∈
L
2
,
where the parameters ρ
i
j
+ (i
∈
I
2),  ρ
i
1
- (i
∈
I
2), ρ
i
k
-
(
i
∈
I
3
), ρ
i
1
+ (i
∈
I
3), ρ
i
k
- (i
∈
I
4), ρ
i
j
+ (i
∈
I
4),   ρ
i
k
s
- (i
∈
I
4), and ρ
l
v
- (l
∈
L
2) can be computed as follows for some given X
- and s
-
i (i
∈
I
2
∪
I
3
∪
I
4) values:
(46)
ρ
i
j
+
=
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
-
e
g
i
j
e
+
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
2
,
ρ
i
1
-
=
μ
i
1
-
α
i
1
-
s
-
i
∏
e
=
1
n
+
m
X
-
e
g
i
1
e
-
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
2
,
ρ
i
k
-
=
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
-
e
g
i
k
e
-
H
i
-
(
X
-
)
+
s
-
i
H
i
+
(
X
-
)
,
i
∈
I
3
,
ρ
i
1
+
=
μ
i
1
+
α
i
1
+
s
-
i
∏
e
=
1
n
+
m
X
-
e
g
i
1
e
+
H
i
-
(
X
-
)
+
s
-
i
H
i
+
(
X
-
)
,
i
∈
I
3
,
ρ
i
k
-
=
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
-
e
g
i
k
e
-
H
i
-
(
X
-
)
,
i
∈
I
4
,
ρ
i
j
+
=
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
-
e
g
i
j
e
+
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
4
,
ρ
i
k
s
-
=
μ
i
k
-
α
i
k
-
s
-
i
∏
e
=
1
n
+
m
X
-
e
g
i
k
e
-
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
4
,
ρ
l
v
-
=
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
-
e
h
l
v
e
-
G
l
-
(
X
-
)
,
l
∈
L
2
.Applying the method above to each denominator of constrains (36), (38)–(40), and (42), we obtain the following optimization problem:
(47)
min
⁡
⁡
1
γ
∏
e
=
1
n
+
m
X
e
-
f
e
+
∑
i
∈
I
2
∪
I
3
∪
I
4
w
i
s
i
,subject to satisfying(48)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
∈
I
1
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
2
(
1
+
θ
)
H
i
-
(
X
)
H
^
2
i
+
(
X
,
s
i
)
≤
1
,
i
∈
I
2
H
i
-
(
X
)
H
i
+
(
X
)
≤
1
,
i
∈
I
3
(
1
+
θ
)
H
i
+
(
X
)
H
^
3
i
-
(
X
,
s
i
)
≤
1
,
i
∈
I
3
H
i
+
(
X
)
H
^
4
i
-
(
X
,
s
i
)
≤
1
,
i
∈
I
4
(
1
+
θ
)
H
i
-
(
X
)
H
^
4
i
+
(
X
,
s
i
)
≤
1
,
i
∈
I
4
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
1
G
l
+
(
X
)
G
^
l
-
(
X
)
≤
1
,
l
∈
L
2
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
m
+
n
θ
≤
s
i
<
1
+
θ
,
i
∈
I
2
∪
I
3
∪
I
4
.
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 X
(
0
), positive constant θ, initial values of auxiliary variables s
i
(
0
) (θ
≤
s
i
(
0
)
<
1
+
θ, i
∈
I
2
∪
I
3
∪
I
4), initial weights w
i
(
0
), and solution accuracy ε
>
0, set iteration counter r
=
0.Step  1.For X
-
=
X
(
r
-
1
) and s
-
i
=
s
i
(
r
-
1
), solve the standard GP (47)-(48) to obtain X
(
r
) and s
i
(
r
) with weighting factors w
i
(
r
-
1
).Step  2. If ∥
X
(
r
)
-
X
(
r
-
1
)
∥
≤
ε and ∥
s
(
r
)
-
s
(
r
-
1
)
∥
≤
ε (s is a vector whose elements are s
i), then stop.Step  3.Update the weighting coefficients w
i
(
r
) with(49)
w
i
(
r
)
=
W
(
w
i
(
r
-
1
)
)
,
i
∈
I
2
∪
I
3
∪
I
4
,
where W is a monotonically increasing function of w
i
(
r
-
1
). Set r
=
r
+
1 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 factorw
i
(
r
).


## 3.1. The GMA Formalism

A biochemical process containingn metabolites is usually modeled as a system of differential equations in which the variation in metabolite concentrations X
i can be represented as the following stoichiometric formalism:
(5)
d
X
i
d
t
=
∑
j
=
1
p
i
+
μ
i
j
+
V
i
j
+
(
X
)
-
∑
k
=
1
p
i
-
μ
i
k
-
V
i
k
-
(
X
)
,
i
=
1,2
,
…
,
n
.
In this representation, both of the fluxes V
i
j
+ and V
i
k
- can be expressed in the following power-law functions, respectively [22, 23, 28]:
(6)
V
i
j
+
(
X
)
=
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
,
j
=
1,2
,
…
,
p
i
+
,
(7)
V
i
k
-
(
X
)
=
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
,
k
=
1,2
,
…
,
p
i
-
,
where the model parameters α
i
j
+
>
0 and α
i
k
-
>
0 are the rate constants for fluxes V
i
j
+ and V
i
k
-, respectively, and g
i
j
e
+
∈
R and g
i
k
e
-
∈
R are the kinetic orders that reflect the direct effects of any given system variable X
e on fluxes V
i
j
+ and V
i
k
-, 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):
(8)
d
X
i
d
t
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
-
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
,
i
=
1,2
,
…
,
n
.Based on (8) the objective function F
(
X
) and constraint (3) can also be written as the following power-law forms, respectively:
(9)
F
(
X
)
=
γ
∏
e
=
1
n
+
m
X
e
f
e
,
∑
u
=
1
q
l
+
(
λ
l
u
+
β
l
u
+
∏
e
=
1
n
+
m
X
e
h
l
u
e
+
)
-
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
≤
0
,
l
=
1,2
,
…
,
d
,
where f
e
∈
R, h
l
u
e
+
∈
R, and h
l
v
e
-
∈
R terms stand for the kinetic orders and γ
>
0, β
l
u
+
>
0, and β
l
v
-
>
0 represent the corresponding rate constants.Now we can rewrite optimization problems (1)–(4) as the following formulation:
(10)
max
⁡
⁡
F
(
X
)
=
γ
∏
e
=
1
n
+
m
X
e
f
e
,subject to satisfying:(11)
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
-
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
=
0
,
i
=
1,2
,
…
,
n
,
(12)
∑
u
=
1
q
l
+
(
λ
l
u
+
β
l
u
+
∏
e
=
1
n
+
m
X
e
h
l
u
e
+
)
-
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
≤
0
,
l
=
1,2
,
…
,
d
,
(13)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
n
+
m
.
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(14)
H
i
+
(
X
)
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
,
i
=
1,2
,
…
,
n
,
H
i
-
(
X
)
=
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
,
i
=
1,2
,
…
,
n
,
G
l
+
(
X
)
=
∑
u
=
1
q
l
+
(
λ
l
u
+
β
l
u
+
∏
e
=
1
n
+
m
X
e
h
l
u
e
+
)
,
l
=
1,2
,
…
,
d
,
G
l
-
(
X
)
=
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
,
l
=
1,2
,
…
,
d
,
where H
i
+
(
X
), H
i
-
(
X
), G
l
+
(
X
), and G
l
-
(
X
) are posynomial functions. Then optimization problem (10)–(13) can be written as the following formulation:
(15)
max
⁡
⁡
F
(
X
)
=
γ
∏
e
=
1
n
+
m
X
e
f
e
,subject to satisfying(16)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
=
1,2
,
…
,
n
,
(17)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
=
1,2
,
…
,
d
,
(18)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
n
+
m
.
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
(19)
min
⁡
1
γ
∏
e
=
1
n
+
m
X
e
-
f
e
+
∑
i
∈
I
2
∪
I
3
∪
I
4
w
i
t
i
,subject to satisfying(20)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
∈
I
1
,
(21)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
2
,
(22)
H
i
+
(
X
)
H
i
-
(
X
)
≥
1
-
t
i
,
i
∈
I
2
,
(23)
H
i
-
(
X
)
H
i
+
(
X
)
≤
1
,
i
∈
I
3
,
(24)
H
i
-
(
X
)
H
i
+
(
X
)
≥
1
-
t
i
,
i
∈
I
3
,
(25)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
4
,
(26)
H
i
+
(
X
)
H
i
-
(
X
)
≥
1
-
t
i
,
i
∈
I
4
,
(27)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
1
,
(28)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
2
,
(29)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
m
+
n
,
(30)
0
≤
t
i
<
1
,
i
∈
I
2
∪
I
3
∪
I
4
,
by relaxing equality constraints (16) with constraints (20)–(26). Here w
i
>
0 are the weighting factors with sufficiently large values, and index sets I
1, I
2, I
3, I
4, L
1, and L
2 are defined, respectively, as
(31)
I
1
=
{
i
∣
i
∈
I
,
p
i
+
=
1
,
p
i
-
=
1
}
,
I
2
=
{
i
∣
i
∈
I
,
p
i
+
≥
2
,
p
i
-
=
1
}
,
I
3
=
{
i
∣
i
∈
I
,
p
i
+
=
1
,
p
i
-
≥
2
}
,
I
4
=
{
i
∣
i
∈
I
,
p
i
+
≥
2
,
p
i
-
≥
2
}
,
L
1
=
{
l
∣
l
∈
L
,
q
l
-
=
1
}
,
L
2
=
{
l
∣
l
∈
L
,
q
l
-
≥
2
}
,
with I
=
{
1,2
,
…
,
n
} and L
=
{
1,2
,
…
,
d
}. It can be easily verified that, as the auxiliary variables t
i decrease from 1 to 0, the feasible region of problem (19)–(30) decreases. When the variables t
i
=
0, the feasible region of relaxed program (19)–(30) is exactly the one of original problem.The variablest
i 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 t
i into a new one in the following representation:
(32)
s
i
=
t
i
+
θ
,
i
∈
I
2
∪
I
3
∪
I
4
,
where θ is a positive constant and θ
≤
s
i
<
1
+
θ. Then we have the following equivalent representation of problem (19)–(30):
(33)
min
⁡
⁡
1
γ
∏
e
=
1
n
+
m
X
e
-
f
e
+
∑
i
∈
I
2
∪
I
3
∪
I
4
w
i
s
i
,subject to satisfying(34)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
∈
I
1
(35)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
2
(36)
(
1
+
θ
)
H
i
-
(
X
)
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
≤
1
,
i
∈
I
2
(37)
H
i
-
(
X
)
H
i
+
(
X
)
≤
1
,
i
∈
I
3
(38)
(
1
+
θ
)
H
i
+
(
X
)
H
i
-
(
X
)
+
s
i
H
i
+
(
X
)
≤
1
,
i
∈
I
3
(39)
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
4
(40)
(
1
+
θ
)
H
i
-
(
X
)
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
≤
1
,
i
∈
I
4
(41)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
1
(42)
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
2
(43)
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
m
+
n
(44)
θ
≤
s
i
<
1
+
θ
,
i
∈
I
2
∪
I
3
∪
I
4
.
It can be observed that optimization problem (33)–(44) is equivalent to original problem (1)–(4) when s
i
=
θ. 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:
(45)
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
+
μ
i
1
-
α
i
1
-
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
-
≥
∏
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
ρ
i
j
+
)
ρ
i
j
+
·
(
μ
i
1
-
α
i
1
-
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
-
ρ
i
1
-
)
ρ
i
1
-
=
H
^
2
i
+
(
X
,
s
i
)
,
i
∈
I
2
,
H
i
-
(
X
)
+
s
i
H
i
+
(
X
)
=
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
+
μ
i
1
+
α
i
1
+
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
+
≥
∏
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
ρ
i
k
-
)
ρ
i
k
-
·
(
μ
i
1
+
α
i
1
+
s
i
∏
e
=
1
n
+
m
X
e
g
i
1
e
+
ρ
i
1
+
)
ρ
i
1
+
=
H
^
3
i
-
(
X
,
s
i
)
,
i
∈
I
3
,
H
i
-
(
X
)
=
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
≥
∏
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
ρ
i
k
-
)
ρ
i
k
-
=
H
^
4
i
-
(
X
,
s
i
)
,
i
∈
I
4
,
H
i
+
(
X
)
+
s
i
H
i
-
(
X
)
=
∑
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
)
+
s
i
∑
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
)
≥
∏
j
=
1
p
i
+
(
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
e
g
i
j
e
+
ρ
i
j
+
)
ρ
i
j
+
·
∏
k
=
1
p
i
-
(
μ
i
k
-
α
i
k
-
s
i
∏
e
=
1
n
+
m
X
e
g
i
k
e
-
ρ
i
k
s
-
)
ρ
i
k
s
-
=
H
^
4
i
+
(
X
,
s
i
)
,
i
∈
I
4
,
G
l
-
(
X
)
=
∑
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
)
≥
∏
v
=
1
q
l
-
(
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
e
h
l
v
e
-
ρ
l
v
-
)
ρ
l
v
-
=
G
^
l
-
(
X
)
,
l
∈
L
2
,
where the parameters ρ
i
j
+ (i
∈
I
2),  ρ
i
1
- (i
∈
I
2), ρ
i
k
-
(
i
∈
I
3
), ρ
i
1
+ (i
∈
I
3), ρ
i
k
- (i
∈
I
4), ρ
i
j
+ (i
∈
I
4),   ρ
i
k
s
- (i
∈
I
4), and ρ
l
v
- (l
∈
L
2) can be computed as follows for some given X
- and s
-
i (i
∈
I
2
∪
I
3
∪
I
4) values:
(46)
ρ
i
j
+
=
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
-
e
g
i
j
e
+
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
2
,
ρ
i
1
-
=
μ
i
1
-
α
i
1
-
s
-
i
∏
e
=
1
n
+
m
X
-
e
g
i
1
e
-
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
2
,
ρ
i
k
-
=
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
-
e
g
i
k
e
-
H
i
-
(
X
-
)
+
s
-
i
H
i
+
(
X
-
)
,
i
∈
I
3
,
ρ
i
1
+
=
μ
i
1
+
α
i
1
+
s
-
i
∏
e
=
1
n
+
m
X
-
e
g
i
1
e
+
H
i
-
(
X
-
)
+
s
-
i
H
i
+
(
X
-
)
,
i
∈
I
3
,
ρ
i
k
-
=
μ
i
k
-
α
i
k
-
∏
e
=
1
n
+
m
X
-
e
g
i
k
e
-
H
i
-
(
X
-
)
,
i
∈
I
4
,
ρ
i
j
+
=
μ
i
j
+
α
i
j
+
∏
e
=
1
n
+
m
X
-
e
g
i
j
e
+
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
4
,
ρ
i
k
s
-
=
μ
i
k
-
α
i
k
-
s
-
i
∏
e
=
1
n
+
m
X
-
e
g
i
k
e
-
H
i
+
(
X
-
)
+
s
-
i
H
i
-
(
X
-
)
,
i
∈
I
4
,
ρ
l
v
-
=
λ
l
v
-
β
l
v
-
∏
e
=
1
n
+
m
X
-
e
h
l
v
e
-
G
l
-
(
X
-
)
,
l
∈
L
2
.Applying the method above to each denominator of constrains (36), (38)–(40), and (42), we obtain the following optimization problem:
(47)
min
⁡
⁡
1
γ
∏
e
=
1
n
+
m
X
e
-
f
e
+
∑
i
∈
I
2
∪
I
3
∪
I
4
w
i
s
i
,subject to satisfying(48)
H
i
+
(
X
)
H
i
-
(
X
)
=
1
,
i
∈
I
1
H
i
+
(
X
)
H
i
-
(
X
)
≤
1
,
i
∈
I
2
(
1
+
θ
)
H
i
-
(
X
)
H
^
2
i
+
(
X
,
s
i
)
≤
1
,
i
∈
I
2
H
i
-
(
X
)
H
i
+
(
X
)
≤
1
,
i
∈
I
3
(
1
+
θ
)
H
i
+
(
X
)
H
^
3
i
-
(
X
,
s
i
)
≤
1
,
i
∈
I
3
H
i
+
(
X
)
H
^
4
i
-
(
X
,
s
i
)
≤
1
,
i
∈
I
4
(
1
+
θ
)
H
i
-
(
X
)
H
^
4
i
+
(
X
,
s
i
)
≤
1
,
i
∈
I
4
G
l
+
(
X
)
G
l
-
(
X
)
≤
1
,
l
∈
L
1
G
l
+
(
X
)
G
^
l
-
(
X
)
≤
1
,
l
∈
L
2
X
e
L
≤
X
e
≤
X
e
U
,
e
=
1,2
,
…
,
m
+
n
θ
≤
s
i
<
1
+
θ
,
i
∈
I
2
∪
I
3
∪
I
4
.
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 X
(
0
), positive constant θ, initial values of auxiliary variables s
i
(
0
) (θ
≤
s
i
(
0
)
<
1
+
θ, i
∈
I
2
∪
I
3
∪
I
4), initial weights w
i
(
0
), and solution accuracy ε
>
0, set iteration counter r
=
0.Step  1.For X
-
=
X
(
r
-
1
) and s
-
i
=
s
i
(
r
-
1
), solve the standard GP (47)-(48) to obtain X
(
r
) and s
i
(
r
) with weighting factors w
i
(
r
-
1
).Step  2. If ∥
X
(
r
)
-
X
(
r
-
1
)
∥
≤
ε and ∥
s
(
r
)
-
s
(
r
-
1
)
∥
≤
ε (s is a vector whose elements are s
i), then stop.Step  3.Update the weighting coefficients w
i
(
r
) with(49)
w
i
(
r
)
=
W
(
w
i
(
r
-
1
)
)
,
i
∈
I
2
∪
I
3
∪
I
4
,
where W is a monotonically increasing function of w
i
(
r
-
1
). Set r
=
r
+
1 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 factorw
i
(
r
).


## 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 inEscherichia coli

We first apply the proposed method to tryptophan biosynthesis inEscherichia coli. A complete description of this metabolic pathway can be found in [33]. The differential equations in dimensionless variables are given as
(50)
d
X
1
d
t
=
V
11
+
-
V
11
-
,
d
X
2
d
t
=
V
21
+
-
V
21
-
,
d
X
3
d
t
=
V
31
+
-
(
V
31
-
+
V
32
-
+
V
33
-
)
.
Here, X
1 is used for mRNA concentration, X
2 is used for enzyme concentration, and X
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]:
(51)
V
11
+
=
0.6403
X
3
-
5.87
×
1
0
-
4
X
5
-
0.8332
,
V
11
-
=
1.0233
X
1
X
4
0.0035
X
11
0.9965
,
V
21
+
=
X
1
,
V
21
-
=
1.4854
X
2
X
4
0.1349
X
12
0.8651
,
V
31
+
=
0.5534
X
2
X
3
-
0.5573
X
6
0.5573
,
V
31
-
=
X
3
X
4
,
V
32
-
=
0.9942
X
3
7.0426
×
1
0
-
4
X
7
,
V
33
-
=
0.8925
X
3
3.5
×
1
0
-
6
X
4
0.9760
X
8
X
9
-
0.0240
X
10
-
3.5
×
1
0
-
6
.Table 1
Optimal solution of Case study 1 obtained by using the proposed approach (s
3
(
0
)
=
1.99).
Variables
Initial steady state
Optimized solution
X
1
0.184654
0.221585
X
2
7.986756
8.908299
X
3
1418.931944
1135.145557
X
4
0.00312
0.005361
X
5
5
4.010578
X
6
2283
5000
X
8
430
1000
F
1.310202
5.169027The reactionV
33
- can be selected as the optimized objective function [5]. Then we have the following steady-state optimization problem:
(52)
max
⁡
⁡
F
=
V
33
-
,subject to satisfying(53)
V
11
+
-
V
11
-
=
0
,
V
21
+
-
V
21
-
=
0
,
V
31
+
-
(
V
31
-
+
V
32
-
+
V
33
-
)
=
0
,
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3
,
0
<
X
4
≤
0.00624
,
(54)
4
≤
X
5
≤
10
,
500
≤
X
6
≤
5000
,
X
7
=
0.0022
X
5
,
0
<
X
8
≤
1000
(
X
9
,
X
10
,
X
11
,
X
12
,
X
13
)
=
(
7.5,0.005,0.9,0.02,0
)
.The following algorithm parameters were chosen in the implementation of the proposed method:w
3
(
0
)
=
1, w
3
(
r
)
=
1
+
r
(
r
≥
1
), θ
=
1, s
3
(
0
)
=
1.95, and ε
=
1
0
-
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: w
3
(
0
)
=
1, w
3
(
r
)
=
1
+
r
(
r
≥
1
), θ
=
1, and ε
=
1
0
-
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 s
3
(
0
) increases from 1.0 to 1.99, the iterations required for the proposed algorithm decrease. This suggests that we should choose a larger s
i
(
0
) value to promote the evolution rate of the proposed approach.Table 2
Comparisons between our proposed algorithm and the method in [14] for Case study 1.
Proposed method
[14]
Initial auxiliaryvariable s
3
(
0
)
Iterations
Objective value
Initial auxiliaryvariable t
3
(
0
)
Iterations
Objective value
1
9
5.169027
0a
—
—
1.1
9
5.169027
0.1
14
5.169027
1.2
5
5.169027
0.2
14
5.169027
1.3
5
5.169027
0.3
14
5.169027
1.4
5
5.169027
0.4
14
5.169027
1.5
5
5.169027
0.5
14
5.169027
1.6
5
5.169027
0.6
14
5.169027
1.7
5
5.169027
0.7
14
5.169027
1.8
5
5.169027
0.8
14
5.169027
1.9
5
5.169027
0.9
14
5.169027
1.99
5
5.169027
0.99
14
5.169027
aThe method in [14] does not allow t
3
(
0
)
=
0 because of the auxiliary variable t
3
>
0 being imposed in the implementation of this algorithm.


### 4.2. Case Study 2: Maximization of Ethanol Production inSaccharomyces cerevisiae

In this case study, we will apply the proposed approach to the production of ethanol bySaccharomyces cerevisiae. The dynamics of this system are described as the following formulations [34]:
(55)
d
X
1
d
t
=
V
in
-
V
HK
,
d
X
2
d
t
=
V
HK
-
V
PFK
-
V
Pol
,
d
X
3
d
t
=
V
PFK
-
V
GAPD
-
0.5
V
Gol
,
d
X
4
d
t
=
2
V
GAPD
-
V
PK
,
d
X
5
d
t
=
2
V
GAPD
+
V
PK
-
V
HK
-
V
Pol
-
V
PFK
-
V
ATPase
,
where X
i represent the following intermediate metabolite concentrations: X
1 is the intracellular glucose concentration, X
2 represents glucose 6-phosphate, X
3 codes for fructose 1, 6-diphosphate, X
4 is phosphoenolpyruvate, and X
5 represents ATP. The indexed quantities V represent the following fluxes: V
in denotes the sugar transport into the cells, V
HK summarizes all hexokinases, V
PFK is the phosphofructokinase reaction, V
GAPD represents glyceraldehyde 3-phosphate dehydrogenase, V
PK represents pyruvate kinase, V
Pol describes glycogen synthetase, V
Gol, the glycerol 3-phosphate dehydrogenase is proportional to V
PK, and V
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]:
(56)
V
in
=
0.8122
X
2
-
0.2344
X
6
,
V
HK
=
2.86
32
X
1
0.7464
X
5
0.0243
X
7
,
V
PFK
=
0.52
32
X
2
0.7318
X
5
-
0.3941
X
8
,
V
Pol
=
0.0009
X
2
8.6107
X
11
,
V
GAPD
=
0.011
X
3
0.6159
X
5
0.1308
X
9
X
14
-
0.6088
,
V
Gol
=
0.04725
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
12
,
V
PK
=
0.0945
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
10
,
V
ATPase
=
X
5
X
13
.Table 3
Optimal solution of Case study 2 obtained by using the proposed approach.
Variables
Initial steady state
Optimized solution
X
1
0.0345
0.038025
X
2
1.0111
1.057148
X
3
9.1437
10.393273
X
4
0.0095
0.010454
X
5
1.1278
1.255394
X
6
19.7
985
X
7
68.5
3147.937162
X
8
31.7
1585
X
9
49.9
2389.028728
X
10
3440
172000
X
13
25.1
1255
F
30.1124
1577.4227The performance index describing the rate of ethanol production is given directly by the flux through the pyruvate kinase,V
PK. Then we have the following optimization problem [3, 36]:
(57)
max
⁡
⁡
F
=
V
PK
,subject to satisfying(58)
V
in
-
V
HK
=
0
,
V
HK
-
V
PFK
-
V
Pol
=
0
,
V
PFK
-
V
GAPD
-
0.5
V
Gol
=
0
,
2
V
GAPD
-
V
PK
=
0
,
2
V
GAPD
+
V
PK
-
V
HK
-
V
Pol
-
V
PFK
-
V
ATPase
=
0
,
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3,4
,
5
,
X
e
(
0
)
≤
X
e
≤
50
X
e
(
0
)
,
e
=
6,7
,
8,9
,
10,1
3
,
V
PK
≤
2
V
in
,
(
X
11
,
X
12
,
X
14
)
=
(
14.31,203,0.042
)
.
This problem has three signomial equality constrains that can be rewritten as follows:
(59)
V
PFK
+
V
Pol
V
HK
=
1
,
V
GAPD
+
0.5
V
Gol
V
PFK
=
1
,
V
HK
+
V
Pol
+
V
PFK
+
V
ATPase
2
V
GAPD
+
V
PK
=
1
.
Thus, three weighting coefficients w
i
(
r
) (r
≥
0 and i
=
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
(60)
max
⁡
⁡
F
=
0.0945
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
10
subject
to
satisfying
:
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3,4
,
5
,
X
e
(
0
)
≤
X
e
≤
50
X
e
(
0
)
,
e
=
6,7
,
8,9
,
10,13
,
min
⁡
⁡
F
=
0.0945
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
10
subject
to
satisfying
:
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3,4
,
5
,
X
e
(
0
)
≤
X
e
≤
50
X
e
(
0
)
,
e
=
6,7
,
8,9
,
10,13
.
The optimal objective values of optimization problems (60) can be easily obtained as follows, respectively:
(61)
F
max
⁡
=
0.0945
(
1.2
X
3
(
0
)
)
0.05
(
1.2
X
4
(
0
)
)
0.533
×
(
0.8
X
5
(
0
)
)
-
0.0822
(
50
X
10
(
0
)
)
=
1731
,
F
min
⁡
=
0.0945
(
0.8
X
3
(
0
)
)
0.05
(
0.8
X
4
(
0
)
)
0.533
×
(
1.2
X
5
(
0
)
)
-
0.0822
(
X
10
(
0
)
)
=
23.0803
.
Then we get an approximate interval [
23.0803,1731
] for objective function F of original problem (58). So we choose the weighting coefficients w
i
(
r
) (r
≥
0) as
(62)
w
i
(
r
)
=
w
-
i
θ
,
i
=
2,3
,
5
,
where w
-
i
∈
[
1
/
1731
,
1
/
23.0803
]
≈
[
0.0006,0.0433
].The following algorithm parameters were firstly used in the implementation of the proposed algorithm:w
i
(
r
)
=
0.0006
(
r
≥
0
), θ
=
1, s
i
(
0
)
=
1.99, and ε
=
1
0
-
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.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: w
i
(
r
)
=
0.0006
(
r
≥
0
), θ
=
1, and ε
=
1
0
-
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.Table 4
Comparisons between the proposed algorithm and the method in [14] for Case study 2.
Proposed method
[14]
Initial auxiliaryvariable s
i
(
0
)
Iterations
Objective value
Initial auxiliaryvariable t
i
(
0
)
Iterations
Objective value
1
5
1577.4227
0b
—
—
1.1
5
1577.4227
0.1
15
1577.422623
1.2
5
1577.4227
0.2
15
1577.422627
1.3
5
1577.4227
0.3
15
1577.422632
1.4
5
1577.4227
0.4
15
1577.422625
1.5
5
1577.4227
0.5
15
1577.422626
1.6
5
1577.4227
0.6
15
1577.422629
1.7
5
1577.4227
0.7
15
1577.422629
1.8
5
1577.4227
0.8
15
1577.422625
1.9
5
1577.4227
0.9
15
1577.422631
1.99
5
1577.4227
0.99
15
1577.422630
bThe method in [14] does not allow t
i
(
0
)
=
0 (i
=
2,3
,
5) because of the auxiliary variables t
i
>
0 being imposed in the implementation of this algorithm.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: w
i
(
r
)
=
0.0006
(
r
≥
0
), s
i
(
0
)
=
θ
+
0.9, and ε
=
1
0
-
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.Table 5
Comparisons of effects of 10 differentθ values on the proposed algorithm.
θ
Iterations
Objective value
0.02
15
1577.4227
0.1
8
1577.4227
0.2
6
1577.4227
0.6
6
1577.4227
0.8
5
1577.4227
1.0
5
1577.4227
1.4
5
1577.4227
2
5
1577.4227
3
5
1577.4227
4
5
1577.4227Next we investigate the influence of weighting coefficientsw
i
(
r
)
(
r
≥
0
) on the performance of our proposed algorithm. Ten computational experiments with different weighting coefficients w
i
(
r
) (w
i
(
r
)
∈
[
0.0006,0.0433
] and r
≥
0) were performed with the following algorithm parameters: θ
=
1, s
i
(
0
)
=
1.99, and ε
=
1
0
-
6. Table 6 compares the influence of different weighting coefficients w
i
(
r
)
(
r
≥
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 w
i
(
r
)
(
r
≥
0
) is chosen as a smaller value. This suggests that we should choose a small w
i
(
r
)
(
r
≥
0
) value in the implementation of the proposed algorithm.Table 6
Comparisons of effects of 10 different weighting coefficients on the proposed algorithm.
Weight coefficient (w
i
(
r
)
,
r
≥
0)
Iterations
Objective value
0.0006
5
1577.4227
0.001
5
1577.4227
0.005
6
1577.4227
0.01
5
1577.4227
0.015
10
1577.4227
0.025
7
1577.4227
0.03
7
1577.4227
0.035
8
1577.4227
0.04
17
1577.4227
0.0433
10
1577.4227


## 4.1. Case Study 1: Tryptophan Biosynthesis inEscherichia coli

We first apply the proposed method to tryptophan biosynthesis inEscherichia coli. A complete description of this metabolic pathway can be found in [33]. The differential equations in dimensionless variables are given as
(50)
d
X
1
d
t
=
V
11
+
-
V
11
-
,
d
X
2
d
t
=
V
21
+
-
V
21
-
,
d
X
3
d
t
=
V
31
+
-
(
V
31
-
+
V
32
-
+
V
33
-
)
.
Here, X
1 is used for mRNA concentration, X
2 is used for enzyme concentration, and X
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]:
(51)
V
11
+
=
0.6403
X
3
-
5.87
×
1
0
-
4
X
5
-
0.8332
,
V
11
-
=
1.0233
X
1
X
4
0.0035
X
11
0.9965
,
V
21
+
=
X
1
,
V
21
-
=
1.4854
X
2
X
4
0.1349
X
12
0.8651
,
V
31
+
=
0.5534
X
2
X
3
-
0.5573
X
6
0.5573
,
V
31
-
=
X
3
X
4
,
V
32
-
=
0.9942
X
3
7.0426
×
1
0
-
4
X
7
,
V
33
-
=
0.8925
X
3
3.5
×
1
0
-
6
X
4
0.9760
X
8
X
9
-
0.0240
X
10
-
3.5
×
1
0
-
6
.Table 1
Optimal solution of Case study 1 obtained by using the proposed approach (s
3
(
0
)
=
1.99).
Variables
Initial steady state
Optimized solution
X
1
0.184654
0.221585
X
2
7.986756
8.908299
X
3
1418.931944
1135.145557
X
4
0.00312
0.005361
X
5
5
4.010578
X
6
2283
5000
X
8
430
1000
F
1.310202
5.169027The reactionV
33
- can be selected as the optimized objective function [5]. Then we have the following steady-state optimization problem:
(52)
max
⁡
⁡
F
=
V
33
-
,subject to satisfying(53)
V
11
+
-
V
11
-
=
0
,
V
21
+
-
V
21
-
=
0
,
V
31
+
-
(
V
31
-
+
V
32
-
+
V
33
-
)
=
0
,
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3
,
0
<
X
4
≤
0.00624
,
(54)
4
≤
X
5
≤
10
,
500
≤
X
6
≤
5000
,
X
7
=
0.0022
X
5
,
0
<
X
8
≤
1000
(
X
9
,
X
10
,
X
11
,
X
12
,
X
13
)
=
(
7.5,0.005,0.9,0.02,0
)
.The following algorithm parameters were chosen in the implementation of the proposed method:w
3
(
0
)
=
1, w
3
(
r
)
=
1
+
r
(
r
≥
1
), θ
=
1, s
3
(
0
)
=
1.95, and ε
=
1
0
-
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: w
3
(
0
)
=
1, w
3
(
r
)
=
1
+
r
(
r
≥
1
), θ
=
1, and ε
=
1
0
-
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 s
3
(
0
) increases from 1.0 to 1.99, the iterations required for the proposed algorithm decrease. This suggests that we should choose a larger s
i
(
0
) value to promote the evolution rate of the proposed approach.Table 2
Comparisons between our proposed algorithm and the method in [14] for Case study 1.
Proposed method
[14]
Initial auxiliaryvariable s
3
(
0
)
Iterations
Objective value
Initial auxiliaryvariable t
3
(
0
)
Iterations
Objective value
1
9
5.169027
0a
—
—
1.1
9
5.169027
0.1
14
5.169027
1.2
5
5.169027
0.2
14
5.169027
1.3
5
5.169027
0.3
14
5.169027
1.4
5
5.169027
0.4
14
5.169027
1.5
5
5.169027
0.5
14
5.169027
1.6
5
5.169027
0.6
14
5.169027
1.7
5
5.169027
0.7
14
5.169027
1.8
5
5.169027
0.8
14
5.169027
1.9
5
5.169027
0.9
14
5.169027
1.99
5
5.169027
0.99
14
5.169027
aThe method in [14] does not allow t
3
(
0
)
=
0 because of the auxiliary variable t
3
>
0 being imposed in the implementation of this algorithm.


## 4.2. Case Study 2: Maximization of Ethanol Production inSaccharomyces cerevisiae

In this case study, we will apply the proposed approach to the production of ethanol bySaccharomyces cerevisiae. The dynamics of this system are described as the following formulations [34]:
(55)
d
X
1
d
t
=
V
in
-
V
HK
,
d
X
2
d
t
=
V
HK
-
V
PFK
-
V
Pol
,
d
X
3
d
t
=
V
PFK
-
V
GAPD
-
0.5
V
Gol
,
d
X
4
d
t
=
2
V
GAPD
-
V
PK
,
d
X
5
d
t
=
2
V
GAPD
+
V
PK
-
V
HK
-
V
Pol
-
V
PFK
-
V
ATPase
,
where X
i represent the following intermediate metabolite concentrations: X
1 is the intracellular glucose concentration, X
2 represents glucose 6-phosphate, X
3 codes for fructose 1, 6-diphosphate, X
4 is phosphoenolpyruvate, and X
5 represents ATP. The indexed quantities V represent the following fluxes: V
in denotes the sugar transport into the cells, V
HK summarizes all hexokinases, V
PFK is the phosphofructokinase reaction, V
GAPD represents glyceraldehyde 3-phosphate dehydrogenase, V
PK represents pyruvate kinase, V
Pol describes glycogen synthetase, V
Gol, the glycerol 3-phosphate dehydrogenase is proportional to V
PK, and V
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]:
(56)
V
in
=
0.8122
X
2
-
0.2344
X
6
,
V
HK
=
2.86
32
X
1
0.7464
X
5
0.0243
X
7
,
V
PFK
=
0.52
32
X
2
0.7318
X
5
-
0.3941
X
8
,
V
Pol
=
0.0009
X
2
8.6107
X
11
,
V
GAPD
=
0.011
X
3
0.6159
X
5
0.1308
X
9
X
14
-
0.6088
,
V
Gol
=
0.04725
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
12
,
V
PK
=
0.0945
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
10
,
V
ATPase
=
X
5
X
13
.Table 3
Optimal solution of Case study 2 obtained by using the proposed approach.
Variables
Initial steady state
Optimized solution
X
1
0.0345
0.038025
X
2
1.0111
1.057148
X
3
9.1437
10.393273
X
4
0.0095
0.010454
X
5
1.1278
1.255394
X
6
19.7
985
X
7
68.5
3147.937162
X
8
31.7
1585
X
9
49.9
2389.028728
X
10
3440
172000
X
13
25.1
1255
F
30.1124
1577.4227The performance index describing the rate of ethanol production is given directly by the flux through the pyruvate kinase,V
PK. Then we have the following optimization problem [3, 36]:
(57)
max
⁡
⁡
F
=
V
PK
,subject to satisfying(58)
V
in
-
V
HK
=
0
,
V
HK
-
V
PFK
-
V
Pol
=
0
,
V
PFK
-
V
GAPD
-
0.5
V
Gol
=
0
,
2
V
GAPD
-
V
PK
=
0
,
2
V
GAPD
+
V
PK
-
V
HK
-
V
Pol
-
V
PFK
-
V
ATPase
=
0
,
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3,4
,
5
,
X
e
(
0
)
≤
X
e
≤
50
X
e
(
0
)
,
e
=
6,7
,
8,9
,
10,1
3
,
V
PK
≤
2
V
in
,
(
X
11
,
X
12
,
X
14
)
=
(
14.31,203,0.042
)
.
This problem has three signomial equality constrains that can be rewritten as follows:
(59)
V
PFK
+
V
Pol
V
HK
=
1
,
V
GAPD
+
0.5
V
Gol
V
PFK
=
1
,
V
HK
+
V
Pol
+
V
PFK
+
V
ATPase
2
V
GAPD
+
V
PK
=
1
.
Thus, three weighting coefficients w
i
(
r
) (r
≥
0 and i
=
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
(60)
max
⁡
⁡
F
=
0.0945
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
10
subject
to
satisfying
:
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3,4
,
5
,
X
e
(
0
)
≤
X
e
≤
50
X
e
(
0
)
,
e
=
6,7
,
8,9
,
10,13
,
min
⁡
⁡
F
=
0.0945
X
3
0.05
X
4
0.533
X
5
-
0.0822
X
10
subject
to
satisfying
:
0.8
X
e
(
0
)
≤
X
e
≤
1.2
X
e
(
0
)
,
e
=
1,2
,
3,4
,
5
,
X
e
(
0
)
≤
X
e
≤
50
X
e
(
0
)
,
e
=
6,7
,
8,9
,
10,13
.
The optimal objective values of optimization problems (60) can be easily obtained as follows, respectively:
(61)
F
max
⁡
=
0.0945
(
1.2
X
3
(
0
)
)
0.05
(
1.2
X
4
(
0
)
)
0.533
×
(
0.8
X
5
(
0
)
)
-
0.0822
(
50
X
10
(
0
)
)
=
1731
,
F
min
⁡
=
0.0945
(
0.8
X
3
(
0
)
)
0.05
(
0.8
X
4
(
0
)
)
0.533
×
(
1.2
X
5
(
0
)
)
-
0.0822
(
X
10
(
0
)
)
=
23.0803
.
Then we get an approximate interval [
23.0803,1731
] for objective function F of original problem (58). So we choose the weighting coefficients w
i
(
r
) (r
≥
0) as
(62)
w
i
(
r
)
=
w
-
i
θ
,
i
=
2,3
,
5
,
where w
-
i
∈
[
1
/
1731
,
1
/
23.0803
]
≈
[
0.0006,0.0433
].The following algorithm parameters were firstly used in the implementation of the proposed algorithm:w
i
(
r
)
=
0.0006
(
r
≥
0
), θ
=
1, s
i
(
0
)
=
1.99, and ε
=
1
0
-
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.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: w
i
(
r
)
=
0.0006
(
r
≥
0
), θ
=
1, and ε
=
1
0
-
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.Table 4
Comparisons between the proposed algorithm and the method in [14] for Case study 2.
Proposed method
[14]
Initial auxiliaryvariable s
i
(
0
)
Iterations
Objective value
Initial auxiliaryvariable t
i
(
0
)
Iterations
Objective value
1
5
1577.4227
0b
—
—
1.1
5
1577.4227
0.1
15
1577.422623
1.2
5
1577.4227
0.2
15
1577.422627
1.3
5
1577.4227
0.3
15
1577.422632
1.4
5
1577.4227
0.4
15
1577.422625
1.5
5
1577.4227
0.5
15
1577.422626
1.6
5
1577.4227
0.6
15
1577.422629
1.7
5
1577.4227
0.7
15
1577.422629
1.8
5
1577.4227
0.8
15
1577.422625
1.9
5
1577.4227
0.9
15
1577.422631
1.99
5
1577.4227
0.99
15
1577.422630
bThe method in [14] does not allow t
i
(
0
)
=
0 (i
=
2,3
,
5) because of the auxiliary variables t
i
>
0 being imposed in the implementation of this algorithm.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: w
i
(
r
)
=
0.0006
(
r
≥
0
), s
i
(
0
)
=
θ
+
0.9, and ε
=
1
0
-
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.Table 5
Comparisons of effects of 10 differentθ values on the proposed algorithm.
θ
Iterations
Objective value
0.02
15
1577.4227
0.1
8
1577.4227
0.2
6
1577.4227
0.6
6
1577.4227
0.8
5
1577.4227
1.0
5
1577.4227
1.4
5
1577.4227
2
5
1577.4227
3
5
1577.4227
4
5
1577.4227Next we investigate the influence of weighting coefficientsw
i
(
r
)
(
r
≥
0
) on the performance of our proposed algorithm. Ten computational experiments with different weighting coefficients w
i
(
r
) (w
i
(
r
)
∈
[
0.0006,0.0433
] and r
≥
0) were performed with the following algorithm parameters: θ
=
1, s
i
(
0
)
=
1.99, and ε
=
1
0
-
6. Table 6 compares the influence of different weighting coefficients w
i
(
r
)
(
r
≥
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 w
i
(
r
)
(
r
≥
0
) is chosen as a smaller value. This suggests that we should choose a small w
i
(
r
)
(
r
≥
0
) value in the implementation of the proposed algorithm.Table 6
Comparisons of effects of 10 different weighting coefficients on the proposed algorithm.
Weight coefficient (w
i
(
r
)
,
r
≥
0)
Iterations
Objective value
0.0006
5
1577.4227
0.001
5
1577.4227
0.005
6
1577.4227
0.01
5
1577.4227
0.015
10
1577.4227
0.025
7
1577.4227
0.03
7
1577.4227
0.035
8
1577.4227
0.04
17
1577.4227
0.0433
10
1577.4227


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



---

*Source: 719496-2014-06-01.xml*