J. Phys. Chem. 1993,97, 8316-8327
8316
Global Analysis of Enzyme Inhibition Kinetics Marc R. Roussel and Simon J. Fraser' Chemical Physics Theory Group, Department of Chemistry and Scarborough College, University of Toronto, Toronto, Ontario, Canada M5S I AI Received: March I. I993
Simple enzyme mechanisms are capable of a surprisingly rich variety of behavior. This paper presents a global analysis of this behavior for several inhibition mechanisms. For instance, in competitive inhibition one can vary the system parameters (typically total inhibitor concentration); then the ordinary differential equations describing this mechanism can display two bifurcations which separate the relaxation behavior into three distinct regimes. In one of these regimes, the usual steady-state treatment does not describe the dynamics, even as a first approximation. These different classes of relaxation kinetics are analyzed in detail, and methods are presented for the extraction of the underlying invariant manifold structure. The global analysis of three related inhibition mechanisms is also discussed.
1. Introduction
= -k,es
This paper examines the flow geometry of the ordinary differential equations (ODEs) of transient and steady-state enzyme inhibition kinetics. While these mechanisms have been treated rigorously by singular perturbation theory,' this method does not provide an obvious geometrical picture of the system's time evolution. And though there is a rudimentary geometry underlying standard kinetics analysis, the associated algebra is formal and lacks rigorous justification. In contrast to previous predictions, we find that simple inhibition mechanisms are dynamically far richer than is commonly supposed-this was missed in earlier studies' which assumed the existence of a welldefined steady state and so inevitably failed to detect new features that occur in the absence of this steady state. While the standard algebraic methods of steady-state kinetics analysis do not detect complex behavior, we show that a geometrical approach clearly demonstrates these features. We begin with an examination of the competitive inhibition mechanism.* We can show that this mechanism exhibits three distinct kinds of relaxation kinetics, depending on the values of the rate constants and total enzyme and inhibitor concentrations. An awareness of these different scenarios is necessary for an understanding of steady-state and transient experiments; for instance, in one experimentallyaccessible case there is no global steady state in the classical Briggs-Haldane3 sense.
+ (k-l + kp)c- k,ei + k-,h
+ k-lc i. = kles - (k-l + k,)c
S = -k,es
p = kpc i = -k,ei
+ k-,h
h = k,ei - k-,h
(1.1)
Because of mass conservation, not all of the concentrations are independent;for instance, e + i. + h = 0, so e c + h is a constant of the motion for this s y ~ t e m .These ~ particular linear conserved quantities are total species concentrations;for species X,we will denote this concentration (chemical mass) by m,.6 The full set of total concentrations for this mechanism is
+
me = e
+c+ h
(1.2) (1 -3)
mi=i+h m, = s
+c+ p
(1.4)
The competitive inhibition mechanism is Using these total concentrations,the six-dimensional system, eqs 1.1, can be reduced to a set of three ODES in 6!/3!3! = 20 ways. Conventionally, the kinetics is expressed in terms of s, i, and c:
+ i ) + k-,c I= -k,i(m, - c - mi + i ) + k-,(mI - i)
S = -k,s(m, ka
E+I=H
(CI)
k-1
- c - mi
i.= k , s ( m , - c - m , + i ) - ( k - , There are six chemical species in this scheme (E, S,C = ES, P, I, and H = EI). We can immediately write a set of autonomous ODEs for the timeevolution oftheconcentrationsof these species;4 let lower case letters represent concentrations, Le., x = [XIwith k = dx/dt, then 0022-3654/93/2097-83 16$04.00/0
+ kp)c
(1.5)
As this is a closed autonomous system,' it can further be reduced to a system of two trajectory equations.* Let x' = dx/ds = k/i, then 0 1993 American Chemical Society
Global Analysis of Enzyme Inhibition Kinetics C c‘=-=
The Journal of Physical Chemistry, Vol. 97, No. 31, 1993 8317
s(m,-c-m,+i)-Kfl +(me - c - mi i) + K s
+ i -i(me - c - mi + i ) + K,(m, - i) i’ ---=a s -s(m, - c - mi + i) + K f i
(a)
2.90
(1.6) c (PM)
where the equilibrium, Michaelis, and other constants are 2.80
KE = k-l/kI
i (PM)
0.00 0.03
0.00
(Y
= k,/k,
(We use these parameter definitions throughout the paper, including section 6,which deals with other mechanisms. Note that the K‘s have concentration units but a is a pure number.) The five rate constants in eqs 1.5 have been replaced by four time-independentparameters in eqs 1.6. By reducing the evolution equations to trajectory equations, we have lost the time scale, and so the number of parameters has decreased by The usual way to restore the time scale is to write the reaction velocity as u = p = kpc
0.00
I
so that if c = c(s) is known exactly, the first-order rate constant
kp implies the time dependence of the system’s evolution. The set of all solutions of the system of ODEs of eqs 1.6 is a trajectory flow” in the phase space r = s X i X c . ~The system of ODEs may be stiff for some values of the rate constants; i.e., the flow it generates naturally separates into modes with very different time scales.12 Steady-state experiments are performed on the time scale of the slowest mode of a stiff system? (Other interesting approaches which use a modal analysis are those of Maas and Pope” and of Palsson et &I.’*) For three-dimensional flows, we have previously identified the dynamical steady state with a line-likeinvariant manifold, and the slowest transient with a two-dimensional surface invariant manifold, 2.16717 (An invariant manifold is a smooth, i.e., infinitely differentiable, set which maps into itself under the flow.18) We also showed that the usual King-Altman appr~ximationl~ is the line-like intersection in r of the nullcline surfaces corresponding to all enzyme complexes.17JO (A nullcline is a surface which is defined by the vanishing of the corresponding time derivative in the set of eqs 1.1 or 1.5.2’) In experiments, the total inhibitor concentration is varied as a means of probing mechanism, i.e., mi is used as a control parameter. Correspondingly, in our theoretical study of the CI mechanism, we use mi as a bifurcation parameter so that, mathematically speaking, the flow changes in a well-defined way as mi passes through a bifurcation value. In fact, the series of trajectory flows in Figure 1 shows how flow structure changes as mi is altered. Figure l a shows trajectories for the CI mechanism at an intermediate value of mi with other constants set so that the flow of eqs 1.6 is stiff or strongly coupled and hence the relaxation structure is clearly resolved.22 Three relaxation modes are discernible: 1. Fast relaxation from an arbitrary point toward the “attracting” surface 2 is one mode. In order to correspond to experiment, this starting point must lie in r+,the physically realizable part of I”. 2. Intermediate rate relaxation within 2 toward the curve .M is a second mode. Here 2 corresponds to the slowest transient. 3. Slow relaxation along the line-like manifold .M (the true dynamical steady state) toward the equilibrium point P23 is the third mode.
(c)
3.2C
i (mM) s (mM)
Figure 1. (a) Relaxation to equilibrium in the tightly coupled case for the CI mechanism (me = 2.90 pM, mi = 2.80 PM, KE = 12.3 pM,KS = 32.3 mM, KI = 1.92 pM, a = 369). Trajectories first relax to the surface Z,then to the slow manifold A, and finally approach the equilibrium point P. Note that Z includes the i axia near P, as demonstrated by the behavior of the trajectory starting from (s,i,c) = (O,mi,m,) at top, back, left corner. (b) Relaxation to equilibrium in the sluggish-inhibitioncase (me = 9.20 pM, mi = 5.00 pM, KE = 913 nM, Ks = 16.3pM, KI= 19.5 pM, (Y = 0.0478). For some values of the rate constants,for small mi the slowestdirectionof approach to the equilibrium Pis the i axis rather than the usual‘hyperbolic” King-Altman trajectory. (c) Relaxation to equilibrium in the fast-inhibition case (all constants as in part b, except mi= 3.20 mM). Note that Z does not include the i axis, cf. part a. This is evident from the trajectory starting from the point (O,m,,m,) at top, back, left corner. This scenario is always obtained for sufficiently large mi.
At fhe end of each relaxation step, the flow velocity decreases abruptly and the trajectory bends sharply.16J7 Note that Z:and M . are special invariant manifolds which “attract” the flow. Figure 1b shows trajectories for the CI mechanism for small values of mi. Although there is a well-defined surface 2, the slow manifold, at least near I ) , is the i axis. This local structure for .Mis not consistent with a Briggs-Haldane (King-Altman) picture of the kinetics3J9 since the slowest mode cannot be approximated by the intersection of the c and h nullclines, i.e., the line where C = h = 0.20 Correspondingly, the CI mechanism itself implies that if both kz and k-2 are sufficiently small, the catalytic step is complete before the inhibition step even begins so there is
Roussel and Fraser
8318 The Journal of Physical Chemistry, Vol. 97, No. 31, 199‘3
essentially no inhibition. However, any CI system which displays the sluggish inhibition shown in Figure l b (with k2 and k-2 both nonzero) can be made to behave like the system of Figure l a by adding a sufficient amount of inhibitor. This crossover occurs at a definite value of mi for fixed values of the other parameters. Figure IC shows trajectories for the CI mechanism for large values of mi. In this case, 2: becomes steeplypitchedcorresponding to fast inhibition, in contrast to small or medium mi where Z is more or less horizontal. This change in relaxation does not affect steady-stateexperimentsbut is important in transient experiments. Chemically speaking, for large mi, the equilibration of i with h is very fast compared with the catalytic step for any value of s; geometrically this means that Z is close to the h nullcline. In terms of phase-space coordinates, inhibition experiments can conveniently be started from at least two different sets of initial conditions: 1. The starting point (s,i, c) = (ms,mi, 0 ) is somewhere along the lower horizontal back edge (c = 0, i = mi) of the phase-space “cube” of Figure 1. In a conventional stopped-flow apparatus, this corresponds to premixing solutions of S and I in one syringe and putting the solution of E in the other. As soon as the contents of the two syringes come into contact, the complexes C = ES and H = E1 begin to form. 2. The starting point is (s, i, c) = (ms,i,, 0 ) where i, is the equilibrium amount of inhibitor. This point is a distance m, along the line i = i, in the (bottom) plane c = 0 of the p h a s e space cube of Figure 1. To set up this state experimentally, E and I are premixed in one syringe, while S is held in the other. Thecomplex C forms when the syringe solutions come in contact. These two sets of initial conditions allow definite regions of I’+ to be explored. Other starting conditions might be harder to realize experimentally. Sections 2 through 5 present a rigorous analysis of the kinetics of competitiveinhibition. In section 2, we examine the parameter ranges displaying each type of relaxation behavior and study the ways in which bifurcationsarise. Insection 3, the range of validity of the classical King-Altman approximation is examined. In section 4, a contraction mapping procedure for solving the functional equation for Jt is discussed.1sJ6~20In section 5, a similar contraction mapping method is used to solve the functional equation for 2.’’ In section6, three relatedinhibition mechanisms are briefly examined. Finally, in section 7,wediscuss our findings.
One can understand the relaxation processes, leading to equilibrium by studying the velocity field v = ( S , i, k ) and its Jacobian matrix J(x) = av/& wherex = (s,i, c). The eigenvalues and eigenvectors of J = J(x) describe the time evolution of small displacements with respect to a representative trajectory passing through x:4 the absolute value of the real part of each eigenvalue gives the corresponding reciprocal relaxation time while the eigenvectors system gives the principal axes in the co-moving frame. The rate equations 1.5 give
[
J= 0 k1e
-k,s -k,(e k1s
+ i ) - k-,
kls + k-1 k2i -k,s - k-, - kp
3
where e = me - c - mi + i. The characteristic polynomial for J is
(J - AI( = -A3 - A2kl[e+ s + K , + a(e + i + K,)] Ak12(a[e(e+ i + s + K , + KI) + K,(i + K I ) + KF]+ e(K, - K E ) ]- k13ae(K,- KE)(e i KI) (2.1)
+ +
Although the roots of this cubic can be expressed analytically,24 the solution cannot be interpreted simply. However, the characteristic polynomial does factor at the equilibrium point P. The
s,
=0
c,
=0
P, = m, i, = ‘f2{mi- me - KI
+ A,’/’)
h, = ‘f2{me+ mi + KI- A,’”) e, = ‘f2{m,- mi - KI
+ Ai/’)
where the discriminant
AI = (me
+ m, + KI)2- 4 m p i
is always positive. At equilibrium, eq 2.1 reduces to
+ Ak,(e, + K,) + k12(Ks- K&,I
IJ, - AI1 = -(A - X,)[A2 where
A, = -k2(i,
+ e,)
- k-,
(2.2)
so that k = $-e,
-K,
f [(e,
- K,)’
+ 4e,KB]’/2) (2.3)
Note that 1x11 C 1x31and that A1 and A, are both real for positive values of the rate constants and total concentrations. Substituting A2 back into the equation J 4 2 = A2u2, we find the eigenvector u2 = (O,l,O)Twhere T indicates transpose,25 so u2 lies along the i axis. The other two eigenvectors are more complicated. We now interpret Figure l a - c in relation to the properties of J. Suppose that A2 is the smallest eigenvalue; then the final (slowest) approach to P would be along u2, Le., along the i axis. This is precisely the situation shown in Figure lb. On the other hand, if A2 were the largest eigenvalue, the fastest (transient) mode of the flow would be parallel to the i axis as in Figure IC. Figure la, the strongly coupled case, corresponds to the ordering bil
2. Local Linear Analysis of the CI Flow
-k,e
full set of equilibrium concentrations is found by setting v = 0 in eqs 1.5 to get
< lA21 < hd.
The linearized picture yields even more information. As ml vanishes, we get lim m
d
A1,3
k, = +me 2
- K , f [(me- K,)’ + 4m&]
’”)
(2.4)
(Note that these are theeigenvaluesobtained near theequilibrium point of the Michaelis-Menten mechanism,lOJ“ the expected result in this limit.) Also lim A, = -k2me - k-, mr-0
In contrast, for large mi A,
-
(2.5)
O(m;’)
Since A2 grows fastest with mi,the fast inhibition case displayed in Figure IC is always obtained for sufficiently large mi. Also, if kzm, or k-2 is large, only fast inhibition is observed, regardless of the value of mi. However, if both k2meand k-2 are sufficiently
Global Analysis of Enzyme Inhibition Kinetics
The Journal of Physical Chemistry, Vol. 97, No. 31, 1993 8319 is almost fixed. At bifurcation, mi = mi2), there is a single eigenvector at P,namely, the i axis, belonging to the doubly degenerate largest eigenvalue. For mi > mj2), as the system approaches P,the eigenvector belonging to the middle eigenvalue is almost fixed and the eigenvector belonging to the largest eigenvalue rotates onto the i axis.
s
(/N
JUU
Figure 2. Surface.8 on which the twosmallest eigenvalues of the Jacobian for CI mechanism are equal (all constants as in Figure lb,c, except that mi = 10.0 pM). Inside 8, the eigenvalues are complex, and the slow manifold is undefined. The slow manifold is the i axis near the equilibrium P (0). On the far side of e, remote from P, the slow manifold is the usual 'hyperbolic" trajectory.
small, then for small mi, (A21 < 1x11 (eqs 2.4 and 2.5), so sluggish inhibition (Figure lb) is observed. Therefore for small k2m, and k-2, mi can be used as a control parameter in a sequence of experiments displaying all three relaxation scenarios. (There is an intermediate case in which the rate constants are such that either tightly coupled kinetics or fast inhibition is possible, but not sluggish inhibition.) Suppose that a given enzyme has rate constants such that all three dynamical regimes can be observed simply by changing mi.26 Imagine increasing mi from zero and consider the global changes in the velocity field v = v(x;mi) in r. At small values of mi, the slowest eigenmode is along the i axis near P,but at very large s the slowest mode is the usual KA-like curve. These two disconnected slow-mode segments, comprising &( CZ), are separated by a folded surface 8,shown in Figure 2, on which the cubic discriminant of eq 2.1 vanishes.24 On 8 the smallest eigenvalues X1 and A2 are degenerate, and inside 8 they are complex conjugates. The eigenvalue degeneracy for J arises by the usual defective-eigensystem mechanism, so that on 8 there is only one eigenvector belonging to the doubly degenerate eigenvalue XI, say. The segment of & coming from large s stops at the boundary 8 n 2;on 2 inside 8 there is no proper steady state because there is no unique smallest eigenvalue-instead, the system's motion is intrinsically two-dimensional and is described locally by the complex eigensystem on 2. Trajectories eventually emerge from this region and approach P along the segment of & lying along the i axis. Let us now increase mi to the first bifurcation value mi1)where X1 = X2 at P (eqs 2.2 and 2.3) and the boundary 8 n 2 just touches P. (Since the eigenvalues at P cannot be complex, P cannot move inside 8.)As mi increases above mjl), the slowest direction of approach to P rotates away from the i axis. However, 8 n 2 still separates the segment of & joined to P from the segment of & at large s. When mi is sufficiently large, these two segments join where & is tangent to 8 n 2; for larger mi, & is globally defined on 2 outside 8. The appearance of small-magnitude complex eigenvalues in sluggish inhibition which breaks & into two pieces is quite different from the appearance of large-magnitude complex eigenvalues associated with transient spiraling about &.17.27-29 Oscillatory transients are not possible for the CI m e ~ h a n i s m . ~ * . ~ ~ J ~ A second bifurcation occurs at large mi = mi2) >> mi'); this corresponds to the motion in mrspace of the fast-relaxation eigenvectors near P . 2 6 (Except at P itself, the corresponding eigenvalues are distinct so the surface 2 is a defined global structure in r.) For mi < mi2),as the system approaches P,the eigenvectors belonging to the middle eigenvalue coalesces with the i axis while the eigenvector belonging to the largest eigenvalue
At both MI1) and mi2) there is a discontinuity in the rate at which the slower mode of a pair moves in mrspace; such bifurcations do not dramatically alter flow structure but are mathematically sharp and lead to experimentally observable changes in relaxation dynamics. Transient relaxation in the strongly coupled and fast inhibition cases (parts a and c of Figure 1) can be described in terms of the eigenbasis of J. In the strongly coupled case, the eigenvalue ordering is lXll < lX2l < lX31 with {ul, u2, us} as the corresponding right (column) unit eigenvectors. (For fast inhibition, the last two eigenvalues and eigenvectors are just relabeled.) Figure 3a shows the time dependence of the eigenbasis decomposition of v along a trajectory in the strongly coupled case.32 Far from 2, v is nearly parallel to u3. Very rapidly, the trajectory approaches 2;v turns and becomes roughly parallel to u2. Finally, v becomes approximately collinear with u1 and approaches P along this vector.33 These eigenvectors remain distinct as they move during the time evolution of the system. In the sluggish inhibition case, relaxation at large s is similar to that just described. However, inside the surface 8 (Figure 2), the eigensystem for the coupled slow modes is complex, so the eigenvalue ordering is lost and the system "drifts" (on 2). On emerging from 8,the system must again relax toward the (unique) slow eigenvector. This process is illustrated in Figure 3b. The geometry of a given relaxation mode dictates its corresponding algebraic form: e.g., the surface 2 can be represented as a scalar function c = c*(s, i) of two independent variables,17 whereas & is a curve (indeed a trajectory)15J6so its equation is a vector function (i = i*(s),c = c*(s)) of a ~inglevariable'~ s that determines all other concentrations in the steady ~ t a t e . 3This ~ important point is often missed in steady-state analysis.
3. Classical King-Altmen Approximation in CI Let us now examine the classical treatment of the CI mechanism. It is generally assumed that after a short time i. becomes small and h comes to equilibrium with e and i. Total enzyme conservation (eq 1.2) provides one further relation. The algebraic equations formally expressing these constraints are
+ c + h = me -k,ei + k-,h = 0 kles - (k-, + kp)c = 0 e
(3.1)
Solving this set of equations for c gives ,C
=
m2
s
+ Ks(l + i/KJ
(3.2)
The subscript KA denotes the classical King-Altman (KA) approximation for the CI mechanism.19 But this approximation, taken literally, cannot represent the steady state since it defines a surface c = c(s, i ) while the steady state, as pointed out previously, is a curve (i = i*(s), c = c*(s)). The problem is that eqs 3.1 do not include total inhibitor conservation. However, combining the set 3.1 with eq 1.3 correctly gives a curve
Roussel and Fraser
8320 The Journal of Physical Chemistry, Vol. 97, No. 31, 1993
1.6. (Steady-state expressions superficially similar to eqs 3.3 have appeared before,)$ but like eq 3.2, they do not represent curves as they should-often the equation obtained defined a point instead!) As eq 3.2 has been used successfully for decades to fit competitiveinhibition data, we now try to give it a precise meaning. Equation 3.3 expanded to first order in me gives
Since iK = mi to @meo),the approximation 3.3 is identical to the classical expression 3.2. The perturbation series in the Appendix shows that eq 3.4 is correct to O(me)for c(s) on &, and to U(me2) this series gives the steady-state formula for c(s) as
+ a ~ s ( m i / ~ I+) [~ss ( +1 ~ ~ / K ~ I I D / (aKr[s + ~ , ( 1 + m , / ~ , ) I ~ ) l m , 2 where c = kp/kl = KS - KE. The coefficient of me2 here is a ~ I (+ 1 ~,/KI)'I
\ \
.-
\
I
_ /,
,
.-(
\
@my') for large mi (but differs from the corresponding coefficient for CK(S)). Therefore, compared with the exact formula for c(s) on A, the relative error in the standard approximation CU(S) (eq 3.2) evaluated at i = mr steadily decreases as mi m for KIfixed. (This agrees with the experimental practice of adding increasinglylarge amounts of inhibitor.) Careful analysis shows that the classical expression is correct for this mechanism under two different conditions: mr >> me, or s >> me. Typically however, me is comparable to other concentrations in vivo36 so that work at physiological concentrations requires better approximations. Furthermore, some texts on enzyme kinetics37 recommend that inhibition experiments be performed in a concentration regime (low s, moderate mi) where the usual approximation to the reaction velocity does not apply. Unless sufficiently low substrate concentrations are used, competitive inhibition might never be observed, and unless values of m,smaller than KI are used, the variation in the rate with inhibitor concentration might not be sufficient to determine the type of inhibition obtained.
-
4. Functional Equation for the Slow Manifold & in CI
The linearized picture described in section 2 can be related to an inverse-iteration method for finding the attracting invariant manifolds A and Z9JoJ5J6,20 For example, the eigenvector at x belonging to the smallesf magnitude eigenvalue of J(x) is reached, from almost any starting vector VO(X), by iterating the inverse matrix equation VH+I(X) = J-~(X)V,,(X).~* A generalized inverse iteration lies at the heart of the numerical method of Maas and Pope for constructing local linearized representations of the attracting invariant manifolds.'3 Similarly, writing the trajectory eq 1.6 for c, i, and s as
i' =fli;s,c) c' = g(c;s,i) where
A(s) = [Ks(me + m, + KI) + K,s]' - 4K&,mi Thesubscript Kdistinguishes this approximation from the related KA function. Although the steady-state approximation (SSA) of eq 3.3 is a curve, it is still not a solution of the trajectory eq
allows formation inversion offwith respect to iand gwith respect to c to produce a solution which is a quadratic in i:
i =J'(i';s,c) c = g-'(c';s,i)
Global Analysis of Enzyme Inhibition Kinetics
The Journal of Physical Chemistry, Vol. 97, No. 31, 1993 8321
The solution space of the ODES has not been changed by this inversion, but (like the algebraic equations of local linearization) the iteratively stable solution of the vector functional eq 4.1 correspondsglobally to the slow time evolution of the system, Le., tothetrajectoryn, withequation(i= i*(s),c=c*(s))?O Starting from (io(s),co(s)), eqs 4.1 can be iterated according to the scheme
to obtain (i*(s),c*(s)). However, the quadratic discriminant of eq 4.1 may be negative in regions of I’ where A2 is well defined. Therefore, working in the phase space s X i X c of the standard dynamical variables is unsatisfactory for finding A iteratively. Instead, we use thes X i X e representation,for which the trajectory equations are j‘=
S
e
+
(e-)
-ei KI(mi- i ) -es KE(m,- e - mi
i
-=a -=
-e(s
s
+ + i) + ai) + Ks(m, - e - mi + i ) + aKI(mi- i ) -es + KE(m,- e - mi + i ) (4.2)
The solution of the pair of trajectory eqs 4.2 is a vector-valued function (i(s),e(s)), the equation for JM being a particular instance. Invertingthese trajectory equationsproduces thevector functional equation
+ aK,mi + KI) (K,- K&(m, - mi + i ) + aKI(mi- i) e= s + K, + ai- e’(s + K E ) i=
w (s-1) Figure 4. Eadie plot of the classical King-Altman approximation (KA, -), steady-state approximation (SSA; - - -), and 4th iterate of the SSA for the CI mechanism (me = 8.90 pM, mi = 8.90 pM,KE = 2.27 pM,Ks=28.0pM,Kr=10.6pM,a=0.758,kl=6.60X 107M-1s-1). The three approximations to the slow manifold merge at large s (small w = v / s ) but are quite different at large s: the SSA and its iterates, along with the slow manifold itself, alwayslie above the KA for this mechanism.
iqes - KE(m,- mi - e ) ]
KEi’+ a(e
(4.3)
in which no radicals appear. This vector functional equation is the basis of an iterative procedure for generating a convergent sequence of approximate formulas for JM.lOJ6.20 Write eq 4.3 as i = F(i’,e;s), e = G(e’,i;s);then iterations converge most rapidly to the fixed-point function ( i = i*(s), e = e*(#)) for A in the ‘chained” form
The other chain, e&l = G(in;in,s),i,+l = F(in;e,+l,s),converges equally rapidly to ( i = i*(s), e = e*(s)).39 In contrast, “parallel” iteration, = F&;e,,s), e,,+l = G(e,;i,,s), converges more slowly.’O.~O This kind of functional mapping converges quickest from an asymptotically good initial approximation like the steady-statetrajectory formula of eq 3.3, represented in s X i X e; this curve lies near at large s.20 The Eadie plot, showing the scalar reaction velocity u = p us “frequency” w = u / s , is a very good way of displayingsteady-state enzyme kinetics data.g Theu X w phase plane (unlike ‘concentration” phase space) depends on the time scale of the problem. Eadie plots have the advantage of emphasizing systematic deviations from KA-like hyperbolic rate laws.9 Figure 4 shows an Eadie plot of the KA approximation (eq 3.2) evaluated at i = mi, the SSA curve of eqs 3.3, and the 4th iterate of the SSA under eqs 4.3 for the CI mechanism for a case where the dynamics is tightly coupled. The numerically computed slow manifold40 is indistinguishablefrom the 4th iterate
The slow manifold is indistinguishablefrom the 4th iterate of the SSA on the scale of this figure.
of the SSA on the scale of this figure. Iteration of the fast inhibition case behaves similarly. As suggestedabove,stabilityof the functionaliterationis related to dynamical stability of the manifold which we are trying to locate.41 In particular, it is necessary for the manifold to be locally “attracting” in order for its explicit formula to be a stable fixed point of the functional equations. However, this condition is not sufficientbecause the functional iteration can diverge slowly at smalls although the underlying dynamics is strongly coupled. Even in these divergent examples, iteration continues to converge at sufficiently large s, whatever the parameter values. This asymptotic iterative stability is connected to the asymptotic dynamical stability of the usual SSA-like slow manifold. We have also found that the iterative scheme only diverges at small s for parameter values very near the bifurcation between sluggish inhibition and strong coupling. Finally, note that when the underlying manifold is globally defined, the divergence is always sufficiently weak that it can be corrected by applying the AitkenSteffensen h2 method.42 Previously, we have used this method to make a convergent sequence from a divergent one by eliminating the linearly growing part;lOJ7conversely,divergenceof theiterative scheme not correctable by the AitkenSteffensen method indicates sluggish inhibition in the CI mechanism, Le., the disappearance or break up of the slow manifold. In eq 4.2, only the expression for e depends on kp, and at kp = 0 this functional equation is identically satisfied by the equilibrium solution
where
+
A,(s) = [KE(m, mi + KI)
+ Kfi]’ - 4K2m,mi
For s 1 0, this solution ( i = i&), e = e&)) at kp = 0 is a line 6 = h ( k p = 0 )of equilibriumpoints. This important result follows from the equations definingthe nullcline intersections. (It implies infinite sets of superstable solutions to eq 4.3 in this case.) Consequently, the Taylor series of eq 4.3 can be written to first
8322 The Journal of Physical Chemistry, Vol. 97, No. 31, 1993
Roussel and Fraser
order in t and its formal iterative solution to O(t) is
where the derivatives in eq 4.4are evaluated a t ( i = is, e = ec) from eq 4 and G&) = G(s;t=O). For large mi, this perturbation formula gives an O(tZ)correction to u = p = eklc. Alternatively, perturbation series solutions in e can be found for eq 1.6 and eq 4.2with i&),ee(s), and&) = c~(s;c=O) as thestartingfunctions (coefficients of €0). The recurrence relations for these series are found by the same method used for the series in me(see Appendix). Note that these formulas are not uniformly applicable because of the breakdown of the usual steady-state approximation for the CI mechanism. Solution of eq 4.2 in the form of asymptotic expansions can be calculated directly. The relation to O(r2) for the fixed-point function c*(s) for .& can be calculated from the asymptotic expansion of (i*(s),e*(s))to this order and gives c*(s)
-
8.00
c (PM)
m.11- Ks( 1 + mi K,); 1+
0.00 10
0 s (PM)
This formula is very good for very large s.43 It follows that CK in eq 3.3 is only correct to O(s-1) and is therefore less accurate at large s than the steady-state result for the Michaelis-MentenHenri mechanism, which is O(S-~).~ Thus, more kinetic information should be accessible from experiments at large s in the CI mechanism since the usual steady-state approximation breaks down sooner.44 5. Functional Equation for the Relaxation Surface Z in CI
In ref 17, we described a method for obtaining the relaxation surface Z of three-dimensional mechanisms. The partial differential equation (PDE) underlying this method is
where e = e(s,i) on the right.45 Inserting the explicit forms for the velocity components, given in eq 4.2, into eq 5.1 produces an equation whose solutions are surfaces containing the flow; i.e., the solutions of eq 5.1 are surfaces everywhere tangent to v. Rearranging this PDE into the form e = $(ae/as,ae/ai;s,i),a scalar functional equation is obtained whose iterative stability properties are similar to those of the functional equation pair 4.3 but whose stable fixed point is Z:46
Figure 5. (a) Approximations to 2 of the CI mechanism for a case where 2 includes the i axis (constants as in Figure 1b). The light surface is the e nullcline; the darker surface is the 4th iterate of the e nullcline under eq 5.2. On the scale of this figure, it is indistinguishablefrom 2. (b) the e nullcline (light) and its 4th iterate (dark) in a case where Z excludes the i axis in the CI mechanism (me = 8.00 pM, mi = 1S O pM, KE = 3.13 pM, Ks = 4.28 pM, KI = 90.9 pM, a = 0.688). Note how 2,which would be indistinguishablefrom c4 on the scale of this figure, is nearly vertical. Figure 5a shows the e nullcline and its 4th iterate under eq 5.2 for the case which we have called sluggish inhibition; the qualitative behavior in the case of tight coupling is identical to this. Convergence of iterates of eq 5.2 to e = e*(s,i),the fixedpoint function representing 2,is rapid.4 Figure 5bshows iterates of the functional equations from eq 5.3 for the fast inhibition case; again, the iterative scheme constructs an excellent approximation to the slow surface with relatively little symbolic labor. Since the two largest eigenvalues never become degenerate in I?+, the iterative mapping obtained from eq 5.2 always converges in the region of interest. At worst, weak divergence has been observed just outside r+. 6. Other Inhibition Mechanisms
We have studied three inhibition mechanisms, similar to the competitiveinhibition mechanism, which are all reducible to threedimensional flows in r. We now summarjze our findings. 1. Competitive Substrate Inhibition (CSI). The CSI mechanism is a modification of the CI mechanism?
The e nullcline, eo = 4(O,O;s,i), i.e.
+ i ) + aKI(mi- i) + K, + ai
KS(me- mi eo(s,i) =
s
(5.3)
isasuitablestartingfunctionfor iterationofeq 5.2. Theenullcline approximates Z in relaxation kinetics in the same way that the SSA approximates .& in steady-state kinetics:47 The e nullcline is very close to Z at large s.
Instead of having a separate inhibitor as in the CI mechanism, the substrate itself can bind nonproductively to the enzyme: presumably, either the substrate has two mutually exclusive binding sites, one of which is inactive, or the substrate can attach itself to the active site in a nonreactive orientation.
The Journal of Physical Chemistry, Vol. 97, No. 31, 1993 8323
Global Analysis of Enzyme Inhibition Kinetics
‘;a’ ’,,
0
0.20
,
.-.,_, ,. 0.30 ,/
I
I
, ‘s* _ ,
,
0.20
0.40
0.30
0.40
t (s)
F v6. LocalJacobianeigenvectordecomposition of the systemvelocity vector v along a trajectory for the CSI mechanism (me= 5.6 pM, KE = 941 nM, Ks = 2.65 pM, KI = 15.7 pM, a = 8.24 X 10-3, kl = 1.70 X
108 M-I 8); the trajectory startingpoint (s,h,c)= ((400,2.00,3.00)pM). (Theinitialrelaxationisnot shown here.) At t = 0.2 s, visnearlycollinear to the eigenvector belonging to the smallest eigenvalue (component indicated by a solid line; the second component is dashed, and the component for the largest eigenvalue is dot-dashed). The eigenvectors of the CSIJacobiancorrespondingto the two smallesteigenvaluesbecome nearly degenerate at t = 0.303 s. After this near collision, these eigenvectorsrotate rapidly into new orientations in roughly 0.01 s; since this is the time scale for relaxation to the slow manifold on Z (0.0040.025 s in this region), v “loses” the slow manifold until the eigenvectors settle down. This is local, mode-mixing behavior. Once the eigenvectors are again moving slowly, v approaches the slow eigenvector again.
This mechanism is three-dimensional: there are five chemical species (E, S,C, P, and H ) but two conservation relations, namely, total substrate and total enzyme. This system can be cast most conveniently in thes X h X c phase space. The functional equations for both 2 and & can be worked out directly from the reduced set of ODES exactly as above. The equilibrium point is trivially found to be (s,,h,,c,) = (O,O,O). The characteristic equation for the eigenvalues of the Jacobian of the CSI mechanism does not factor, even at the equilibrium point. The functional equations for & diverge for some values of the rate constants and total enzyme concentration, and we find numerically that there appears to be the same sort of bifurcation behavior as found in the CI mechanism. Figure 6 shows a (Jacobian) eigenvector decomposition of the velocity along a trajectory; a comparison of this figure with Figure 3b suggests that the two flows are related. Again, we can apply chemical intuition to the dynamics. Suppose that the equilibration involving complex C is fast compared to that involving complex H. Then h remains roughly constant in time through most of the course of reaction; the asymptotic dynamics is dominated by the conversion of S to P through C. However, in the final stages of relaxation H is a reservoir for S and influences the dynamics near equilibrium. The two eigenvectors belonging to the smallest eigenvalues of the Jacobian become nearly collinear just before the apparent bifurcation; at this point, they both point roughly in the (1,0,0) direction, i.e., along the s axis. During the following period, the two smallest eigenvectors rotate rapidly into completely new orientations; as suggested by the qualitative analysis presented above, the eigenvector belonging to the smallest eigenvalue acquires a large component in the h direction. This mechanism is in one way significantly different from the CI model: in the CSI mechanism, the two smallest eigenvalues never cross and the corresponding eigenvectors never become fully degenerate; there is no bifurcation analogous to that encountered in the CI mechanism. The velocity vector comes close to the middle eigenvector due to the short period during which the eigenvectors rotate compared to the time scale for
Figure 7. Differencebetween the KA velocity and the velocity along an integrated trajectory for the CSI mechanism (all constants and starting point as in Figure 6). In cases where there is either mode mixing or a bifurcationin the slowmanifold structure,a plot of the differencebetween the King-Altman velocity and the experimentally-determinedvelocity will show a characteristic kink. This kink should be experimentally resolvable;in this case, the maximum amplitude of the difference is over 20% of the reaction velocity.
relaxation to the slow manifold on the slow surface. We call this mode-mixing behavior. In normal cases (where the rates in the two branches are more nearly comparable), there is no rapid rotation of the eigenvectors; the “normal modes”I4 cannot be understood as being associated with separate branches of the mechanism (except perhaps a t infinite s). Because the eigenvectors change so rapidly, the velocity vector cannot closely follow the direction of the slow eigenvector in the mixing region. In other words, there can only be transient behavior in the mixing region and thus no slow (steady-state) manifold. The measures which we have so far discussed for detecting the breakup of the slow manifold are probably not experimentally practical. It is in principle possible to obtain the eigenvalues and eigenvectors of the Jacobian from perturbation experiments,49 but producing a plot such as Figure 6 from experimental data is probably unrealistically difficult. Similarly, it is difficult to imagine a straightforward way to use the divergence of the functional iteration method to detect bifurcations or mixing. In an earlier paper, we showed that & in a normal chemical relaxation scenariostays near the set of nullcline intersections of the flow.20.50 The KA approximation is a nullcline intersections1 and is constructible from data obtained a t large s . ~It follows that in cases where the slow manifold is a smooth curve, the difference Av = v - VK = k,(c - CK) should also be a smooth function of time; this is indeed what we find from numerical experiments. However, where the slow manifold becomes undefined, the reaction velocity should veer away from the KA curve; this phenomenon is illustrated by Figure 7. 2. Uncompetitive Inhibition (UI).In the UI mechanism, the inhibitor acts by forming a complex with the enzyme-substrate complex C rather than with the enzyme E:2
Due to total enzyme, substrate, and inhibitor conservation, this mechanism is also reducible to a three-dimensional system; we choose s, i, and c as the independent variables. In the corresponding phase space, the equilibrium point is (sq,i,,cq) = (0,miS)).
8324
Roussel and Fraser
The Journal of Physical Chemistry, Vol. 97, No. 31, 1993
We can show that the UI mechanism displays exactly the same types of relaxation as the CSI mechanism: for most values of the rate constants, the usual relaxation sequence, I’+(Cl?) 2 At P is followed. However, when the substrate is converted sufficiently quickly to product in the catalytic step, the complex B acts as a reservoir of substrate until very late in the reaction, at which point the character of the relaxation changes abruptly. The slow manifold then becomes undefined for small s through the defective eigensystem mechanism described for CI. 3. Uncompetitive Substrate Inhibition (USI).In the presence of large amounts of substrate, the productive complex C = ES can combine with a second substrate molecule to form a nonproductive complex B = ES2.2,3’ This US1 mechanism2 is dynamically interesting because catalysis is slowed at small [SI and effectively stopped at very large [SI. Cornish-Bowden remarks that although this effect is not likely to occur to any great extent at physiological concentrations it is a useful diagnostic tool for distinguishing between possible reaction path~ays,3~8 and so we examine it carefully. The US1 mechanism should be contrasted with the previously discussed CSI mechanism in which both the productive complex C and the nonproductive complex H can form between E and a single S molecule. The US1 scheme is
--
-
kp
kl
E+S+C+P+E k-i k2
(UW
C+S*B k-2
This mechanism reduces to a three-dimensional system due to total substrate and total enzyme conservation. The evolution equations are
S = -k,s(m, - c - b)
+ k-,c - k2cs + k-,b
b = k2cs - k-,b t = k,s(m, - c - b) - (k-,
+ kp)c - k2cs + k-,b
(6.1)
and the functional equations for At, in thes X b X c representation, are
+ arcs] + arcs b’(s + aKI) + aK, c’[s(m,- b) - aK#] + s(m, - b) + arK# C= + KE- + s + K, + CIS b=
b’[s(me - c) - K f i
C’(S
CYS)
w (s-1)
Figure 8. &die plot of the KA approximation (-) and its 4th iterate for the US1 mechanism (m, = 3.90 pM, under iteration of q 6.2 KE = 1.13 pM, Ks = 5.13 pM, KI = 2.02 mM, CY = 0.01 11, kl = 8.00 X lo8 M-l s-l). Note that the reaction velocity approaches 0 both for (.e.)
larges (smallw)andnear theequilibriumpoint(wapproachingitsmaximal value) for this mechanism. The KA approximationunderestimatesthe ‘maximumvalue of w obtained by about 4596, the 4th iterate returns this intercept within 1.5%. In general, any measure of the accuracy of representationof the slow manifold is improvedgeometricallyby functional iteration from order to order. of the KA line with the s nullcline surface. Consequently, these equilibria occur at (O,O,O) as well as a t the “points” lim-,(s,bK(s),cK(s)), but only the origin is stable. The characteristic equation of the Jacobian for this mechanism is
X 3 + X 2 k l [ a r ( c + ~ + K I ) + s + K , + m , - c - b ]+Xk,’X {ar[(c + K,)(s + K,) + K# + s2] + (me-C - b)[a(2s + KI) + K, - K E ] )+ kI3a(Ks- KE)[KAme-c - b ) - cs] = 0 Unlike the other systems studied here, the coefficients of this equation are not positive everywhere within I’+: the constant term vanishes on the surface
me - c - b - cs/K, = 0 At such points, the smallest eigenvalue vanishes. At the equilibrium (O,O,O), this term is positive and so the eigenvalues are all negative; on the other hand, since b and c are bounded by me, this term always becomes negative as s QD so that the smallest eigenvalue crosses zero exactly once on At. Thevelocity on At vanishes as s m, so on Jt a t large s the reaction velocity u = p first increases a s s decreases, then reaches a maximum near the maximum of c*(s) (where the smallest eigenvalue is zero), and finally decreases toward equilibrium. The eigenvalues of the Jacobian near the equilibrium point can be obtained analytically for this mechanism:
-
(6.2)
At corresponds to the stable fixed point (b*(s),c*(s)) of eq 6.2. The line-like KA approximation for this mechanism is the intersection of the b and c nullclines, which gives
A2
A1,3
CK
=
m2
s(l
+ s/K,) + K,
-
When S is in excess, the nonproductive complex B removes the productive complex C through the reaction C S B so that as s m the (scalar) reaction velocityp = kfi*(s)approaches zero. This can be seen in Figure 8, which shows Eadie plots of both the KA approximation and its 4th iterate under the functional eqs 6.2. The system 6.1 has “equilibrium” points a t the intersections of the three nullcline surfaces, Le., the intersection
-
+
-
kl
= -k-2
= -{-me- K, f [(me - K,)’ 2
+ 4m&]’/’)
The eigenvaluesA1 and A3 are identical to those for the equilibrium of the Michaelis-Menten mechanismlOJk because, near equilibrium, the system 6.1 reduces to the two-dimensional MichaelisMenten system and a decoupled (kinetically insignificant) inhibition step. Correspondingly, the eigenvectors belonging to AI and A3 both have b components of 0. This behavior follows from the “quadratic” dependence of b on s implied by the form of the inhibition step. At the beginning of the reaction, when free E and S are mixed, C = ES must form before B = ES2 in the initial (fast) relaxation to 2. Toward the end of the reaction,
Global Analysis of Enzyme Inhibition Kinetics
The Journal of Physical Chemistry, Vol. 97, No. 31, 1993 8325
when s is small, b(s) = O(s2) on At as long as k-2 is not too small, and this leads to the loss of any inhibition effect. Since A1 and do not depend on k-2, there can also be bifurcations of the flow near the equilibrium point. It turns out that the first bifurcation (corresponding to the crossing of X1 and A2 at the equilibrium point) leads to behavior of the mixing type rather than to the formation of a degenerate surface as in the CI mechanism. In other words, the eigenvalues only become degenerate at the boundary of the physical region and not in its interior; it is a rapid change in the orientation of the eigenvectors that is responsible for the loss of the slow manifold. Finally,note that transient oscillationsanalogous to those found in the extended Michaelis-Menten (EMM) mechanism can be observed in USI. In the EMM mechanism, the free enzyme, ES, and EP complexes form a Wangle” of reacting enzyme species which allows the two largest eigenvalues to collide and become ~ o m p l e x . l ~ JIn~ -this ~ ~regime, relaxation to the slow manifold is oscillatory for some physically realizable range of concentrations. In the US1 mechanism, the appearance of oscillatory transients is associated with a curved degenerate surface under which the eigenvalues are complex. (This situation should be compared to the quasi-planar flow encountered in the CI mechanism, where the two smallest eigenvalues can become complex.) At the bifurcation value of k-2, this surface touches the equilibriumpoint; for values of the constants above this critical value of k-2, the slow manifold must puncture this surface to reach the equilibrium point. This is not the only way in which this second bifurcation can occur, however. The degenerate surface is sometimes very flat and (e.g., in the s X b X c phase space) nearly vertical; in such cases, it just touches the equilibrium point and does not pass into I’+ as the critical value of k-2 is passed. This behavior is similar to that of the second bifurcation at mi2)in the CI mechanism. In the US1 mechanism, the iterative stability of the functional equations can be reduced to the study of a planar map near equilibrium, so that when At exists locally it takes the form
+ O(s3) c = ys + O(s2)
b = Os2
(6.3)
Substituting eqs 6.3 into eqs 6.2 yields a mapping for 0 and y:
To lowest order, y (and thus c) is decoupled from B (b) near equilibrium. Since the fixed point (B*,y*) of eqs 6.4 is invariant, Le., @* = f ( @ * , y * )and y* = g(y*), we can represent the (small) displacement vector from (B*,y*) by
and consider the dynamics of this vector under the linearization of eqs 6.4 at (B*,y*), Le., the trajectories of the planar map
-
This system can easily be solved.52 One can show that 0 < ag/ ay* < 1 so that, under the map, by,, 0 geometrically fast with multiplier (eigenvalue) ag/ay*. Therefore, if 0 < af/aB* < 1, under the map bg,,- 0 with multiplier max(af/a#?*,ag/ay*), and
the fixed point (j3*,y*)is stable. In fact, we can show
so iteration of the local model 6.4 only converges when X1/X2 < 1/2. This is a necessary condition for convergence of the full iterative scheme based on eqs 6.2. Consequently, for 1/2 < X,/X2 < 1, JM exists but the iterative scheme does not converge; however, numerical work suggests that Aitken’s method removes this divergence. 7. Discussion
Throughout a series of papers (refs 9,lS-17, and 20),we have examined chemical kinetics from a global geometricalperspective. Generally speaking, there are several reasons for wanting a correct and complete description of the kinetics of any given process: 1. Natural processes may not always operate in experimentally convenient parameter ranges, and dynamically important effects may be missed or extreme effects overemphasized by working outside the. physiological range.36 2. Correct interpretation of experimental results requires knowledge of the range of validity of whatever approximation is being used. In principle, for very accurate theory and results, all rate constants can be determined. 3. New experimental design for a system requires an understanding of the overall dynamical behavior for a postulated mechanism. 4. Geometrical description in r provides a simple but rigorous way to teach steady-state (or transient) kinetics. The last point deserves comment. Steady-state kinetics is invariably presented in a self-contradictory way: On the one hand steady-state formulas are used to describe motion across nullclines where one or more velocity components vanishes; on the other hand, the same formulas are used to describe motion along the slow manifold in integrated rate law formulas. These two situations are mathematically quite distinct, and endless confusion (and indeed published discussion) arises from failing to distinguish them clearly. In the phasespace description, this distinction is obvious. Specifically, this paper has focused on two dynamical aspects of enzyme inhibition kinetics: the bifurcation analysis of the phase-space flow and the invariant manifold structure associated with the relaxation modes of these enzyme mechanisms. This two-pronged approach gave us a geometrical picture of the global dynamics of a number of inhibition mechanisms. The case of competitive inhibition was studied in detail. Here the existence of a dynamical steady state for all substrate concentrations was found to depend on total inhibitor concentration, which appears as a bifurcation parameter in the analysis and is also a natural control parameter in experiments. At any particular value of the system parameters, the dynamics of this system reduces to the problem of extracting time-independent curves and surfaces in r since the flow is “steady”, A complete analysis of this kind is necessary to avoid the ambiguities of the standard description (section 3). First, the dimensionalityand geometry of the objects associated with the transients and steady state are usually not considered or are patently confused.19.35 Second, the important bifurcations associated with these dynamical structures have not been discussed in connectionwith ordinary batch kinetics. Third, oversimplificationof the standard analysisobscures the dynamical richness that is easily apparent in the phase-space picture. In view of these deficiencies, we believe that the standard analysis should be applied cautiously or with some additional evidence for a well-defined dynamical steady state. The principal methods which we have developed for extracting invariant manifolds (either At or 2) rely on inversion of the equations describing the flow (trajectory equations or PDE for trajectory sheets), Inversion has also been used either explicitly
Roussel and Fraser
8326 The Journal of Physical Chemistry, Vol. 97, No. 31, 1993 or implicitly in other methods for solving kinetic equations. For instance, Poland constructs convergent series representing chemical ev0lution:5~noting that the inverse function t ( s ) of the time dependence of concentration s ( t ) normally contains no movable singularities,Poland showsthat series solutions for t ( s ) can readily be obtained starting from any point of I’. However, this does not selectively represent the slow manifold or relaxation behavior. In comparison, from the fixed-point formula we can directly obtain expressions for the time to traverse any segment of More obviously, the method of Maas and Pope” is a local numerical analog of our global analytic method. The MaasPope method can easily be extended to higher-dimensionalsystems sinceit is relativelyeasy and efficient to treat linear reconstructions of local tangent spaces. While it is clear how to extract the 1and (n - 1)-dimensionalmanifolds, we have not worked with the intermediate equations in this hierarchy-although they are easy to write forma1ly.*gs4 These two methods are in some respect complementary. The iterative methods we have used in this series of p a p e r ~ ~ J ~ 1can ~ J be O thought of as functional equation analogs of iteration methods for the solutionof algebraic or trancendental equations.55 However, the spectrum associated with the linearization of these functional equations about their fixed point is continuous, not discrete as in the simpler examples; convergence difficulties can sometimes be traced to this difference. Stability and uniqueness of the solution of our method are ensured by the Banach-Cacciopoli fixed-point theorem56when the linearization about the fixed point is a contraction mapping. By solving an “inverse” problem, we contrive convergence to “slow” (smalleigenvalue) structures in the system. Finally, we mention the work of Constantin et al., who have developed methods for extracting low-dimensional manifolds generated by PDEsS7 Their method uses a different approach. Rather than constructing the manifolds directly, they construct the boundary of the manifold andobtain thisobject by integration. Their method is intended for infinite dimensional systems. In contrast, our method automatically determines the asymptotic behavior of a system-no boundary condition need be imposed, since the iterative scheme is superstable asymptotically.
Acknowledgment. This work has been funded by the Natural Sciencesand Engineering Research Council of Canada (NSERC); M.R.R. thanks NSERC for the award of a postgraduate scholarship. We thank Dr. John Bunting for many very helpful discussions. Appendix: Perturbation Series in me for At in the CI Mechanism There are many ways to solve functional equations, e.g., iteration and perturbation expansion techniques.10JO Here, we obtain a perturbation solution to eq 1.6 with me as a small parameter that is the phase-space analog to the outer solution of singular perturbation theory.’ The perturbation solutiongiven justifies the expansion in eqs 3.4 for & for the CI mechanism. Consider the pair of ODES of eqs 1.6 as a coupled set of functional equations, and let us seek a solution of the form
OD
Wesubstitutetheseseriesintoeqs 1.6andrearrange. Thesolution depends on the choice of zeroth-order solution, and for A we choose
xo = 0
(A.2) (A.3)
Using this, the coefficients of me are
This first-order solution justifies the expansion in eqs 3 for &. For j 1 2 and the initial condition eqs A.3 and A.2, the perturbation series specializes to
r
i- I
I
k=l
The terms for j = 2 and higher can be directly evaluated from this expression. For j = 2, the coefficientsare O(mp),with mi/KI fixed, so the correction to the standard approximation is not necessarily small even if mi (and thus KI for mi/& fixed) is laree--ma must also be small! The functional equations have other perturbative solutions of the form A.l, e.g., “inner solutions”, but we have forced an expansion for &,Le., the outer solutionof thesingular perturbation problem,’ by choosing the initial conditions A.3 and A.2. A few trials with the extended perturbation series show that it converges to the equation for & even for moderate me.
References and Notes (1) (a) Heineken, F. G.; Tsuchiya, H. M.; Aris, R. Math. Biosci. 1967, 1,95-113. (b) Rubinow, S.I.; Lebowitz, J. L. J. Am. Chem. Soc. 1970,92, 3888-3893. (c) Lin, C. C.; Segel,L. A. MathematicsAppliedtoDeterministic Problems in the Natural Sciences; Macmillan: New York, 1974; Chapters 9,lO. (d) Murray, J. D.Lectures on Nonlinear-Differential-EquationModels in Biolonv; Clarendon: Oxford. 1977: ChaDter 1. (2) &e, e.g.: Ainsworth, S. Steady-Sthe Enzyme Kinetics; University Park Baltimore. 1977. (3) Briggs, 6. E.;Haldane, J. B. S. Biochem. J. 1925, 19, 338-339. (4) Phase space(s) will be denoted generically by r and the coordinates of a particular phase space will he indicated as the corresponding Cartesian product. Each r has the usual Euclidean metric. See, e.g.: Tikhonov,A. N.; Vasil’eva, A. B.; Svwhnikov, A. G. Differential Equations; Springer-Verlag: Berlin, 1985. (5) In cases where the total concentrations and other conservedquantities are difficult to determine by inspection, an algebraic method can be used.See, for example: (a) Holmes, M. H.; Bell, J. J. Comput. Chem. 1991,12, 12231231. (b) Schuster, S.; HBfer, T. J. Chem. Soc., Faraday Trans. 1991,87, 2561-2566. (6) We use this notation because eo, io, etc. will denote the iterates of functional equations. (7) Equation 1.4 does not appear in the reduced system 1.5 because P is formed irreversibly. (8) In the general case, the entire hierarchy of these equations can be derived by methods formallyequivalentto thoseusedin center manifold theory. See: Carr, J. Applications of Centre Manifold Theory;Appl. Math. Sciencur Vol. 35; Springer: New York, 1981. (9) Roussel, M. R.; Fraser, S. J. J. Phys. Chem. 1991,95,8762-8770. (10) Roussel, M. R. MSc. Thesis, University of Toronto, 1990. (11) Mees, A. I.; Rapp, P. E. J . Math. Biol. 1978, 5, 99-114. (12) Curtiss, C. F.; Hirschfelder, J. F. Proc. Natl. Acad. Sci. U.S.A. 1952, 38, 235-243. (13) Maas, U. A.; Pope, S. B. Combusr. Flame 1992,88, 239-264. (14) (a) Palsson, B. 0.;Lightfoot, E.N. J. Theor. Biol. 1984,111,273302. (b) Palsson, B. 0.;Jamier, R.; Lightfoot, E. N. J . Theor. Biol. 1984, 111, 303-321. (c) Palsson, B. 0. Chem. Eng. Sci. 1987, 42, 447458. (15) Fraser, S.J. J . Chem. Phys. 1988,88,47324738. (16) Nguyen, A. H.; Fraser, S.J. J. Chem. Phys. 1989, 91, 186193.
Global Analysis of Enzyme Inhibition Kinetics (17) Roussel, M. R.; Fraser, S . J. J. Chem. Phys. 1991, 94, 71067113. (18) Wiggins, S.Introduction to Applied Nonlinear Dynamical Sysrems and Chaos; Springer-Verlag: New York, 1990. (19) King, E. L.; Altman, C. J. Phys. Chem. 1956, 60, 1375-1378. (20) Roussel, M. R.; Fraser, S . J. J. Chem. Phys. 1990,93, 1072-1081. (21) Hubbard, J. H.; West, B. H. Differential Equations: A Dynamical Systems Approach, Part I; Springer-Verlag: New York, 1990; Chapter 1. (22) Throughout this work, we have used values of the rate constants and total species concentrations consistent with those observed in vivo. Typical ranges for the rate constants were obtained from Fersht: Fersht, A. Enzyme Structure and Mechanism, 2nd ed.; Freeman: New York, 1985; Chapter 4, Part 2, section B. Physiologically realistic total concentration values were obtained from ref 36a. (23) Strictly the euuilibrium point P is the only attractor of the . speaking, _ system 1.5. (24) Pachner, J. Handbook of Numerical Analysis Applications; McGraw-Hill: New York, 1984; Chapter 6. (25) Right eigenvectors of J are column vector; v is compared with the right eigenvectors of J in this representation. (26) The phase-space flow corresponds to a vector field V(X) C r parametrized by mi. Then as m, moves continuously in the space 'l X m,, bifurcationsat particular mlvaluescorrspond to mathematicallysharpchanges in v(x;m,). (27) Ouellet, L.; Laidler, K. J. Can. J. Chem. 1956, 34, 146150. (28) Darvey, I. G. J. Theor. Eiol. 1968,19, 215-231. (29) Ryde-Pettersson, U. Eur. J. Eiochem. 1989, 186, 145-148. (30) Hearon, J. Z . Ann. N.Y. Acad. Sci. 1963, 108, 36-68. (31) Ryde-Pettersson, U. Eur. J. Eiochem. 1990, 194, 431-436. (32) The velocity comppnents must be scaled for plotting because the velocity norm Ilv(r)ll=2!( + P + t 2 ) 1 / 2 decreases exponentially fast with time. Thevelocityvectorv(t) ISfirst projectedonto thebasisofright unit eigenvectors (u,(t)k = 1,2,3)giving v = Zaph Then the scaled components a,(r)/llv(r)ll are plotted in the figure. (33) The slow manifold JU is very close to the line in r o n which the phase flow moves most slowly, Le., the one-dimensional set of smallest magnitude vectors. However, this 'slow velocity line" is not quite tangent to the flow and so is not a trajectory. (34) In the irreversible, single-substrate enzyme systems studied here, substrate concentration s is unbounded, so the equation for JU in terms of s has no vertical asymptotes anywhere in the physical regime (see ref 10). However, all other concentrations are bounded, so an attempt to represent JU as a function of, e.g., inhibitor concentration, fails. However, in reversible mechanisms all concentrations are bounded, Le., the physical domain in r is bounded. (35) Henderson, P. J. F. Eiochem. J . 1973, 135, 101-107. (36) (a) Srere, P. A. Science 1967,158,936937. (b) Sols,A.; Marco, R. Curr. Topics Cell. Reg. 1970,2,227-273. (c) Laurent, M.; Kellershohn, N. J . Mol. Bi01. 1984, 174, 543-555. (37) (a) Cornish-Bowden, A. Fundamentals of Enzyme Kinetics; Butterworths: London, 1979, section 5.9. (b) Segel, I. H. Enzyme Kinetics; Wiley: New York, 1975; pp 46, 106. (38) See, e.&: Wilkinson, J. H. The Algebraic Eigenualue Problem; Clarendon: Oxford, 1965. (39) Note that the vector function ( i - i*(s), e e*(s)) is annihilated by the substantial derivative operator d/dr = k a/& i a/ai + e a p e restricted to JU and so behaves formally like a first integral of the motion on A. (40) Since JU is the penultimate 'attractor" for the flow, its position can be determined numerically by using pairs of trajectories.I0 The same method can be used to locate A in (41) Attempts to construct JU in a region where it does not exist will diverge. We have also studied perturbation expansions in m, (discussedin the
+
The Journal of Physical Chemistry, Vol. 97,No. 31, I993 8327 Appendix) and in kp (similar to the expansion discussed in refs 20 and 10). Theadvantageofiterationis that itsconvergenceis not limited by thesmallness of some parameter, as discuued in ref 10,so (provided the functional equations are nonquadratic, see section 4) divergence will usually be indicative of dynamical anomalies in the system. (42) (a) Aitken, A. C. Proc. R . Soc. Edinburgh 1926,46,289-305. (b) Steffensen, J. F. Skand. Aktuarietidskrift 1933, 16, 64-72. (43) The asymptotic expansion of eq 4 to O(r2) is also formally q m ; ) . Therefore, q 4 also agrees with the asymptotic expansion to O(s2) of the perturbation series for c(s) (see Appendix) up to the term x2(s)mt. (44) (a) Alberty, R. A. J. Am. Chem. Soc. 1953, 75, 1928-1932. (b) I . 1957,66,34P. (c) Alberty, R. A. J.Am. Chem. Soc. Dalziel, K. Biochem. . 1958, 80, 1777-1782. These papers give reciprocal steady-state velocity expressions that are asymptotic expansions to second order in two reciprocal concentrations. Such expressions may easily be rewritten as asymptotic expansions for the steady-state reaction velocity u. However, these are still two-variable expansions and so represent surfaces, not lines like the slow manifold. (45) In this case, there is no difficulty in using thes X i X c representation. In either representation, the functional equation for the relaxation surface is well behaved. (46) All the conserved total concentrations are used in obtaining the equation for the two-dimensional Z and the related two-dimensionalnullcline surfaces. (47) The SSA is the line-like intersections of e, c, i, and h nullclines for any choice of independent dynamical variables in eqs 1.1 because, in the reduced three-dimensionalsystem, any two of these nullclinesurfaccs intersect in a curve and, from eqs 1.1, this implies the intersection of the remaining two nullclines with the same curve. (48) This means that the surface function e - e*(s,i) is annihilated by the substantial derivative operator39 restricted to Z and so behaves like a first integral of the motion on 2. (49) Hiromi, K. Kinetics of Fast Enzyme Reactions; Wiley: New York, 1979. (50) The "funnels" introduced in ref 21 are related to the notion of the nullcline prism as a trapping region in ref 20. (51) In this mechanism, because there is no distinct inhibitor, the classical KA approximation is identical to the geometrically correct SSA; Le., the distinction discussed in section 3 does not exist. (52) When af/a@y#,ag/**, the evolution of eq 6.5 is found in the usual way by projecting the inibal state {6@0,670) onto the eigenvectorsof the matrix in eq 6.5. However, if the diagonal elements of this matrix are equal, their value is the degenerate eigenvalueof the matrix. The eigensystemb defective in this case with only one eigenvector, so the time n appears linearly in the evolution of the system. (53) Poland, D.J. Comput. Chem. 1990, 11, 382-395. (54) This problem reduces to presenting a set of equations for the structure ofinterest. We havedemonstratedthetechnique for the twoobvious"extreme" reductions from n dimensions to 1 dimension and from n to n - 1 dimensions. Intermediate reductions are not unique, even after using all conserved total concentrations and separating the remaining variables into dependent and independent sets. (55) Thomas, R.; Richelle, J.; d'Ari, R. Bull. Classe Sci. Acad. R . Eelg. 1987, 73,62-83. (56) Saaty, T. L. Modern Nonlinear Equations; Dover: New York, 1981; Chapter 1. Hale, J. Ordinary Differenrial Equations; Wiley: New York, 1969. (57) Constantin, P.; Foias, C.; Nicolaenko, E.; Temam, R. Inzepal Manifolds and Inertial Manifolds for Dissipative Parrial Differenrial Equations; Springer-Verlag: New York, 1989.
