Nonlinear Metabolic Control Analysis

metabolic engineering 1, 7587 (1999) article no. MT980108

Nonlinear Metabolic Control Analysis Vassily Hatzimanikatis* Computational Biology Group, Dupont Central Research 6 Development, Experimental Station, Wilmington, DE 19880-0328

Mechanistic models have been traditionally used for the analysis of metabolic systems. This kind of model requires knowledge of the metabolic stoichiometry, of the kinetic characteristics of the system, i.e., rate expressions and kinetic constants, and of several in vivo variables, such as metabolite concentrations and metabolic fluxes, for model validation and, in some cases, for parameter identification. Consequently, the construction of such models requires numerous data and experiments specifically designed to generate these classes of data. Furthermore, when modelers use kinetic parameter values determined by in vitro experiments the results are always subject to criticism due to the uncertainty of how representative these values are of the corresponding in vivo parameters. Alternatively, if one simply tries to study a model system's dependence on parameter values, then one is faced with the large number of the parameters that make conceptual understanding difficult. In contrast to mechanistic models, metabolic control analysis (MCA), introduced in the biology community 25 years ago, proposes a quantitative description of metabolic systems which is based on information derived from a smaller number of experiments, and uses simplified kinetic parameters (Kacser and Burns, 1973; Cornish-Bowden and Cardenas, 1990; Fell 1992; Fell 1997). Within the MCA framework, the responses of metabolic systems to the manipulation of metabolic parameters, such as enzyme expression levels, can be expressed in terms of control coefficients (the percentage change in the observable metabolic output per unit percent change in the enzyme expression level). The kinetic properties of the metabolic processes can be described in terms of elasticities (the percentage change in the process for a unit percentage change in the associated metabolic variable, such as metabolite concentration, when the rest of the metabolic parameters and variables remain invariant.) Advances in recombinant DNA technology and analytical biochemistry techniques enable sophisticated manipulation of metabolic parameters and precise measurements of intracellular variables that provide the information for calculating MCA parameters (i.e., elasticities and control coefficients). These parameters can be used for the construction of (log)linear MCA models that describe the linear correlation between the logarithmic deviations of the metabolic variables and the logarithmic deviations of the

Mathematical description of metabolic systems allows the calculation of the expected responses of metabolism to genetic modifications and the identification of the most promising targets for metabolic engineering. Metabolic control analysis (MCA) provides such a description in the form of quantitative indices (elasticities and control coefficients). These indices are determined by perturbation experiments around a reference steady state and, therefore, the predictive power of MCA is limited to small changes in the metabolic parameters. The modeling framework introduced here allows accurate description of the metabolic responses over wide range of changes in the metabolic parameters. The framework requires information about the MCA indices at the reference state and the corresponding values of the metabolic reaction rates, and employs simplifying assumptions about the reaction mechanisms. It is shown that knowledge of the intracellular metabolite concentrations is not necessary for the application of the framework. The performance of the methodology is illustrated using three elementary metabolic systems that display highly nonlinear responses to the modification in their parameters: an unbranched pathway, an interconvertible enzyme system, and a branched pathway subject to feedback inhibition.  1999 Academic Press

INTRODUCTION The formulation of mathematical models and the development of mathematical frameworks for the analysis of these models have been proven to be invaluable approaches in science and engineering by providing insight and guidance for the design of the systems these models represent. Within the first demarcation of metabolic engineering science, mathematical methods have been identified as an essential component of the emerging field as they ``are the only way that net consequences of simultaneous, coupled, and often counteracting processes can be simulated and evaluated consistently and quantitatively'' (Bailey, 1991) Mathematical methodologies have been successfully applied for the identification of promising targets for metabolic engineering and for the explanation of unexpected experimental results. * Current address: Cargill Corn Milling Division, 2301 Crosby Road, Wayzata, MN 55391-2397. Fax: 612-742-2381. E-mail: vassilyinnocent. com.

75

1096-717699 30.00 Copyright  1999 by Academic Press All rights of reproduction in any form reserved.

metabolic engineering 1, 7587 (1999) article no. MT980108

Vassily Hatzimanikatis

metabolic parameters. (Log)linear models can adequately simulate and predict the responses of a metabolic system to spatiotemporal variations in the metabolic parameters around a reference steady state (Hatzimanikatis and Bailey, 1997.) Analysis of a (log)linear model can provide guidance for engineering of the metabolic system it represents by identifying the metabolic parameters that should be manipulated in order to achieve the desired goal. In metabolic engineering practice, the metabolic parameters are subject to manipulation within wide ranges dictated by the desirable objective or by the available methods. The outcome of such large changes cannot be predicted by (log)linear models since these models are approximations around a specific reference steady state and their predictive power is limited to infinitesimal changes in their parameters. Therefore, (log)linear models based on MCA information around a specific steady state will most likely lead to false expectations. Recognizing these limitations several researchers have addressed this problem (Hofer 6 Heinrich, 1993; Small 6 Kacser, 1993; Kacser 6 Acerenza, 1993; Small 6 Kacser, 1994). We present here a novel framework which addresses the limitations of (log)linear models discussed above. The key concept of the framework arises from the realization that the parameters of the (log)linear models, namely the elasticities, depend on the metabolic parameters and on the concentrations of the metabolites, which also only depend on the metabolic parameters. We show that, for the most common enzymatic rate expressions, the changes in the elasticities associated with the changes in the metabolic parameters depend only on the values of the elastisticities at the reference state. Therefore, the relatively simple knowledge of the elasticities around the reference state and some assumptions on the enzymatic rate expressions are sufficient information to predict the nonlinear response of a metabolic system to large changes in its parameters.

additional examples are presented in the next section, one of which features enzyme inhibition. The mass balance for the first metabolite can be written as follows: dX 1 =V 1 &V 2 . dt

For a time-invariant concentration of the external substrate, the system will reach a steady state. Linearizing around this state, the mass balance becomes d(X 1 &X 1, 0 ) V 1 = dt S

}

(S&S 0 )& 0

V 2 X 1

} (X &X 1

1, 0

),

(2)

0

where the subscript ``0'' denotes the steady-state values of the corresponding variables. After the appropriate scaling is introduced, the linearized mass balance can be written as d(X 1 &X 1, 0 )X 1, 0 dt =

V 1, 0 S 0 V 1 X 1, 0 V 1, 0 S

}

(S&S 0 ) V 2, 0 X 1, 0 V 2 & S0 X 1, 0 V 2, 0 X 1 0

}

(X 1 &X 1, 0 ) . X 1, 0 0 (3)

From the MCA definitions we can identify the elasticities in equation (3): = VX21 #= 21 #

X 1, 0 V 2  ln V 2 = V 2, 0 X 1 0  ln X 1

}

}

(4a) 0

and ? VS 1 #? 1S #

S 0 V 1  ln V 1 = . V 1, 0 S 0  ln S 0

}

}

(4b)

These are the elasticities of the corresponding enzymes with respect to metabolites and to metabolic parameters, respectively. The approximate linear equation (3) can be transformed into an approximate nonlinear equation in the original variables employing the Taylor expansion:

NONLINEAR METABOLIC CONTROL ANALYSIS An unbranched enzymatic reaction pathway (Fig. 1) will be used here as an example that can help us review the basics of MCA and demonstrate the methodology we call nonlinear metabolic control analysis (NMCA). The kinetics of the reaction steps are assumed to be MichaelisMenten and they are presented in Appendix I. The external substrate, S, is considered to be the independently manipulated metabolic parameter. This particular model system has been chosen for its simplicity to illustrate the ability of NMCA to predict nonlinear responses of metabolic systems (e.g. changes in metabolite concentrations) to large changes in metabolic parameters (e.g. external substrate). Two

ln

Xi X i, 0

\ + =ln

X i, 0 1 + (X &X i, 0 )+(higher order terms). X i, 0 X i, 0 i

\ + =0

(5)

Defining logarithmic deviations of the metabolite concentrations and the metabolic parameters as new variables: 76

File: DISTL2 010802 . By:JB . Date:03:09:98 . Time:11:02 LOP8M. V8.B. Page 01:01 Codes: 5415 Signs: 4102 . Length: 57 pic 1 pts, 240 mm

(1)

metabolic engineering 1, 7587 (1999) article no. MT980108

Nonlinear Metabolic Control Analysis

X1

\X + S q =ln \S + ,

z 1 =ln

The control coefficient of the flux through each enzyme in the unbranched pathway with respect to external substrate can be calculated from Eqs. (10) and (11) as

(6a)

1, 0

(6b)

s

0

C VS j #

we can express for the mass balances as a (log)linear model (Hatzimanikatis 6 Bailey, 1997): dz 1 V 1, 0 1 V ? q s & 2, 0 = 21 z 1 = dt X 1, 0 S X 1, 0

(7a)

i=2, 3, 4, (7b)

For a change in the concentration of the external substrate, the new steady state can be estimated by setting the right-hand side of the mass balance equations (Eqs. (7a) and (7b)) equal to zero, and solving for z i as a function of q S : zi =

Xi S ? 1S ?1 q  ln = i+1 ln . i+1 S =i X i, 0 =i S0

\ +

\ +

(8)

The concentration control coefficient of a metabolite X i with respect to a metabolic parameter, p, is defined within MCA as C ip #

d ln X i dz i # d ln p dq p

(9)

= i+1 = i

and quantifies the logarithmic deviation of the concentration of the metabolite X i for a unit change in the logarithmic deviation of the parameter p. From Eqs. (8) and (9), the concentration control coefficients of the metabolites in the unbranched pathway with respect to the concentration of the external substrate (or, response coefficient (Kacser and Burns, 1973)) can be calculated as ? 1S C iS = i+1 i=1, ..., 4. =i

K1 . K 1 +S

(14)

An incremental change in the external substrate concentration will result in changes in the metabolite concentrations. For a subsequent change in the external substrate the elasticities will have values different than the ones used in the previous change and they must be updated. The following iterative algorithm should be considered

(10)

ln (11)

\

? 1S | k X i | k+1 S | k+1 = i+1 ln Xi | k =i |k S |k

+

\

 ln(X i | k+1 )&ln(X i | k )=

where n is the number of the metabolites in the metabolic system. 77

File: DISTL2 010803 . By:JB . Date:03:09:98 . Time:11:02 LOP8M. V8.B. Page 01:01 Codes: 4919 Signs: 3360 . Length: 57 pic 1 pts, 240 mm

(13)

? 1S =

n

d ln V j = : = j C i +? pj , d ln p i=1 i p

K i+1 K i+1 +X i

and

The control coefficient of a metabolic flux with respect to a metabolic parameter is similarly defined (Heinrich 6 Rapoport, 1974; Reder, 1998) C Vp j #

(12)

One can observe from Eqs. (10) and (12) that the control coefficients do not depend on the reference concentration of the metabolites or on the reference value of the metabolic fluxes. Further, the control coefficients for any metabolic system without conserved metabolites are independent of the value of the reference concentration of the metabolites and they solely depend on the reference value of the metabolic fluxes, the stoichiometry and the values of the elasticities (Appendix III). The unbranched pathway is a special case for which the control coefficients are also independent of the values of the metabolic fluxes due to the simple stoichiometry that requires the steady-state fluxes through each enzyme to be equal. Therefore, knowledge of the elasticities and the distribution of the fluxes in a metabolic network suffice to calculate the control coefficients of any metabolic system. The change in the metabolite concentrations in response to changes in the external substrate can be calculated using Eq. 8. However, the value of the elasticities depends on the concentration of the corresponding metabolites. For the MichaelisMenten kinetics assumed here (Appendix I), and based on the definitions for the elasticities (Eqs. (4a) and (4b)), we have the following functional forms:

and dz i V i, 0 i V i+1, 0 i+1 = i&1 z i&1 & = zi , = dt X i, 0 X i, 0 i

d ln V j =? 1S . d ln S

+

? 1S | k [ln(S | k+1 )&ln(S | k )] = i+1 |k i (15)

metabolic engineering 1, 7587 (1999) article no. MT980108

= i+1 | k+1 == i+1 | k+ i i

Vassily Hatzimanikatis

= i+1 | k &= i+1 | k&1 i i (X i | k &X i | k&1 ) X i | k &X i | k&1

are a function of the elasticities only. Similar results can be obtained for other enzyme rate expressions (Hatzimanikatis, manuscript in preparation), as well as for typical expressions for enzyme regulation (see subsequent section.) Moreover, the above results allow us to formulate Eqs. (18a)(18c) as a system of ordinary differential equations. Using the definitions from Eqs. (6a) and (6b), the system can be written as follows

(Xi | k &Xi | k&1 )  0

=========O = i+1 | k+1 i $= i+1 | k+ i

= i+1 i X i

}

(X i | k &X i | k&1 )

(16)

k

and

? 1S | k+1 =? 1S | k +

? 1S | k &? 1S | k-1 (S | k &S | k&1 ) S | k &S | k-1

(S | k &S | k&1 )  0

========O ? 1S | k+1 $? 1S | k +

? 1S S

}

(S | k &S | k&1 ), k

(17)

d ln X i =

? =

d= i+1 = i

d ln S O dz i =

1 S i+1 i

? =

dq

d? 1S =

 ln = i+1 i d ln X i  ln X i

 ln ? 1S d ln S.  ln S

 ln ? 1S S ? 1S S = =& =? 1S &1.  ln S ? 1S S K 1 +S

(18b)

(18c)

(19a) (19b)

The above equations prove an interesting result: for a MichaelisMenten rate expression, the derivatives of the natural logarithm of the elasticities with respect to the natural logarithm of the metabolite concentration

(20c)

FIG. 1. 78

File: 741J 010804 . By:XX . Date:03:09:98 . Time:09:40 LOP8M. V8.B. Page 01:01 Codes: 5307 Signs: 3077 . Length: 57 pic 1 pts, 240 mm

du 1S =? 1S &1, dq S

(21a) (22a)

and the system of Eqs. (20a)(20c) can be treated as a system of ordinary differential equations with q s the independent variable, z i, w i+1 , and u 1S the dependent variables, i i+1 1 and = i | 0 and ? S | 0 the independent parameters. Integration of the system of Eqs. (20a)(20c) using zero as the initial condition for every independent variable, we obtain the percent change in the metabolite concentrations as a function of the percentage change in the manipulated metabolic parameter. The only information required for such integration is the one about the value of the elasticities at the reference steady state. The model described by the set of ordinary differential Eqs. (20a)(20c) will be called an NMCA model. We will apply the NMCA methodology to the unbranched pathway (Fig. 1) assuming the following, arbitrarily chosen, values for the reference elasticities: = i+1 | 0 =[0.5, 0.4118, 0.3333, 0.2308] and ? 1S | 0 =[0.4286]. i The kinetic parameters of the enzyme rate expressions that give these same reference elasticities are given in Appendix I. The results from the integration of the NMCA model

From Eqs. (13) and (14) we can calculate the partial derivatives for the elasticities: X i = i+1 Xi  ln = i+1 i i = i+1 =& == i+1 &1 i  ln X i =i X i K i+1 +X i

(20b)

? 1S =? 1S | 0 exp(u 1S )

(18a)

d? 1 S ? 1 S ? 1S S O 1S = 1 S S ? S ? S S S O d ln ? 1S =

dw i+1 ? 1S dw i+1 i == i+1 &1 O i =(= i+1 &1) i+1 i i dz i dq S =i

= i+1 == i+1 | 0 exp(w i+1 ) i i i

= i+1 d= i+1 X i = i+1 dX i i i dX i O i+1 = i+1 i X i =i =i X i X i O d ln = i+1 = i

(20a)

=ln(= i+1 = i+1 | 0 ) and u 1S =ln(? 1S ? 1S | 0 ), with where w i+1 i i i i+1 1 = i | 0 and ? S | 0 being the reference values for the elasticities. From these definitions, we can write the elasticities as functions of the new variables

where the subscript ``k'' denotes the value of the corresponding variable during the k th change in the external substrate concentration. For infinitesimally small changes in the external substrate, Eqs. (15)(17) can be written in differential form: 1 S i+1 i

? 1S dz i = i+1 dq s = i

Linear biochemical pathway.

metabolic engineering 1, 7587 (1999) article no. MT980108

Nonlinear Metabolic Control Analysis

FIG. 2. Steady-state responses of the concentrations of the intermediate metabolites, relative to their reference concentration, for relative changes in the concentration of the external substrate for the linear pathway (Fig. 1). Solid lines correspond to the NMCA model, dashed lines correspond to the MCA model ((log)linear model with constant elasticities (Eq. III.5)), and open circles correspond to the nonlinear model. (A) Responses of the metabolites X 1 , X 2 and X 3 . (B) Responses of the metabolite X 4 .

had the information about the maximum enzyme activity, V max. The parameters of the corresponding nonlinear model (Appendix I) suggest the possibility of the accumulation of metabolite X 4 since the V max for the enzyme that transforms X 4 is the smallest in the pathway.

(Eqs. (20a)(20c)) are presented in Figs. 2a and 2b, and demonstrate the striking performance of the NMCA model. The model captures the accumulation of metabolite X 4 as the concentration of the external metabolite approaches a critical value. This accumulation would be expected if we 79

File: 741J 010805 . By:XX . Date:03:09:98 . Time:09:40 LOP8M. V8.B. Page 01:01 Codes: 1801 Signs: 1154 . Length: 57 pic 1 pts, 240 mm

metabolic engineering 1, 7587 (1999) article no. MT980108

Vassily Hatzimanikatis

In many biological reaction networks featuring unbranched pathways, the reactions are reversible and the enzyme can act as a catalyst of the reverse reaction. The method can also be applied to the kinetics of reversible enzymatic reactions as discussed in Appendix II. The NMCA model can be also used in order to estimate the elasticities of the reference system since they are treated as parameters in a system of ordinary differential equations. There are many robust optimization algorithms that allow the estimation of the parameters of such a system from data about the dependent variables as functions of the independent variable. These data can be obtained from experimental measurements of the metabolic fluxes and of some of the metabolite concentrations for different values of the manipulated parameters. This type of experiments is relatively simple since they are steady-state experiments and the data required should be in the form of relative values. Moreover, advanced optimization methods allow parameter estimation in the presence of experimental uncertainties (ACSL, Mitchell and Gauthier, Concord, USA).

FIG. 3.

The NMCA model for this system (Eq. 23) can be derived following the same procedure as for the unbranched pathway. Considering the maximum activity of the protein kinase, V max, 1 , as the manipulated parameter, the NMCA model can be written as follows: dz ? 11 = 2 dq = C* &= 1C = CC* dw 2C* ? 11 =(= 2C* &1) 2 dq = C* &= 1C = CC*

EXAMPLE APPLICATIONS OF NMCA We will consider here two examples with highly nonlinear responses to the manipulation of the metabolic parameters: an interconvertible enzyme system, characteristic of signal transduction pathways, and a branched pathway.

dw 1C ? 11 =(= 1C &1) = CC* 2 dq = C* &= 1C = C C* dw CC* ? 11 , =(1&= CC* ) 2 dq = C* &= 1C = C C*

1. Interconvertible Enzyme System The elementary interconvertible enzyme system is presented in Fig. 3. It consists of a protein C that is phosphorylated into an active form, C*, by a protein kinase (enzyme E 1 ) and dephosphorylated by a phosphatase (enzyme E 2 ). These systems are very common to signal transduction pathways which switch abruptly from the inactive to active state when the activity of the protein kinase increases beyond a threshold value, relative to the activity of the phosphatase, under the action of an effector. Goldbeter and Koshland (1981) have introduced a modeling framework for these systems that considers Michaelis-Menten kinetics for the two enzymes: dC* C C* =V max, 1 &V max, 2 . dt K 1 +C K 2 +C*

Interconvertible enzyme system.

(24b) (24c) (24d)

where z and q are the logarithmic deviations of C* and V max, 1 , respectively, and w 1C and w 2C* are the logarithmic deviations of the elasticities, = 1C and = 2C* , respectively. We also need to introduce a new variable, = CC* , which arises from the conservation of the total enzyme and is defined as =C C* #

 ln C  ln(C T &C*) C* = =&  ln C*  ln C* C

(25)

which is the ratio of the amount of the active form of the interconvertible enzyme over the inactive form. In order to integrate the NMCA model we simply have to know the reference values of the elasticities (Eqs. (21a) and (21b)) as well as the reference value of the ratio C*C. We should emphasize here that we do not need to estimate the concentration of the interconvertible enzyme but the ratio of the concentrations. This ratio can be estimated from the amount of the enzyme in the various forms from relatively simple experiments whereas estimation of the absolute concentrations would be much more difficult due to uncertainties in the cell volume and the possible compartmentalization. Compartmentalization is a very common

(23)

They showed that if the total amount of the interconvertible enzyme, C T =C+C*, is time-invariant, and the Michaelis constants K 1 and K 2 are much smaller than C T , the interconvertible enzyme system is capable of generating threshold responses to changes in the relative value of the V max's (Fig. 4b). 80

File: 741J 010806 . By:XX . Date:03:09:98 . Time:09:40 LOP8M. V8.B. Page 01:01 Codes: 5200 Signs: 3984 . Length: 57 pic 1 pts, 240 mm

(24a)

metabolic engineering 1, 7587 (1999) article no. MT980108

Nonlinear Metabolic Control Analysis

FIG. 4. Steady-state responses of the interconvertible enzyme system (Fig. 3) to relative changes of the protein kinase (E 1 ). (A). Solid lines correspond to the NMCA model, dashed lines correspond to the MCA model ((log)linear model with constant elasticities (Eq. (III.5)), and open circles correspond to the nonlinear model. (B) Responses predicted from the NMCA model for different values of the reference elasticities [= 1C , = 2C* , = CC* ]: [0.0099, 0.9010, 0.0011] (solid line), [0.09, 0.9, 0.0011] (dashed line), [0.09, 0.80, 0.0011] (dotted line), and [0.009, 0.80, 0.0011] (dasheddotted line).

81

File: 741J 010807 . By:XX . Date:03:09:98 . Time:09:40 LOP8M. V8.B. Page 01:01 Codes: 1288 Signs: 690 . Length: 57 pic 1 pts, 240 mm

metabolic engineering 1, 7587 (1999) article no. MT980108

Vassily Hatzimanikatis

phenomenon in signal transduction pathways and in most of the cases the volume of the compartment where the pathway operates as well as the compartment itself are unknown. Next we will study the interconvertible enzyme system using an NMCA model. We will consider reference elasticities with values [= 1C , = 2C* , = CC* ]=[0.0099, 0.9010, 0.0011] and ? 11 =1. The results of the integration are presented in Fig. 4a. The parameters of the nonlinear system that would give the same values for these reference elasticities are [V max, 1 , V max, 2 , K 1 , K 2 ]=[0.1, 1, 0.01, 0.01]. The performance of the NMCA model is excellent in capturing the highly nonlinear response of the system and predicting the threshold value of the kinase activity. Assuming different reference values for the elasticities, we can generate different responses (Fig. 4b). This observation suggests that we can estimate the reference values for the elasticities using experimental data from the steady-state responses of the system to manipulations of the activities of the protein kinase (E 1 ) andor phosphatase (E 2 ).

FIG. 5.

Branched biochemical pathway. Dashed lines denote inhibition.

= VS =

K(I) K(I)+S

= VI = &

K(I) = K = &= VS = K I , K(I)+S I

=K I =

Branched pathways are very common in metabolic systems. In many cases they feature regulatory characteristics such as feedback inhibition. The pathway considered here is presented in Fig. 5 and consists of three intermediate metabolites, an external substrate, and five enzymes. Two of the enzymes are feedback inhibited by their products. There are various kinetic expressions that can be derived based on the mode of action of the inhibitor on an enzyme. The most common mode of action is the competitive inhibition where the product of the enzyme, due to its structural similarity with the substrate, binds to the catalytic site and thus inhibits substrate binding. Assuming Michaelis Menten kinetics for the action of the enzyme on the substrate and competitive inhibition for the action of the inhibitor on the enzyme, we can derive the following rate expression S

\

K 0 1+

I +S Ki

+

=V max

S , K(I)+S

d ln K(I) I = d ln I K i +I

(29)

is the elasticity of the Michaelis constant with respect to V K inhibitor. When the elasticities = V S and = I are known, = I can K V V be calculated from Eq. (28): = I =&= I = S . In order to formulate a NMCA model we have to introduce two equations that will describe the dependency of the logarithmic deviation of the elasticities on the logarithmic deviations of the substrate and the inhibitor. These descriptions can be derived from Eqs. (27)(29):

(26)

 ln = VS S =& == VS &1  ln S K(I)+S

(30a)

S  ln = VS = = K =(1&= VS ) = KI  ln I K(I)+S I

(30b)

Ki  ln = KI = =1&= }I .  ln I K i +I

(30c)

The above equations lead to an interesting conclusion: the logarithmic deviations of the elasticities with respect to the logarithmic deviations of the substrate and the inhibitor are functions of the elasticities only, when competitive inhibition is considered. The NMCA model for the branched pathway studied here can be formulated based on the Eqs. (27)(30c), considering the external substrate as the manipulated variable

where S and I are the concentrations of the substrate and the inhibitor, respectively, and K(I) denotes the dependency of the apparent Michaelis constant on the concentration of the inhibitor. From Eq. (26) the elasticities of the enzyme with respect to the substrate and to the inhibitor can be derived 82

File: 741J 010808 . By:XX . Date:03:09:98 . Time:09:40 LOP8M. V8.B. Page 01:01 Codes: 4780 Signs: 3512 . Length: 57 pic 1 pts, 240 mm

(28)

where

2. BRANCHED PATHWAY

V=V max

(27)

metabolic engineering 1, 7587 (1999) article no. MT980108

Nonlinear Metabolic Control Analysis

4 2 K2 1 (V 1 +V 2 ) ? 0S (= 32 += 11 = K dz 1 2 )(= 3 += 1 = 3 ) = 1 3 4 2 K2 2 4 3 dq s V 1 = 1 = 2 (= 3 += 1 = 3 )+V 2 = 1 = 3 (= 2 += 11 = K2 )

(31.a)

dz 2 dz 1 =1 = 3 11 K1 dq S = 2 += 1 = 2 dq S

(31.b)

dz 3 = 21 dz 1 = += 21 = K3 2 dq S = 43 dq S

(31.c)

1

dz dw 11 dz 2 =(= 11 &1) 1 +(1&= 11 ) = K2 dq S dq S dq S 1

dw K2 1 dz 2 =(1&= K2 1 ) dq S dq S dz 3 dz 1 dw 21 2 =(= 21 &1) +(1&= 21 ) = K 3 dq S dq S dq S

(31.d) (31.e) (31.f)

dz dw K3 2 =(1&= K3 2 ) 3 dq S dq S

(31.g)

dz 2 dw 32 =(= 32 &1) dq S dq S

(31.h)

dz 3 dw 43 =(= 43 &1) dq S dq S

(31.i)

du 0S =(? 01 &1) dq S

(31.k)

dz 1 dz 2 dr 1 == 11 &= 11 = K2 1 dq S dq S dq S

(31.1)

dz dz dr 2 == 21 1 &= K3 2 3 , dq S dq S dq S

(31.m)

FIG. 6. Steady-state responses of the difference (V 1 &V 2 ) in the branched pathway (Fig. 5) for relative changes in the concentration of the external substrate. Solid lines correspond to the NMCA model, dashed lines correspond to the MCA model ((log)linear model with constant elasticities and rates used in the scaling of the stoichiometricmatrix (Eq. (III.5)), and open circles correspond to the nonlinear model.

and

[V 1 , V 2 ]=[0.0294, 0.0003]. The corresponding nonlinear system is presented in Appendix III. The inhibition present in the system increases the nonlinearity of the system. However, the NMCA model performs extremely well in capturing the steady-state responses of the system as shown in Fig. 6, where the response of the difference (V 1 &V 2 ) has been studied. Using exactly the same information available for a MCA model, and by simply assuming the form of the rate expressions, the NMCA model can predict responses that the MCA model would have failed to capture even approximately.

where z i and q S are the logarithmic deviations of the metabolite X i and the external substrate, respectively, w and u are the logarithmic deviations of the corresponding elasticities, and r j are the logarithmic deviations of the reaction rate V j . At steady state, V 0 =V 1 +V 2 , V 1 =V 3 , and V 2 =V 4 , and therefore, calculation of only two of the rates is adequate for the description of the system. In order to integrate the NMCA model, the only information required is the value of the reference elasticities and the value of the rates at the reference state. We studied the performance of the NMCA model considering the reference state as described by the following elasticity and reaction rate values:

SUMMARY The development of a novel framework for modeling and analysis of biochemical reaction networks has been introduced. Based on information provided from MCA and assuming common rate expressions for the reaction steps in the network, an alternative model for the system has been

[= 11 , = 12 , = 21 , = 23 , = 32 , = 43 , ? 0S ] =[0.9706, &0.0027, 0.9997, &10 &6, 0.9804, 0.99980.9901] 83

File: 741J 010809 . By:XX . Date:03:09:98 . Time:09:41 LOP8M. V8.B. Page 01:01 Codes: 4498 Signs: 2600 . Length: 57 pic 1 pts, 240 mm

metabolic engineering 1, 7587 (1999) article no. MT980108

Vassily Hatzimanikatis

with [V max, 2 , V max, 3 , V max, 4 , V max, 5 ]=[2, 1.7, 1.5, 1.3], [K 2 , K 3 , K 4 , K 5 ]=[0.3333, 0.6667, 0.6429, 0.1875], and

constructed. The model, called an NMCA model, has been set up in the form of a set of ordinary differential equations which, when solved simultaneously, predicts the responses of the metabolic system to large changes in its parameters. Moreover, the NMCA model allows the calculation of elasticities and, consecutively, control coefficients as continuous functions of the metabolic parameters. Therefore, it offers the capability of predicting the changes in the distribution of the ``control'' throughout a metabolic system associated with large changes in its parameters. The NMCA model has been tested using three examples that display nonlinear responses and the results demonstrated its excellent performance. The rate expressions considered in the examples were limited to MichaelisMenten kinetics for single substrate and competitive inhibition. However, the methodology is applicable to various classes of rate expressions for single and multiple substrates, and for various types of enzyme regulation (Hatzimanikatis, manuscript in preparation.) One of the most known and cited works on metabolic control analysis for large changes in the metabolic parameters is the work by Small and Kacser (1993). There are significant differences between the work presented here and their work. Small and Kacser considered linear kinetics whereas the NMCA framework can be used for systems featuring certain forms of nonlinear kinetics. Their method requires only knowledge of the control coefficients and provides an estimate of the responses of the reaction rates and the concentration of the metabolites of the system. The NMCA methodology requires knowledge of the elasticities and, in the case of reversible reactions, the thermodynamic properties of the system. In addition, it also allows accurate prediction of the responses of the reaction rates and the concentration of the metabolites of the system for large changes in the metabolic parameters. Quantitative metabolic engineering is critically dependent on MCA quantities such as elasticities, reaction rates, and control coefficients. NMCA offers a new methodology for estimating this information from simple, steady-stateresponse experiments.

V=1.75

S 1.5+S

for the biotransformation of the external substrate. At the reference steady state we considered S=2, which resulted in [X 1 , X 2 , X 3 , X 4 ]=[0.3333, 0.9524, 1.2858, 0.6250] for the steady-state metabolite concentrations. APPENDIX II Consider an unbranched metabolic pathway consisting of n steps X 0  S 1  S 2  S 3 } } } S n&1  X n . The nonlinear rate expression for the transformation of Si to S j can be written as follows vi =

(V max, i K m, i )(S i &S j K ij ) , 1+S i K m, i +S j K m, j

(II.1)

where v i is the net rate (positive for the formation of S j ), V max, i is the maximum rate of S j formation, K m, i , and K m, j are the K m constants for S i and S j , respectively, and K ij is the equilibrium constant for the reaction. The above rate expression can also be written as the difference between the forward and backward reaction rates v i =v if &v bi =

f V max, i S i K m, i 1+S i K m, i +S j K m, j

&

V bmax, i S j K m, j , 1+S i K m, i +S j K m, j

(II.2)

where V fmax, i and V bmax, i are the maximum forward and backward rates, respectively (note that V fmax, i =V max, i in Eq. (II.1)). The elasticities of the forward and backward reaction rates with respect to S i and S j can be written as:

APPENDIX I

= bi #

 ln v bi S i K m, i =&  ln S i 1+S i K m, i +S j K m, j

(II.3a)

The rate expressions of the enzymes in the linear pathway (Fig. 1) were considered to follow MichaelisMenten kinetics

= if #

 ln v if 1+S j K m, j = =1+= bi  ln S i 1+S i K m, i +S j K m, j

(II.3b)

= jf #

 ln v if S j K m, j =&  ln S j 1+S i K m, j +S j K m, j

(II.3c)

= bj #

 ln v bi 1+S i K m, i = =1+= fj . (II.3d)  ln S j 1+S i K m, i +S j K m, j

V=V max, i

Xi K i +X i

i=2, 3, 4, 5,

(I.1) 84

File: DISTL2 010810 . By:JB . Date:03:09:98 . Time:11:02 LOP8M. V8.B. Page 01:01 Codes: 5559 Signs: 3892 . Length: 57 pic 1 pts, 240 mm

metabolic engineering 1, 7587 (1999) article no. MT980108

Nonlinear Metabolic Control Analysis

The elasticities of the net reaction rate, v I , with respect to Si and S j can be written as functions of the above elasticities (Eqs. (II.3a)(II.3d)), and the equilibrium constant, K ij : = vi #

 ln v i = b K ij = b i += bi  ln S i = i K ij &= jf

(II.4a)

= vj #

&= jf  ln v i = b += fj .  ln S j = i K ij &= jf

(II.4b)

The above equations (II.4a)(II.4b) and (II.5a)(II.5h) suggest that knowledge of the thermodynamic properties of the system (K ij ) and an estimate of the values of the reference elasticities allow the application of the NMCA framework on biochemical systems that feature reversible reactions. Alternatively, the employment of the above system of equations can be used to estimate reference elasticities from experimental data as discussed in the main body of the article.

The equilibrium constant, K ij , can be estimated from the thermodynamic properties of the system, and, if the values of elasticities = vi and = vj can be estimated, the solution of the system of Eqs. (II.4a) and (II.4b) will provide with the values of the elasticities = if , = bi, = fj , and = bj . In order to apply the NMCA framework to this system we have to show that the derivatives of the logarithmic deviations of the elasticities with respect to the logarithmic deviations of the metabolite concentrations are functions of the elasticities. It can be shown that

APPENDIX III The general formulation of the (log)linear model has been presented before (Hatzimanikatis 6 Bailey, 1997). We will consider here a metabolic system consisting of n metabolites with no conserved species among them and m enzymatically catalyzed reactions. For simplicity, we consider that those m reactions are the only processes in the system by which concentrations of metabolites change. Under the assumption of negligible concentration gradients in the cell, the mass balances for the metabolites may be written as

d ln = vi K ij = bi (1+2= bi ) = d ln S i (K ij = bi &= fj )(K ij = bi &= fj &1) &

(K ij = bi &= fj ) = bi (K ij = bi &= fj &1)

dx =Nv(x; p), dt (II.5a)

where x is the n-dimensional metabolite concentration vector, v is the m-dimensional reaction rate vector, p is the s-dimensional manipulated metabolic parameter vector (e.g., external substrate concentration, enzyme concentration), and N is the stoichiometric matrix. Let x 0 be a steady state of the system for a given set of the parameter values, p 0, and let x 0 be non-zero with positive elements. The linearized system of the mass balance equations around this steady state can be written as follows:

K ij = bi (1+2= jf ) d ln = vi =& d ln S j (K ij = bi &= fj )(K ij = bi &= fj &1) + d ln = vj d ln S i

=&

d ln S j

=

(K ij = bi &= fj &1)

(K ij = bi &= fj )(K ij +K ij = bi &= fj &1) (K ij = bi &= fj )(1+= bi ) (K ij +K ij = bi &= fj &1)

(K ij = bi &= fj ) = fj (K ij +K ij = bi &= fj &1)

b i

d ln = =&= fj d ln S j d ln S i d ln = fj d ln S j

d(x&x 0 ) dt v =N x

(K ij = bi &= fj )(K ij +K ij = bi &= fj &1)

d ln = bi =1+= bi d ln S i

d ln =

(II.5c)

K ij = fj (1+2= fj )

+

f j

(II.5b)

K ij = fj (1+2= bi )

+ d ln = vi

(K ij = bi &= fj )(1+= fj )

(III.1)

}

(x&x 0 )+N x0 , p0

v p

(II.5e) X &1 0

d(x&x 0 ) v &1 =X &1 0 NV 0 V 0 dt x

(II.5f) &1 +X &1 0 NV 0 V 0

(II.5g)

=1+= fj .

(II.5h)

(III.2)

}

X 0 X &1 0 (x&x 0 ) x0 , p0

v p

}

P 0 P &1 0 (p&p 0 ), x 0 , p0

(III.3) where X 0 is a diagonal n_n matrix with X 0, ii , P 0 is a diagonal s_s matrix with P 0, ii =P i, 0 , and V 0 is a diagonal 85

File: DISTL2 010811 . By:JB . Date:03:09:98 . Time:11:02 LOP8M. V8.B. Page 01:01 Codes: 6065 Signs: 3339 . Length: 57 pic 1 pts, 240 mm

(p&p 0 ). x 0 , p0

The above linearized equation can be scaled with respect to the reference values of the metabolite concentrations and the metabolic parameter values

(II.5d)

=&= bi

}

metabolic engineering 1, 7587 (1999) article no. MT980108

Vassily Hatzimanikatis

s_s matrix with V 0, ii =V i, 0 . The approximation introduced in Eq. (5) suggests that X &1 0 (x&x 0 )$z and P &1 (p&p )$q with z =ln(x x ) and q i =ln( p i p i, 0 ). 0 i i i, 0 0 Moreover, we can define the matrices

=#V &1 0

v x

}

and

X0

6#V &1 0

x 0 , p0

v p

}

species in the metabolite pool, the concentration control coefficients and the flux control coefficients can be calculated without explicit knowledge of the metabolite concentrations. The additional knowledge of the concentration of the metabolites is required for the examination of the local stability characteristics of the system by calculating the eigenvalues of the matrix X &1 0 NV 0 =. In the case of a system with conserved pools of metabolites it can be shown that the additional knowledge of the relative distribution of the metabolites in the conserved pool is required (see ``Interconvertible Enzyme'' subsection in the Example Applications).

P0 x0 , p0

with elements X

=ji # V i, 0

j, 0

V j X i

}

and

6 jk #

x 0 , p0

P k, 0 V j V j, 0 P k

}

x0 , p0

APPENDIX IV

which by definition (Eqs. (5.a) and (5.b)) are the elasticities of the enzyme catalyzing reaction step j with respect to metabolite X i and to parameter P j , respectively. Employing the above defintions we can write for the (log)linear model the following equation: dz &1 =X &1 0 NV 0 =z+X 0 NV 0 6q. dt

The branched pathway considered is presented in Fig. 5. The external substrate, S, was the manipulated variable. The kinetics used are the following (numbered after the numbering in Fig. 5):

(III.4)

V 0 =3

The steady-state responses of the metabolic system to changes in the metabolic parameters (or: the logarithmic deviations of the concentration of the metabolites for given values of the logarithmic deviations in the metabolic parameters) can be calculated from the steady-state solution of the (log)linear model z= &(NV 0 =) &1 NV 0 6q O z=C zq q,

V1=

V2=

V 4 =2 (III.6)

dV j d ln V j = dq k d ln P k

(IV.3)

+

X2 5+X 2

(IV.4)

X3 . 10+X 3

(IV.5)

For the reference state we set S=0.01 which results in the steady-state metabolite concentrations: [X 1 , X, X 3 ]= [0.0003, 0.999, 0.0017]. (III.7)

ACKNOWLEDGMENTS

and can be calculated from the elasticity matrices and the concentration control coefficient matrix: C vq = =C zq +6.

(IV.2)

+

X1 5 1 1+ +X 1 X3

V 3 =1.5

The flux control coefficient matrix is defined as (C vq ) jk =

X1 1 0.01 1+ +X 1 X2

\

(III.5)

where C is, by the MCA definition, the concentration control coefficient matrix with elements: dz i d ln X i (C ) = = . dq k d ln P k

(IV.1)

\

z q

z q ik

S 1+S

The author is grateful to Dr. S. Wilkinson for his suggestions, that have significantly contributed to this work, and to Drs W. D. Smith, W. Foster, and J. S. Schwaber.

(III.8) REFERENCES

Equations (III.5) and (III.8) suggest that if the stoichiometry and the elasticities of the metabolic system described by Eq. (III.1) are known, and there are no conserved

Bailey, J. E. (1991). Toward a science of metabolic engineering. Science 252, 16681675. 86

File: DISTL2 010812 . By:JB . Date:03:09:98 . Time:11:02 LOP8M. V8.B. Page 01:01 Codes: 5079 Signs: 2885 . Length: 57 pic 1 pts, 240 mm

metabolic engineering 1, 7587 (1999) article no. MT980108

Nonlinear Metabolic Control Analysis

Cornish-Bowden, A., and Cardenas, M. L. (Eds.) (1990). ``Control of Metabolic Processes,'' NATO ASI Series A: Life Sciences Vol. 190, Plenum Press, New York. Fell, D. A. (1992). Metabolic control analysis: A survey of its theoretical and experimental development. Biochem. J. 286, 313330. Fell, D. (1997). ``Understanding the Control of Metabolism,'' Portland Press, London. Goldbeter, A., and Koshland, D. E. (1981). An amplified sensitivity arising from covalent modification in biological systems. Proc. Natl. Acad. Sci. USA 78(11), 68406844. Hatzimanikatis, V., and Bailey, J. E. (1997). Effects of spatiotemporal variations on metabolic control: Approximate analysis using (log)linear kinetic models. Biotechol. Bioeng. 54(2), 91104. Heinrich, R., and Rapoport, T. A. (1974). A linear steady state Tetatment of Ezeymatic chains. Eur. J. Biochem. 42, 8995.

Hofer, T., and Heinrich, R. (1993), A second order approach to metabolic control analysis. J. Theor. Biol. 164, 85102. Kacser, H., and Acerenza, L. (1993). A universal method for achieving increases in metabolite production. Eur. J. Biochem. 216, 361367. Kacser, H., and Burns, J. A. (1973). The control of flux. Symp. Soc. Exp. Biol. 27, 65104. Reder, C. (1988). Metabolic control theory: A structural approach. J. Theor. Biol. 135, 175201. Small, J. R., and Kacser, H. (1993). Response of metabolic systems to large changes in enzyme activities and effectors. 1. The linear treatment of unbranched chains. Eur. J. Biochem. 213, 613624. Small, J. R., and Kacser, H. (1993). Response of metabolic systems to large changes in enzyme activities and effectors. 2. The linear treatment of branched pathways and metabolite concentrations. Assement of the general non-linear case. Eur. J. Biochem. 213, 625640.

87

File: DISTL2 010813 . By:JB . Date:03:09:98 . Time:11:02 LOP8M. V8.B. Page 01:01 Codes: 5253 Signs: 1905 . Length: 57 pic 1 pts, 240 mm