- Email: [email protected]
MCA Has More to Say
J. theor. Biol. (1996) 182, 233–242
MCA Has More to Say V H†‡ J E. B†§ † Institute of Biotechnology, ETH-Zurich, CH-8093 Zurich, Switzerland and the ‡ Department of Chemical Engineering, California Institute of Technology, Pasadena, CA-91125, U.S.A.
The recent development of two mathematical frameworks for the analysis and design of metabolic systems is reviewed here. Both frameworks employ a (log)linear kinetic metabolic model that is constructed making direct, explicit use of the individual enzyme kinetic parameters which are a foundation of Metabolic Control Analysis (MCA) information. The first framework allows the description of the dynamic responses of metabolic systems subject to fluctuations in their parameters. The second framework considers the problem of optimizing the regulatory structure of metabolic networks. Both frameworks are powerful in enabling greater insights into metabolism based on MCA quantities. 7 1996 Academic Press Limited
understanding of the determinants of metabolic fluxes (Kacser & Burns, 1973). Since its introduction 30 years ago, MCA has received much theoretical and experimental attention [an excellent survey is found in Fell (1992)]. Many theoretical and experimental methodologies have been developed that allow determination of quantitative indices defined within MCA, such as control coefficients and elasticities. The description of metabolic systems by MCA quantities is commonly used because of limitations in the available information about these systems. Furthermore, when a nonlinear model is available, it can be linearized and well-studied within the same framework. One of the major limitations of the standard formulation of MCA has been its restriction to a system in a steady state. The MCA framework has not been used to describe transient phenomena that arise from fluctuations in metabolic parameters and inputs (such as enzyme expression levels, and the concentrations of external substrates and independent effectors). However, as will be shown here, this limitation has been overcome recently by the development of a (log)linear kinetic model of metabolic systems based explicitly on MCA information (i.e., elasticities and control coefficients).
1. Introduction Metabolism is an extremely complicated, multicomponent, dynamic system. In addition to numerous different chemical reactions and intermolecular interactions, many coupled through common cofactors, metabolism is complicated by an overlay of multiple control systems operating at both genetic and protein levels. Understanding the genetic bases for metabolism, and how metabolism is influenced by environment, is central to biology, medicine, and biotechnology in many respects (Bailey, 1991). Metabolic control analysis (MCA), a sensitivity analysis framework, is one of the most developed methods for the quantitative description of metabolism and microbial physiology that allows the analysis and study of the responses of metabolic systems to changes in their parameters (Kacser & Burns, 1973; Heinrich & Rapoport, 1974; Kell & Westerhoff, 1986; Cornish-Bowden & Cardenas, 1990; Rutgers et al., 1991; Fell, 1992). Kacser and Burns’ formulation and lucid explanation of MCA paved the way for an ongoing dialogue between biology and mathematics, aimed at quantitative § Author to whom correspondence should be addressed. 0022–5193/96/190233 + 10 $25.00/0
233
7 1996 Academic Press Limited
234
. . .
This model can be used to simulate dynamic responses of the system subject to fluctuations in its inputs and parameters. The major advantages of the (log)linear model are the specification of the system dynamics explicitly using the same parameters as employed in MCA, and the linearity (in terms of logarithms) of the model. As discussed in more detail later, after examining two applications of the (log)linear representation, neither of these features are shared by the well-known Biochemical Systems Theory (BST) (Savageau, 1976), which also considers the dynamics of metabolic systems. Moreover, the formulation of this (log)linear model allowed the development of a novel mathematical framework for determining changes in the structure and strength of metabolic regulation which should be considered in order to optimize a particular metabolic process. A mixed-integer linear programming (MILP) formulation is proposed for the general case of optimizing performance of metabolic processes which can be adequately approximated using a (log)linear model. The solution of the MILP formulation provides information on which enzymes should be present at different levels, the extent of such changes needed, and the accompanying modifications in the regulatory structure that will optimize the metabolic process. The results of an application of this method are also presented, illustrating the power of the method and the additional results which can be obtained based on MCA quantities.
2. A (Log)Linear Kinetic Model of Metabolic Reaction Networks Prediction of the changes in metabolic fluxes which will occur after a change in any parameter affecting metabolism requires a kinetic representation relating fluxes to metabolic concentrations and system parameters. In general such models are not available, and the kinetic model must therefore be developed from measurements of fluxes and metabolite concentrations. A variety of formulations are possible. One attractive choice is reviewed next (Hatzimanikatis & Bailey, 1995). Consider a metabolic system consisting of n metabolites and m enzymatically-catalysed reactions. We are interested in studying how modifications of the metabolic parameters and of the properties of the enzymes that catalyze these reactions affect the time response characteristics of metabolic functions of the system, such as metabolic concentrations, fluxes, and the specific growth rate. Assuming that any concentration gradients in the cell are of negligible
significance, mass balances on the metabolites of the system may be written as: dx = f(v(x;p(t)),x;p(t),t) dt
(1)
where x is the n-dimensional metabolite concentration vector, f is a function determined by the mass balances, v is the m-dimensional reaction rate vector, and p is the s-dimensional manipulated parameter vector (e.g., enzyme concentrations, external nutrient or substrate concentrations). In addition to metabolic reaction rates, the mass balance equations also include terms that account for other processes by which concentrations of metabolites change (such as the dilution brought about by increases in the biomass volume (Fredrickson, 1976) and transport through the cell envelope). In addition, we consider the r-dimensional vector of metabolic outputs, h, which may be written in general in the form: h = h(v(x;p(t)),x;p(t),t).
(2)
In eqn (2) h is a function of the rates of interest, of the metabolic concentration, and of the parameters. Let x0 be a steady-state—or one of the steady states in the case of steady-state multiplicity—that corresponds to the parameter vector values p0 , and let x0 be non-zero with positive elements. Linearization around this steady state will result in the linear system: d(x − x0 ) 1f = dt 1x
b
x0 ,p0
(x − x0 ) +
1f 1p
b
x0 ,p0
(p(t) − p0 ) (3)
for the mass balances, and
h = h(v0 (x0 ;p0 ),x0 ;p0 ,0) +
1h 1x
b
+
x0 ,p0
1h 1p
(x − x0 )
b
x0 ,p0
(p(t) − p0 )
(4)
for the metabolic outputs. In order to describe the dependence of the metabolic system explicitly on the reaction rates and on the rest of the metabolic processes, the following final representation is adopted, after the appropriate scaling of the linearized system [eqns (3) and (4)] and transformation of variables and parameters
(Hatzimanikatis & Bailey, 1995; Hatzimanikatis et al., 1995): dz = Nez + kz + NPq + Lq dt
(5)
w = Jez + Hz + JPq + Uq
(6)
zi = ln(xi /xi,0 ) qk = ln(pk /pk,0 )
(7)
wl = ln(hl /hl,0 ) and N, J, k, L, H, U, e, and P, are matrices, defined as:
6
vj,0 1fi xi,0 1vj
0 1 7
6
vj,0 1hl hl,0 1vj
J = jl,j =jl,j =
x0 ,p0
0 1 7
6
xk,0 1fi xi,0 1xk
6
pk,0 1fi xi,0 1pk
6
xi,0 1hl hl,0 1xi
6
pk,0 1hl hl,0 1pk
k = ki,k =ki,k =
L = li,k =li,k =
H = hl,i =hl,i =
U = ul,k =ul,k =
6
e = ej,i =ej,i =
6
,
0 1 7 x0 ,p0
0 1 7 x0 ,p0
,
,
0 1 7
,
x0 ,p0
0 1 7 x0 ,p0
,
0 1 7
xi,0 1vj vj,0 1xi
P = pj,k =pj,k =
,
x0 ,p0
x0 ,p0
,
0 1 7
pk,0 1vj vj,0 1pk
x0 ,p0
(8)
At steady state, the solution of (5) and (6) yields: w = Cq
where C = −(Je + H)(Ne + k)−1(NP + L) + JP + U (10) with
where, z, q, and w, are the logarithmic deviations of the metabolic concentrations, the enzyme levels, and the metabolic outputs, respectively:
N = ni,j =ni,j =
235
(9)
6
C = cl,k \ cl,k =
0 17
pk,0 dhl hl,0 dpk
p0
The (log)linear kinetic metabolic model described above depends on the same information as that employed within the framework of metabolic control analysis (MCA) (Reder, 1988; Cornish-Bowden & Cardenas, 1990; Schlosser & Bailey, 1990). Matrices e and P are the elasticity matrices with respect to metabolites and to parameters, respectively. Matrix N is a generalized scaled stoichiometric matrix for the mass balances, and matrix J is the scaled stoichiometric matrix for the outputs, i.e. the linear relations between the reaction rates and the metabolic outputs. Matrices k, H, L, and U, could be defined as the elasticities with respect to metabolites (k and H) and to parameters (L and U) for metabolic processes that are independent of the reaction rates, v. The matrix C is the control coefficient matrix of the metabolic functions h with respect to parameters p. Equation (10) defines the control coefficients in a very general and systematic fashion independently of the stoichiometry of the metabolic network and the complexity of the metabolic outputs. The coefficient cl,k in the l-th row and k-th column of C represents the fractional change of metabolic function hl , in response to a fractional change in metabolic parameter or input, pk , which is the mathematical definition of the control coefficient of hl with respect to pk . Experimental determination of the parameters for this (log)linear system has been the subject of several studies (Fell, 1992; Cornish-Bowden & Cardenas, 1990; Schlosser & Bailey, 1990; Schlosser et al., 1993) and in many cases such parameters are the only available kinetic information about a metabolic system. We should notice here that the final representation [eqns (5) and (6)] allows the explicit description of the system with quantities that are characteristic of each enzyme. Therefore, we could study the effects of modifications of the catalytic properties of enzyme i with respect to its substrate (or regulatory effector) j, by changing the value of the corresponding elasticity eij . One of the main advantages of the (log)linear model is its ability to describe the time responses of the system to transient changes in the metabolic parameters or to temporal perturbations in the
. . .
236
concentration of the metabolites. The linearity of eqn (5) allows us to solve explicitly for z as a function of time (Seinfeld & Lapidus, 1974): z(t) = e (Ne + k)(t − t0 )z0 +
g
×(NP + L)q(t)dt
(11)
z0 = z(t0 ) For the time dependence of the logarithmic deviations of the metabolic outputs we have from eqns (6) and (11)
0
w(t) = (Je − H) e (Ne + k)(t − t0 )z0 +
g
t
e (Ne + k)(t − t)
t0
1
×(NP + L)q(t)dt + (JP + U)q(t)
(12)
The linearity of the model permits fast numerical calculation of the integral in eqns (11) and (12). For the time-dependent concentration of the metabolites and for the metabolic outputs we simply have to calculate the transformations: xi (t) = x0,i e zi (t)
(13)
hl (t) = h0,l e wl (t)
(14)
and
which can also be written as
k=1
0 1 0 1 xk (t) x0,k
alk
pj (t) p0,j
s
P
j=1
blj
(15)
where
0
n
1
(16)
1
(17)
alk = s jli eik + hlk i=1
and
0
n
blj = s jli pil + ulj i=1
Equations (11–17) account for the time dependence of metabolites and metabolic outputs based explicitly on standard MCA quantities and thus substantially extend the scope of metabolic analysis and design which can be achieved using these quantities. : A simple linear pathway (Fig. 1) was studied as a test of the performance of the (log)linear metabolic 2.1.
2 S2
3 S3
4 S4
P
F. 1. Linear pathway with Michaelis-Menten kinetics used for assessing the (log)linear dynamic model.
e (Ne + k)(t − t)
where
n
1 S1
t
t0
hl = h0,l P
0 S0
model in approximating a given nonlinear metabolic system (which here is the full set of unsteady-state nonlinear material balances on the components in this pathway) (Hatzimanikatis & Bailey, 1995). In the nonlinear system, Michaelis-Menten kinetics were considered for every reaction rate and the parameters used are presented in Appendix A. The (log)linear model parameters were derived as described above. We studied the dynamic responses of the flux through step 4, to the final product P, for step changes [Fig. 2(a)] and for sinusoidal variation [Figs 2(b), (c) and (d)] of the external substrate, s0 . The time response characteristics of the flux through the pathway are in good agreement between the original nonlinear model and the approximate (log)linear model. The observed deviations are due to the approximations inherent in the MCA quantities used to construct the (log)linear model. The final steady-state differences for step processes [Fig. 2(a)] are the differences between the MCA calculations of new steady states and the responses of the original system. The question of the sensitivity of a system either to changes in its parameters or to fluctuations of the concentrations of intracellular metabolites has been addressed and studied using the simple horizontal model of the system; i.e., the system described by eqns (3) and (4) (Kahn & Westerhoff, 1993). However, the (log)linear model is in much better agreement with the original system model than is a simple linear model, especially with respect to response to periodic inputs [Figs 2(b), (c) and (d)]. The excellent performance of the (log)linear model exhibited here under temporal variation of system parameters suggests that MCA information can be used directly in order to analyse the behavior of metabolic systems subject to fluctuations in their parameters. In particular, metabolic systems, like any nonlinear dynamical system, can exhibit a shift in time-averaged properties (such as metabolic fluxes and control coefficients) in response to fluctuations with zero-time-average value (Bailey, 1974, 1977). The (log)linear dynamic model presented here provides a powerful, convenient vehicle for analysis of such phenomena (Hatzimanikatis & Bailey, 1995). Good agreement between response characteristics of the (log)linear model and of the nonlinear model has also been observed for various more complicated and strongly nonlinear metabolic systems [such as a
(a)
1.20
237
(c) II
1.10
1.000
I 1.00 III
0.980
IV
0.960
Flux to product
0.90 0.80 0.70 (b)
(d)
1.000
1.000 0.990 0.995 0.980 0.990 0.970 0.985
0
20
40 Time
60
0
20.0
40.0
60.0
Time
F. 2. Dynamic responses of the flux through the linear pathway. Solid lines correspond to the (log)linear model, dashed lines correspond to the nonlinear model, and dotted-dashed lines correspond to the linear model. Vertical axes correspond to the flux to product P through reaction step 4. A. Responses to step changes of the external substrate, S0 . I: +20%; II: +50%; III: −20%; IV: −50%. B–D. Responses to sinusoidal variation of the external substrate with different amplitudes and frequencies: B: s0 = 2(1 − 0.2 sin(t)). C: s0 = 2(1 − 0.5 sin(t)). D: s0 = 2(1 − 0.5 sin(t/0.1)).
branched pathway and yeast glycolysis (Hatzimanikatis & Bailey, 1995)]. This agreement suggests that, in the absence of a nonlinear model, and when MCA quantities can be experimentally determined or estimated, the (log)linear model can be used as a good first approximation for analysing the dynamic response characteristics of metabolic systems. To summarize, the (log)linear model is an attractive representation of metabolic pathways since it can be easily constructed using MCA information, it has an unsteady-state analytical solution which makes computation easier, and it can describe very well the dynamic response characteristics of metabolic networks. Moreover, the analytical solution of the (log)linear model allows the derivation of indices such as the regulatory strength and the homeostatic strength (Kahn & Westerhoff, 1993) by simple manipulations of eqns (11–15) (Hatzimanikatis & Bailey, in preparation). 3. Optimization of Regulatory Architectures: Analysis and Synthesis Problems Metabolic reaction networks in living organisms have evolved through natural selection and are therefore configured to maximize the probability of survival of the species, at least insofar as evolution has
progressed. While explicit formulation of the objective function for natural metabolism is non-trivial, prior investigators have proposed maximization of growth rates or most efficient utilization of cellular energetic and chemical resources as the objective function for evolution of natural metabolism (Savageau, 1976; Ramkrishna, 1983; Heinrich et al., 1987; Marr, 1991; Schuster & Heinrich, 1991). However, biotechnological objectives require the identification of a different configuration of fluxes which directs raw materials to products efficiently at high rates and in the presence of high concentrations of product (Bailey, 1991). Therefore, a productionoriented evolution is needed to achieve these goals. Optimization approaches to the analysis of biochemical systems for improvement of biotechnological objectives have been followed in various studies (Majewski & Domach, 1990; Regan et al., 1992; Voit, 1992; Varma & Palsson, 1994) for the identification of the optimum manipulation of the external inputs to the system (such as independent effectors and external substrates) and of enzyme expression levels. However, these studies do not address the problem of optimizing the regulatory structure of the metabolic network. Optimization of a metabolic function with respect to changes in the amounts of enzymes in the pathway and in
. . .
238
modifications in regulatory characteristics of those enzymes can be undertaken using the kinetic description provided by the (log)linear metabolic model (Hatzimanikatis et al., 1995). In this optimization context, the model defines constraints in the forms of steady-state mass balances for the intracellular metabolites and of (log)linear expressions of metabolic functions that should be constrained within certain values to force the optimization search to stay in physiologically reasonable territory. When MCA information is available for a metabolic system, the information concerning the regulatory structure is included in the elasticity matrix e, which can be written as a sum of two matrices: e = e s + er
(18)
where the elements in matrix es correspond to the substrate elasticities of the enzymes, that is, the sensitivities of enzyme activities with respect to concentrations of their substrates, and the elements of matrix er correspond to the regulatory elasticities of the enzymes, that is, the sensitivities of enzyme activities with respect to concentrations of regulatory metabolites. In this representation, the substrates themselves can also be considered as regulatory metabolites (in cases such as substrate inhibition):
6
es = ej,is =
0 1
xis ,0 1vj vj,0 1xis
x0 ,p0
=xis is a substrate for reaction j
7
and
6
er = ej,ir =
0 1
xir , 0 1vj vj,0 1xir
x0 ,p0
=xir is a regulator for reaction j
7
Changes of the elements in matrix er from non-zero values to zero or vice versa correspond to modifications in the regulatory structure of the system. The problem of optimizing the performance of a metabolic network by modifying its regulatory structure has been defined in two different contexts, called analysis and synthesis problems (Hatzimanikatis et al., 1995). In the context of a given regulatory structure, which can be modified by
deleting certain control interactions, the analysis problem was defined as follows: Which of the existent regulatory loops should be inactivated, and what associated changes should be made in the manipulated variables (e.g., enzyme expression levels, environmental conditions, effectors external to the system) in order to optimize the performance of the metabolic network? In a more extended context a regulatory superstructure can be postulated, in which a set of alternative regulatory elasticities for each enzyme, corresponding to different kinds of regulation by each metabolite, is embedded. Then, the synthesis problem considers the following equations: What kind of regulation (i.e., activation or inhibition, by which metabolite and of what strength) should be assigned to each enzyme in the network, and what associated changes should be made in the manipulated parameters (e.g., enzyme expression levels, environmental conditions, effectors external to the system), in order to optimize the performance of the metabolic network? Depending on the measure of ‘‘the performance of the metabolic network’’ (e.g., the metabolic function to be optimized), the answer(s) to the above questions can suggest possible evolutionary criteria that gave rise to naturally existing regulatory structures, or, in a biotechnological context, provide potential useful quantitative guidance for promising targets for protein engineering to achieve a preferred flux distribution. Inactivation or activation of a regulatory loop is equivalent to elimination or introduction of non-zero terms in the er matrix in eqn (18). Moreover, the synthesis problem will typically be subject to some constraints such as the possible number of regulatory actions on each enzyme and the requirement that an enzyme cannot be activated and inhibited by the same metabolite. These kind of decisions, present in both the analysis and synthesis problems, are discrete decisions. Biochemical analysis criteria or experimental limitations might also constrain the maximum number of continuously adjustable manipulated parameters, imposing another set of discrete constraints. The mixed discrete and continuous nature of the problem and the linear (in terms of logarithms) description of the system led to the formulation of the analysis and synthesis problems as MILP problems (Hatzimanikatis et al., 1995). A multitude of questions can be addressed within this framework by
239
PEP G6P + PEP
1
2
DAHP ATP
CHR
3
PHP
4
PHE
4ATP
ADP G6P 6
5
4ADP TRP
PHE
9
Biomass
7
TYR
TRP
8 TYR F. 3. The aromatic amino acid synthesis pathway. Solid arrows indicate reactions and dashed arrows indicate feedback inhibition loops. Chemical species: G6P = glucose-6-phosphate; PEP = phosphoenolpyruvate; ATP = adenosine triphosphate; ADP = adenosine diphosphate; DAHP = 3-deoxy-D-arabino-heptulosonate-7-phosphate; CHR = chorismate; PHP = prephenate; PHE = phenylalanine; TYR = tyrosine; TRP = tryptophan.
postulating the appropriate continuous and discrete constraints, while using MCA information to characterize quantitative interactions in the metabolic system. : This novel optimization formulation was applied to a prototype mathematical model of bacterial amino acid production (Fig. 3) (Schlosser & Bailey, 1990; Hatzimanikatis et al., 1995) for the maximization of phenylalanine selectivity, defined as:
regulatory structures with only two loops, none of which exists in the original pathway, that could increase selectivity by 92% (Fig. 5) (Hatzimanikatis et al., 1995).
3.1.
Sphe =
V4 V4 + V5 + V6
(19)
where the numbers refer to the reaction steps in Fig. 3, while keeping the specific growth rate of the bacterium unaffected. The number of regulatory loops is equal to eight and the number of the enzymes that can be manipulated is equal to six. Therefore, there are 2(8 + 6) = 16 384 alternative solution structures. The optimization calculations indicated that, from the eight feedback inhibitory loops in the original regulatory structure of the pathway, inactivation of at least three loops and overexpression of three enzymes will increase phenylalanine selectivity by 42% without significant effect on the growth rate. The optimal regulatory structures are presented in Fig. 4. Moreover, postulation of a regulatory superstructure and solution of the synthesis problem suggested novel
4. Conclusions The methods presented above show that, after 25 years, MCA has indeed more to say. MCA quantities can be used for the construction of (log)linear kinetic models of metabolism. These kinetic models can simulate the dynamic response characteristics of metabolic systems using simple computational methods, and they have been shown to be a very good first approximation of the original metabolic systems they describe. Moreover, the formulation of a mathematical framework that employs mixed-integer programming optimization algorithms and (log)linear kinetic models allows solution, in a general and systematic way, of the problem of optimizing the regulatory structure of a network (while also considering, simultaneously, beneficial changes in the activities of enzymes in the network). Discussion of the relationship between MCA and BST has been a major concern in theoretical biology. This brief exposition of features and applications of the MCA-based (log)linear model would be incomplete without some comment about the relationship between this (log)linear model and the canonical S-system models used in BST.
. . .
240
1
2 DAHP
3 CHR 6 TRP
4 PHP
PHE
DAHP
CHR
PHP
TRP
TYR
CHR
PHP
TRP
TYR
CHR
PHP
TRP
TYR
CHR
PHP
TRP
TYR
PHE
5 TYR (e)
(a)
DAHP
CHR
PHP
TRP
TYR
PHE
(b)
DAHP
PHE
(f)
DAHP
CHR
PHP
TRP
TYR
PHE
(c)
DAHP
PHE
(g)
DAHP
CHR
PHP
TRP
TYR
PHE
(d)
DAHP
PHE
(h)
F. 4. The eight best solutions for the analysis problem. Solid arrows indicate reactions, dashed arrows indicate inhibitory loops, and thick solid arrows indicate enzyme overexpression for the respective reaction. In solutions b–h the reaction numbering has been omitted for clarity.
BST involves a power-law approximation of the nonlinear kinetic expressions appearing in metabolite mass balances, then proceeds to local stability analysis, numerical simulations, and other results. The (log)linear model is a different approximation to the same original nonlinear metabolite mass balances, which is reached by a different route. Here the original nonlinear equations are linearized in a standard way, and then these approximate linear
equations [eqns (3) and (4)] are transformed into approximate nonlinear eqns in the original variables [see eqns (5–7)] making use of the following Taylor expansion of ln (xi /xi,0 ): ln
0 1 0 1
01
xi xi 1 = ln + xi,0 xi,0 x = x xi x = x i i ,0 i i,0 zxcxv 0 × (xi − xi,0 ) + (higher order terms) (20)
1
DAHP
2
CHR
3
+
–
6 TRP
4
PHP
PHE
241
DAHP
CHR +
5 TYR
PHP
PHE
– TRP
TYR
(d)
(a)
DAHP
CHR
PHP +
TRP
PHE
DAHP
–
CHR
PHP
+
TYR
TRP
(b)
PHE –
TYR
(e)
DAHP
CHR
PHP
PHE
DAHP
–
+ TRP
TYR
(c)
CHR +
PHP
PHE
– TRP
TYR
(f)
F. 5. The six best solutions for the synthesis problem that feature activation and inhibition. Solid and thick arrows as in Fig. 4. Dashed arrows indicate regulation, ( + ) indicates activation and ( − ) inhibition.
to write the deviation in xi (that is xi − xi,0 ) in terms of the logarithmic deviation of xi (ln(xi /xi,0 )). In spite of this unusual route to a simpler, approximate nonlinear representation of the original nonlinear system the (log)linear approach has been shown here, and in other studies, to provide high quality approximations in numerous demanding situations so far examined. Based on analytical considerations and computational experience, the following comparisons can be made between the S-system canonical and (log)linear models: (1) both give the same steady-states; (2) both give the same local stability characteristics at the reference steady state; (3) the dynamic responses of the two representations are in general different; (4) the BST dynamic representation is nonlinear [see, for example, Irvine and Savageau’s (1990) discussion of pertinent numerical methods]
while the (log)linear representation is linear and can be analytically integrated [see eqns (11) and (12)]; (5) the (log)linear model uses explicitly, the same parameters as used in MCA. S-systems models use different parameters, which can be derived from MCA parameters. Metabolism is extremely complicated and only partly understood and characterized (Bailey, 1991). Consequently analysing and engineering metabolic networks can benefit from all of the different viewpoints and tools, such as MCA and BST, which can be brought to bear on these difficult nonlinear problems. The (log)linear model and application frameworks summarized here complement and build upon the prior contributions of metabolic control analysis. This research was supported by the Swiss Priority Program in Biotechnology (SPP Biotech).
. . .
242 REFERENCES
B, J. E. (1974). Periodic Operation of Chemical Reactors: A Review. Chem. Eng. Comm. 1, 111–124. B, J. E. (1977). Periodic Phenomena. In: Chemical Reactor Theory. A Review. (Lapidus, L. & Amundson, N. R., eds) pp. 758–813. Englewood Cliffs, NJ: Prentice-Hall. B, J. E. (1991). Toward a Science of Metabolic Engineering, Science 252, 1668–1675. C-B, A. & C, M. L. (Eds.) (1990). Control of Metabolic Processes, NATO ASI Series A: Life Sciences Vol. 190, New York: Plenum Press. F, D. A. (1992). Metabolic Control Analysis—A Survey of its Theoretical and Experimental Development. Biochem. J. 152, 313–330. F, A. G. (1976). Formulation of Structured Growth Models. Biotechnol. Bioeng. 28, 1481–1486. H, V. & B, J. E. (1996). Effects of Spatiotemporal Variations on Metabolic Control: Approximate Analysis Using (Log)Linear Kinetic Models. Biotechnol. Bioeng. (in press). H, V., F, C. A. & B, J. E. (1996). Analysis and Design of Metabolic Reaction Networks via Mixed-Integer Linear Optimization. AIChE J. 42, 1277–1292. H, R., H¨, H.-G. & S, S. (1987). A Theoretical Approach to the Evolution and Structural Design of Enzymatic Networks. Bull. Math. Biol. 49, 539–595. H, R. & R, T. A. (1974). A linear steady-state treatment of enzymatic chains. Eur. J. Biochem. 42, 89–95. I, D. H. & S, M. A. (1990). Efficient Solution of Nonlinear Ordinary Differential Equations Expressed in S-system Canonical Form. SIAM J. Numer. Anal. 27(3), 704–735. K, H. & B, J. A. (1973). The control of flux. Symp. Soc. Exp. Biol. 27, 65–104. K, D. & W, H. V. (1993). Regulation and Homeostatis in Metabolic Control Theory: Interplay between fluctuations of variables and parameter changes. In: Modern Trends in Biothermokinetics (Schuster, S. et al., eds) pp. 199–204. New York: Plenum Press. K, D. B. & W, H. V. (1986). Metabolic control theory: its role in microbiology and biotechnology. FEMS Microbiol. Rev. 39, 305–320. M, R. A. & D, M. M. (1990). Simple ConstraintedOptimization View of Acetate Overflow in E. coli. Biotechnol. Bioeng. 35, 732–738. M, A. G. (1991). Growth Rate of Escherichia coli. Microbiol. Rev. 55, 316–333. R, D. (1983). A Cybernetic Perspective of Microbial Growth. ACS Symp. Ser. 207, 161–178. R, C. (1988). Metabolic Control Theory: A structural approach. J. theor. Biol. 135, 175–201. R, L., B, I. D. L., W, J. M. & D, P. (1992). Simulation and optimization of metabolic pathways. Comp. Chem. Eng. 16, 237–247. R, M., V D K. & W, H. V. (1991). Control and Thermodynamics of Microbial Growth: Rational tools for Bioengineering. Crit. Rev. Biotech. 11, 367–395.
S, M. A. (1976). Biochemical systems analysis: A study of function and design in molecular biology. Reading, MA: Addison-Wesley. S, P. M. & B, J. E. (1990). An Integrated Modelling-Experimental Strategy for the Analysis of Metabolic Pathways. Math. Biosci. 100, 87–114. S, P, M., H, T. & B, J. E. (1993). Determining Metabolic Sensitivity Coefficients Directly from Experimental Data. Biotechnol. Bioeng. 41, 1027–1038. S, S. & H, R. (1991). Minimization of Intermediate Concentrations as a Suggested Optimality Principle for Biochemical Networks I. Theoretical Analysis. J. Math. Biol. 29, 425–442. S, J. H. & L, L. (1974). Process modeling, estimation and identification. Englewood Cliffs, NJ: Prentice-Hall. V, A. & P, B. O. (1994). Metabolic Flux Balancing: Basic Concepts, Scientific and Practical Use. Bio/Technology 12, 994–998. V, E. O. (1992). Optimization of Integrated Biochemical Systems. Biotechnol. Bioeng. 40, 572–582.
APPENDIX A The linear pathway considered is presented in Fig. 1. The external substrate, s0 , was the manipulated variable. The kinetics used are the following (numbered according to the numbering of the steps in Fig. 1): v0 = 1.75
s0 1.5 + s0
s1 v1 = 2 0.3333 + s1 v2 = 1.7
s2 0.6667 + s2
v3 = 1.5
s3 0.6429 + s3
v4 = 1.3
s4 0.1875 + s4
For the reference steady state we set s0 = 2 which results in the steady-state metabolite concentrations: [s1 ,s2 ,s3 ,s4 ] = [0.3333,0.9524,1.2858,0.6250].
MCA Has More to Say V H†‡ J E. B†§ † Institute of Biotechnology, ETH-Zurich, CH-8093 Zurich, Switzerland and the ‡ Department of Chemical Engineering, California Institute of Technology, Pasadena, CA-91125, U.S.A.
The recent development of two mathematical frameworks for the analysis and design of metabolic systems is reviewed here. Both frameworks employ a (log)linear kinetic metabolic model that is constructed making direct, explicit use of the individual enzyme kinetic parameters which are a foundation of Metabolic Control Analysis (MCA) information. The first framework allows the description of the dynamic responses of metabolic systems subject to fluctuations in their parameters. The second framework considers the problem of optimizing the regulatory structure of metabolic networks. Both frameworks are powerful in enabling greater insights into metabolism based on MCA quantities. 7 1996 Academic Press Limited
understanding of the determinants of metabolic fluxes (Kacser & Burns, 1973). Since its introduction 30 years ago, MCA has received much theoretical and experimental attention [an excellent survey is found in Fell (1992)]. Many theoretical and experimental methodologies have been developed that allow determination of quantitative indices defined within MCA, such as control coefficients and elasticities. The description of metabolic systems by MCA quantities is commonly used because of limitations in the available information about these systems. Furthermore, when a nonlinear model is available, it can be linearized and well-studied within the same framework. One of the major limitations of the standard formulation of MCA has been its restriction to a system in a steady state. The MCA framework has not been used to describe transient phenomena that arise from fluctuations in metabolic parameters and inputs (such as enzyme expression levels, and the concentrations of external substrates and independent effectors). However, as will be shown here, this limitation has been overcome recently by the development of a (log)linear kinetic model of metabolic systems based explicitly on MCA information (i.e., elasticities and control coefficients).
1. Introduction Metabolism is an extremely complicated, multicomponent, dynamic system. In addition to numerous different chemical reactions and intermolecular interactions, many coupled through common cofactors, metabolism is complicated by an overlay of multiple control systems operating at both genetic and protein levels. Understanding the genetic bases for metabolism, and how metabolism is influenced by environment, is central to biology, medicine, and biotechnology in many respects (Bailey, 1991). Metabolic control analysis (MCA), a sensitivity analysis framework, is one of the most developed methods for the quantitative description of metabolism and microbial physiology that allows the analysis and study of the responses of metabolic systems to changes in their parameters (Kacser & Burns, 1973; Heinrich & Rapoport, 1974; Kell & Westerhoff, 1986; Cornish-Bowden & Cardenas, 1990; Rutgers et al., 1991; Fell, 1992). Kacser and Burns’ formulation and lucid explanation of MCA paved the way for an ongoing dialogue between biology and mathematics, aimed at quantitative § Author to whom correspondence should be addressed. 0022–5193/96/190233 + 10 $25.00/0
233
7 1996 Academic Press Limited
234
. . .
This model can be used to simulate dynamic responses of the system subject to fluctuations in its inputs and parameters. The major advantages of the (log)linear model are the specification of the system dynamics explicitly using the same parameters as employed in MCA, and the linearity (in terms of logarithms) of the model. As discussed in more detail later, after examining two applications of the (log)linear representation, neither of these features are shared by the well-known Biochemical Systems Theory (BST) (Savageau, 1976), which also considers the dynamics of metabolic systems. Moreover, the formulation of this (log)linear model allowed the development of a novel mathematical framework for determining changes in the structure and strength of metabolic regulation which should be considered in order to optimize a particular metabolic process. A mixed-integer linear programming (MILP) formulation is proposed for the general case of optimizing performance of metabolic processes which can be adequately approximated using a (log)linear model. The solution of the MILP formulation provides information on which enzymes should be present at different levels, the extent of such changes needed, and the accompanying modifications in the regulatory structure that will optimize the metabolic process. The results of an application of this method are also presented, illustrating the power of the method and the additional results which can be obtained based on MCA quantities.
2. A (Log)Linear Kinetic Model of Metabolic Reaction Networks Prediction of the changes in metabolic fluxes which will occur after a change in any parameter affecting metabolism requires a kinetic representation relating fluxes to metabolic concentrations and system parameters. In general such models are not available, and the kinetic model must therefore be developed from measurements of fluxes and metabolite concentrations. A variety of formulations are possible. One attractive choice is reviewed next (Hatzimanikatis & Bailey, 1995). Consider a metabolic system consisting of n metabolites and m enzymatically-catalysed reactions. We are interested in studying how modifications of the metabolic parameters and of the properties of the enzymes that catalyze these reactions affect the time response characteristics of metabolic functions of the system, such as metabolic concentrations, fluxes, and the specific growth rate. Assuming that any concentration gradients in the cell are of negligible
significance, mass balances on the metabolites of the system may be written as: dx = f(v(x;p(t)),x;p(t),t) dt
(1)
where x is the n-dimensional metabolite concentration vector, f is a function determined by the mass balances, v is the m-dimensional reaction rate vector, and p is the s-dimensional manipulated parameter vector (e.g., enzyme concentrations, external nutrient or substrate concentrations). In addition to metabolic reaction rates, the mass balance equations also include terms that account for other processes by which concentrations of metabolites change (such as the dilution brought about by increases in the biomass volume (Fredrickson, 1976) and transport through the cell envelope). In addition, we consider the r-dimensional vector of metabolic outputs, h, which may be written in general in the form: h = h(v(x;p(t)),x;p(t),t).
(2)
In eqn (2) h is a function of the rates of interest, of the metabolic concentration, and of the parameters. Let x0 be a steady-state—or one of the steady states in the case of steady-state multiplicity—that corresponds to the parameter vector values p0 , and let x0 be non-zero with positive elements. Linearization around this steady state will result in the linear system: d(x − x0 ) 1f = dt 1x
b
x0 ,p0
(x − x0 ) +
1f 1p
b
x0 ,p0
(p(t) − p0 ) (3)
for the mass balances, and
h = h(v0 (x0 ;p0 ),x0 ;p0 ,0) +
1h 1x
b
+
x0 ,p0
1h 1p
(x − x0 )
b
x0 ,p0
(p(t) − p0 )
(4)
for the metabolic outputs. In order to describe the dependence of the metabolic system explicitly on the reaction rates and on the rest of the metabolic processes, the following final representation is adopted, after the appropriate scaling of the linearized system [eqns (3) and (4)] and transformation of variables and parameters
(Hatzimanikatis & Bailey, 1995; Hatzimanikatis et al., 1995): dz = Nez + kz + NPq + Lq dt
(5)
w = Jez + Hz + JPq + Uq
(6)
zi = ln(xi /xi,0 ) qk = ln(pk /pk,0 )
(7)
wl = ln(hl /hl,0 ) and N, J, k, L, H, U, e, and P, are matrices, defined as:
6
vj,0 1fi xi,0 1vj
0 1 7
6
vj,0 1hl hl,0 1vj
J = jl,j =jl,j =
x0 ,p0
0 1 7
6
xk,0 1fi xi,0 1xk
6
pk,0 1fi xi,0 1pk
6
xi,0 1hl hl,0 1xi
6
pk,0 1hl hl,0 1pk
k = ki,k =ki,k =
L = li,k =li,k =
H = hl,i =hl,i =
U = ul,k =ul,k =
6
e = ej,i =ej,i =
6
,
0 1 7 x0 ,p0
0 1 7 x0 ,p0
,
,
0 1 7
,
x0 ,p0
0 1 7 x0 ,p0
,
0 1 7
xi,0 1vj vj,0 1xi
P = pj,k =pj,k =
,
x0 ,p0
x0 ,p0
,
0 1 7
pk,0 1vj vj,0 1pk
x0 ,p0
(8)
At steady state, the solution of (5) and (6) yields: w = Cq
where C = −(Je + H)(Ne + k)−1(NP + L) + JP + U (10) with
where, z, q, and w, are the logarithmic deviations of the metabolic concentrations, the enzyme levels, and the metabolic outputs, respectively:
N = ni,j =ni,j =
235
(9)
6
C = cl,k \ cl,k =
0 17
pk,0 dhl hl,0 dpk
p0
The (log)linear kinetic metabolic model described above depends on the same information as that employed within the framework of metabolic control analysis (MCA) (Reder, 1988; Cornish-Bowden & Cardenas, 1990; Schlosser & Bailey, 1990). Matrices e and P are the elasticity matrices with respect to metabolites and to parameters, respectively. Matrix N is a generalized scaled stoichiometric matrix for the mass balances, and matrix J is the scaled stoichiometric matrix for the outputs, i.e. the linear relations between the reaction rates and the metabolic outputs. Matrices k, H, L, and U, could be defined as the elasticities with respect to metabolites (k and H) and to parameters (L and U) for metabolic processes that are independent of the reaction rates, v. The matrix C is the control coefficient matrix of the metabolic functions h with respect to parameters p. Equation (10) defines the control coefficients in a very general and systematic fashion independently of the stoichiometry of the metabolic network and the complexity of the metabolic outputs. The coefficient cl,k in the l-th row and k-th column of C represents the fractional change of metabolic function hl , in response to a fractional change in metabolic parameter or input, pk , which is the mathematical definition of the control coefficient of hl with respect to pk . Experimental determination of the parameters for this (log)linear system has been the subject of several studies (Fell, 1992; Cornish-Bowden & Cardenas, 1990; Schlosser & Bailey, 1990; Schlosser et al., 1993) and in many cases such parameters are the only available kinetic information about a metabolic system. We should notice here that the final representation [eqns (5) and (6)] allows the explicit description of the system with quantities that are characteristic of each enzyme. Therefore, we could study the effects of modifications of the catalytic properties of enzyme i with respect to its substrate (or regulatory effector) j, by changing the value of the corresponding elasticity eij . One of the main advantages of the (log)linear model is its ability to describe the time responses of the system to transient changes in the metabolic parameters or to temporal perturbations in the
. . .
236
concentration of the metabolites. The linearity of eqn (5) allows us to solve explicitly for z as a function of time (Seinfeld & Lapidus, 1974): z(t) = e (Ne + k)(t − t0 )z0 +
g
×(NP + L)q(t)dt
(11)
z0 = z(t0 ) For the time dependence of the logarithmic deviations of the metabolic outputs we have from eqns (6) and (11)
0
w(t) = (Je − H) e (Ne + k)(t − t0 )z0 +
g
t
e (Ne + k)(t − t)
t0
1
×(NP + L)q(t)dt + (JP + U)q(t)
(12)
The linearity of the model permits fast numerical calculation of the integral in eqns (11) and (12). For the time-dependent concentration of the metabolites and for the metabolic outputs we simply have to calculate the transformations: xi (t) = x0,i e zi (t)
(13)
hl (t) = h0,l e wl (t)
(14)
and
which can also be written as
k=1
0 1 0 1 xk (t) x0,k
alk
pj (t) p0,j
s
P
j=1
blj
(15)
where
0
n
1
(16)
1
(17)
alk = s jli eik + hlk i=1
and
0
n
blj = s jli pil + ulj i=1
Equations (11–17) account for the time dependence of metabolites and metabolic outputs based explicitly on standard MCA quantities and thus substantially extend the scope of metabolic analysis and design which can be achieved using these quantities. : A simple linear pathway (Fig. 1) was studied as a test of the performance of the (log)linear metabolic 2.1.
2 S2
3 S3
4 S4
P
F. 1. Linear pathway with Michaelis-Menten kinetics used for assessing the (log)linear dynamic model.
e (Ne + k)(t − t)
where
n
1 S1
t
t0
hl = h0,l P
0 S0
model in approximating a given nonlinear metabolic system (which here is the full set of unsteady-state nonlinear material balances on the components in this pathway) (Hatzimanikatis & Bailey, 1995). In the nonlinear system, Michaelis-Menten kinetics were considered for every reaction rate and the parameters used are presented in Appendix A. The (log)linear model parameters were derived as described above. We studied the dynamic responses of the flux through step 4, to the final product P, for step changes [Fig. 2(a)] and for sinusoidal variation [Figs 2(b), (c) and (d)] of the external substrate, s0 . The time response characteristics of the flux through the pathway are in good agreement between the original nonlinear model and the approximate (log)linear model. The observed deviations are due to the approximations inherent in the MCA quantities used to construct the (log)linear model. The final steady-state differences for step processes [Fig. 2(a)] are the differences between the MCA calculations of new steady states and the responses of the original system. The question of the sensitivity of a system either to changes in its parameters or to fluctuations of the concentrations of intracellular metabolites has been addressed and studied using the simple horizontal model of the system; i.e., the system described by eqns (3) and (4) (Kahn & Westerhoff, 1993). However, the (log)linear model is in much better agreement with the original system model than is a simple linear model, especially with respect to response to periodic inputs [Figs 2(b), (c) and (d)]. The excellent performance of the (log)linear model exhibited here under temporal variation of system parameters suggests that MCA information can be used directly in order to analyse the behavior of metabolic systems subject to fluctuations in their parameters. In particular, metabolic systems, like any nonlinear dynamical system, can exhibit a shift in time-averaged properties (such as metabolic fluxes and control coefficients) in response to fluctuations with zero-time-average value (Bailey, 1974, 1977). The (log)linear dynamic model presented here provides a powerful, convenient vehicle for analysis of such phenomena (Hatzimanikatis & Bailey, 1995). Good agreement between response characteristics of the (log)linear model and of the nonlinear model has also been observed for various more complicated and strongly nonlinear metabolic systems [such as a
(a)
1.20
237
(c) II
1.10
1.000
I 1.00 III
0.980
IV
0.960
Flux to product
0.90 0.80 0.70 (b)
(d)
1.000
1.000 0.990 0.995 0.980 0.990 0.970 0.985
0
20
40 Time
60
0
20.0
40.0
60.0
Time
F. 2. Dynamic responses of the flux through the linear pathway. Solid lines correspond to the (log)linear model, dashed lines correspond to the nonlinear model, and dotted-dashed lines correspond to the linear model. Vertical axes correspond to the flux to product P through reaction step 4. A. Responses to step changes of the external substrate, S0 . I: +20%; II: +50%; III: −20%; IV: −50%. B–D. Responses to sinusoidal variation of the external substrate with different amplitudes and frequencies: B: s0 = 2(1 − 0.2 sin(t)). C: s0 = 2(1 − 0.5 sin(t)). D: s0 = 2(1 − 0.5 sin(t/0.1)).
branched pathway and yeast glycolysis (Hatzimanikatis & Bailey, 1995)]. This agreement suggests that, in the absence of a nonlinear model, and when MCA quantities can be experimentally determined or estimated, the (log)linear model can be used as a good first approximation for analysing the dynamic response characteristics of metabolic systems. To summarize, the (log)linear model is an attractive representation of metabolic pathways since it can be easily constructed using MCA information, it has an unsteady-state analytical solution which makes computation easier, and it can describe very well the dynamic response characteristics of metabolic networks. Moreover, the analytical solution of the (log)linear model allows the derivation of indices such as the regulatory strength and the homeostatic strength (Kahn & Westerhoff, 1993) by simple manipulations of eqns (11–15) (Hatzimanikatis & Bailey, in preparation). 3. Optimization of Regulatory Architectures: Analysis and Synthesis Problems Metabolic reaction networks in living organisms have evolved through natural selection and are therefore configured to maximize the probability of survival of the species, at least insofar as evolution has
progressed. While explicit formulation of the objective function for natural metabolism is non-trivial, prior investigators have proposed maximization of growth rates or most efficient utilization of cellular energetic and chemical resources as the objective function for evolution of natural metabolism (Savageau, 1976; Ramkrishna, 1983; Heinrich et al., 1987; Marr, 1991; Schuster & Heinrich, 1991). However, biotechnological objectives require the identification of a different configuration of fluxes which directs raw materials to products efficiently at high rates and in the presence of high concentrations of product (Bailey, 1991). Therefore, a productionoriented evolution is needed to achieve these goals. Optimization approaches to the analysis of biochemical systems for improvement of biotechnological objectives have been followed in various studies (Majewski & Domach, 1990; Regan et al., 1992; Voit, 1992; Varma & Palsson, 1994) for the identification of the optimum manipulation of the external inputs to the system (such as independent effectors and external substrates) and of enzyme expression levels. However, these studies do not address the problem of optimizing the regulatory structure of the metabolic network. Optimization of a metabolic function with respect to changes in the amounts of enzymes in the pathway and in
. . .
238
modifications in regulatory characteristics of those enzymes can be undertaken using the kinetic description provided by the (log)linear metabolic model (Hatzimanikatis et al., 1995). In this optimization context, the model defines constraints in the forms of steady-state mass balances for the intracellular metabolites and of (log)linear expressions of metabolic functions that should be constrained within certain values to force the optimization search to stay in physiologically reasonable territory. When MCA information is available for a metabolic system, the information concerning the regulatory structure is included in the elasticity matrix e, which can be written as a sum of two matrices: e = e s + er
(18)
where the elements in matrix es correspond to the substrate elasticities of the enzymes, that is, the sensitivities of enzyme activities with respect to concentrations of their substrates, and the elements of matrix er correspond to the regulatory elasticities of the enzymes, that is, the sensitivities of enzyme activities with respect to concentrations of regulatory metabolites. In this representation, the substrates themselves can also be considered as regulatory metabolites (in cases such as substrate inhibition):
6
es = ej,is =
0 1
xis ,0 1vj vj,0 1xis
x0 ,p0
=xis is a substrate for reaction j
7
and
6
er = ej,ir =
0 1
xir , 0 1vj vj,0 1xir
x0 ,p0
=xir is a regulator for reaction j
7
Changes of the elements in matrix er from non-zero values to zero or vice versa correspond to modifications in the regulatory structure of the system. The problem of optimizing the performance of a metabolic network by modifying its regulatory structure has been defined in two different contexts, called analysis and synthesis problems (Hatzimanikatis et al., 1995). In the context of a given regulatory structure, which can be modified by
deleting certain control interactions, the analysis problem was defined as follows: Which of the existent regulatory loops should be inactivated, and what associated changes should be made in the manipulated variables (e.g., enzyme expression levels, environmental conditions, effectors external to the system) in order to optimize the performance of the metabolic network? In a more extended context a regulatory superstructure can be postulated, in which a set of alternative regulatory elasticities for each enzyme, corresponding to different kinds of regulation by each metabolite, is embedded. Then, the synthesis problem considers the following equations: What kind of regulation (i.e., activation or inhibition, by which metabolite and of what strength) should be assigned to each enzyme in the network, and what associated changes should be made in the manipulated parameters (e.g., enzyme expression levels, environmental conditions, effectors external to the system), in order to optimize the performance of the metabolic network? Depending on the measure of ‘‘the performance of the metabolic network’’ (e.g., the metabolic function to be optimized), the answer(s) to the above questions can suggest possible evolutionary criteria that gave rise to naturally existing regulatory structures, or, in a biotechnological context, provide potential useful quantitative guidance for promising targets for protein engineering to achieve a preferred flux distribution. Inactivation or activation of a regulatory loop is equivalent to elimination or introduction of non-zero terms in the er matrix in eqn (18). Moreover, the synthesis problem will typically be subject to some constraints such as the possible number of regulatory actions on each enzyme and the requirement that an enzyme cannot be activated and inhibited by the same metabolite. These kind of decisions, present in both the analysis and synthesis problems, are discrete decisions. Biochemical analysis criteria or experimental limitations might also constrain the maximum number of continuously adjustable manipulated parameters, imposing another set of discrete constraints. The mixed discrete and continuous nature of the problem and the linear (in terms of logarithms) description of the system led to the formulation of the analysis and synthesis problems as MILP problems (Hatzimanikatis et al., 1995). A multitude of questions can be addressed within this framework by
239
PEP G6P + PEP
1
2
DAHP ATP
CHR
3
PHP
4
PHE
4ATP
ADP G6P 6
5
4ADP TRP
PHE
9
Biomass
7
TYR
TRP
8 TYR F. 3. The aromatic amino acid synthesis pathway. Solid arrows indicate reactions and dashed arrows indicate feedback inhibition loops. Chemical species: G6P = glucose-6-phosphate; PEP = phosphoenolpyruvate; ATP = adenosine triphosphate; ADP = adenosine diphosphate; DAHP = 3-deoxy-D-arabino-heptulosonate-7-phosphate; CHR = chorismate; PHP = prephenate; PHE = phenylalanine; TYR = tyrosine; TRP = tryptophan.
postulating the appropriate continuous and discrete constraints, while using MCA information to characterize quantitative interactions in the metabolic system. : This novel optimization formulation was applied to a prototype mathematical model of bacterial amino acid production (Fig. 3) (Schlosser & Bailey, 1990; Hatzimanikatis et al., 1995) for the maximization of phenylalanine selectivity, defined as:
regulatory structures with only two loops, none of which exists in the original pathway, that could increase selectivity by 92% (Fig. 5) (Hatzimanikatis et al., 1995).
3.1.
Sphe =
V4 V4 + V5 + V6
(19)
where the numbers refer to the reaction steps in Fig. 3, while keeping the specific growth rate of the bacterium unaffected. The number of regulatory loops is equal to eight and the number of the enzymes that can be manipulated is equal to six. Therefore, there are 2(8 + 6) = 16 384 alternative solution structures. The optimization calculations indicated that, from the eight feedback inhibitory loops in the original regulatory structure of the pathway, inactivation of at least three loops and overexpression of three enzymes will increase phenylalanine selectivity by 42% without significant effect on the growth rate. The optimal regulatory structures are presented in Fig. 4. Moreover, postulation of a regulatory superstructure and solution of the synthesis problem suggested novel
4. Conclusions The methods presented above show that, after 25 years, MCA has indeed more to say. MCA quantities can be used for the construction of (log)linear kinetic models of metabolism. These kinetic models can simulate the dynamic response characteristics of metabolic systems using simple computational methods, and they have been shown to be a very good first approximation of the original metabolic systems they describe. Moreover, the formulation of a mathematical framework that employs mixed-integer programming optimization algorithms and (log)linear kinetic models allows solution, in a general and systematic way, of the problem of optimizing the regulatory structure of a network (while also considering, simultaneously, beneficial changes in the activities of enzymes in the network). Discussion of the relationship between MCA and BST has been a major concern in theoretical biology. This brief exposition of features and applications of the MCA-based (log)linear model would be incomplete without some comment about the relationship between this (log)linear model and the canonical S-system models used in BST.
. . .
240
1
2 DAHP
3 CHR 6 TRP
4 PHP
PHE
DAHP
CHR
PHP
TRP
TYR
CHR
PHP
TRP
TYR
CHR
PHP
TRP
TYR
CHR
PHP
TRP
TYR
PHE
5 TYR (e)
(a)
DAHP
CHR
PHP
TRP
TYR
PHE
(b)
DAHP
PHE
(f)
DAHP
CHR
PHP
TRP
TYR
PHE
(c)
DAHP
PHE
(g)
DAHP
CHR
PHP
TRP
TYR
PHE
(d)
DAHP
PHE
(h)
F. 4. The eight best solutions for the analysis problem. Solid arrows indicate reactions, dashed arrows indicate inhibitory loops, and thick solid arrows indicate enzyme overexpression for the respective reaction. In solutions b–h the reaction numbering has been omitted for clarity.
BST involves a power-law approximation of the nonlinear kinetic expressions appearing in metabolite mass balances, then proceeds to local stability analysis, numerical simulations, and other results. The (log)linear model is a different approximation to the same original nonlinear metabolite mass balances, which is reached by a different route. Here the original nonlinear equations are linearized in a standard way, and then these approximate linear
equations [eqns (3) and (4)] are transformed into approximate nonlinear eqns in the original variables [see eqns (5–7)] making use of the following Taylor expansion of ln (xi /xi,0 ): ln
0 1 0 1
01
xi xi 1 = ln + xi,0 xi,0 x = x xi x = x i i ,0 i i,0 zxcxv 0 × (xi − xi,0 ) + (higher order terms) (20)
1
DAHP
2
CHR
3
+
–
6 TRP
4
PHP
PHE
241
DAHP
CHR +
5 TYR
PHP
PHE
– TRP
TYR
(d)
(a)
DAHP
CHR
PHP +
TRP
PHE
DAHP
–
CHR
PHP
+
TYR
TRP
(b)
PHE –
TYR
(e)
DAHP
CHR
PHP
PHE
DAHP
–
+ TRP
TYR
(c)
CHR +
PHP
PHE
– TRP
TYR
(f)
F. 5. The six best solutions for the synthesis problem that feature activation and inhibition. Solid and thick arrows as in Fig. 4. Dashed arrows indicate regulation, ( + ) indicates activation and ( − ) inhibition.
to write the deviation in xi (that is xi − xi,0 ) in terms of the logarithmic deviation of xi (ln(xi /xi,0 )). In spite of this unusual route to a simpler, approximate nonlinear representation of the original nonlinear system the (log)linear approach has been shown here, and in other studies, to provide high quality approximations in numerous demanding situations so far examined. Based on analytical considerations and computational experience, the following comparisons can be made between the S-system canonical and (log)linear models: (1) both give the same steady-states; (2) both give the same local stability characteristics at the reference steady state; (3) the dynamic responses of the two representations are in general different; (4) the BST dynamic representation is nonlinear [see, for example, Irvine and Savageau’s (1990) discussion of pertinent numerical methods]
while the (log)linear representation is linear and can be analytically integrated [see eqns (11) and (12)]; (5) the (log)linear model uses explicitly, the same parameters as used in MCA. S-systems models use different parameters, which can be derived from MCA parameters. Metabolism is extremely complicated and only partly understood and characterized (Bailey, 1991). Consequently analysing and engineering metabolic networks can benefit from all of the different viewpoints and tools, such as MCA and BST, which can be brought to bear on these difficult nonlinear problems. The (log)linear model and application frameworks summarized here complement and build upon the prior contributions of metabolic control analysis. This research was supported by the Swiss Priority Program in Biotechnology (SPP Biotech).
. . .
242 REFERENCES
B, J. E. (1974). Periodic Operation of Chemical Reactors: A Review. Chem. Eng. Comm. 1, 111–124. B, J. E. (1977). Periodic Phenomena. In: Chemical Reactor Theory. A Review. (Lapidus, L. & Amundson, N. R., eds) pp. 758–813. Englewood Cliffs, NJ: Prentice-Hall. B, J. E. (1991). Toward a Science of Metabolic Engineering, Science 252, 1668–1675. C-B, A. & C, M. L. (Eds.) (1990). Control of Metabolic Processes, NATO ASI Series A: Life Sciences Vol. 190, New York: Plenum Press. F, D. A. (1992). Metabolic Control Analysis—A Survey of its Theoretical and Experimental Development. Biochem. J. 152, 313–330. F, A. G. (1976). Formulation of Structured Growth Models. Biotechnol. Bioeng. 28, 1481–1486. H, V. & B, J. E. (1996). Effects of Spatiotemporal Variations on Metabolic Control: Approximate Analysis Using (Log)Linear Kinetic Models. Biotechnol. Bioeng. (in press). H, V., F, C. A. & B, J. E. (1996). Analysis and Design of Metabolic Reaction Networks via Mixed-Integer Linear Optimization. AIChE J. 42, 1277–1292. H, R., H¨, H.-G. & S, S. (1987). A Theoretical Approach to the Evolution and Structural Design of Enzymatic Networks. Bull. Math. Biol. 49, 539–595. H, R. & R, T. A. (1974). A linear steady-state treatment of enzymatic chains. Eur. J. Biochem. 42, 89–95. I, D. H. & S, M. A. (1990). Efficient Solution of Nonlinear Ordinary Differential Equations Expressed in S-system Canonical Form. SIAM J. Numer. Anal. 27(3), 704–735. K, H. & B, J. A. (1973). The control of flux. Symp. Soc. Exp. Biol. 27, 65–104. K, D. & W, H. V. (1993). Regulation and Homeostatis in Metabolic Control Theory: Interplay between fluctuations of variables and parameter changes. In: Modern Trends in Biothermokinetics (Schuster, S. et al., eds) pp. 199–204. New York: Plenum Press. K, D. B. & W, H. V. (1986). Metabolic control theory: its role in microbiology and biotechnology. FEMS Microbiol. Rev. 39, 305–320. M, R. A. & D, M. M. (1990). Simple ConstraintedOptimization View of Acetate Overflow in E. coli. Biotechnol. Bioeng. 35, 732–738. M, A. G. (1991). Growth Rate of Escherichia coli. Microbiol. Rev. 55, 316–333. R, D. (1983). A Cybernetic Perspective of Microbial Growth. ACS Symp. Ser. 207, 161–178. R, C. (1988). Metabolic Control Theory: A structural approach. J. theor. Biol. 135, 175–201. R, L., B, I. D. L., W, J. M. & D, P. (1992). Simulation and optimization of metabolic pathways. Comp. Chem. Eng. 16, 237–247. R, M., V D K. & W, H. V. (1991). Control and Thermodynamics of Microbial Growth: Rational tools for Bioengineering. Crit. Rev. Biotech. 11, 367–395.
S, M. A. (1976). Biochemical systems analysis: A study of function and design in molecular biology. Reading, MA: Addison-Wesley. S, P. M. & B, J. E. (1990). An Integrated Modelling-Experimental Strategy for the Analysis of Metabolic Pathways. Math. Biosci. 100, 87–114. S, P, M., H, T. & B, J. E. (1993). Determining Metabolic Sensitivity Coefficients Directly from Experimental Data. Biotechnol. Bioeng. 41, 1027–1038. S, S. & H, R. (1991). Minimization of Intermediate Concentrations as a Suggested Optimality Principle for Biochemical Networks I. Theoretical Analysis. J. Math. Biol. 29, 425–442. S, J. H. & L, L. (1974). Process modeling, estimation and identification. Englewood Cliffs, NJ: Prentice-Hall. V, A. & P, B. O. (1994). Metabolic Flux Balancing: Basic Concepts, Scientific and Practical Use. Bio/Technology 12, 994–998. V, E. O. (1992). Optimization of Integrated Biochemical Systems. Biotechnol. Bioeng. 40, 572–582.
APPENDIX A The linear pathway considered is presented in Fig. 1. The external substrate, s0 , was the manipulated variable. The kinetics used are the following (numbered according to the numbering of the steps in Fig. 1): v0 = 1.75
s0 1.5 + s0
s1 v1 = 2 0.3333 + s1 v2 = 1.7
s2 0.6667 + s2
v3 = 1.5
s3 0.6429 + s3
v4 = 1.3
s4 0.1875 + s4
For the reference steady state we set s0 = 2 which results in the steady-state metabolite concentrations: [s1 ,s2 ,s3 ,s4 ] = [0.3333,0.9524,1.2858,0.6250].