Journal of Applied Mathematics

Volume 2014 (2014), Article ID 719496, 10 pages

http://dx.doi.org/10.1155/2014/719496

Research Article

An Improved Geometric Programming Approach for Optimization of Biochemical Systems

Gongxian Xu and Lei Wang

Department of Mathematics, Bohai University, Jinzhou 121013, China

Received 31 January 2014; Accepted 18 May 2014; Published 1 June 2014

Academic Editor: Chein-Shan Liu

Copyright © 2014 Gongxian Xu and Lei Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [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. 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]:

Table 1: Optimal solution of Case study 1 obtained by using the proposed approach ().

Variables
Initial steady state
Optimized solution

0.184654
0.221585

7.986756
8.908299

1418.931944
1135.145557

0.00312
0.005361

5
4.010578

2283
5000

430
1000

1.310202
5.169027

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.

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

Proposed method
[14]

Initial auxiliary 

variable
Iterations
Objective value
Initial auxiliary 

variable
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

The method in [14] does not allow because of the auxiliary variable being imposed in the implementation of this algorithm.

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]:

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

Variables
Initial steady state
Optimized solution

0.0345
0.038025

1.0111
1.057148

9.1437
10.393273

0.0095
0.010454

1.1278
1.255394

19.7
985

68.5
3147.937162

31.7
1585

49.9
2389.028728

3440
172000

25.1
1255

30.1124
1577.4227

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.

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

Proposed method
[14]

Initial auxiliary 

variable
Iterations
Objective value
Initial auxiliary 

variable
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

The method in [14] does not allow () because of the auxiliary variables 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: , , 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.

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

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.

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

Weight coefficient ()
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.

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

References

E. O. Voit, “Optimization in integrated biochemical systems,” Biotechnology and Bioengineering, vol. 40, no. 5, pp. 572–582, 1992.

N. V. Torres, E. O. Voit, and C. González-Alcón, “Optimization of nonlinear biotechnological processes with linear programming: application to citric acid production by Aspergillus niger,” Biotechnology and Bioengineering, vol. 49, no. 3, pp. 247–258, 1996.

N. V. Torres, E. O. Voit, C. González-Alcón, and F. Rodríguez, “An indirect optimization method for biochemical systems: description of method and application to the maximization of the rate of ethanol, glycerol and carbohydrate production in Saccharomyces cerevisiae,” Biotechnology and Bioengineering, vol. 55, no. 5, pp. 758–772, 1997.

V. Hatzimanikatis, C. A. Floudas, and J. E. Bailey, “Analysis and design of metabolic reaction networks via mixed-integer linear optimization,” AIChE Journal, vol. 42, no. 5, pp. 1277–1292, 1996.

A. Marín-Sanguino and N. V. Torres, “Optimization of tryptophan production in bacteria. Design of a strategy for genetic manipulation of the tryptophan operon for tryptophan flux maximization,” Biotechnology Progress, vol. 16, no. 2, pp. 133–145, 2000.

A. Marín-Sanguino, E. O. Voit, K. C. Gonzalez-Alcon, and N. V. Torres, “Optimization of biotechnological systems through geometric programming,” Theoretical Biology and Medical Modelling, vol. 4, article 38, 16 pages, 2007.

Y.-J. Chang and N. V. Sahinidis, “Optimization of metabolic pathways under stability considerations,” Computers & Chemical Engineering, vol. 29, no. 3, pp. 467–479, 2005.

P. K. Polisetty, E. P. Gatzke, and E. O. Voit, “Yield optimization of regulated metabolic systems using deterministic branch-and-reduce methods,” Biotechnology and Bioengineering, vol. 99, no. 5, pp. 1154–1169, 2008.

C. Pozo, G. Guillén-Gosálbez, A. Sorribas, and L. Jiménez, “Outer approximation-based algorithm for biotechnology studies in systems biology,” Computers & Chemical Engineering, vol. 34, no. 10, pp. 1719–1730, 2010.

C. Pozo, G. Guillén-Gosálbez, A. Sorribas, and L. Jiménez, “A spatial branch-and-bound framework for the global optimization of kinetic models of metabolic networks,” Industrial and Engineering Chemistry Research, vol. 50, no. 9, pp. 5225–5238, 2011.

A. Sorribas, C. Pozo, E. Vilaprinyo, G. Guillén-Gosálbez, L. Jiménez, and R. Alves, “Optimization and evolution in metabolic pathways: global optimization techniques in Generalized Mass Action models,” Journal of Biotechnology, vol. 149, no. 3, pp. 141–153, 2010.

G. Xu, C. Shao, and Z. Xiu, “A modified iterative IOM approach for optimization of biochemical systems,” Computers & Chemical Engineering, vol. 32, no. 7, pp. 1546–1568, 2008.

G. Xu, “Bi-objective optimization of biochemical systems by linear programming,” Applied Mathematics and Computation, vol. 218, no. 14, pp. 7562–7572, 2012.

G. Xu, “Steady-state optimization of biochemical systems through geometric programming,” European Journal of Operational Research, vol. 225, no. 1, pp. 12–20, 2013.

J. Vera, C. González-Alcón, A. Marín-Sanguino, and N. V. Torres, “Optimization of biochemical systems through mathematical programming: methods and applications,” Computers & Operations Research, vol. 37, no. 8, pp. 1427–1438, 2010.

J. O. Sendin, J. Vera, N. V. Torres, and J. Banga, “Model based optimization of biochemical systems using multiple objectives: a comparison of several solution strategies,” Mathematical and Computer Modelling of Dynamical Systems, vol. 12, no. 5, pp. 469–487, 2006.

J. O. Sendin, O. Exler, and J. R. Banga, “Multi-objective mixed integer strategy for the optimisation of biological networks,” IET Systems Biology, vol. 4, no. 3, pp. 236–248, 2010.

A. Kostromins, I. Mozga, and E. Stalidzans, “ConvAn: a convergence analyzing tool for optimization of biochemical networks,” BioSystems, vol. 108, no. 1–3, pp. 73–77, 2012.

P.-K. Liu and F.-S. Wang, “Hybrid differential evolution including geometric mean mutation for optimization of biochemical systems,” Journal of the Taiwan Institute of Chemical Engineers, vol. 41, no. 1, pp. 65–72, 2010.

G. M. de Hijas-Liste, E. Klipp, E. Balsa-Canto, and J. R. Banga, “Global dynamic optimization approach to predict activation in metabolic pathways,” BMC Systems Biology, vol. 8, article 1, 15 pages, 2014.

A. R. Zomorrodi, P. F. Suthers, S. Ranganathan, and C. D. Maranas, “Mathematical optimization applications in metabolic networks,” Metabolic Engineering, vol. 14, no. 6, pp. 672–686, 2012.

M. A. Savageau, “Biochemical systems analysis: I. Some mathematical properties of the rate law for the component enzymatic reactions,” Journal of Theoretical Biology, vol. 25, no. 3, pp. 365–369, 1969.

M. A. Savageau, “Biochemical systems analysis. II. The steady-state solutions for an n-pool system using a power-law approximation,” Journal of Theoretical Biology, vol. 25, no. 3, pp. 370–379, 1969.

M. A. Savageau, “Biochemical systems analysis: III. Dynamic solutions using a power-law approximation,” Journal of Theoretical Biology, vol. 26, no. 2, pp. 215–226, 1970.

M. A. Savageau, Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology, Addison-Wesley, Reading, Mass, USA, 1976.

M. A. Savageau, E. O. Voit, and D. H. Irvine, “Biochemical systems theory and metabolic control theory: 1. Fundamental similarities and differences,” Mathematical Biosciences, vol. 86, no. 2, pp. 127–145, 1987.

M. A. Savageau, E. O. Voit, and D. H. Irvine, “Biochemical systems theory and metabolic control theory: 2. The role of summation and connectivity relationships,” Mathematical Biosciences, vol. 86, no. 2, pp. 147–169, 1987.

E. O. Voit, Computational Analysis of Biochemical Systems. A Practical Guide for Biochemists and Molecular Biologists, Cambridge University Press, Cambridge, UK, 2000.

S. Boyd, S.-J. Kim, L. Vandenberghe, and A. Hassibi, “A tutorial on geometric programming,” Optimization and Engineering, vol. 8, no. 1, pp. 67–127, 2007.

M. Chiang, “Geometric programming for communication systems,” Foundations and Trends in Communications and Information Theory, vol. 2, no. 1-2, pp. 1–156, 2005.

A. Mutapcic, K. Koh, S. J. Kim, and S. Boyd, “GGPLAB: a simple Matlab toolbox for geometric programming,” 2006, http://www.stanford.edu/~boyd/ggplab/

MOSEK Modeling Manual, MOSEK ApS, 2013, http://www.mosek.com/resources/doc/

Z. Xiu, A.-P. Zeng, and W.-D. Deckwer, “Model analysis concerning the effects of growth rate and intracellular tryptophan level on the stability and dynamics of tryptophan biosynthesis in bacteria,” Journal of Biotechnology, vol. 58, no. 2, pp. 125–140, 1997.

J. L. Galazzo and J. E. Bailey, “Fermentation pathway kinetics and metabolic flux control in suspended and immobilized Saccharomyces cerevisiae,” Enzyme and Microbial Technology, vol. 12, no. 3, pp. 162–172, 1990.

R. Curto, A. Sorribas, and M. Cascante, “Comparative characterization of the fermentation pathway of Saccharomyces cerevisiae using biochemical systems theory and metabolic control analysis: model definition and nomenclature,” Mathematical Biosciences, vol. 130, no. 1, pp. 25–50, 1995.

N. V. Torres and E. O. Voit, Pathway Analysis and Optimization in Metabolic Engineering, Cambridge University Press, Cambridge, UK, 2002.