- Email: [email protected]
Control analysis of metabolic systems involving quasi-equilibrium reactions
Biochimica et Biophysica Acta 1379 Ž1998. 337–352
Control analysis of metabolic systems involving quasi-equilibrium reactions Boris N. Kholodenko a
a,b
, Stefan Schuster c , Jordi Garcia a , Hans V. Westerhoff d , Marta Cascante a,)
Departament de Bioquımica i Fisiologia, Facultat de Quımica, UniÕersitat de Barcelona, Martı´ i Franques, ´ ´ ` 1, 08028 Barcelona, Spain b A.N. Belozersky Institute of Physico-Chemical Biology, Moscow State UniÕersity, 119899 Moscow, Russian Federation c Max Delbruck ¨ Center for Molecular Medicine, 13125 Berlin-Buch, Germany d BioCentrum, Department of Microbial Physiology, Faculty of Biology, Free UniÕersity, De Boelelaan 1087, NL-1081 HV Amsterdam, Netherlands Received 17 April 1997; revised 3 September 1997; accepted 5 September 1997
Abstract Reactions for which the rates are extremely sensitive to changes in the concentrations of variable metabolite concentrations contribute little to the control of biochemical reaction networks. Yet they do interfere with the calculation of the system’s behaviour, both in terms of numerical integration of the rate equations and in terms of the analysis of metabolic control. We here present a way to solve this problem systematically for systems with time hierarchies. We identify the fast reactions and fast metabolites, group them apart from the other Ž‘‘slow’’. reactions and metabolites, and then apply the appropriate quasi-equilibrium condition for the fast subsystem. This then makes it possible to eliminate the fast reactions and their elasticity coefficients from the calculations, allowing the calculation of the control coefficients of the slow reactions in terms of the elasticity coefficients of the slow reactions. As expected, the elasticity coefficients of the fast reactions drop out of the calculations, and they are irrelevant for control at the time resolution of the steady state of the slow reactions. The analysis, when applied iteratively, is expected to be particularly valuable for the control analysis of living cells, where a time hierarchy exists, the fastest being at the level of enzyme kinetics and the slowest at gene expression. q 1998 Elsevier Science B.V. Keywords: Control coefficient; Metabolic control analysis; Quasi-equilibrium approximation; Time hierarchy
1. Introduction The pronounced hierarchy of the characteristic times of enzyme reactions is a fundamental property of cellular metabolism w1,2x. Due to differences in the enzyme activities it is often possible to subdivide the reactions into
)
Corresponding author. Fax: q34 3 402 1219.
0304-4165r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 0 4 - 4 1 6 5 Ž 9 7 . 0 0 1 1 4 - 1
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
338
two classes: Ž very. fast reactions and slow reactions. Experimental measurements often show that at the operating state the ratios of substrate and product concentrations Žmass action ratios. of the fast reactions are close to the corresponding equilibrium constants Žsee, e.g., the experimental determination of metabolite concentrations of glycolysis in erythrocytes w3x.. The sensitivities of relevant variables in steady states of chemical and biochemical reaction systems Ž such as steady-state fluxes and concentrations of metabolites. to changes in the activities of particular reactions can be studied by the use of a theoretical framework known as metabolic control analysis Žfor a review, see w4x.. A general way of defining control coefficients is w5–8x: CiJk s
E ln < Jk
i , j,k s 1,2,3 . . .
Ž 1.1 .
Here, Õi , Jk and x j denote the velocity of reaction i, the flux through reaction k at a stationary state under consideration, and the concentration of a substance X j , respectively. pi is a parameter affecting only reaction i directly, i.e.,
E ln < Õi < E pi
/ 0,
E ln < Õi < E pk
s 0 for any k / i.
Ž 1.2 .
The derivatives of the logarithmic rates ln < Õi < with respect to the parameters in Eqs. Ž1.1. and Ž1.2. are taken as if reaction i proceeded in isolation from the remaining network. When there is a one-to-one correspondence between enzymes and reactions and each reaction rate is proportional to the concentration of the corresponding enzyme, definition Ž 1.1. is equivalent to the definition of control coefficients given by Kacser and Burns w9x in terms of enzyme concentrations. The kinetic properties of the individual enzymes are expressed in terms of their elasticity coefficients Ž e ., which describe the sensitivities of the reaction rates to changes in the metabolite concentrations w5,9x:
´ ki s d ln < Õi
Ž 1.3 .
In contrast to the definition of the control coefficients Ž Eq. Ž 1.1.. , no evolution of the entire pathway to a steady state is presupposed in Eq. Ž 1.3. . Comparing the effects of small changes in the activities of Žvery. fast reactions and of slow reactions on the steady-state variables one can usually neglect the effects of the former in comparison to the latter. Therefore, as has often been stated, e.g., w4,5,10x, the control coefficients of those reactions which have mass action ratios close to the equilibrium constants are virtually equal to zero. This is in line with a fundamental property of fast reactions: the elasticity coefficients of these reactions with respect to the substrates and products assume very high positive and negative values, respectively, since at the system’s steady states these reactions are sufficiently close to thermodynamic equilibrium Ž although the flux through these reactions is non-zero. w11,12x. This feature hampers the determination of control coefficients Ž quantitative indicators of ‘‘global’’ control properties of a pathway. in terms of the elasticity coefficients Žreflecting local kinetic properties of an enzyme reaction. . It is therefore useful to elaborate methods for determining control properties using the elasticities of the slow reactions only. The quasi-equilibrium approximations have been known for a long time Ž see, e.g., w13x, and references therein.. These approximations Žalso referred to as rapid-equilibrium approximations. have been employed in the modelling of the dynamics of various metabolic systems Ž see, e.g., w1,14–19x. . Although a metabolic pathway as a whole may be far from equilibrium, the quasi-equilibrium approximation allows one to consider a subsystem comprising fast reactions to be in a Ž quasi-. equilibrium state. This makes it possible to avoid a kinetic analysis of the fast reactions, i.e., to take into consideration only their equilibrium constants instead of considering their rates.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
339
In the general case the conditions under which one can reduce the number of dynamic variables in the system using the quasi-steady-state approximation were formulated for the first time by Russian mathematician Tikhonov w20–22x. The way of using this general formalism to validate the rapid-equilibrium approximation in the case of biochemical pathways was elaborated by several authors, see, e.g., w23–26x. Until now, the control properties of metabolic pathways with quasi-equilibrium steps have not yet been studied in a systematic way. The goal of this paper is to develop, in the framework of quasi-equilibrium approximation, a general method for determining the control coefficients in terms of the elasticity coefficients of the slow reactions only. In Section 2, the method is illustrated by way of an example. Section 3 is aimed at presenting the rapid-equilibrium approximation in a matrix notation so as to provide the basis for the control analysis set out in Section 4. As will become clear from below, when the system involves conservation relations saying that linear combinations of concentrations remain constant in time, these relations have to be taken into account in the control analysis of networks with rapid-equilibrium subsystems. This is particularly complicated when such relations involve substances participating in fast reactions. To begin with a clear presentation which shows the essential ideas and to avoid an unnecessarily complicated analysis, we here exclude conservation relations of this type. We allow for conservation relationships among the substances participating in slow reactions only. The general case of systems with arbitrary conserved-moiety structures will be analysed in a sequel paper.
2. A simple example We shall begin with illustrating the principles of the method using the exemplifying reaction system shown in Scheme 1. This is an unbranched metabolic pathway consisting of three sequential monomolecular reactions. They may represent, e.g., the reactions catalysed by hexokinase, phosphoglucoisomerase, and phosphofructokinase involved in glycolysis. The symbols Õ s and Õ f stand for the rates of ‘‘slow’’ and ‘‘fast’’ reactions, respectively. We here assume that direct interactions between the enzymes are absent w27x. The concentrations of the initial substrate, S, and of the final product, P, are taken to be constant. The differential equations for the time evolution of the concentrations x 1 and x 2 read d x 1rd t s Õ 1s y Õ 2f ,
Ž 2.1a .
d x 2rd t s Õ 2f y Õ 3s .
Ž 2.1b .
At any steady state of the system the time derivatives are equal to zero. Therefore, all three reactions carry the same flux Ž J1 ., and one may write s f s Õ 1
Ž 2.2 .
Here and below, the subscript st.st. specifies steady-state conditions. System Ž 2.2. consisting of three equations allows one to determine the concentrations x 1 and x 2 and the flux J1 through the pathway of Scheme 1 at steady state, provided that the rate laws are known. The control properties of this metabolic pathway are described by the control coefficients Ž CiJ1 and Cix k . of the enzymes over the pathway flux and metabolite concentrations. The control coefficients of the enzymes depend
Scheme 1.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
340
on their kinetic properties and on pathway structure. One way of expressing this dependence is by using the summation and connectivity theorems of metabolic control analysis w6,9,12,28,29x. They can be written in a concise and elegant form of a matrix equation w6,30–33x. For the pathway of Scheme 1 this matrix equation can be written as follows:
C1J1 C1x 1 C1x 2
1 C3J1 C3x 1 P 1 C3x 2 1
C2J1 C2x 1 C2x 2
0
s
y´ 1Õ1 f
y´ 1Õ 2 y´
Õ 3s 1
s
y´ 2Õ1
1 f y´ 2Õ 2 s 0 s 0 y´ Õ3 2
0
0 1 0
0 0 . 1
0
Ž 2.3 .
The products of the rows of the control matrix ŽC. and the first column of the elasticity matrix ŽE. give all the summation theorems for the flux and concentration control coefficients for the pathway of Scheme 1. The products of the rows of matrix C and the next two columns of matrix E give all the connectivity theorems. Writing these theorems in the form of a matrix equation Ž Eq. Ž 2.3.. allows one to calculate the control coefficients by inversion of the elasticity matrix E. Now we assume that the activity of enzyme 2 is so high that reaction 2 can be considered to be, after a very short transient period, at a quasi-equilibrium state Ž which is actually the case for phosphoglucoisomerase in the glycolytic pathway w3x, for example.. For the time evolution of the metabolites ŽEqs. Ž2.1a. and Ž2.1b.. this implies that the ratio of the concentrations of x 2 and x 1 approaches the equilibrium constant of reaction 2: x 2rx 1 s K 2eq .
Ž 2.4 .
According to this equation only one concentration out of the two concentrations x 1 and x 2 can now be considered independent. Without loss of generality, we shall choose x 2 as the independent variable. Summing Eqs. Ž2.1a. and Ž2.1b. , one arrives at an equation describing the time evolution of a so-called ‘‘slow variable’’ or ‘‘pool variable’’ w1,2x. Only slow enzyme rates enter the new equation: d Ž x 1 q x 2 . rd t s d yrd t s Õ 1s y Õ 3s .
Ž 2.5 .
The pool variable, y s x 1 q x 2 , is the sum of the concentrations of the metabolites participating in the fast reaction. Since the requirements of Tikhonov’s theorem are fulfilled Ž see, e.g., w24,26x. one can substitute Eq. Ž2.4. into Eq. Ž2.5.. The obtained reduced equation system consists of one differential equation only and describes the time evolution of the independent concentration Ž x 2 .: d x 2rd t s
K 2eq 1 q K 2eq
P Ž Õ 1s y Õ 3s . .
Ž 2.6 .
At steady state the rates of the slow reactions 1 and 3 are equal Ž cf. Eq. Ž 2.2.. . The control coefficients of the enzymes of the slow reactions can be determined by differentiation of their rates in Eq. Ž2.2. Žthis method of calculating the control coefficients goes back to the work of Crabtree and Newsholme w34x. . Expressing in the rates Õis < st. st. in Eq. Ž2.2., the dependent concentration x 1 via the independent concentration x 2 according to Eq. Ž2.4., one obtains, by the chain rule of differentiation,
E ln Jj E pi
s
ž
E ln Õjs E ln x 1 E ln x 1 E ln x 2
q
E ln Õjs E ln x 2
/
E ln x 2 E pi
q
E ln Õjs E pi
,
i , j s 1,3.
Ž 2.7 .
Here and in the following, for simplicity’s sake, we omit the bars indicating the absolute value of reaction rates. From Eq. Ž2.4. it follows that d ln x 1rd ln x 2 s 1. Therefore, dividing Eq. Ž2.7. by E ln ÕisrE pi gives, due to definition Ž 1.1. and Eq. Ž1.2., CiJj y Ž ´ 1j q ´ 2j . P Cix 2 s d ji ,
i s 1,3, j s 1,3,
Ž 2.8 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
341
where d i j is the Kronecker symbol:
d ji s 1 if i s j, otherwise d ji s 0.
Ž 2.9 .
The four equations contained in Eq. Ž2.8. Žfor all i, j . can be written in the form of the following matrix equation:
1
y Ž ´ 11 q ´ 21 .
1
y Ž ´ 13 q ´ 23 .
0ž P
C1J1 C1x 2
C3J1 1 s 0 C3x 2
/
ž
0 . 1
/
Ž 2.10 .
To compare this equation to Eq. Ž2.3. it is convenient to permute the matrices E and C Žthis is always possible since their product is the identity matrix. :
ž
C1J1 C1x 2
1 C3J1 x2 P C3 1
/
y Ž ´ 11 q ´ 21 . y Ž ´ 13 q ´ 23 .
0
s
ž
1 0
0 . 1
/
Ž 2.11 .
This equation contains the summation and connectivity theorems for the control coefficients calculated in the framework of quasi-equilibrium approximation. The products of the rows of matrix C and the first column of matrix E show that the summation theorems derived for the original system still apply provided that the control coefficients of the fast reaction 2 are assumed to be equal to zero. The products of the rows of matrix C and the second column of matrix E represent the ‘‘new’’ connectivity relationships which correspond to simultaneous equal relative changes in the concentrations x 1 and x 2 so that their equilibrium ratio remains constant Ž cf. the perturbation method in w12x.. Note that the elasticity coefficients of the fast reaction which tend to plus and minus infinity, respectively, do not enter the new relationships. Importantly, one can regard Eqs. Ž 2.10. and Ž2.11. as corresponding to a simplified pathway involving reactions 1 and 3 only. As has been derived earlier by Fell and Sauro w35x, the effective elasticities of reactions 1 and 3 are the sums of the elasticities of these reactions with respect to x 1 and x 2 .
3. General analysis: Fundamentals We shall now consider a metabolic pathway of any complexity, involving n reactions and m metabolites. We shall assume that due to differences in the activities of enzymes or in the rate constants of uncatalysed reactions, the rates Ž Õif . of some very fast reactions Ž f in number. greatly exceed the rates Ž Õjs . of the n y f remaining, slow reactions Ž for any metabolite concentrations except a small region, see below.: < Õjs < < < Õif < for all i , j Ž i s 1, . . . , f ; j s f q 1, . . . ,n . .
Ž 3.1 .
Note that the numbering in Scheme 1 differs from this general numbering for the sake of clarity. As the reaction rates depend on concentrations, one should indicate in which region of the concentration space inequality Ž3.1. is to hold. Since the fast reactions are assumed to relax very rapidly to quasi-equilibrium states, where Õif f 0, we restrict the validity of relation Ž3.1. to a domain which does not include the region defined by Õif s 0 Žthe ‘‘slow submanifold’’. and a small neighbourhood around it. This neighbourhood includes the region where the steady-state Eq. Ž 2.2. holds. In other words, due to large differences in the kinetic parameters, some reactions have much larger rates than others except in situations where the concentrations adapt to these differences in such a way that a steady state is nevertheless possible. Let q be the number of metabolites which participate in the fast reactions and possibly in the slow reactions. Their concentrations will be designated by x if , i s 1, . . . , q. The remaining m y q metabolites Ž m is the number
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
342
of all metabolites. participate only in the slow reactions, and their concentrations will be designated by x is , i s q q 1, . . . , m. The number q of metabolites x if is usually greater than the number f of fast reactions. With these notations one may write the differential equations for the time evolution of metabolites: d x if
f
s
dt d x is
n
Ý
Ni fj P Õjf q
js1
Ý
Ni fsj P Õjs ,
i s 1,2, . . . ,q,
Ž 3.2a .
jsfq1
n
s
dt
Ý
Nisj P Õjs ,
i s q q 1,q q 2, . . . ,m.
Ž 3.2b .
jsfq1
Here Ni fj is the stoichiometric coefficient of metabolite X if in the jth fast reaction Ž j s 1, . . . , f .. We have if X if is a substrate, and Ni fj ) 0 if X is a product of the reaction. Ni fsj is the stoichiometric coefficient of metabolite X if in the jth slow reaction Ž j s f q 1, . . . , n.. Nisj is the stoichiometric coefficient of metabolite X is in the jth slow reaction Ž j s f q 1, . . . , n.. Below we shall use matrix notation due to its conciseness. We gather the reaction rates Õjs and Õif in the column vectors U and W, respectively. In matrix form Eqs. Ž3.2a. and Ž3.2b. read: Ni fj - 0
ž
f dQrdt W sNP s N dJrdt U 0
/ ž / ž
N fs P W . U Ns
/ž /
Ž 3.3 .
s s Here Q s Ž x 1f , x 2f , . . . , x qf . T and J s Ž x qq1 , x qq2 , . . . , x ms . T are the vectors of metabolite concentrations. The superscript T denotes transposition of vectors or matrices. N is the stoichiometric matrix of the entire pathway. It has the dimension m = n and it is composed by the q = f matrix N f of the fast reactions and by the m = Ž n y f . matrix ŽN fs N s . T of the slow reactions. When the conditions for applicability of Tikhonov’s theorem are fulfilled, then the fast variables obey the following steady-state equations w20,21x:
N f P W ( Q ) s 0.
Ž 3.4 .
Note that besides the usual continuity requirements, the conditions for applicability of Tikhonov’s theorem state that the fast system should have an asymptotically stable steady-state solution with the initial concentration values lying in its domain of attraction w20,21x, cf. w22,43x. Let r denote the rank of matrix N f. Then f y r columns of this matrix can be expressed in terms of r linearly independent columns. Accordingly, there is a full-rank, f = Ž f y r . matrix that fulfills the equation N f P ff s 0.
Ž 3.5 .
The columns of the matrix ff form a basis of the kernel of the matrix N f. In many reaction systems, the matrix N f is full-rank Ž r s f . Žsee the above example and the models in w16,18,24x.. In this case, ff s 0, and Eq. Ž 3.4. implies W s 0 Žstrictly speaking, this equation as well as Eq. Ž3.4. hold in a fast time scale only.. This means that the fast subsystem attains a quasi-equilibrium state after some initial transient period. The mass action ratios of the fast reactions are then nearly equal to the equilibrium constants of these reactions. This holds true, in particular, at any steady state of the pathway. In other words, the affinities Ž i.e., negative Gibbs free energy changes. of those reactions are close to zero: q
Ý Nifj P ln x if y ln K jeq s 0, is1
j s 1, . . . , f .
Ž 3.6 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
343
Using matrix notation one may write T
Ž N f . P ln Q y ln K eq s 0.
Ž 3.7 .
Here K eq s Ž K 1eq , K 2eq , . . . , K feq . T is the vector of equilibrium constants of the fast reactions. Logarithms of vectors are defined via the logarithms of their components. Eq. Ž 3.7. represents a linear system of f equations in the q logarithmic concentrations ln x if. When the rank, r, of matrix N f is less than f, the quasi-steady-state condition Ž3.4. may nevertheless entail the quasi-equilibrium condition W s 0, if the values of K ieq obey some restrictions. For the usual mass action kinetics, Volpert and Khudiaev w36x and Feinberg w37x derived a necessary and sufficient condition for the identity of steady state and equilibrium. Schuster and Schuster w26x proved this condition to be valid also for a generalized mass-action kinetics, which comprises all reversible enzymatic rate laws. To illustrate this condition we note that a linear dependence between the columns of matrix N f often implies the existence of a cycle in the scheme of fast reactions. Alternatively, it can represent the existence of a ‘‘tree’’ connecting source and sink metabolites w37x. At equilibrium the net flux through any cycle vanishes. This implies the existence of a constraint on the possible values of the equilibrium constants of reactions within a cycle. For a simple cycle in which at any steady state all rates must be equal, this constraint is the well-known ‘‘detailed balance’’ relationship or Wegscheider condition w38,39x:
Ł K ieq s Ł kqi rkyi s 1, i
Ž 3.8 .
i
y Ž . where kq i and k i are the forward and backward rate constants, respectively. The product in Eq. 3.8 is taken over all the reactions participating in the cycle. Note that the concentrations of ligands which are assumed to be constant must be included into the equilibrium constants and rate constants. Consider, e.g., the reaction system shown in Scheme 2. The Wegscheider condition for the cycle composed of the three fast reactions reads
k1 k 2 k 3 ky1 ky2 ky3
s 1.
Ž 3.9 .
To generalize the condition given by Eq. Ž 3.8. for a pathway of arbitrary stoichiometry one can premultiply Eq. Ž3.7. by the matrix Ž ff . T. Then, using Eq. Ž3.5. one arrives at the condition for the equivalence of steady state and equilibrium state of the fast subsystem w26,36,37x:
Ž ff .
T
P ln K eq s 0.
Ž 3.10 .
Under this condition, r linearly independent equations of Eq. Ž3.7. allow one to express r Ždependent. concentrations of metabolites participating in the fast reactions in terms of the remaining q y r concentrations chosen as independent concentrations. To do so one may choose a submatrix, NRf , of N f by using r linearly independent columns Ž without loss of generality we shall consider that the reactions are renumbered such that
Scheme 2.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
344
the first r columns are linearly independent.. Then one may replace the f equations in Eq. Ž3.7. by r linearly independent equations only:
Ž NRf .
T
P ln Q s ln K eq R .
Ž 3.11 .
Here K eq R is the reduced vector of the equilibrium constants of those reactions which correspond to the reduced matrix NRf . To construct the system of equilibrium relationships Ž3.11. , one should cancel one reaction for each cycle and one reaction for each branch of a ‘‘tree’’ Ž such reactions whose deletion does not decrease the rank of the stoichiometry matrix of the fast subsystem.. When the equilibrium state Ž defined by Eq. Ž 3.7.. exists, it is globally asymptotically stable for thermodynamic reasons w26,40,41x. This implies that there cannot exist any other stable steady state for the same moiety conservation quantities and that the stability hypothesis in Tikhonov’s theorem is justified. It is convenient to express the reduced q = r matrix NRf as the product of a q = r link matrix, Lf , w6x and a square r = r matrix, Nrr, which consists of r linearly independent rows of the matrix NRf : NRf s Lf P Nrr ,
Lf s
Ir
ž / Lf0
.
Ž 3.12 .
The Ž q y r . = r matrix Lf0 expresses the dependent q y r rows of the matrix into NRf its r independent rows. The link matrices for the whole matrix N f and for the reduced matrix NRf are identical because NRf results from N f by cancelling linearly dependent columns. When r - q, there exist conservation relations in the fast subsystem which hold in a fast timescale only. It is important to note that this does not contradict our above assumption that conservation relations in the entire system Ž if any. only comprise metabolites participating in slow reactions. In the fast subsystem brought about by cancelling all slow reactions, generally additional conservation relations occur. Substituting Eq. Ž3.12. into Eq. Ž 3.11. and choosing the first r concentrations of the metabolites participating in the fast reactions as the dependent concentrations, Qdep s Ž x 1f , x 2f , . . . , x rf . T, and the last q y r concentrations f f as the independent concentrations, Q ind s Ž x rq1 , x rq2 , . . . , x qf . T one may write after some rearrangement: T
ln Qdep s y Ž Lf0 . P ln Q ind q Ž Nrr .
ž
T y1
/
P ln K eq R .
Ž 3.13 .
Eq. Ž3.13. expresses the first r concentrations of metabolites participating in the fast reactions in terms of the equilibrium constants of these reactions and of the last q y r concentrations chosen as the independent concentrations. Let h be the number of independent conservation relations in the slow subsystem. These relations can be expressed by J s Ls J ind q const,
Ž 3.14 .
where Ls is the link matrix of the slow subsystem w6x. This matrix links the stoichiometry matrix N s to a submatrix of linearly independent rows. J ind is the vector of ‘‘slow’’ metabolites independent with respect to conservation relations.
4. Control analysis in the framework of quasi-equilibrium approximation As we wish to carry out a steady-state control analysis, we consider not only the fast reactions of the network to attain a steady state, but also the slow reactions to do so. As mentioned in Section 1 we shall assume that all moiety conservation constraints valid in the original network Ž if any. only involve metabolites not participating in fast reactions.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
345
At any steady state all rates of slow reactions can be expressed via a set of linearly independent slow rates. Let g denote the number of these independent slow rates. This number is given by the number of slow reactions Ž n y f . minus the number of independent steady-state equations with respect to metabolite concentrations. The latter number is given by the number of metabolites Ž m. minus the number of independent quasi-equilibrium conditions Ž r . minus the number of independent conservation relations in the slow subsystem Ž h.. Therefore, we have g s n y f y m q r q h.
Ž 4.1 .
As for Scheme 2, e.g., n s 5, f s 3, m s 3, r s 2 Žwhen the generalized Wegscheider condition is fulfilled in the cycle., and h s 0. This gives, correctly, g s 1. The two slow reactions actually carry the same flux at steady state. We shall designate by JR the Žreduced. vector consisting of g linearly independent steady-state rates of slow reactions T
s s s JR s Ž Õ fq1 ,Õ fq2 , . . . ,Õ fqg . < st .st .
Ž 4.2 .
Žwithout loss of generality it is assumed that the slow reactions are renumbered such that at steady state the first g slow rates are linearly independent. . The vector U of slow rates can be expressed in terms of the reduced vector JR as follows Žcf. w6,8,33x, and Eq. Ž2.2..: J s U
Ž 4.3 .
Here the vector of the slow steady-state rates is designated by J, f is an Ž n y f . = g matrix, which can be considered as the null-space matrix of the reduced reaction system. Its first g rows coincide with the g-dimensional identity matrix fs
Ig
ž / f0
.
Ž 4.4 .
We now calculate the g = Ž n y f . matrix C UJR containing the control coefficients over the fluxes of the independent slow reactions. Using Eq. Ž 4.3. and definition Ž 1.1., we obtain
˜ UJR . C UJ s fC
Ž 4.5 .
˜ relates to the matrix f of Eq. Ž4.3. as follows: The Ž n y f . = g matrix f ˜ s Ž d ln Urd ln JR . < st .st .s Ž diag U . < st .st .P f P Ž diag JR . . f y1
Ž 4.6 .
Here and below Ždiag F. designates Žfor any vector F. the square matrix with diagonal elements equal to the components Fi and all off-diagonal elements equal to zero. We employ the usual notation Ž d ln Ard ln B. for the matrix Žd ln A ird ln Bj . consisting of the derivatives of each component of the vector ln A with respect to each component of the vector ln B. On the other hand, the matrix of control coefficients, C UJ , of the slow reactions can be calculated by differentiation of the steady-state rates of the slow reactions considered as functions of the metabolite concentrations and the parameters pi Žsee Eq. Ž2.7. and w8x.:
E ln J Ep
s
ž
E ln U
E ln Qdep
E ln Qdep E ln Q ind
E ln U q
E lnQ ind
/
E ln Q ind Ep
E ln U q
E ln J E ln J ind
E ln J E ln J ind
Ep
E ln U q
Ep
.
Ž 4.7 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
346
Since the system is to contain no conservation relations involving the components of Q, there are no dependencies between the vectors Q and J and no dependencies between the components of Qdep and Q ind other than the quasi-equilibrium relations. Owing to Eq. Ž3.13., we have
Ž d ln Qdeprd ln Q ind . s y Ž Lf0 .
T
.
Ž 4.8 .
Therefore, y Ž Lf0 . Ž d ln Qrd ln Q ind . s I qyr
T
0
˜ r. sL
Ž 4.9 .
˜ r in analogy to the link matrix L Žsee Eq. Ž3.12... Most importantly, We designate this q = Ž q y r . matrix L this matrix does not depend on the steady-state concentrations, nor even on the equilibrium constants of the fast reactions. Dividing Eq. Ž4.7. by E ln UrE p gives T
C UJ s y´ QUdep P Ž Lf0 . q ´ QUind P C QU ind q ´ U J
ž
/
E ln J E ln J ind
P C JU ind q I nyf .
Ž 4.10 .
ind Ž q y r . = Ž n y f . matrix of the control coefficients over the independent Here C Q and C J U U are the Ž . concentrations Q ind of metabolites participating in fast reactions and the Ž m y q . = Ž n y f . matrix of the control coefficients over the concentrations of metabolites participating only in slow reactions, respectively. By Eq. Ž3.14. and the rules for differentiation of logarithms, we conclude
E ln J E ln J ind
s Ž diag J .
y1
EJ E J ind
y1 Ž diag J ind . s Ž diag J . Ls Ž diag J ind . .
Ž 4.11 .
The control coefficients over the vector J can therefore be expressed by the control coefficients over the concentrations J ind independent with respect to conservation relations as follows: CJ U s diag Ž J .
y1
P Ls P diag Ž J ind . P C JU ind .
Ž 4.12 .
After defining the Ž m y h. = Ž m y r y h. matrix,
˜s L
y Ž Lf0 .
T
0
I qyr 0
0 diag Ž J .
y1
Ls diag Ž J ind .
0
Ž 4.13 .
Eq. Ž4.10. can be written as
˜ P C UX ind q I nyf . C UJ s ´ xU P L
Ž 4.14 .
The Ž m y r y h.-dimensional vector X ind is generated from the vector X by cancelling the r concentrations Qdep dependent due to the quasi-equilibrium conditions and the h concentrations J dep dependent due to the conservation relations in the slow subsystem. ´Ux is the Ž n y f . = Ž m y h. matrix of the elasticities of the slow reaction rates ŽU. with respect to all the pathway concentrations except for J dep , ´Ux s ´UQdep < ´UQ ind < ´UJ ind .
ž
/
Ž 4.15 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
347
C UX ind is the Ž m y r y h. = Ž n y f . matrix of the control coefficients over all the independent concentrations of the reduced system. This matrix is defined as d ln Q ind d ln U C UX ind s
ž / ž /
dp
y1
dp
d ln J ind d ln U dp
y1
dp
0
s
C QU ind
ž / C JU ind
.
Ž 4.16 .
Combining Eqs. Ž4.5. and Ž 4.14. one arrives at:
ž
˜ < y ´UX P L ˜ P f
/
C UJR
ž / C UX ind
s E P C U s I nyf .
Ž 4.17 .
By inverting the matrix E, one can express the control coefficients of the slow reactions in terms of the elasticities of the slow reactions, without knowledge of the elasticities of the fast reactions. On the other hand, by inverting the matrix C U , one can calculate the elasticity coefficients of the slow reactions in terms of their ˜ and ye Ux P L ˜ have dimensions Ž n y f . = g and Ž n y f . = Ž m y r y h., and control coefficients. The matrices f JR X ind the matrices C U and C U have dimensions g = Ž n y f . and Ž m y r y h. = Ž n y f ., respectively. Therefore, Eq. Ž4.1. guarantees that the two matrices on the left-hand side of Eq. Ž4.17. are square. From these general formulas, we can readily derive the results for the example in Scheme 1. We have m s q s 2, n s 3, f s 1. The stoichiometry matrix of the fast reactions reads N f s Žy1 1. T with the submatrix ˜ s Ž1 Nrf s y1. Then, Q s Ž x 1, x 2 . T, Q ind s x 2 . The link matrix reads Lf s Ž 1 y1. T, hence Lf0 s y1 and L T T ˜ ˜ ˜ 1. . Since the pathway is unbranched, f s Ž1 1. . Substituting the matrices f and L into Eq. Ž4.17. one can again derive Eq. Ž2.10..
5. An illustration We shall illustrate our method applying it to the metabolic pathway shown in Scheme 3. The system involves six internal metabolites Ž m s 6. and nine reactions Ž n s 9., three of which are considered fast Ž f s 3. . There are five metabolites involved in fast reactions Ž q s 5. . If we write T
X s Ž X 1f , X 2f , X 3f , X 4f , X5f , X6s . , T
v s Ž Õ 1f ,Õ 2f ,Õ 3f ,Õ4s ,Õ5s ,Õ6s ,Õ 7s ,Õ 8s ,Õ 9s . ,
Scheme 3.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
348
the dynamics of the system is described by the equation dX dt
s N P v,
Ž 5.1 .
where N is the stoichometric matrix
y1 0 0 Ns 1 0 0
0 y1 0 1 0 0
0 0 1 0 y1 0
0 0 0 0 0 1
1 0 0 0 0 y1
1 0 y1 0 0 0
0 0 0 y1 0 0
0 1 0 0 0 y1
0 0 0 . 0 1 y1
0
Ž 5.2 .
For the fast subsystem:
y1 0 Nfs 0 1 0
0 y1 0 1 0
0 0 1 . 0 y1
0
Ž 5.3 .
As the rank of N f is 3; NRf s N f , and the equation N f s Lf P Nrf gives us 1 0 f Ls 0 y1 0
0 1 0 y1 0
0 0 1 , 0 y1
0
Nrf s
y1 0 0
0 y1 0
0 0 . 1
0
Ž 5.4 .
For this case, n s 9, f s 3, m s 6, r s 3, h s 0, so the number of independent fluxes g is 3. We choose J4 , J5 and J6 : J 7 s J4 ; J8 s J4 y J5 y J6 ; J9 s J6 ; JR s Ž J4 , J5 , J6 . T and U s Ž J4 , J5 , J6 , J 7 , J8 , J9 . T. 1 0 ˜s 0 f 1 J4rJ8 0
0 1 0 0 yJ5rJ8 0
0 0 1 , 0 yJ6rJ8 1
0
Ž 5.5 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
°´ ´
4 4 1 2
¶
´ 34´ 44´ 54´ s4
1 1 ˜s 0 L 1 0 0
´ 15´ 25 ´ 35´ 45´ 55´ s5 ´UX s
´ 16´ 26 ´ 36´ 46´ 56´ s6 ´ 17´ 27 ´ 37 ´ 47´ 57´ s7 9 9 1 2
0 0 1 0 1 0
0 0 0 . 0 0 1
0
,
´ 18´ 28 ´ 38´ 48´ 58´ s8
¢´ ´
349
ß
´ 39 ´ 49´ 59´ s9
Ž 5.6 .
˜ . and inverting the result we can obtain the matrix Multiplying the elasticity matrix Ž ´Ux . by the link matrix ŽL of control coefficients:
°C
¶
J4 4
C5J4
C6J4
C 7J4
C8J4
C9J4
C4J5
C5J5
C6J5
C 7J5
C8J5
C9J5
C4J6
C5J6
C6J6
C 7J6
C8J6
C9J6
C4x 4
f
C5x 4
f
C6x 4
f
C 7x 4
f
C8x 4
f
C9x 4
C4x 5
f
C5x 5
f
C6x 5
f
C 7x 5
f
C8x 5
f
C9x 5
x 6s 4
C5x 6
s
C6x 6
s
C 7x 6
s
C8x 6
s
C9x 6
¢C
f
f
ß
s
°1
0
0
y Ž ´ 14 q ´ 24 q ´ 44 . y Ž ´ 34 q ´ 54 . y ´ s4
0
1
0
y Ž ´ 15 q ´ 25 q ´ 45 . y Ž ´ 35 q ´ 55 . y ´ s5
0
0
1
y Ž ´ 16 q ´ 26 q ´ 46 . y Ž ´ 36 q ´ 56 . y ´ s6
1
0
0
y Ž ´ 17 q ´ 27 q ´ 47 . y Ž ´ 37 q ´ 57 . y ´ s7
J4rJ8
yJ5rJ8
yJ6rJ8
y Ž ´ 18 q ´ 28 q ´48 . y Ž ´ 38 q ´ 58 . y ´ s8
0
1
y Ž ´ 19 q ´ 29 q ´49 . y Ž ´ 39 q ´ 59 . y ´ s9
s
¢0
¶
y1
.
Ž 5.7 .
ß
As in the simple example dealt with in Section 2, the effective elasticities of slow reactions are sums of elasticities with respect to the metabolites participating in fast reactions. One of these sums runs over the metabolites X 1f , X 2f and X 4f forming one component of the fast subsystem. The other sum runs over X 3f and X5f , i.e., over another component of this subsystem.
6. Discussion In the above presentation, we have shown how the control coefficients of slow reactions over the fluxes through these reactions can be calculated in terms of the elasticities of the slow reactions and vice versa. So we have solved what is sometimes called the direct and inverse problems of MCA w33,42x for the situation that some reactions subsist in quasi-equilibrium. The method presented is of considerable help in the control analysis of metabolic systems characterized by time hierarchy because it allows one to draw conclusions about the control properties even if the kinetic properties of the fast reactions are imperfectly known. Despite the recent availability of powerful computers, kinetic and control analysis of large reaction networks is often hampered or completely impeded by this incompleteness in knowledge of parameters. Our analysis concerns quasi-equilibrium reactions, i.e., reactions which still carry a non-zero flux although their mass-action ratios nearly equal their equilibrium constants. A related problem is to calculate control
350
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
coefficients for true equilibrium reactions for which the net flux is zero when the whole system is at steady state. For the case of adenylate kinase this problem was tackled by Fell and Sauro w35x. Eqs. Ž4.13. and Ž 4.17. show that neither the kinetic parameters, nor the elasticities, nor the equilibrium constants of the fast reactions need to be known. Since the matrix Ž d ln X deprd ln X ind . is a matrix of constant coefficients Žit depends neither on metabolite concentrations nor on equilibrium constants. we need to know only network stoichiometry and the elasticity coefficients of the slow reactions at the operating steady state. Note that f . the derivative ŽdX frdX ind , which would be needed for calculating non-normalized control coefficients does depend on the equilibrium constants. This is an advantage of normalized Ž i.e., logarithmic. control coefficients. For further advantages and drawbacks of normalization see w7x. One can see that upon derivation of Eq. Ž4.17. for determining the control properties we avoid to consider ‘‘slow variables’’, also called ‘‘pool variables’’ w2x Ž these were considered in the example in Section 2 only. . Using ‘‘pool variables’’ is suitable if one wishes to consider the time evolution of the system in the framework of the quasi-equilibrium approximation w24,43x. Moreover, it may be helpful to use these variables when pathways with arbitrary conservation relations are studied. The partitioning of concentrations into independent and dependent ones is of importance for introducing any non-linear constraints among the concentrations of metabolites Ž cf. w29,44x.. For expressing control coefficients in terms of elasticities, it is then necessary to calculate the matrix Ž d ln X deprd ln X ind . of derivatives of dependent concentrations with respect to the independent ones. For the study of in vivo situations, it is worthwhile combining the analysis presented here with other extensions of MCA, such as top-down or modular approaches Ž see, e.g., w45,46x.. In the analysis of metabolite channelling w47,48x the elimination of quasi-equilibrium steps can also be of great help.
Acknowledgements This work was supported by the Research grant of Comision ´ Interministerial de Ciencia y Tecnologıa, ´ Spain ŽCICYT. PM 96-0099, by Programa Catedra Ž FBBV. , Spain by the Netherlands Organization for Scientific ´ Research and by the E.U. grant B104-CT96-0662.
References w1x R. Heinrich, S.M. Rapoport, T.A. Rapoport, Metabolic regulation and mathematical models, Prog. Biophys. Mol. Biol. 32 Ž1977. 1–82. w2x J.G. Reich, E.E. Selkov, Energy Metabolism of the Cell, Academic Press, London, 1981. w3x S. Minakami, H. Yoshikawa, Studies on erythrocyte glycolysis. II. Free energy changes and rate limiting steps in erythrocyte glycolysis, J. Biochem. ŽTokyo. 59 Ž1969. 139–144. w4x D.A. Fell, Metabolic control analysis: A survey of its theoretical and experimental development, Biochem. J. 286 Ž1992. 313–330. w5x R. Heinrich, T.A. Rapoport, A linear steady-state treatment of enzymatic chains. General properties, control and effector strength, Eur. J. Biochem. 42 Ž1974. 89–95. w6x C. Reder, Metabolic control theory. A structural approach, J. Theor. Biol. 135 Ž1988. 175–201. w7x S. Schuster, R. Heinrich, The definitions of metabolic control analysis revisited, BioSystems 27 Ž1992. 1–15. w8x B.N. Kholodenko, D. Molenaar, S. Schuster, R. Heinrich, H.V. Westerhoff, Defining control coefficients in non-ideal metabolic pathways, Biophys. Chem. 56 Ž1995. 215–226. w9x H. Kacser, J.A. Burns, The control of flux, Symp. Soc. Exp. Biol. 27 Ž1973. 65–104 Župdated reprint by: H. Kacser, J.A. Burns, D.A. Fell, Biochem. Soc. Trans. 23 Ž1995. 341–366.. w10x J. Delgado, J.C. Liao, Control of metabolic pathways by time-scale separation, BioSystems 36 Ž1995. 55–70. w11x H.V. Westerhoff, K. van Dam, Thermodynamics and Control of Biological Free-energy Transduction, Elsevier, Amsterdam, 1987. w12x B.N. Kholodenko, Control of molecular transformations in multienzyme systems: quantitative theory of metabolic regulation, Molek. Biol. ŽMoscow. 22 Ž1988. 1238–1256 Žin Russian.; English trans. in: Molec. Biol. ŽMoscow. 22 Ž1998. 990–1005.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
351
w13x A.N. Zhabotinsky, Autonomous Oscillations in Concentrations Žin Russian., Nauka, Moscow, 1974. w14x E.E. Selkov, Stabilization of energy charge, generation of oscillations and multiple steady states in energy metabolism as a result of purely stoichiometric regulation, Eur. J. Biochem. 59 Ž1975. 151–157. w15x B.N. Kholodenko, A study of the dynamic behavior of the glycolytic system, Biofizika 24 Ž1979. 640–645 Žin Russian.; English trans. in: Biophysics ŽUSSR. 24 Ž1979. 659–665. w16x B.N. Kholodenko, A.V. Stepanchicova, L.I. Ehrlich, F.I. Ataullakhanov, A.M. Zhabotinsky, A.V. Pichugin, The role of the 2,3-diphosphoglycerate bypass in the regulation of energetic metabolism in human erythrocytes Žin Russian., Izv. Acad. Sci. USSR ŽBiol.. 2 Ž1985. 196–205. w17x H.-G. Holzhutter, G. Jacobasch, A. Bisdorff, Mathematical modelling of metabolic pathways affected by an enzyme deficiency, Eur. ¨ J. Biochem. 149 Ž1985. 101–111. w18x R. Schuster, H.G. Holzhutter, G. Jacobasch, Interrelations between glycolysis and the hexose monophosphate shunt in erythrocytes ¨ as studied on the basis of a mathematical model, BioSystems 22 Ž1988. 19–36. w19x G. Pettersson, U. Ryde-Pettersson, Dependence of the Calvin cycle activity on kinetic parameters for the interaction of non-equilibrium cycle enzymes with their substrates, Eur. J. Biochem. 186 Ž1989. 683–686. w20x A.N. Tikhonov, On the dependence of the solutions of differential equations on a small parameter Žin Russian., Mat. Sborn. 22 Ž1948. 193–204. w21x A.N. Tikhonov, Systems of differential equations containing small parameters in their derivatives Žin Russian., Mat. Sborn. 31 Ž1952. 575–586. w22x W. Klonowski, Simplifying principles for chemical and enzyme reaction kinetics, Biophys. Chem. 18 Ž1983. 73–87. w23x V.M. Vasilev, A.I. Volpert, S.I. Khudiaev, On the method of quasi-equilibrium concentrations for chemical kinetic equations Žin Russian., Zhurn. Vychysl. Mat. Matemat. Fiz. 13 Ž1973. 683–697. w24x M. Schauer, R. Heinrich, Quasi-steady-state approximation in the mathematical modeling of biochemical reaction networks, Math. Biosci. 65 Ž1983. 155–171. w25x J.C. Liao, E.N. Lightfoot Jr., Lumping analysis of biochemical reaction systems with time scale separation, Biotechnol. Bioeng. 31 Ž1988. 869–879. w26x S. Schuster, R. Schuster, A generalization of Wegscheider’s condition. Implications for properties of steady states and for quasi-steady-state approximation, J. Math. Chem. 3 Ž1989. 25–42. w27x B.N. Kholodenko, H.V. Westerhoff, The macroworld versus the microworld of biochemical regulation and control, Trends Biochem. Sci. 20 Ž1995. 52–54. w28x H.V. Westerhoff, Y.-D. Chen, How do enzyme activities control metabolite concentrations? An additional theorem in the theory of metabolic control, Eur. J. Biochem. 142 Ž1984. 425–430. w29x B.N. Kholodenko, Metabolic Control Theory. New relationships for determining control coefficients of enzymes and response coefficients of system variables, J. Nonlin. Biol. 1 Ž1990. 107–126. w30x H.M. Sauro, J.R. Small, D.A. Fell, Metabolic control and its analysis. Extensions to the theory and matrix method, Eur. J. Biochem. 165 Ž1987. 215–221. w31x H.V. Westerhoff, D.B. Kell, Matrix method for determining steps most rate-limiting to metabolic fluxes in biotechnological processes, Biotechnol. Bioeng. 30 Ž1987. 101–107. w32x M. Cascante, R. Franco, E.I. Canela, Use of implicit methods from general sensitivity theory to develop a systematic approach to metabolic control. I. Unbranched pathways, Math. Biosci. 94 Ž1989. 271–288; II. Complex systems, Math. Biosci. 94 Ž1989. 289–309 w33x H.V. Westerhoff, J.-H.S. Hofmeyr, B.N. Kholodenko, Getting to the inside of cells using metabolic control analysis, Biophys. Chem. 50 Ž1994. 273–283. w34x B. Crabtree, E.A. Newsholme, The derivation and interpretation of control coefficients, Biochem. J. 247 Ž1987. 113–120. w35x D.A. Fell, H.M. Sauro, Metabolic control and its analysis. Additional relationships between elasticities and control coefficients, Eur. J. Biochem. 148 Ž1985. 555–561. w36x A.I. Volpert, S.I. Khudiaev, Analysis in Classes of Continuous Functions and the Equations of Mathematical Physics Žin Russian., Nauka, Moscow, 1975 w37x M. Feinberg, Necessary and sufficient conditions for detailed balancing in mass action systems of arbitrary complexity, Chem. Eng. Sci. 44 Ž1989. 1819–1827. ¨ w38x R. Wegscheider, Uber simultane Gleichgewichte und die Beziehungen zwischen Thermodynamik und Reaktionskinetik homogener Systeme, Z. Phys. Chem. 39 Ž1902. 257–303. w39x J.Z. Hearon, The kinetics of linear systems with special reference to periodic reactions, Bull. Math. Biophys. 15 Ž1953. 121–141. w40x Ya.B. Zeldovich, Proof of a unique solution to the mass action law Žin Russian., Zh. Fiz. Khim. 11 Ž1938. 685–687. w41x F. Horn, R. Jackson, General mass action kinetics, Arch. Rational Mech. Anal. 47 Ž1972. 81–116. w42x B.N. Kholodenko, H.M. Sauro, H.V. Westerhoff, M. Cascante, Coenzyme cycles and metabolic control analysis: the determination of the elasticity coefficients from the generalized connectivity theorem, Biochem. Mol. Biol. Int. 35 Ž1995. 615–625.
352 w43x w44x w45x w46x
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
R. Heinrich, S. Schuster, The Regulation of Cellular Systems, Chapman and Hall, New York, 1996. S. Schuster, Control analysis in terms of generalized variables characterizing metabolic systems, J. Theor. Biol. 182 Ž1996. 259–268. M.D. Brand, Top down metabolic control analysis, J. Theor. Biol. 182 Ž1996. 351–360. S. Schuster, D. Kahn, H.V. Westerhoff, Modular analysis of the control of complex metabolic pathways, Biophys. Chem. 48 Ž1993. 1–17. w47x B.N. Kholodenko, H.V. Westerhoff, Metabolic channelling and control of the flux, FEBS Lett. 320 Ž1993. 71–74. w48x B.N. Kholodenko, M. Cascante, H.V. Westerhoff, Control theory of metabolic channelling, Mol. Cell. Biochem. 143 Ž1995. 151–168.
Control analysis of metabolic systems involving quasi-equilibrium reactions Boris N. Kholodenko a
a,b
, Stefan Schuster c , Jordi Garcia a , Hans V. Westerhoff d , Marta Cascante a,)
Departament de Bioquımica i Fisiologia, Facultat de Quımica, UniÕersitat de Barcelona, Martı´ i Franques, ´ ´ ` 1, 08028 Barcelona, Spain b A.N. Belozersky Institute of Physico-Chemical Biology, Moscow State UniÕersity, 119899 Moscow, Russian Federation c Max Delbruck ¨ Center for Molecular Medicine, 13125 Berlin-Buch, Germany d BioCentrum, Department of Microbial Physiology, Faculty of Biology, Free UniÕersity, De Boelelaan 1087, NL-1081 HV Amsterdam, Netherlands Received 17 April 1997; revised 3 September 1997; accepted 5 September 1997
Abstract Reactions for which the rates are extremely sensitive to changes in the concentrations of variable metabolite concentrations contribute little to the control of biochemical reaction networks. Yet they do interfere with the calculation of the system’s behaviour, both in terms of numerical integration of the rate equations and in terms of the analysis of metabolic control. We here present a way to solve this problem systematically for systems with time hierarchies. We identify the fast reactions and fast metabolites, group them apart from the other Ž‘‘slow’’. reactions and metabolites, and then apply the appropriate quasi-equilibrium condition for the fast subsystem. This then makes it possible to eliminate the fast reactions and their elasticity coefficients from the calculations, allowing the calculation of the control coefficients of the slow reactions in terms of the elasticity coefficients of the slow reactions. As expected, the elasticity coefficients of the fast reactions drop out of the calculations, and they are irrelevant for control at the time resolution of the steady state of the slow reactions. The analysis, when applied iteratively, is expected to be particularly valuable for the control analysis of living cells, where a time hierarchy exists, the fastest being at the level of enzyme kinetics and the slowest at gene expression. q 1998 Elsevier Science B.V. Keywords: Control coefficient; Metabolic control analysis; Quasi-equilibrium approximation; Time hierarchy
1. Introduction The pronounced hierarchy of the characteristic times of enzyme reactions is a fundamental property of cellular metabolism w1,2x. Due to differences in the enzyme activities it is often possible to subdivide the reactions into
)
Corresponding author. Fax: q34 3 402 1219.
0304-4165r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 0 4 - 4 1 6 5 Ž 9 7 . 0 0 1 1 4 - 1
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
338
two classes: Ž very. fast reactions and slow reactions. Experimental measurements often show that at the operating state the ratios of substrate and product concentrations Žmass action ratios. of the fast reactions are close to the corresponding equilibrium constants Žsee, e.g., the experimental determination of metabolite concentrations of glycolysis in erythrocytes w3x.. The sensitivities of relevant variables in steady states of chemical and biochemical reaction systems Ž such as steady-state fluxes and concentrations of metabolites. to changes in the activities of particular reactions can be studied by the use of a theoretical framework known as metabolic control analysis Žfor a review, see w4x.. A general way of defining control coefficients is w5–8x: CiJk s
E ln < Jk
i , j,k s 1,2,3 . . .
Ž 1.1 .
Here, Õi , Jk and x j denote the velocity of reaction i, the flux through reaction k at a stationary state under consideration, and the concentration of a substance X j , respectively. pi is a parameter affecting only reaction i directly, i.e.,
E ln < Õi < E pi
/ 0,
E ln < Õi < E pk
s 0 for any k / i.
Ž 1.2 .
The derivatives of the logarithmic rates ln < Õi < with respect to the parameters in Eqs. Ž1.1. and Ž1.2. are taken as if reaction i proceeded in isolation from the remaining network. When there is a one-to-one correspondence between enzymes and reactions and each reaction rate is proportional to the concentration of the corresponding enzyme, definition Ž 1.1. is equivalent to the definition of control coefficients given by Kacser and Burns w9x in terms of enzyme concentrations. The kinetic properties of the individual enzymes are expressed in terms of their elasticity coefficients Ž e ., which describe the sensitivities of the reaction rates to changes in the metabolite concentrations w5,9x:
´ ki s d ln < Õi
Ž 1.3 .
In contrast to the definition of the control coefficients Ž Eq. Ž 1.1.. , no evolution of the entire pathway to a steady state is presupposed in Eq. Ž 1.3. . Comparing the effects of small changes in the activities of Žvery. fast reactions and of slow reactions on the steady-state variables one can usually neglect the effects of the former in comparison to the latter. Therefore, as has often been stated, e.g., w4,5,10x, the control coefficients of those reactions which have mass action ratios close to the equilibrium constants are virtually equal to zero. This is in line with a fundamental property of fast reactions: the elasticity coefficients of these reactions with respect to the substrates and products assume very high positive and negative values, respectively, since at the system’s steady states these reactions are sufficiently close to thermodynamic equilibrium Ž although the flux through these reactions is non-zero. w11,12x. This feature hampers the determination of control coefficients Ž quantitative indicators of ‘‘global’’ control properties of a pathway. in terms of the elasticity coefficients Žreflecting local kinetic properties of an enzyme reaction. . It is therefore useful to elaborate methods for determining control properties using the elasticities of the slow reactions only. The quasi-equilibrium approximations have been known for a long time Ž see, e.g., w13x, and references therein.. These approximations Žalso referred to as rapid-equilibrium approximations. have been employed in the modelling of the dynamics of various metabolic systems Ž see, e.g., w1,14–19x. . Although a metabolic pathway as a whole may be far from equilibrium, the quasi-equilibrium approximation allows one to consider a subsystem comprising fast reactions to be in a Ž quasi-. equilibrium state. This makes it possible to avoid a kinetic analysis of the fast reactions, i.e., to take into consideration only their equilibrium constants instead of considering their rates.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
339
In the general case the conditions under which one can reduce the number of dynamic variables in the system using the quasi-steady-state approximation were formulated for the first time by Russian mathematician Tikhonov w20–22x. The way of using this general formalism to validate the rapid-equilibrium approximation in the case of biochemical pathways was elaborated by several authors, see, e.g., w23–26x. Until now, the control properties of metabolic pathways with quasi-equilibrium steps have not yet been studied in a systematic way. The goal of this paper is to develop, in the framework of quasi-equilibrium approximation, a general method for determining the control coefficients in terms of the elasticity coefficients of the slow reactions only. In Section 2, the method is illustrated by way of an example. Section 3 is aimed at presenting the rapid-equilibrium approximation in a matrix notation so as to provide the basis for the control analysis set out in Section 4. As will become clear from below, when the system involves conservation relations saying that linear combinations of concentrations remain constant in time, these relations have to be taken into account in the control analysis of networks with rapid-equilibrium subsystems. This is particularly complicated when such relations involve substances participating in fast reactions. To begin with a clear presentation which shows the essential ideas and to avoid an unnecessarily complicated analysis, we here exclude conservation relations of this type. We allow for conservation relationships among the substances participating in slow reactions only. The general case of systems with arbitrary conserved-moiety structures will be analysed in a sequel paper.
2. A simple example We shall begin with illustrating the principles of the method using the exemplifying reaction system shown in Scheme 1. This is an unbranched metabolic pathway consisting of three sequential monomolecular reactions. They may represent, e.g., the reactions catalysed by hexokinase, phosphoglucoisomerase, and phosphofructokinase involved in glycolysis. The symbols Õ s and Õ f stand for the rates of ‘‘slow’’ and ‘‘fast’’ reactions, respectively. We here assume that direct interactions between the enzymes are absent w27x. The concentrations of the initial substrate, S, and of the final product, P, are taken to be constant. The differential equations for the time evolution of the concentrations x 1 and x 2 read d x 1rd t s Õ 1s y Õ 2f ,
Ž 2.1a .
d x 2rd t s Õ 2f y Õ 3s .
Ž 2.1b .
At any steady state of the system the time derivatives are equal to zero. Therefore, all three reactions carry the same flux Ž J1 ., and one may write s f s Õ 1
Ž 2.2 .
Here and below, the subscript st.st. specifies steady-state conditions. System Ž 2.2. consisting of three equations allows one to determine the concentrations x 1 and x 2 and the flux J1 through the pathway of Scheme 1 at steady state, provided that the rate laws are known. The control properties of this metabolic pathway are described by the control coefficients Ž CiJ1 and Cix k . of the enzymes over the pathway flux and metabolite concentrations. The control coefficients of the enzymes depend
Scheme 1.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
340
on their kinetic properties and on pathway structure. One way of expressing this dependence is by using the summation and connectivity theorems of metabolic control analysis w6,9,12,28,29x. They can be written in a concise and elegant form of a matrix equation w6,30–33x. For the pathway of Scheme 1 this matrix equation can be written as follows:
C1J1 C1x 1 C1x 2
1 C3J1 C3x 1 P 1 C3x 2 1
C2J1 C2x 1 C2x 2
0
s
y´ 1Õ1 f
y´ 1Õ 2 y´
Õ 3s 1
s
y´ 2Õ1
1 f y´ 2Õ 2 s 0 s 0 y´ Õ3 2
0
0 1 0
0 0 . 1
0
Ž 2.3 .
The products of the rows of the control matrix ŽC. and the first column of the elasticity matrix ŽE. give all the summation theorems for the flux and concentration control coefficients for the pathway of Scheme 1. The products of the rows of matrix C and the next two columns of matrix E give all the connectivity theorems. Writing these theorems in the form of a matrix equation Ž Eq. Ž 2.3.. allows one to calculate the control coefficients by inversion of the elasticity matrix E. Now we assume that the activity of enzyme 2 is so high that reaction 2 can be considered to be, after a very short transient period, at a quasi-equilibrium state Ž which is actually the case for phosphoglucoisomerase in the glycolytic pathway w3x, for example.. For the time evolution of the metabolites ŽEqs. Ž2.1a. and Ž2.1b.. this implies that the ratio of the concentrations of x 2 and x 1 approaches the equilibrium constant of reaction 2: x 2rx 1 s K 2eq .
Ž 2.4 .
According to this equation only one concentration out of the two concentrations x 1 and x 2 can now be considered independent. Without loss of generality, we shall choose x 2 as the independent variable. Summing Eqs. Ž2.1a. and Ž2.1b. , one arrives at an equation describing the time evolution of a so-called ‘‘slow variable’’ or ‘‘pool variable’’ w1,2x. Only slow enzyme rates enter the new equation: d Ž x 1 q x 2 . rd t s d yrd t s Õ 1s y Õ 3s .
Ž 2.5 .
The pool variable, y s x 1 q x 2 , is the sum of the concentrations of the metabolites participating in the fast reaction. Since the requirements of Tikhonov’s theorem are fulfilled Ž see, e.g., w24,26x. one can substitute Eq. Ž2.4. into Eq. Ž2.5.. The obtained reduced equation system consists of one differential equation only and describes the time evolution of the independent concentration Ž x 2 .: d x 2rd t s
K 2eq 1 q K 2eq
P Ž Õ 1s y Õ 3s . .
Ž 2.6 .
At steady state the rates of the slow reactions 1 and 3 are equal Ž cf. Eq. Ž 2.2.. . The control coefficients of the enzymes of the slow reactions can be determined by differentiation of their rates in Eq. Ž2.2. Žthis method of calculating the control coefficients goes back to the work of Crabtree and Newsholme w34x. . Expressing in the rates Õis < st. st. in Eq. Ž2.2., the dependent concentration x 1 via the independent concentration x 2 according to Eq. Ž2.4., one obtains, by the chain rule of differentiation,
E ln Jj E pi
s
ž
E ln Õjs E ln x 1 E ln x 1 E ln x 2
q
E ln Õjs E ln x 2
/
E ln x 2 E pi
q
E ln Õjs E pi
,
i , j s 1,3.
Ž 2.7 .
Here and in the following, for simplicity’s sake, we omit the bars indicating the absolute value of reaction rates. From Eq. Ž2.4. it follows that d ln x 1rd ln x 2 s 1. Therefore, dividing Eq. Ž2.7. by E ln ÕisrE pi gives, due to definition Ž 1.1. and Eq. Ž1.2., CiJj y Ž ´ 1j q ´ 2j . P Cix 2 s d ji ,
i s 1,3, j s 1,3,
Ž 2.8 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
341
where d i j is the Kronecker symbol:
d ji s 1 if i s j, otherwise d ji s 0.
Ž 2.9 .
The four equations contained in Eq. Ž2.8. Žfor all i, j . can be written in the form of the following matrix equation:
1
y Ž ´ 11 q ´ 21 .
1
y Ž ´ 13 q ´ 23 .
0ž P
C1J1 C1x 2
C3J1 1 s 0 C3x 2
/
ž
0 . 1
/
Ž 2.10 .
To compare this equation to Eq. Ž2.3. it is convenient to permute the matrices E and C Žthis is always possible since their product is the identity matrix. :
ž
C1J1 C1x 2
1 C3J1 x2 P C3 1
/
y Ž ´ 11 q ´ 21 . y Ž ´ 13 q ´ 23 .
0
s
ž
1 0
0 . 1
/
Ž 2.11 .
This equation contains the summation and connectivity theorems for the control coefficients calculated in the framework of quasi-equilibrium approximation. The products of the rows of matrix C and the first column of matrix E show that the summation theorems derived for the original system still apply provided that the control coefficients of the fast reaction 2 are assumed to be equal to zero. The products of the rows of matrix C and the second column of matrix E represent the ‘‘new’’ connectivity relationships which correspond to simultaneous equal relative changes in the concentrations x 1 and x 2 so that their equilibrium ratio remains constant Ž cf. the perturbation method in w12x.. Note that the elasticity coefficients of the fast reaction which tend to plus and minus infinity, respectively, do not enter the new relationships. Importantly, one can regard Eqs. Ž 2.10. and Ž2.11. as corresponding to a simplified pathway involving reactions 1 and 3 only. As has been derived earlier by Fell and Sauro w35x, the effective elasticities of reactions 1 and 3 are the sums of the elasticities of these reactions with respect to x 1 and x 2 .
3. General analysis: Fundamentals We shall now consider a metabolic pathway of any complexity, involving n reactions and m metabolites. We shall assume that due to differences in the activities of enzymes or in the rate constants of uncatalysed reactions, the rates Ž Õif . of some very fast reactions Ž f in number. greatly exceed the rates Ž Õjs . of the n y f remaining, slow reactions Ž for any metabolite concentrations except a small region, see below.: < Õjs < < < Õif < for all i , j Ž i s 1, . . . , f ; j s f q 1, . . . ,n . .
Ž 3.1 .
Note that the numbering in Scheme 1 differs from this general numbering for the sake of clarity. As the reaction rates depend on concentrations, one should indicate in which region of the concentration space inequality Ž3.1. is to hold. Since the fast reactions are assumed to relax very rapidly to quasi-equilibrium states, where Õif f 0, we restrict the validity of relation Ž3.1. to a domain which does not include the region defined by Õif s 0 Žthe ‘‘slow submanifold’’. and a small neighbourhood around it. This neighbourhood includes the region where the steady-state Eq. Ž 2.2. holds. In other words, due to large differences in the kinetic parameters, some reactions have much larger rates than others except in situations where the concentrations adapt to these differences in such a way that a steady state is nevertheless possible. Let q be the number of metabolites which participate in the fast reactions and possibly in the slow reactions. Their concentrations will be designated by x if , i s 1, . . . , q. The remaining m y q metabolites Ž m is the number
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
342
of all metabolites. participate only in the slow reactions, and their concentrations will be designated by x is , i s q q 1, . . . , m. The number q of metabolites x if is usually greater than the number f of fast reactions. With these notations one may write the differential equations for the time evolution of metabolites: d x if
f
s
dt d x is
n
Ý
Ni fj P Õjf q
js1
Ý
Ni fsj P Õjs ,
i s 1,2, . . . ,q,
Ž 3.2a .
jsfq1
n
s
dt
Ý
Nisj P Õjs ,
i s q q 1,q q 2, . . . ,m.
Ž 3.2b .
jsfq1
Here Ni fj is the stoichiometric coefficient of metabolite X if in the jth fast reaction Ž j s 1, . . . , f .. We have if X if is a substrate, and Ni fj ) 0 if X is a product of the reaction. Ni fsj is the stoichiometric coefficient of metabolite X if in the jth slow reaction Ž j s f q 1, . . . , n.. Nisj is the stoichiometric coefficient of metabolite X is in the jth slow reaction Ž j s f q 1, . . . , n.. Below we shall use matrix notation due to its conciseness. We gather the reaction rates Õjs and Õif in the column vectors U and W, respectively. In matrix form Eqs. Ž3.2a. and Ž3.2b. read: Ni fj - 0
ž
f dQrdt W sNP s N dJrdt U 0
/ ž / ž
N fs P W . U Ns
/ž /
Ž 3.3 .
s s Here Q s Ž x 1f , x 2f , . . . , x qf . T and J s Ž x qq1 , x qq2 , . . . , x ms . T are the vectors of metabolite concentrations. The superscript T denotes transposition of vectors or matrices. N is the stoichiometric matrix of the entire pathway. It has the dimension m = n and it is composed by the q = f matrix N f of the fast reactions and by the m = Ž n y f . matrix ŽN fs N s . T of the slow reactions. When the conditions for applicability of Tikhonov’s theorem are fulfilled, then the fast variables obey the following steady-state equations w20,21x:
N f P W ( Q ) s 0.
Ž 3.4 .
Note that besides the usual continuity requirements, the conditions for applicability of Tikhonov’s theorem state that the fast system should have an asymptotically stable steady-state solution with the initial concentration values lying in its domain of attraction w20,21x, cf. w22,43x. Let r denote the rank of matrix N f. Then f y r columns of this matrix can be expressed in terms of r linearly independent columns. Accordingly, there is a full-rank, f = Ž f y r . matrix that fulfills the equation N f P ff s 0.
Ž 3.5 .
The columns of the matrix ff form a basis of the kernel of the matrix N f. In many reaction systems, the matrix N f is full-rank Ž r s f . Žsee the above example and the models in w16,18,24x.. In this case, ff s 0, and Eq. Ž 3.4. implies W s 0 Žstrictly speaking, this equation as well as Eq. Ž3.4. hold in a fast time scale only.. This means that the fast subsystem attains a quasi-equilibrium state after some initial transient period. The mass action ratios of the fast reactions are then nearly equal to the equilibrium constants of these reactions. This holds true, in particular, at any steady state of the pathway. In other words, the affinities Ž i.e., negative Gibbs free energy changes. of those reactions are close to zero: q
Ý Nifj P ln x if y ln K jeq s 0, is1
j s 1, . . . , f .
Ž 3.6 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
343
Using matrix notation one may write T
Ž N f . P ln Q y ln K eq s 0.
Ž 3.7 .
Here K eq s Ž K 1eq , K 2eq , . . . , K feq . T is the vector of equilibrium constants of the fast reactions. Logarithms of vectors are defined via the logarithms of their components. Eq. Ž 3.7. represents a linear system of f equations in the q logarithmic concentrations ln x if. When the rank, r, of matrix N f is less than f, the quasi-steady-state condition Ž3.4. may nevertheless entail the quasi-equilibrium condition W s 0, if the values of K ieq obey some restrictions. For the usual mass action kinetics, Volpert and Khudiaev w36x and Feinberg w37x derived a necessary and sufficient condition for the identity of steady state and equilibrium. Schuster and Schuster w26x proved this condition to be valid also for a generalized mass-action kinetics, which comprises all reversible enzymatic rate laws. To illustrate this condition we note that a linear dependence between the columns of matrix N f often implies the existence of a cycle in the scheme of fast reactions. Alternatively, it can represent the existence of a ‘‘tree’’ connecting source and sink metabolites w37x. At equilibrium the net flux through any cycle vanishes. This implies the existence of a constraint on the possible values of the equilibrium constants of reactions within a cycle. For a simple cycle in which at any steady state all rates must be equal, this constraint is the well-known ‘‘detailed balance’’ relationship or Wegscheider condition w38,39x:
Ł K ieq s Ł kqi rkyi s 1, i
Ž 3.8 .
i
y Ž . where kq i and k i are the forward and backward rate constants, respectively. The product in Eq. 3.8 is taken over all the reactions participating in the cycle. Note that the concentrations of ligands which are assumed to be constant must be included into the equilibrium constants and rate constants. Consider, e.g., the reaction system shown in Scheme 2. The Wegscheider condition for the cycle composed of the three fast reactions reads
k1 k 2 k 3 ky1 ky2 ky3
s 1.
Ž 3.9 .
To generalize the condition given by Eq. Ž 3.8. for a pathway of arbitrary stoichiometry one can premultiply Eq. Ž3.7. by the matrix Ž ff . T. Then, using Eq. Ž3.5. one arrives at the condition for the equivalence of steady state and equilibrium state of the fast subsystem w26,36,37x:
Ž ff .
T
P ln K eq s 0.
Ž 3.10 .
Under this condition, r linearly independent equations of Eq. Ž3.7. allow one to express r Ždependent. concentrations of metabolites participating in the fast reactions in terms of the remaining q y r concentrations chosen as independent concentrations. To do so one may choose a submatrix, NRf , of N f by using r linearly independent columns Ž without loss of generality we shall consider that the reactions are renumbered such that
Scheme 2.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
344
the first r columns are linearly independent.. Then one may replace the f equations in Eq. Ž3.7. by r linearly independent equations only:
Ž NRf .
T
P ln Q s ln K eq R .
Ž 3.11 .
Here K eq R is the reduced vector of the equilibrium constants of those reactions which correspond to the reduced matrix NRf . To construct the system of equilibrium relationships Ž3.11. , one should cancel one reaction for each cycle and one reaction for each branch of a ‘‘tree’’ Ž such reactions whose deletion does not decrease the rank of the stoichiometry matrix of the fast subsystem.. When the equilibrium state Ž defined by Eq. Ž 3.7.. exists, it is globally asymptotically stable for thermodynamic reasons w26,40,41x. This implies that there cannot exist any other stable steady state for the same moiety conservation quantities and that the stability hypothesis in Tikhonov’s theorem is justified. It is convenient to express the reduced q = r matrix NRf as the product of a q = r link matrix, Lf , w6x and a square r = r matrix, Nrr, which consists of r linearly independent rows of the matrix NRf : NRf s Lf P Nrr ,
Lf s
Ir
ž / Lf0
.
Ž 3.12 .
The Ž q y r . = r matrix Lf0 expresses the dependent q y r rows of the matrix into NRf its r independent rows. The link matrices for the whole matrix N f and for the reduced matrix NRf are identical because NRf results from N f by cancelling linearly dependent columns. When r - q, there exist conservation relations in the fast subsystem which hold in a fast timescale only. It is important to note that this does not contradict our above assumption that conservation relations in the entire system Ž if any. only comprise metabolites participating in slow reactions. In the fast subsystem brought about by cancelling all slow reactions, generally additional conservation relations occur. Substituting Eq. Ž3.12. into Eq. Ž 3.11. and choosing the first r concentrations of the metabolites participating in the fast reactions as the dependent concentrations, Qdep s Ž x 1f , x 2f , . . . , x rf . T, and the last q y r concentrations f f as the independent concentrations, Q ind s Ž x rq1 , x rq2 , . . . , x qf . T one may write after some rearrangement: T
ln Qdep s y Ž Lf0 . P ln Q ind q Ž Nrr .
ž
T y1
/
P ln K eq R .
Ž 3.13 .
Eq. Ž3.13. expresses the first r concentrations of metabolites participating in the fast reactions in terms of the equilibrium constants of these reactions and of the last q y r concentrations chosen as the independent concentrations. Let h be the number of independent conservation relations in the slow subsystem. These relations can be expressed by J s Ls J ind q const,
Ž 3.14 .
where Ls is the link matrix of the slow subsystem w6x. This matrix links the stoichiometry matrix N s to a submatrix of linearly independent rows. J ind is the vector of ‘‘slow’’ metabolites independent with respect to conservation relations.
4. Control analysis in the framework of quasi-equilibrium approximation As we wish to carry out a steady-state control analysis, we consider not only the fast reactions of the network to attain a steady state, but also the slow reactions to do so. As mentioned in Section 1 we shall assume that all moiety conservation constraints valid in the original network Ž if any. only involve metabolites not participating in fast reactions.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
345
At any steady state all rates of slow reactions can be expressed via a set of linearly independent slow rates. Let g denote the number of these independent slow rates. This number is given by the number of slow reactions Ž n y f . minus the number of independent steady-state equations with respect to metabolite concentrations. The latter number is given by the number of metabolites Ž m. minus the number of independent quasi-equilibrium conditions Ž r . minus the number of independent conservation relations in the slow subsystem Ž h.. Therefore, we have g s n y f y m q r q h.
Ž 4.1 .
As for Scheme 2, e.g., n s 5, f s 3, m s 3, r s 2 Žwhen the generalized Wegscheider condition is fulfilled in the cycle., and h s 0. This gives, correctly, g s 1. The two slow reactions actually carry the same flux at steady state. We shall designate by JR the Žreduced. vector consisting of g linearly independent steady-state rates of slow reactions T
s s s JR s Ž Õ fq1 ,Õ fq2 , . . . ,Õ fqg . < st .st .
Ž 4.2 .
Žwithout loss of generality it is assumed that the slow reactions are renumbered such that at steady state the first g slow rates are linearly independent. . The vector U of slow rates can be expressed in terms of the reduced vector JR as follows Žcf. w6,8,33x, and Eq. Ž2.2..: J s U
Ž 4.3 .
Here the vector of the slow steady-state rates is designated by J, f is an Ž n y f . = g matrix, which can be considered as the null-space matrix of the reduced reaction system. Its first g rows coincide with the g-dimensional identity matrix fs
Ig
ž / f0
.
Ž 4.4 .
We now calculate the g = Ž n y f . matrix C UJR containing the control coefficients over the fluxes of the independent slow reactions. Using Eq. Ž 4.3. and definition Ž 1.1., we obtain
˜ UJR . C UJ s fC
Ž 4.5 .
˜ relates to the matrix f of Eq. Ž4.3. as follows: The Ž n y f . = g matrix f ˜ s Ž d ln Urd ln JR . < st .st .s Ž diag U . < st .st .P f P Ž diag JR . . f y1
Ž 4.6 .
Here and below Ždiag F. designates Žfor any vector F. the square matrix with diagonal elements equal to the components Fi and all off-diagonal elements equal to zero. We employ the usual notation Ž d ln Ard ln B. for the matrix Žd ln A ird ln Bj . consisting of the derivatives of each component of the vector ln A with respect to each component of the vector ln B. On the other hand, the matrix of control coefficients, C UJ , of the slow reactions can be calculated by differentiation of the steady-state rates of the slow reactions considered as functions of the metabolite concentrations and the parameters pi Žsee Eq. Ž2.7. and w8x.:
E ln J Ep
s
ž
E ln U
E ln Qdep
E ln Qdep E ln Q ind
E ln U q
E lnQ ind
/
E ln Q ind Ep
E ln U q
E ln J E ln J ind
E ln J E ln J ind
Ep
E ln U q
Ep
.
Ž 4.7 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
346
Since the system is to contain no conservation relations involving the components of Q, there are no dependencies between the vectors Q and J and no dependencies between the components of Qdep and Q ind other than the quasi-equilibrium relations. Owing to Eq. Ž3.13., we have
Ž d ln Qdeprd ln Q ind . s y Ž Lf0 .
T
.
Ž 4.8 .
Therefore, y Ž Lf0 . Ž d ln Qrd ln Q ind . s I qyr
T
0
˜ r. sL
Ž 4.9 .
˜ r in analogy to the link matrix L Žsee Eq. Ž3.12... Most importantly, We designate this q = Ž q y r . matrix L this matrix does not depend on the steady-state concentrations, nor even on the equilibrium constants of the fast reactions. Dividing Eq. Ž4.7. by E ln UrE p gives T
C UJ s y´ QUdep P Ž Lf0 . q ´ QUind P C QU ind q ´ U J
ž
/
E ln J E ln J ind
P C JU ind q I nyf .
Ž 4.10 .
ind Ž q y r . = Ž n y f . matrix of the control coefficients over the independent Here C Q and C J U U are the Ž . concentrations Q ind of metabolites participating in fast reactions and the Ž m y q . = Ž n y f . matrix of the control coefficients over the concentrations of metabolites participating only in slow reactions, respectively. By Eq. Ž3.14. and the rules for differentiation of logarithms, we conclude
E ln J E ln J ind
s Ž diag J .
y1
EJ E J ind
y1 Ž diag J ind . s Ž diag J . Ls Ž diag J ind . .
Ž 4.11 .
The control coefficients over the vector J can therefore be expressed by the control coefficients over the concentrations J ind independent with respect to conservation relations as follows: CJ U s diag Ž J .
y1
P Ls P diag Ž J ind . P C JU ind .
Ž 4.12 .
After defining the Ž m y h. = Ž m y r y h. matrix,
˜s L
y Ž Lf0 .
T
0
I qyr 0
0 diag Ž J .
y1
Ls diag Ž J ind .
0
Ž 4.13 .
Eq. Ž4.10. can be written as
˜ P C UX ind q I nyf . C UJ s ´ xU P L
Ž 4.14 .
The Ž m y r y h.-dimensional vector X ind is generated from the vector X by cancelling the r concentrations Qdep dependent due to the quasi-equilibrium conditions and the h concentrations J dep dependent due to the conservation relations in the slow subsystem. ´Ux is the Ž n y f . = Ž m y h. matrix of the elasticities of the slow reaction rates ŽU. with respect to all the pathway concentrations except for J dep , ´Ux s ´UQdep < ´UQ ind < ´UJ ind .
ž
/
Ž 4.15 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
347
C UX ind is the Ž m y r y h. = Ž n y f . matrix of the control coefficients over all the independent concentrations of the reduced system. This matrix is defined as d ln Q ind d ln U C UX ind s
ž / ž /
dp
y1
dp
d ln J ind d ln U dp
y1
dp
0
s
C QU ind
ž / C JU ind
.
Ž 4.16 .
Combining Eqs. Ž4.5. and Ž 4.14. one arrives at:
ž
˜ < y ´UX P L ˜ P f
/
C UJR
ž / C UX ind
s E P C U s I nyf .
Ž 4.17 .
By inverting the matrix E, one can express the control coefficients of the slow reactions in terms of the elasticities of the slow reactions, without knowledge of the elasticities of the fast reactions. On the other hand, by inverting the matrix C U , one can calculate the elasticity coefficients of the slow reactions in terms of their ˜ and ye Ux P L ˜ have dimensions Ž n y f . = g and Ž n y f . = Ž m y r y h., and control coefficients. The matrices f JR X ind the matrices C U and C U have dimensions g = Ž n y f . and Ž m y r y h. = Ž n y f ., respectively. Therefore, Eq. Ž4.1. guarantees that the two matrices on the left-hand side of Eq. Ž4.17. are square. From these general formulas, we can readily derive the results for the example in Scheme 1. We have m s q s 2, n s 3, f s 1. The stoichiometry matrix of the fast reactions reads N f s Žy1 1. T with the submatrix ˜ s Ž1 Nrf s y1. Then, Q s Ž x 1, x 2 . T, Q ind s x 2 . The link matrix reads Lf s Ž 1 y1. T, hence Lf0 s y1 and L T T ˜ ˜ ˜ 1. . Since the pathway is unbranched, f s Ž1 1. . Substituting the matrices f and L into Eq. Ž4.17. one can again derive Eq. Ž2.10..
5. An illustration We shall illustrate our method applying it to the metabolic pathway shown in Scheme 3. The system involves six internal metabolites Ž m s 6. and nine reactions Ž n s 9., three of which are considered fast Ž f s 3. . There are five metabolites involved in fast reactions Ž q s 5. . If we write T
X s Ž X 1f , X 2f , X 3f , X 4f , X5f , X6s . , T
v s Ž Õ 1f ,Õ 2f ,Õ 3f ,Õ4s ,Õ5s ,Õ6s ,Õ 7s ,Õ 8s ,Õ 9s . ,
Scheme 3.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
348
the dynamics of the system is described by the equation dX dt
s N P v,
Ž 5.1 .
where N is the stoichometric matrix
y1 0 0 Ns 1 0 0
0 y1 0 1 0 0
0 0 1 0 y1 0
0 0 0 0 0 1
1 0 0 0 0 y1
1 0 y1 0 0 0
0 0 0 y1 0 0
0 1 0 0 0 y1
0 0 0 . 0 1 y1
0
Ž 5.2 .
For the fast subsystem:
y1 0 Nfs 0 1 0
0 y1 0 1 0
0 0 1 . 0 y1
0
Ž 5.3 .
As the rank of N f is 3; NRf s N f , and the equation N f s Lf P Nrf gives us 1 0 f Ls 0 y1 0
0 1 0 y1 0
0 0 1 , 0 y1
0
Nrf s
y1 0 0
0 y1 0
0 0 . 1
0
Ž 5.4 .
For this case, n s 9, f s 3, m s 6, r s 3, h s 0, so the number of independent fluxes g is 3. We choose J4 , J5 and J6 : J 7 s J4 ; J8 s J4 y J5 y J6 ; J9 s J6 ; JR s Ž J4 , J5 , J6 . T and U s Ž J4 , J5 , J6 , J 7 , J8 , J9 . T. 1 0 ˜s 0 f 1 J4rJ8 0
0 1 0 0 yJ5rJ8 0
0 0 1 , 0 yJ6rJ8 1
0
Ž 5.5 .
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
°´ ´
4 4 1 2
¶
´ 34´ 44´ 54´ s4
1 1 ˜s 0 L 1 0 0
´ 15´ 25 ´ 35´ 45´ 55´ s5 ´UX s
´ 16´ 26 ´ 36´ 46´ 56´ s6 ´ 17´ 27 ´ 37 ´ 47´ 57´ s7 9 9 1 2
0 0 1 0 1 0
0 0 0 . 0 0 1
0
,
´ 18´ 28 ´ 38´ 48´ 58´ s8
¢´ ´
349
ß
´ 39 ´ 49´ 59´ s9
Ž 5.6 .
˜ . and inverting the result we can obtain the matrix Multiplying the elasticity matrix Ž ´Ux . by the link matrix ŽL of control coefficients:
°C
¶
J4 4
C5J4
C6J4
C 7J4
C8J4
C9J4
C4J5
C5J5
C6J5
C 7J5
C8J5
C9J5
C4J6
C5J6
C6J6
C 7J6
C8J6
C9J6
C4x 4
f
C5x 4
f
C6x 4
f
C 7x 4
f
C8x 4
f
C9x 4
C4x 5
f
C5x 5
f
C6x 5
f
C 7x 5
f
C8x 5
f
C9x 5
x 6s 4
C5x 6
s
C6x 6
s
C 7x 6
s
C8x 6
s
C9x 6
¢C
f
f
ß
s
°1
0
0
y Ž ´ 14 q ´ 24 q ´ 44 . y Ž ´ 34 q ´ 54 . y ´ s4
0
1
0
y Ž ´ 15 q ´ 25 q ´ 45 . y Ž ´ 35 q ´ 55 . y ´ s5
0
0
1
y Ž ´ 16 q ´ 26 q ´ 46 . y Ž ´ 36 q ´ 56 . y ´ s6
1
0
0
y Ž ´ 17 q ´ 27 q ´ 47 . y Ž ´ 37 q ´ 57 . y ´ s7
J4rJ8
yJ5rJ8
yJ6rJ8
y Ž ´ 18 q ´ 28 q ´48 . y Ž ´ 38 q ´ 58 . y ´ s8
0
1
y Ž ´ 19 q ´ 29 q ´49 . y Ž ´ 39 q ´ 59 . y ´ s9
s
¢0
¶
y1
.
Ž 5.7 .
ß
As in the simple example dealt with in Section 2, the effective elasticities of slow reactions are sums of elasticities with respect to the metabolites participating in fast reactions. One of these sums runs over the metabolites X 1f , X 2f and X 4f forming one component of the fast subsystem. The other sum runs over X 3f and X5f , i.e., over another component of this subsystem.
6. Discussion In the above presentation, we have shown how the control coefficients of slow reactions over the fluxes through these reactions can be calculated in terms of the elasticities of the slow reactions and vice versa. So we have solved what is sometimes called the direct and inverse problems of MCA w33,42x for the situation that some reactions subsist in quasi-equilibrium. The method presented is of considerable help in the control analysis of metabolic systems characterized by time hierarchy because it allows one to draw conclusions about the control properties even if the kinetic properties of the fast reactions are imperfectly known. Despite the recent availability of powerful computers, kinetic and control analysis of large reaction networks is often hampered or completely impeded by this incompleteness in knowledge of parameters. Our analysis concerns quasi-equilibrium reactions, i.e., reactions which still carry a non-zero flux although their mass-action ratios nearly equal their equilibrium constants. A related problem is to calculate control
350
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
coefficients for true equilibrium reactions for which the net flux is zero when the whole system is at steady state. For the case of adenylate kinase this problem was tackled by Fell and Sauro w35x. Eqs. Ž4.13. and Ž 4.17. show that neither the kinetic parameters, nor the elasticities, nor the equilibrium constants of the fast reactions need to be known. Since the matrix Ž d ln X deprd ln X ind . is a matrix of constant coefficients Žit depends neither on metabolite concentrations nor on equilibrium constants. we need to know only network stoichiometry and the elasticity coefficients of the slow reactions at the operating steady state. Note that f . the derivative ŽdX frdX ind , which would be needed for calculating non-normalized control coefficients does depend on the equilibrium constants. This is an advantage of normalized Ž i.e., logarithmic. control coefficients. For further advantages and drawbacks of normalization see w7x. One can see that upon derivation of Eq. Ž4.17. for determining the control properties we avoid to consider ‘‘slow variables’’, also called ‘‘pool variables’’ w2x Ž these were considered in the example in Section 2 only. . Using ‘‘pool variables’’ is suitable if one wishes to consider the time evolution of the system in the framework of the quasi-equilibrium approximation w24,43x. Moreover, it may be helpful to use these variables when pathways with arbitrary conservation relations are studied. The partitioning of concentrations into independent and dependent ones is of importance for introducing any non-linear constraints among the concentrations of metabolites Ž cf. w29,44x.. For expressing control coefficients in terms of elasticities, it is then necessary to calculate the matrix Ž d ln X deprd ln X ind . of derivatives of dependent concentrations with respect to the independent ones. For the study of in vivo situations, it is worthwhile combining the analysis presented here with other extensions of MCA, such as top-down or modular approaches Ž see, e.g., w45,46x.. In the analysis of metabolite channelling w47,48x the elimination of quasi-equilibrium steps can also be of great help.
Acknowledgements This work was supported by the Research grant of Comision ´ Interministerial de Ciencia y Tecnologıa, ´ Spain ŽCICYT. PM 96-0099, by Programa Catedra Ž FBBV. , Spain by the Netherlands Organization for Scientific ´ Research and by the E.U. grant B104-CT96-0662.
References w1x R. Heinrich, S.M. Rapoport, T.A. Rapoport, Metabolic regulation and mathematical models, Prog. Biophys. Mol. Biol. 32 Ž1977. 1–82. w2x J.G. Reich, E.E. Selkov, Energy Metabolism of the Cell, Academic Press, London, 1981. w3x S. Minakami, H. Yoshikawa, Studies on erythrocyte glycolysis. II. Free energy changes and rate limiting steps in erythrocyte glycolysis, J. Biochem. ŽTokyo. 59 Ž1969. 139–144. w4x D.A. Fell, Metabolic control analysis: A survey of its theoretical and experimental development, Biochem. J. 286 Ž1992. 313–330. w5x R. Heinrich, T.A. Rapoport, A linear steady-state treatment of enzymatic chains. General properties, control and effector strength, Eur. J. Biochem. 42 Ž1974. 89–95. w6x C. Reder, Metabolic control theory. A structural approach, J. Theor. Biol. 135 Ž1988. 175–201. w7x S. Schuster, R. Heinrich, The definitions of metabolic control analysis revisited, BioSystems 27 Ž1992. 1–15. w8x B.N. Kholodenko, D. Molenaar, S. Schuster, R. Heinrich, H.V. Westerhoff, Defining control coefficients in non-ideal metabolic pathways, Biophys. Chem. 56 Ž1995. 215–226. w9x H. Kacser, J.A. Burns, The control of flux, Symp. Soc. Exp. Biol. 27 Ž1973. 65–104 Župdated reprint by: H. Kacser, J.A. Burns, D.A. Fell, Biochem. Soc. Trans. 23 Ž1995. 341–366.. w10x J. Delgado, J.C. Liao, Control of metabolic pathways by time-scale separation, BioSystems 36 Ž1995. 55–70. w11x H.V. Westerhoff, K. van Dam, Thermodynamics and Control of Biological Free-energy Transduction, Elsevier, Amsterdam, 1987. w12x B.N. Kholodenko, Control of molecular transformations in multienzyme systems: quantitative theory of metabolic regulation, Molek. Biol. ŽMoscow. 22 Ž1988. 1238–1256 Žin Russian.; English trans. in: Molec. Biol. ŽMoscow. 22 Ž1998. 990–1005.
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
351
w13x A.N. Zhabotinsky, Autonomous Oscillations in Concentrations Žin Russian., Nauka, Moscow, 1974. w14x E.E. Selkov, Stabilization of energy charge, generation of oscillations and multiple steady states in energy metabolism as a result of purely stoichiometric regulation, Eur. J. Biochem. 59 Ž1975. 151–157. w15x B.N. Kholodenko, A study of the dynamic behavior of the glycolytic system, Biofizika 24 Ž1979. 640–645 Žin Russian.; English trans. in: Biophysics ŽUSSR. 24 Ž1979. 659–665. w16x B.N. Kholodenko, A.V. Stepanchicova, L.I. Ehrlich, F.I. Ataullakhanov, A.M. Zhabotinsky, A.V. Pichugin, The role of the 2,3-diphosphoglycerate bypass in the regulation of energetic metabolism in human erythrocytes Žin Russian., Izv. Acad. Sci. USSR ŽBiol.. 2 Ž1985. 196–205. w17x H.-G. Holzhutter, G. Jacobasch, A. Bisdorff, Mathematical modelling of metabolic pathways affected by an enzyme deficiency, Eur. ¨ J. Biochem. 149 Ž1985. 101–111. w18x R. Schuster, H.G. Holzhutter, G. Jacobasch, Interrelations between glycolysis and the hexose monophosphate shunt in erythrocytes ¨ as studied on the basis of a mathematical model, BioSystems 22 Ž1988. 19–36. w19x G. Pettersson, U. Ryde-Pettersson, Dependence of the Calvin cycle activity on kinetic parameters for the interaction of non-equilibrium cycle enzymes with their substrates, Eur. J. Biochem. 186 Ž1989. 683–686. w20x A.N. Tikhonov, On the dependence of the solutions of differential equations on a small parameter Žin Russian., Mat. Sborn. 22 Ž1948. 193–204. w21x A.N. Tikhonov, Systems of differential equations containing small parameters in their derivatives Žin Russian., Mat. Sborn. 31 Ž1952. 575–586. w22x W. Klonowski, Simplifying principles for chemical and enzyme reaction kinetics, Biophys. Chem. 18 Ž1983. 73–87. w23x V.M. Vasilev, A.I. Volpert, S.I. Khudiaev, On the method of quasi-equilibrium concentrations for chemical kinetic equations Žin Russian., Zhurn. Vychysl. Mat. Matemat. Fiz. 13 Ž1973. 683–697. w24x M. Schauer, R. Heinrich, Quasi-steady-state approximation in the mathematical modeling of biochemical reaction networks, Math. Biosci. 65 Ž1983. 155–171. w25x J.C. Liao, E.N. Lightfoot Jr., Lumping analysis of biochemical reaction systems with time scale separation, Biotechnol. Bioeng. 31 Ž1988. 869–879. w26x S. Schuster, R. Schuster, A generalization of Wegscheider’s condition. Implications for properties of steady states and for quasi-steady-state approximation, J. Math. Chem. 3 Ž1989. 25–42. w27x B.N. Kholodenko, H.V. Westerhoff, The macroworld versus the microworld of biochemical regulation and control, Trends Biochem. Sci. 20 Ž1995. 52–54. w28x H.V. Westerhoff, Y.-D. Chen, How do enzyme activities control metabolite concentrations? An additional theorem in the theory of metabolic control, Eur. J. Biochem. 142 Ž1984. 425–430. w29x B.N. Kholodenko, Metabolic Control Theory. New relationships for determining control coefficients of enzymes and response coefficients of system variables, J. Nonlin. Biol. 1 Ž1990. 107–126. w30x H.M. Sauro, J.R. Small, D.A. Fell, Metabolic control and its analysis. Extensions to the theory and matrix method, Eur. J. Biochem. 165 Ž1987. 215–221. w31x H.V. Westerhoff, D.B. Kell, Matrix method for determining steps most rate-limiting to metabolic fluxes in biotechnological processes, Biotechnol. Bioeng. 30 Ž1987. 101–107. w32x M. Cascante, R. Franco, E.I. Canela, Use of implicit methods from general sensitivity theory to develop a systematic approach to metabolic control. I. Unbranched pathways, Math. Biosci. 94 Ž1989. 271–288; II. Complex systems, Math. Biosci. 94 Ž1989. 289–309 w33x H.V. Westerhoff, J.-H.S. Hofmeyr, B.N. Kholodenko, Getting to the inside of cells using metabolic control analysis, Biophys. Chem. 50 Ž1994. 273–283. w34x B. Crabtree, E.A. Newsholme, The derivation and interpretation of control coefficients, Biochem. J. 247 Ž1987. 113–120. w35x D.A. Fell, H.M. Sauro, Metabolic control and its analysis. Additional relationships between elasticities and control coefficients, Eur. J. Biochem. 148 Ž1985. 555–561. w36x A.I. Volpert, S.I. Khudiaev, Analysis in Classes of Continuous Functions and the Equations of Mathematical Physics Žin Russian., Nauka, Moscow, 1975 w37x M. Feinberg, Necessary and sufficient conditions for detailed balancing in mass action systems of arbitrary complexity, Chem. Eng. Sci. 44 Ž1989. 1819–1827. ¨ w38x R. Wegscheider, Uber simultane Gleichgewichte und die Beziehungen zwischen Thermodynamik und Reaktionskinetik homogener Systeme, Z. Phys. Chem. 39 Ž1902. 257–303. w39x J.Z. Hearon, The kinetics of linear systems with special reference to periodic reactions, Bull. Math. Biophys. 15 Ž1953. 121–141. w40x Ya.B. Zeldovich, Proof of a unique solution to the mass action law Žin Russian., Zh. Fiz. Khim. 11 Ž1938. 685–687. w41x F. Horn, R. Jackson, General mass action kinetics, Arch. Rational Mech. Anal. 47 Ž1972. 81–116. w42x B.N. Kholodenko, H.M. Sauro, H.V. Westerhoff, M. Cascante, Coenzyme cycles and metabolic control analysis: the determination of the elasticity coefficients from the generalized connectivity theorem, Biochem. Mol. Biol. Int. 35 Ž1995. 615–625.
352 w43x w44x w45x w46x
B.N. Kholodenko et al.r Biochimica et Biophysica Acta 1379 (1998) 337–352
R. Heinrich, S. Schuster, The Regulation of Cellular Systems, Chapman and Hall, New York, 1996. S. Schuster, Control analysis in terms of generalized variables characterizing metabolic systems, J. Theor. Biol. 182 Ž1996. 259–268. M.D. Brand, Top down metabolic control analysis, J. Theor. Biol. 182 Ž1996. 351–360. S. Schuster, D. Kahn, H.V. Westerhoff, Modular analysis of the control of complex metabolic pathways, Biophys. Chem. 48 Ž1993. 1–17. w47x B.N. Kholodenko, H.V. Westerhoff, Metabolic channelling and control of the flux, FEBS Lett. 320 Ž1993. 71–74. w48x B.N. Kholodenko, M. Cascante, H.V. Westerhoff, Control theory of metabolic channelling, Mol. Cell. Biochem. 143 Ž1995. 151–168.