A Global Optimization Approach for Metabolic Flux

Quantification
Carlos A. M. Riascos1, Andreas K. Gombert1 and Jose M. Pinto1,2*
Department of Chemical Engineering, University of Sao Paulo, So Paulo (SP), Brazil,
05508-900
2
Department of Chemical and Biological Sciences and Engineering, Polytechnic
University, Brooklyn (NY), USA, 11201

Abstract
Metabolic flux analysis (MFA) is an important tool in metabolic engineering that allows
characterizing the organisms in terms of in vivo fluxes. Non-linear MFA quantifies
fluxes by minimizing the error between measurements and model predictions. It
generates a non-convex non-linear programming (NLP) model with bilinear constraints
that result from metabolite and isotope balances. In this work, we propose a global
optimization approach that relies on a spatial branch and bound search algorithm to
solve the NLP. An example that estimates fluxes in the central metabolism of S.
cerevisiae is developed. The results from evolutionary methods and local NLP solvers
generate local optima with central pathways that are different from the global solution,
which predicts values that are closer to the measured ones.
Keywords: metabolic flux analysis, global optimization, spatial branch and bound

1. Introduction
MFA involves the quantification of biochemical reaction rates in a metabolic network
and requires the combination of experimental data, mathematical modeling and a flux
quantification procedure. Measurement includes the concentration of biomass and of
extra-cellular metabolites, the determination of the biomass composition, and the 13Clabeling pattern in intracellular compounds. The latter type of measurement allows flux
estimation in more complex metabolic networks. Quantification estimates the fluxes
that minimize the error between measurements and model predictions.
For error minimization, several optimization methods have been applied, such as the
simplex method, a simulated annealing algorithm and a Newton-like algorithm
(Wiechert, 2001); an evolutionary strategy (Gombert et al., 2001); a genetic algorithm
(Zhao and Shimizu, 2003); and sequential quadratic programming (Wiechert et al.,
2001). All the optimization-based methods rely on non-linear programming (NLP)
techniques, but do not guarantee that the global minimum is found, due to the nonconvexity of the feasible region. Recently, there has been a significant effort to develop
deterministic approaches for global optimization (Quesada and Grossmann, 1995;
Esposito and Floudas, 1998; Smith and Pantelides, 1999).
*

Author to whom correspondence should be addressed: [email protected]

The objective of this work is to develop a global optimization technique that relies on a
spatial branch and bound search algorithm. Fluxes in the central metabolism of S.
cerevisiae have been efficiently estimated, based on a stoichiometric model and GC-MS
experimental data previously reported (Gombert et al., 2001).

2. NLP Model for Metabolic Flux Quantification Problem

The mathematical modeling of a MFA problem, using fractional labeling balances,
generates a non-linear optimization problem (NLP). The constraints of the NLP are
derived from the metabolite balances, the labeling balances, the definition of the
measurable labelings, and the variable bounds. The modeling is detailed in Riascos et al.
(2004) and relies on the following assumptions:
Bi-directional fluxes are divided into net and exchange fluxes (Gombert et al., 2001).
Pseudo-steady state (no accumulation) for intracellular metabolite balances.
Fractional labeling balances for each carbon atom in the metabolic network.
The summed fractional labeling (SFL) for a metabolite fragment measured by GC-MS
is a linear function of the fractional labelings in the intracellular metabolites.
The NLP objective is to minimize the weighted squared error between the
measurements and the model predictions. The measured variables are a set of the
metabolic fluxes, determined from extra-cellular fluxes and biomass composition, and
the summed fractional labelings. Thus the following NLP is generated:
2

min

s.t.

ri ri

iI i

SFL
SFL
m
m

m
mM

r = 0

(2)

g(ri lj,k) = 0
SFLm =

j J m k K j ,m

(1)

l j ,k

j, k Kj

(3)

mM

(4)

where G is the stoichiometric matrix, r is the vector of metabolic fluxes ri

(Stephanopoulos et al., 1998), lj,k is the fractional labeling in carbon atom k of
metabolite j, Kj is the set of carbon atoms in metabolite j, M is the set of measured
fragments, Jm is the set of intracellular metabolites that generate the fragment m, Kj,m is
the set of carbon atoms from metabolite j which are incorporated in fragment m,

r,
i SFLm are the fitted values for variables; ri , SFLm are the measured values; i, m are

their standard deviations; and I, M are the sets of measured fluxes and SFLs.

3. Global Optimum Search

In order to solve the metabolic flux quantification problem to global optimality, we
have adapted a spatial branch and bound algorithm (Esposito and Floudas, 1998). The
global search is developed by the generation of lower (LB) and upper bounds (UB) to
the global solution. The convergence to the global optimum is achieved through
successive division of the search region (branching) that allows the generation of tighter
bounds for the solution (bounding). The lower bound is obtained from the approximated

convex solutions (lb) in the search sub-regions, and the upper bound is the smallest of
the non-convex solutions (ub) for the original problem in the search sub-regions.
To develop an approximate convex model, only the constraints from labeling balances
(Eq. 3) must be linearized. The implementation of the global optimization strategy was
developed in the GAMS environment (Brooke et al., 2000), employing CONOPT2
(Drud, 1996) and is presented in Fig. 1.
3.1 Linear approximation of labeling balances
The linear approximation of fractional labeling balances is based on the method
originally developed by McCormick (1976). Each bilinear term, ri lj,k, is replaced by an
auxiliary variable, rli,j,k, and the new variables are constrained by:
(5)

ri L l j ,k + l Lj , k ri ri L l Lj , k rli , j , k 0
ri U l j ,k + l Uj ,k ri ri U l Uj , k rli , j , k 0
ri U l j ,k l Lj ,k ri + ri U l Lj , k + rli , j , k 0
ri L l j ,k l Uj ,k ri + ri L l Uj , k + rli , j ,k 0

where superscripts L and U denote the lower and upper bounds of variables.
3.2 Measured and non-measured variable bounds
Due to the structure of the McCormick estimators, bounds for measured and nonmeasured variables strongly affect convergence to the global optimum. The measured

Set bounds on variables and parameters

Compute UB1
Tighten initial parameter bounds
Generate LB1
Yes
Stop
G.O. = UB

(1-) UB LB

Compute lb

(1-) UB > lb
Yes
Tighten the bounds

No
Branch the search region
Evaluation of the two generated subregions (detailed on the RHS)
Select region to branch
LB = min lb

Figure 1: Flowsheet of the global optimization procedure.

Compute new lb
Compute new ub
if ub < UB UB = ub

No

variable bounds are obtained from their measured values ( x ) and a parameter (),
which reflects the uncertainty in the measurements, as follows:
xL = (1 - ) x

xU = (1 + ) x

(6)

For non-measured variables, an initial set of bounds is obtained from network

observation and prior knowledge. In the case study the ranges are: 0 rexch 1500 for
exchange fluxes, 0 rnet 150 for net fluxes, and 1.1 lj,k 100 for labelings.
For global optimization convergence, a procedure for tightening parameter bounds is
necessary; it consists of a set of min/max problems that compute the parameter bounds
as the limiting values that satisfy the convex constraints.
3.3 Branch and bound procedure
In the branching steps, the region containing the lower bound (the smallest of the
convex solutions) is divided by bisecting along a selected variable. The branching
procedure must continue while there are still open nodes; a region is considered an open
node if its convex solution (lb) is smaller than the current stopping value (lb < (1- )
UB), where is the relative tolerance for the global convergence.
There are several alternatives to select the branching variable (Esposito and Floudas,
1998; Smith and Pantelides, 1999; Quesada and Grossmann, 1995). In this work, the
selected branching variable is the one that presents the largest difference between its
values in the convex and non-convex problem. To obtain faster convergence, a bound
tightening step for the new sub-regions into each branching was added, as suggested by
Esposito and Floudas (1998).

4. Case Study
We have applied the described methodology to quantify fluxes in the central carbon
metabolism of aerobically growing Saccharomyces cerevisiae cells, which were
harvested from a chemostat at steady state, with a specific growth rate of 0.1 h-1
(Gombert et al., 2001). The network of metabolic fluxes is shown in Fig. 2.
Initial variable bounds employed an uncertainty parameter ( ) value of 0.2. Table 1
shows the local search results from randomly generated initial guesses and the best
feasible result obtained by an evolutionary strategy, thus confirming the existence of
local minima. Moreover, in Fig. 2 the fluxes in PP and EMP pathways show different
values in the locally and globally minimum solutions in the metabolic network.

Table 1: Results from the local NLP and the evolutionary algorithm.
Method
NLP
Evolutionary method
(Gombert et al., 2001)

Obj. function Summed quadratic errors

Fluxes
Labelings
22.74
0.86
21.88
201.85
1.02
200.83
26.31

1.70

24.61

Central metabolic fluxes

EMP
PP
TCA
34.37
43.74
62.14
47.17
29.85
64.23
33.64

44.16

59.89

The global convergence considers a 10% convergence tolerance ( = 0.1). The largest
reduction in the feasibility gap was observed by branching on non-measured net fluxes
and summed fractional labelings (SFL). The selected branching variable is the one that
presents the largest difference between its values in the convex and non-convex
problem. Global convergence required several branching operations that included the
bound tightening procedure for each new region (see Fig. 1). Overall, 132 branching
steps and 10.4 CPU hours (in a 2 GHz Pentium IV PC) are required to achieve
convergence. Note that the large CPU time corresponds to the compilation and solution
of all problems at each node of the search tree.
The obtained global solution corresponds to the first local solution in Table 1. The
estimated flux distributions by global optimization and by an evolutionary algorithm
(Gombert et al., 2001) are presented in Fig. 2. The central metabolic fluxes are similar
for both approaches, but the quadratic error on the variables (objective function) is
smaller in the global approach. Moreover, the obtained fluxes for the transport of
oxaloacetate across the mitochondrial membrane, as well as for the Acetyl-CoA
generation ways, are not the same in the global optimization and evolutionary solutions.
Glucose

Proposed Global solution

Evolutionary solution
Measurement (when available)

100

PP pathway

GLY

60.8
60.3

2.9 (0.0)
2.9 (0.1)

4.7

SER

4.7

128.1
127.1

OAA

30.1
27.3

123.2
122.1

21.2
18.3

41.5
35.8

0.0
0.1
0.0
0.0
0.1
0.0

4.0
4.0
4.07

PYR
21.6
45.1

67.5
45.7

0.0
0.1
0.0

5.0
5.0
4.99

PEP

0.8 (0.1)
0.8 (0.1)

THR

0.8
0.8
0.82

G3P

ACA

4.4
4.4

22.3
22.9
22.0

Biomass

F6P

1.8
1.8
1.81

3.7
3.7
3.66

5.3
5.4
5.3

34.4 (1420)
33.6 (1000)

12.0
12.1

Biomass

21.9
22.2
22.2

G6P

26.5
26.7

5.3
5.4
5.36

3.6
3.6
3.63

43.7
44.2

ACA

65.0
39.9
62.1
59.9

OAA
MAL

FUM

22.3
45.7

AKG

AcCoA

0.0
22.8

2.8
2.8
2.8

ICI
62.1
59.9

51.6
49.1

Acetate

AcCoA

TCA cycle

51.6 (262)
49.1(237)

ACE

Ethanol

12.3
12.5
12.4

PYR

9.8
6.7

22.3
45.8

Glycerol

10.5
10.8
9.62

Biomass

mitochondrion

Figure 2: Results for the metabolic network (reversible reactions - flux in parenthesis).

It is important to note that tighter initial bounds and consequently a smaller feasibility
gap, which are obtained employing a smaller uncertainty parameter value (), might
reduce the number of branching steps and the computational effort, but the current
solution might be excluded from the search region. The largest reduction on the gap
occurs when bounds are tightened; consequently, in spite of the fact that this step
consumes 58% of the computational effort (CPU time), it is critical to achieve global
convergence. A more extensive analysis is developed in Riascos et al. (2004).

5. Conclusions
When measurements by GC-MS and fractional labeling balances are employed, the
resulting bilinear constraints add a special difficulty in MFA problems. This difficulty is
due to the fact that the resulting optimization model not only becomes non-convex but
also presents multiple local minima, each of which corresponds to a different
distribution in the metabolic network.
A spatial branch and bound method was satisfactorily applied to non-linear MFA. The
summed fractional labelings (SFLs) and the non-measured net fluxes have been
identified as the most suitable sets of branching variables. The greater reduction on the
feasibility gap is due to the bound tightening procedure, which is indispensable to
achieve global convergence. The results have shown that global optimization is
effective for flux quantification, despite the higher computational effort.
Acknowledgments
The authors would like to acknowledge financial support from CAPES and
PADCT/CNPq under grant 62.0239/97 QEQ.
References
Brooke, A., D. Kendrick, A.A. Meeraus and R. Raman, 2000, GAMS A User's Guide
(Release 2.50). The Scientific Press, Redwood City.
Drud, A., 1996, CONOPT: A system for large scale nonlinear optimization. Reference
manual. ARKI Consulting and Development A/S, Bagsvaerd.
Esposito, W.R. and C.A. Floudas, 1998, Comp. Chem. Engng. 22, S213-S220.
Gombert, A.K., M.M. dos Santos, B. Christensen and J. Nielsen, 2001, J. Bacteriology
183 (4), 1441-1451.
McCormick, G.P., 1976, Mathematical Programming 10, 147-175.
Quesada, I. and I.E. Grossmann, 1995, Comp. Chem. Engng. 19 (12), 1219-1242.
Riascos, C.A.M., A.K. Gombert and J.M. Pinto, 2004, Comp. Chem. Engng. (in press).
Smith, E.M.B. and C.C. Pantelides, 1999, Comp. Chem. Engng. 23, 457-478.
Stephanopoulos, G.N., A.A. Aristidou and J. Nielsen, 1998, Metabolic Engineering
Principles and Methodologies. Academic Press, San Diego.
Wiechert, W., 2001, Met. Engng. 3, 195-206.
Wiechert, W., M. Mllney, S. Petersen and A.A. de Graaf, 2001, Met. Engng. 3, 265-283.
Zhao, J. and K. Shimizu, 2003, J. Biotechnology 101, 101-117.