, I Phys, , Chem. 1991, 95, 8762-8770
8762
size, which limited the ultimate time duration for the diffusion process, In this sense, the present study is a complement of the earlier related work of Flanders and Fricke.I5 In their work, the total reaction fractions between the first-order and second-order reactions were obtained as the reactions are completed. They also showed that Jaffe's method was inadequate in the cases of high-recombination situations. The present study shows that the reactive species diffuse more slowly in the approximate method than in a real system. Jaffe proposed the number density function to be N(t)*(r,t), in which N(t) is a diffusion-independent term. Apparently, a global average account of the recombination reaction underestimates its contribution in the high-density region but overestimates its contribution in the low-density region. Nevertheless, the overall approximate result still overestimates the recombination contribution. In the approximate model, the reactive species tend to be more localized than the true distribution and the total number is always less than the real values. In most pulsed-laser devices the fundamental spatial mode (TEMW)was usually mixed with some minor higher order modes. Instead of going into a full analysis by Hermitffiaussian functions, a more practical treatment is approximately fitting the real intensity profile with the fundamental Gaussian intensity dis-
tribution. If, say, 10% of the intensity c o m a from the TEMlo mode, the initial local density may have more than a 10% deviation locally from the average Gaussian distribution. However, according to numerical calculations, the error in total number would be less than 1.5% under the recombinationdominant regime. The higher order mode may contribute more errors. Nevertheless, the contribution of the higher order mode to the intensity is usually less important. In the diffusive motion regime, the error would even be much less than the above figure. Up to now we have been focusing our attention on the single-photon dissociation process in generating the reactive radicals. Two-photon or multiphoton dissociation of molecules has been under intensive investigation ever since the introduction of high-power lasers. The present calculations could be easily extended to multiphoton dissociation cases. Since multiphoton dissociation tends to localize the reactive radicals more than the single-photon process under the same laser spatial intensity, the approximate solutions would show larger deviations from the exact ones. The general patterns of the reaction-diffusion behavior are expected to be the same.
Acknowledgment. This work is financially supported by the National Science Council, Republic of China.
Accurab Stmdy-State Aqwoxhalions: Implkatlons tor Klnetka Experknents and Mechanism Marc R. Rorrssel and Simon J. Fraser' Department of Chemistry and Scarborough College, University of Toronto, Toronto, Ontario MSS 1Al. Canada (Received: May 13, 1991)
A theoretical analysis of accurate steady-state experiments is presented. In principle, mechanisms can be distinguished and all rate constants found independently by using this approach. The method is illustrated by comparing two enzymscatalysis mechanisms originally proposed by Henri: the usual enzyme-complex model, and a bimolecular model. It is well-known that both mechanisms give the hyperbolic rate law in the standard kinetic analysis; therefore, it is believed that steady-state kinetics is unable to distinguish between them. However, we show that the hyperbolic rate law is only the lowest order approximation in a systematic procedure for getting more accurate formulas for the reaction velocity. In higher order approximations the two model mechanisms are (mathematically) distinct: At any substrate concentration, the enzymecomplex mechanism has a larger velocity than the hyperbolic law, and the bimolecular mechanism a smaller velocity. This diffenmx is most obvious in Woolf or Eadie-Hofstee plots where the hyperbolic law is linear, but the two mechanisms give curves with distinctive characteristics. In general, these systematic methods provide very accurate steady-state velocity expressions, and they can be applied to accurate experimental results to yield more kinetic information.
native bimolecular mechanism is (B), where the (nonproductive
1. Introduction
In Henri's pioneering research on enzyme catalysis he remarked that either a productive or a nonproductive enzyme-complex mechanism could explain the observed kinetics, because in the mathematical approximation he used, both mechanisms gave the same reaction velocity expression.' In this paper we will show that such cases of "aincidmt" steady-state kinetics can be resolved by accurate mathematical analysis; this improved description should greatly increase the power and usefulness of steady-state kinetics studies. The first model Henri proposed is the standard Henri-Michaelis-Menten2 mechanism (H), where E is enzyme, S substrate, k
k2
E+S+H-E+P k-i
H = ES enzyme-substrate complex, and P product. The alter-
enzyme-subtrate) complex ES is now denoted by B. Here, in the first step, we have made complex formation the forward process, and complex decomposition the reverse process, so that (H) and (B) correspond as closely as possible. While the bimolecular mechanism (B) is not nowadays seriously considered as a probable enzyme reaction pathway, mechanisms (H) and (B) are alternatives in other important elementary processes such as chargetransfer-complex reactions: where the detected complex may not be productive. Since both mechanisms give the hyperbolic rate law in the standard steady-state kinetics analysis, it is commonly thought that steady-state kinetics cannot distinguish between them. However, in this paper we show that these cases can be differentiated by an accurate theoretical steady-state
( I ) (a) Henri, V. Compr. Rend. Acad. Sci. (Paris) 1902, 135, 916-919.
(b) Henri, V. Lds GMrales de rAction des Diastawq Hennann: Paris, 1903. (2) Michaelis, L.: Menten, M. L. Biuchem. Z. 1913,49, 333-369.
0022-3654/91/209S-8762$02.50/0
~
~~
(3) See e.&: Fukuzumi, S.; Koshi, S. J . Am. Chem. Suc. 1982, 104,
7599-7609.
0 1991 American Chemical Society
Accurate Steady-State Approximations
treatment, and in principle this analysis allows all the specific rate constants to be determined by accurate experiment. In fact, the following two-step mechanistic models have all been studied by steady-state methods: simple enzyme ~atalysis,'.~*~ exciplex formation,5 nonproductive complex formation in bimolecular processes,'*3and the Lindemann mechanism.6 The equations of motion for these models cannot be solved explicitly, but, characteristically, they have a crossouer regime where the order of reaction changes. Our analysis provides an exact description of this regime, implying the complete determination of the rate constants in these cases. This procedure can be generalized to multistep reactions. Because of the practical importance of steady-state kinetics,' the dynamical steady state, Le., the precise slow time evolution of the system, has been treated by a number of mathematical techniques: perturbation theory? singular perturbation theory? series expansion,1° topological reaction network and functional equation methods.I3 Each of these treatments provides a generalized description of the constraints on differential rates or gives accurate formulas for the evolving steady state of the reaction. There have been several extensive kinetics studies in which many or all of the rate constants have been extracted from transient or steady-state experiments. From a reading of the sample list of papers in ref 14, it can be appreciated that the work required is both careful and elaborate. The theory presented here could provide similar detailed information on rate constants from far fewer, accurate steady-state measurements. To understand any steady-state method fully, we must describe the geometrical structures (line manifolds, nullclines, etc.) that underlie the evolution of chemical rate equations in the phase space r of species concentration ~ a r i a b l e s . ' ~ JSuch ~ a picture gives (4) (a) Wong, J. T.-F. Kinetics of Enzyme Mechanisms; Academic: New York, 1975. (b) Segel, 1. H. Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady-State Enzyme Systems; Wiley: New York, 1975. (c) Laidler, K. J.; Bunting, P. S. The Chemical Kinetics of Enzyme Action; Clarendon: Oxford, 1973. (d) For a detailed discussion and (mode) analysis of mechanism (H) see: Palsson, B. 0. Chem. Eng. Sci. 1987, 42, 447-458. ( 5 ) Gordon, M., Ware, W. R., Ed. The Exciplex; Academic: New York, 1975. (6) (a) Lindemann, F. A. Trans. Faraday SOC.1922, 17, 598-599. (b) Christiansen, J. A. Ph.D. Thesis, University of Copenhagen, 1921. (c) Hin-
schelwood, C. N . The Kinetics of Chemical Change; Clarendon: Oxford, 1955. (7)
See e.g.: (a) Benson, S . W. The Foundations of Chemical Kinetics; McGraw-Hill: New York, 1960. (b) Laidler, K. J. Chemical Kinetics; Harper and Row: New York, 1987. (8) Hirschfelder, J. 0.J . Chem. Phys. 1957,26, 271-273. This is an early application of time-dependent perturbation theory to steady-state chemical kinetics. (9) (a) Heineken, F. G.; Tsuchiya, H. M.; Ark, R. Marh. Biosci. 1967, I , 95-1 13. As an example of the singular perturbation method, this paper gives a very detailed analysis of mechanism (H). (b) Bowen, J. R.; Acrivos, A.; Oppenheim, A. K. Chem. Eng. Sci. 1963, 18, 177-187. (c) Richardson, W.; Volk, L.; Lau, K. H.; Lin, S.H.; Eyring, H. Proc. Natl. Acad. Sci. U S A . 1973, 70, 1588-1592.
(IO) Poland, D. J . Comp. Chem. 1990, 1 1 , 382-395. (11) King, E. L.; Altman, C. J . Phys. Chem. 1956,60, 1375-1378. (12) (a) Clarke, B. L. Ado. Chem. Phys. 1980,43, 1-215. (b) Stucki, J. W. frog. Biophys. Mol. Biol. 1978,33,99-187. (c) Helfferich, F. G. J. Phys. Lightfoot, E. N. J . Theor. Chem. 1989,93,6676-6681. (d) Palsson, B. 0.; Biol. 1984, I l l , 273-302. (e) Chou, K.-C. J. Biol. Chem. 1989, 264, 12074-1 2079. (13) (a) Fraser, S . J. J . Chem. Phys. 1988,88.4732-4738. (b) Nguyen, A. H.; Fraser, S.J. J . Chem. Phys. 1989, 91, 186-193. (c) Roussel, M. R.; Fraser, S.J. J. Chem. Phys. 1990,93, 1072-1081. (d) Roussel, M. R.; Fraser, S.J. J . Chem. Phys. 1991, 94, 7106-71 13. (14) See e.&: (a) Chance, B. Arch. Biochem. 1949, 22, 224-252. (b) Bailey, J. M.; French, D. J . Biol. Chem. 1957, 226, 1-181. (c) Alberty. R. A,; Peirce, W. H. J . Am. Chem. Soc. 1957, 79, 15261530. (d) Theorell, H.; McKinley-McKee, J. S . Acta Chem. Scand. 1961, 15, 1797-1810. (e) Hammes, G.G.; Kochavi, D. J . Am. Chem. SOC.1%2,84, 2073-2079. (0 Eigen, M.; Hammes, G. G.Ado. Enzymol. 1963,25, 1-37. (g) Bloomfield, V.; Peller, L.; Alberty, R. A. J . Am. Chem. SOC.1972,84, 4375-4381. (h) Chandra Sekhar, V.; Plapp, 8. V. Biochemistry 1990, 29,4289-4295. ( I 5 ) A description of geometrical flow structure in the phase space is given in e.g.: (a) Arnold, V. 1. Ordinary Dijjerential Equations; MIT Cambridge, 1978. (b) Kaplan, W. Ordinary Dvferential Equations; Addison-Wesley: London, 1958.
The Journal of Physical Chemistry. Vol. 95, No. 22, 1991 8763 an unequivocal meaning to exact or approximate steady-state solution^.'^ In F, almost all initial states during transient decay are "attracted" by a unique trajectory, the linelike slow manifold JU, which then guides them toward the eventual equilibrium of the system.'* All steady-state methods try to locate At, and we remark that knowing the position of JU accurately allows all rate constants to be extracted independentlyfrom steady-state data. Motion along JU describes reaction kinetics not only in the usual "steady-state" regime but also closer to equilibrium, where there is often a change in the reaction order. This paper describes a functional equation method for locating A;"this procedure is related to the backward differentiation method for integrating stiff differential equationsI9 first described by Curtiss and Hirschfelder.20 Local (linearized) versions of this procedure are used in Gear's method for integrating stiff differential equations.I" However a kineticist frequently requires an explicit, global formula for the (steady-state) reaction velocity, and this paper illustrates a general procedure for deriving such formulas. From now on we will refer generically to the entire slow time evolution of any reaction (which takes place very close to JU) as the steady state, and we now embark on a description of the phase-space geometry of the steady state. For d chemical species, the instantaneous state of reaction can be represented as a point, x ( t ) . in the positive part (xj 1 Ob = 1, 2, ..., 4 of phase space r = I'd.z1 The reaction as a whole is a flow, in I', parametrized by a set of rate constants (k). Since d first-order differential equations govern this flow, Le. the velocity vector v is uniquely determined by x, thereby fixing the entire evolutionary history on the trajectory passing through it. Since there are linear conservation laws for eq 1.1 resulting from stoichiometry, the flow is restricted to a Euclidean subset of r of dimension less than d e.g., for a single-step reaction this subset is the positive real line, and for a two-step reaction it is the positive quadrant of the phase plane. From now on eq 1.1 will be written by using these linear conservation laws to reduce dimension, so that d will now denote this reduced dimension of r = rd, the corresponding "smaller" phase space. In chemical reaction models, eq 1.1 is often a stiff system, Le., the various decay modes have well-separated relaxation times. The transient decay in such systems has a simple but important geometry in Each relationfi(x; (k))= constant derived from eq 1.1 defines a (d - 1)dimensional surface. The nullcline surfaces f = 0, have a special role in stiff systems because their intersections Iorm a (nested) geometrical hierarchy of sets through which the (trajectory) flow cascades. Eventually all trajectories come to' lie close to the set of line intersections of these surfaces, which form the edges of a curved, prism-shaped region in r. The edges of the prism pinch together at the point equilibrium of the flow I':13b9d
(16) A description of manifold theory and differential topology is given in: (a) Carr, J. Applications of Centre Manifold Theory, Appl. Marh. Sciences 35; Springer: New York, 1981. (b) Marsden, J. E.; McCracken, M. The Hop/ Bifurcation and Its Applications, Appl. Marh. Sciences 19; Springer: New York, 1976. (d) Chillingworth, D. R. J. Differential Topology with a View to Applications; Pitman: San Francisco, 1976. (17) The idea that an understanding of the steady state can be achieved through the elimination of time as a variable appears in the discussion section of Wong, J. T.-F. J . Am. Chem. Soc. 1%5,87, 1788-1793. (18) The flows considered in this paper, to which the SSA is commonly applied, are very simple, having just one stable fixed point. Mechanism (H) and (B) belong to this class. Note that these irreversible mechanisms have an unstable fixed point at infinity, but if product formation is reversible, then this fixed point moves into the (unphysical) finite part of the phase plane r2, discussed in the Introduction. We do not consider the dynamically richer flows with limit cycles or strange attractors. (19) (a) Gear, C. W. Numerical Initial Value Problems in Ordinary Difjerenrial Equations; Prentice-Hall: Englewood Cliffs, NJ, 197 1; see especially Chapter 11. (b) Lambert, J. D. Computational Methods in Ordinary Dflerential Equations; Wiley: New York, 1973; see the discussion in Chapter 8. (20) Curtiss, C. F.; Hirschfelder, J. 0. Proc. Narl. Acad. Sci. U.S.A. 1952, 38, 235-243. (21) Phase space r d is given the usual Euclidean metric, is., the {x,) form
an orthogonal d-dimensional coordinate system and distance d p , y ) between points x and y of r d is defined by d(x.y) = [X,".,(x, - yj)']'/ .
8764 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 where v = 0.22 The steady-state approximation (SSA)11*23 corresponds to one of the prism-edge lines and the equilibrium approximation (EA)'J to another. Strictly speaking, these approximations imply constrained motion along the appropriate line." In reality however, it is not S A , EA, or any other line intersection in 'I that attracts every trajectory, but the slow manifold lying clm by (and, we argue,u entirely within the prism). In stiff systems, the true motion along M is imitated by the constrained motion along the nearby SSA and EA lines. Evidently, the steady-state kinetics of the system 1.1 is exactly solved by locating A. While SSA and EA are not uniformly accurate at all concentrations, particularly in the later stages of a reaction, the fact that M is a line means that all concentrations can be written explicitly in terms of a single (variable) concentration. in normal steady-state experiments, the measured (scalar) reaction velocity, u say, is just p, the rate-of-product-formation component of the velocity vector v4 along A. Determining u requires only one series of steady-state concentration mcasurantnts but gives limited kinetic information from the usual SSA or EA formulas. However, by use of better dynamical-steady-state formulas, suitable data from a single accurate experiment could be used to discriminate between alternative mechanisms and also determine all rate constants for the correct mechanism. In section 2 we set up a rigorous mathematical framework for systematically solving the steady-state kinetics of the two models. For planar systems such as (H) and (B), iteration or perturbation method^'^*^^ converge very rapidly to the exact solution. Throughout we deliberately set up the mathematical formulas in terms of experimentally accessible parameters and express the small-[SI and larga[S] (asymptotic) steady-state behavior using these parameters. In section 3 we show how mechanisms (B) and (H) may be distinguished by combining theory with accurate steady-state experiments. Of the various graphical methodsz6 for linearizing the hyperbolic law, the Woolf or Eadie-Hofstee plots2' seem to be best because they emphasize the important crossover regime.28 Calculating reaction velocity u as a function of the reciprocal time (frequency) w = u/ [SI,displays the crucial fact that globally, in the accurate steady-state treatment, (B) always lies below the hyperbolic law because it has a smaller reaction velocity, and (H) always lies above the hyperbolic law because it has a greater reaction velocity. It is worth commenting that "lag" in the transient kinetics29allows (H)and (B)to be distinguished,Mbut (22) Here we are considering differential flow for closed chrmicar sptcms,
so that there is a unique stabk equilibrium. For a proof see: Shear D. J.
Theor. Biol. 1967, 16,212-228. (23) (a) Chapman, D. L.; Underhill, L. K. J. Chem. Soc., (Transacrions) 1913,103,496-508. (b) Bodmstein, M.2.Phys. Chem. 1913,85,329-397. (c) Briggs, G. E.; Haldane, J. B. S.Biochem. J. 1925, 19, 338-339. (d) Walter, C. Steady-Srare Applicarions in Enzyme Kinetics; Ronald: New York, 1975. (e) Fromm. H. J. Inltlol Rate Enzyme Kinerics; Springer: New York, 1975. (0 Gutfreund. H. Enzymes: Physical Principles; Wiley: New York, 1972. (24) On A, all the velocity components x, arc small but nomanishing since they are stiffly coupled. Therefore & cannot intersect a (d 1)-dimensional nullcline surface in r,,; indeed, A would have to be o r t h o p a l to some x, wis at such a crossing, and we remark that .4t is the only trajectory that lies everywhere close to all the nullcline line intersections. Since we always find that & lies within the prism both near equilibrium (v = 0) and asymptotically far from equilibrium (x -), and since the nullclines arc surfaces in ,'I At must therefore lie globally within the prism. (25) Applications of perturbation theory to dynamical problems are given in e.g.: (a) Giacaglia. G. E. 0.Perrurbation Methods in Non-linear Systems, Appl. Math. Sci. Ser. 8; Springer: New York, 1972. (b) Nayfeh, A. H. Perturbation Methods; Wile : New York, 1973. (c) Nayfeh, A. H. Inrro. ducrion ro Pcrrur6ation T e c h p e s ; Wiley: New York, 1981 (26) For a comparison of t h a c methods, see: Cornish-Bowden, A. Principles of Enzyme Kinetics; Butterworths: London, 1976. (27) (a) Eadie, G. S.J. Biol. Chem. 1942, 146, 85-93. (b) Hofstee, B. H. J. Science 1952,116,329-331. (c) Hofstee, B. H. J. Narun 1959, 184, 12961298. This method was wed by Wodf (unpubkhcd), but see: Haldane, J. B. S.;Sta, K. 0.Allgetnek Chltnie der Enzyme 11%Vetlag von Tkdor Steinkopff: D d e n , 1932. Haldane, J. B. S . Nature 1957, 179, 832-832. (28) The "advantage" of the Woolf or EaditHofstae repmentations is that the physically relevant parts of the velocity manifolds uH and vg appears in a closed, finite (compact) part of the plane o X D . Departures from the hyperbolic law must then appear as curvature in Y, which is easy to see and analyze.
-
-
Roussel and Fraser 1.
0.
0.
c 0.
0.
p\
0 0 0 0
'\\\
\ , v
0 2
\\
\I
0 6
0 4
0 8
1 0
S
Figure 1. Trajectory flow in the s X h plane for mechanism (H). Rate constants are kl = 8.00, k-, = 1.00, k2 = 1.00, with eo = 1.00. For these
[k},the SSA nullcline h&) for mechanism (I$ is identical with the EA nullcline b&) for mechanism (B) with the lk}values given in Figure 2.
The flow crosses the S = 0 nullcline (EA) vertically, and the k = 0 nullcline (SSA) horizontally, attracted by slow manifold AH,lying between EA and SSA. At large s, AHis closer to SSA. The third iterate of cq 2.5, starting from hs(s), coincides with AHon the scale of this figure. Note how close AHis to & in Figure 3, and note the differencar in the slope of the fast transients.
the values of the rate constants that make this lag easy to detect also allow these mechanisms to be resolved by accurate steady-state
analysis. For both (H) and (B) we give accurate representations of u as a function of o and discuss how all rate constants may be found from accurate steady-state experiments, Deviations from the hyperbolic law in enzyme mechanisms are well doc~mented,3~3~ but if the reaction velocity increases m~notonically~~ as a function of substrate concentration in such cases,our analysis suggests that mechanism (H) might still apply and could then be used to extract all the rate constants. In section 4 we discuss the implications of this analysis and its applications to other kinetics problems.
2. Geometrical Picture of Two Models of Enzyme Catalysis In this section we derive the equations for the slow time evolution of (H) and (B) that locate the corresponding slow manifolds and obtain the standard approximations (SSA and EA) as limiting cases. It is important to discuss the velocity and manifold formulas at large substrate concentrations, since they determine the asymptotic forms for the Michaelis constants for the two mechanisms. This allows the kinetics to be compared correctly when [SI is very large and the enzyme nearly saturated. Throughout this paper we use lower case italic letters for concentrations, so that for species X,the concentration [XI is written x. Mechanism H. Using the conservation of enzyme, e + h = eo, (H) can be reduced exactly to the following system of ordinary differential equations: = kl(eo- h)s - (k-' + kz)h (2.1) j. = -kl(eo - h)s klh (2.2)
+
p = kzh (2.3) Only the first two equations are coupled nontrivially because in (29) See e.&, articles in: (a) Kustin, K., Ed. Merhods Enzymol: Fasr Reactions 1969, 16. (b) Bernasconi, C . F. Relaxation Kinetics; Academic: New York, 1976. (30) See: Hiromi, K.;Kinetics of Fasr Enzyme Reacrions: Theory and Pracrlce: Wiley: New York, 1979; pp 41-42, for discussion of differentiation between mcchanhms (H) and (B) through the obeervation of lag. (31) Hill, C.M.;Waight, R. D.; Bardsley, W. 0. Mol. Cell. Biochem. 1977, 15, 173-178. (32) Hanes, C. S.;Bronskill, P. M.;Gurr, P. A.; Wong, J. T.-F. Can. J . Biochem. 1972,50, 1385-1413. (33) See e.g.: (a) Hyalop, J. M.Real Variable;Interscience: New York, 1960, pp 60-62. (b) Apostol, T. Mathematical Analysis: A Modern A p proach ro Advanced Calculus; Addison-Wesley: London, 1963; p 78.
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8765
Accurate Steady-State Approximations model (H) product formation is irreversible. Total substrate, s + h p = so, is conserved but does not appear in the reduced equations of motion. This means that the evolution of the system can be represented in the phase plane r2(H) = s X h. The experimentally measured variable p appears in the substrate conservation law, so knowing p and s determines h. Figure 1 shows trajectories for (H)in s X h at kl = 8.00, k-l = 1.00, k2 = 1-00, and eo = 1.OO. It is evident that all trajectories reach a unique trajectory, the slow manifold A H . The slow time evolution of the is located. system is precisely defined once We now follow the method of center manifold theoryI6 to obtain a single equation for the phase-space trajectories. .First write h = h(s) and differentiate with respect to time to get h = (dh/ds)S; introducing eqs 2.1 and 2.2 into this relation then gives the trajectory or one-dimensional invariant manifold differential equation
+
dh A’= ds
-6 = ki(eo - h)s - (k-1 + k2)h -kl(eo - h)s + k-lh
S
(2.4)
The coupled eqs 2.1 and 2.2 define theflow in s X h, and the single eq 2.4 defines the phase-plane trajectories of this flow. The equilibrium (stagnation point h = S = 0) of this system is at h = s = 0. Equation 2.4, can be rearranged to give the functional equation
h(s) =
eos
s
+ K E + c/(I + h’(s)) (2.6) (2.7)
is another rate constant ratio, which (for slow decay at least) is a small parameter. Equation 2.5 defines a relation between h on the left and its derivative h’on the right; it can be solved by direct iteraction from a smooth (infinitely differentiable) starting function, or by perturbation theory. In fact, for such functional equations, the network analysis version of the SSA”J2 or EA” provides good starting functions for more accurate solutions. Also note that eq 2.5 contains only two independent rate parameters, KE and t, and defines a purely geometrical structure-a curved line (invariant manifold) in s X h. However, (H) contains three independent rate constants. The third (required) relation is implied by writing the reaction velocity in terms of the solution of eq 2.5, so restoring a time scale into the problem. It is easy to extract equations for EA and SSA from eq 2.5 and interpret their meaning. Suppose t is fixed in eq 2.5; then lim h(s) = /I&)
h*--
e,,s/(s
+ KE)
(2.8)
gives the EA curve, so h&) C s X h, the vertical flow nullcline, is the image of any vertical line under the functional eq 2.5. Since S = 0, or equivalently h’ = 0 3 , everywhere on hE, all trajectories The EA reaction cross it vertically, as shown in Figure (34) There are some important conceptual differences between eq 2.4 and 2.5, although they are algebraically equivalent. The trajectory equation 2.4 has movable, branch-point singularities on the S = 0 nullcline, where the trajectories turn over producing a local square-root dependence of h on s. (Similar singularities appear in ds/dh on the h = 0 nullcline.) This means that AHfor example cannot be found by extrapolation procedures. In contrast, the functional equation 2.5 operates on the space of functions that are at least once differentiable (C’),and usually infinitely differentiable (C)for s E [0, m). Of all trajectories (solutions of eq 2.4) AHis the only bounded, irerorioely sroble solution lying in the function space implied by eq 2.5. All other solutions of eq 2.4, i.e. trajectories h = h(s), are umrobleflxed points of eq 2.5. Those solutions of eq 2.5 (trajectories) that cross EA (h’ a) diverge at finite s; those solutions (trajectories) that cross SSA (h’= 0), at some value of s, have increasingly negative derivative at larger s. so that lim, h’(s) = -1; therefore these solutions of q 2.5 diverge as s m. For a discussion of movable singularities see e.g.: Davis, H. T. Introduction to Nonlineor Dif/crential and Integral Equotiom; Dover: New York, 1962; pp 8, 185.
-
h ‘-0
gives SSA, where Ks = KE
+
t
= (k-1
+ k2)/kl
(2.1 1)
the SSA-Michaelis constant. All trajectories cross hs C s X h horizontally as shown in Figure 1, because h = 0 (or equivalently h’ = 0) everywhere along this line, and hs(s) is the image of any horizontal line under the functional equation 2.5. The SSA reaction velocity is U S O ) = k2hS.W (2.12) Rigorously, the SSA is constrained motion along hs combined with the differential eqs 2.2 and 2.3, which together imply s + p = To obtain the constraint, set the right-hand side of eq 2.1 to zero and differentiate with respect to time to get kl(eo - h)S - ( k l s &l + k2)h = 0; then, using eq 2.1 again, set h = hs(s) to relate h and S. In Figure 1 the nullcline edges hs and hE pinch together at 6 = S = 0, as described in the Introduction. Finally note that lim h(s) = h&) (2.13)
+
where KE is the EA-Michaelis constant for (H), given by KE = k - l / k l
velocity for the formation of p is UE(s) = k2hE(s) (2.9) (Note that from now on we will understand v to mean p.) Rigorously speaking, the EA is constrained motion along hE,together with the differential eqs 2.1 and 2.3, which together imply h + p = constant.’jB To obtain the constraint, set the right-hand side of eq 2.2 to zero and differentiate with respect to time to get kl(eo - h)S - ( k l s + k-# = 0; then, using eq 2.2 again, set h = hE(s) to relate h and S. Again, if e is fixed in eq 2.5, then (2.10) lim h(s) = hs(s) = e $ / @ + Ks)
-
f-0
so that for vanishing decay rate ( k , = 0), hE(s) is a line of equilibrium points in the s X h plane. Notice that KE and KS have units of concentration, but the unimolecular “fission” ratio of the first-order decay rate constants for complex H, t/KE = k 2 / k - , , and the decay probability of H to product, c / K s = k , / ( k - , k2), are pure numbers and so are invariant under any scale transformation (change of units). In the equations of motion for (H) there arefour parameters: eo, k l , k-l, and k,, so we can choose a concentration scale and a time scale in any way we like; thus the scaled dynamics for (H) has only two parameters. To put all variables in dimensionless (scaled) form,9a.bJ2dit is convenient to measure concentrations in units of either eo or the asymptotically correct Michaelis constant, which will turn out to be Ks for (H). For completeness, however, we will retain eo and the Michaelis constants in all formulas. From the analysis given below, all rate constants can be found within a factor depending on the purity of the enzyme. Equation 2.5 can be written in the form h,+l(s) =fis,h,,’(s)) and solved iteratively. This procedure is convergent: From a suitable, smooth starting function ho(s) we will eventually reach a fixed-point (function) h*(s) = limp- h,(s), which represents 0, eq 2.5 becomes the slow manifold A H in s X h. As t iteratively superstable, and for any choice of (finite) h’(s) on the right, converges in one step to h* = hE, which is the line of equilibrium states for k2 = 0, i.e., t = 0. Equation 2.5 can also be solved by perturbation expansion. Using A&) as starting function, iteration and perturbation solutions are completely e q ~ i v a l e n t . 1In ~ ~s x h, A Hlies between hE and hs because the flow is into the region bounded by these lines. The third iterate of eq 2.5, starting from h&), is indistinguishable from A H(found by numerical integration) on the scale of Figure I . What is of fundamental importance experimentally is that the slow time evolution of the system can be represented by progress along AH,and the exact reaction velocity u = v* is bounded between uE in eq 2.9 and us in eq 2.12, so (2.14) uE(s) > u*(s) = d = kzh*(s) > u~(s)
+
-
8766 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991
Roussel and Fraser
2 shows trajectories i n s X b a t 1,= 4.49, kl = 1.12, f 2 = 4.00, and eo = 1.00. Dividing eq 2.19 by 2.20 gives the trajectory (invariant manifold) differential equation
MB
n
Rearranging eq 2.22 gives the functional equation e@[l + (1 + ?)6'(s)] b(s) = s +k E + [(I + 2)s + kE]b'(s)
(2.23)
where - 0 0
0 2
0 6
0 4
-
0 8
kE 1-1/ f , 2 = i2/iI
1 0
It is easy to see that the nullclines and slow manifold converge at large s, because lim h&) = lim h*(s) = lim hs(s) = eo (2.15) r-
e"
Note that 2 is a pure number, being the ratio of two second-order rate constants, so it is not changed by scale transformations. Also, 2/( 1 + 2), the decay probability to product P for each reactive E-S collision is a pure number and similarly invariant. Letting 6' 0 or 2 0 in eq 2.23 gives EA, with equation
- -
bE(S)
v-
However, it can also be shown, by substitution in the functional equation 2.5,that at very large s u*(s) = Us($)
+ 0(s-3)
(2.17)
= e@/($
+ $1
which is the (b = 0)-nullcline. Letting 6'with equation
(2.26) 03
formally gives
bS(4 = w / ( s + R S )
(2.27)
which is the (S = 0)-nullcline, with
ks = kE/(1 + 2)
e-
and the approximate velocities also converge to the true maximum velocity V,,, because V,,, = lim u&) = lim u*(s) = lim us(s) = kleo (2.16) e-
(2.25)
S
Figure 2. Trajectory flow in the s X, b plane for mechanism (B). Rate constants are kl = 4.49, k-, 1.12, k2 = 4.00, with eo = 1.00. For these If),the E A nullcline b&) for mechanism (B)is identical with the SSA nullcline hs(s) for mechanism (H)at the (k)values given in Figure 1. The flow crosses the b = 0 nullcline (EA) horizontally, and the J = 0 nullcline lying between E A and (SSA) vertically, attracted by slow manifold AB, SSA. At large s, ABis closer to EA. The third iterate of eq 2.23, starting from bE(s),coincides with ABon the scale of this figure. Compare the position of the slow manifold and the direction of the transient decay with Figure I .
Y-
(2.24)
(2.28)
Notice how the labeling S and E has been reversed in (B)compared with (H), because of the switching of position of the enzyme complex in the two mechanisms. The solution of eq 2.23 gives the slow manifold ."ICB for mechanism (B),which attracts all other trajectories as is evident from Figure 2. The third iterate of eq 2.23,starting from b&), is indistinguishable from AB(fgund by numerical integration) on the scale of Figure 2. The manifolds are ordered 6S(s) 1 b*(s) 1 be($)
but u*(s) = uE(s)
+
(2.18)
so that us is asymptotically a far better description of steady-state behavior than 0E.35 This can be understood geometrically from Figure 1, since the flow is horizontal on h* and hs but vertical on hE, as s a. (Note that the 0 notation used throughout this paper is explained in ref 33.) The velocity u* contains exact information about the slow time evolution of mechanism (H). In section 3, we show that by accurately measuring s, p, eo, and so and knowing the function u*, all the rate constants { k )can be extracted by suitable fitting procedures from steady-state data. Finally, we note that using h = h*(s) in eq 2.2 gives S = + (kls+ kl)h*(s),which exactly expresses the slow time evolution of the system on JCH. We comment on this 'center-manifold" differential equation16ain the Discussion. Mechanism (B). The above analysis can be also applied to mechanism (B). Using the conservation of enzyme e + b = eo, the system of differential equations for (B) is b = -Lib + i l ( e 0- b)s (2.19) s = R 1 b- (1,+ R2)(eo- b)s (2.20) (2.21) p = R2(e0- b)s
-
Only eqs 2.19 and 2.20are coupled nontrivially. (Velocities and rate constants in mechanism (B) are labeled with a hat?) Figure (35) From the asymptotic analysis in ref 13b for (H)it is immediately apparent that lu*(s) - us(s)l is at most qr2) at large s, but with more care we can show that coefficient of f2vanishes identically and our result sharpens to O(s-)) decay as s -.
-
(2.29)
In contrast to (H) however, the ordering of the reaction velocities for (B)(the rate of product formation) is the reverse of that of the manifold^:^' For a known amount of enzyme eo, we have DE(s)
i2(e0- ~ E ( s ) ) s = 1&e@/(s + kE) 2p*(s)= kz(eo- b,*(s))s 1 4(s) = k2(eo- bs(s))s = k2kSe@/(s + $) (2.30)
While the manifolds converge in the phase plane r2(B)= s X b for this mechanism, i.e. lim b&) = lim b*(s) = lim bE(s) = eo -
e-
P"
(2.31)
s-.-
the velocities are not asymptotically the same, because lim DE(^) = i2kEeo
(2.32)
e-
but
(36) It might seean uausurl to BIlociptc the s nullclime with the S A for mechanism (By,but the 6 nullcline clearly cormponda to the EA, as we have indicated. For (B), twa6@a could be considered PI intermsdiptes: E and S. It is tempting to assign the SSA 20 the e nullcline, s h the catalyst E is regenerated, but enzyme canrervrtion dictate8 that db/dr = -de/dt 80 that the b and e nuflclw both comrpond to the EA. We must therefore urign the S A to the s nullcline. (37) This property is a trivial c o n a c q u e ~of~the ~ ~ nu 'fdd ord ring and qf the form of e velocity exprwion: If 92 b*, then -&h 5 and k$(eO - &) Ip$(eo - b*). The rest of the inequality 2.30 follows from eqs 2.21 and 2.29 in the same way.
-E&*
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8767
Accurate Steady-State Approximations
r
It is therefore important to find which approximation is asymp totically correct. In fact, it can shown, again by substitution in the functional eq 2.23, that at very large s
b*(s) = bs(s) + O(s-I)
(2.34)
+
but
b*(S) = 6e(S)
+ O(S-')
(2.35)
so that, asymptotically, EA is more accurate t h y SSA for mechanism (B), and the true maximum velocity, V,,,,,, is given by V,,,,,, = lim DE(s) i2&0 (2.36)
Like (H), this agrees with the geometry of the flow a t large s, which has the same horizontal direction on 6* and & but is vertical on bS,as shown in Figure 2. This fact is important in determining the relationships between the rate constants for the two mechanisms.
3. Dynamical Comparison of Mechanisms H and B We are now ready to compare the velocity (vector) fields of (H)and (B) so that the rate constants can be related to experimental data. Since Michaelis constants and reaction velocities at very large or very small substrate concentration are experimentally relatively easy to extract, we have used these quantities to parametrize our equations. Parametrization of the Hyperbolic Rate Law. We can fix two relations between {k}and (&I if we can pin down the steady-state behavior at very large values of s, where the enzyme is saturated. The asymptotic analysis of the previous section shows that for (H), SSA is more accurate than EA, but for (B), EA is more accurate than SSA. Therefore comparing reaction velocities as s m, we get
-
k2 =
(3.1)
which, in terms of rate constants is
k2 = i2i-,/il (3.2) Comparing the asymptotic behavior of the SSA manifold for (H) with the EA manifold for (B), Le., equating the asymptotically correct Michaelis constants, we get Ks = KE
(3.3)
which in terms of rate constants is
(k-l
+ k 2 ) / k l = f-,/L1
(3.4)
This large-s analysis therefor: fixes two of the three possible relations between the (k}and {kJsets of rate constants rigorously in terms of experimentally determinable parameters. Parametrization of Linear Decay at Small Substrate Concentration The third relation is obtained by comparing the manifolds a t very small values of s (close to equilibrium). Note that, having pinned down the hyperbolic law as we did in the previous subsection, we cannot compare velocities a t very small s since this leads to an inconsistency. At very small s most of the enzyme is free and the reaction dynamics is linear, i.e., the reactant concentrations decay exponentially with time. The planar flows of (H) and (B) can therefore be expressed as a linear combination of a fast and a slow decay mode.4dJ0J2dJ3b(In the irreversible examples considered here, the fast decay mode lies outside the physical domain of r2,but for reversible product formation it is possible in principle to prepare the system in the fast mode, i.e., on the fast manifold. Because flow is locally parallel to this fast manifold, it can be found by Taylor series expansion of the corresponding trajectory differential equation about the equilibrium point.) Of more practical importance, we note that the slow decay mode is just the local, linear representation of the slow manifold. In the neighborhood of the origin (equilibrium), h* can be written as h*(s) = as + O(s2)
(3.5)
n
0.2
O
'
0.of 0.0
V I
I
I
I
0.2
0.4
0.6
0 8
J
1.0
s
Figure 3. Common nullcline and the slow manifolds for mechanisms (H) and (B)shown in the phase plane r2. In thislfigure,the (k]rate constants for (H) are those used in Figure 1 and the (k)rate constants for (B)are those used in Figure 2. For these rate constants, (H) and (B)have the same nullcline be = hs, and the flow crosses this nullcline horizontally everywhere. The (B) slow manifold AB= 6*(s), lies above the (H) slow = h*(s), which lies above the common nullcline b&), = manifold h&); cf. eq 3.9. The slow manifolds always lie below the corresponding nullclines for vertical flow, which are not shown here. Note the proximity between the slow manifolds, and that h* and b* coincide as s 0 and
as s
-
-
a.
which leads to a quadratic in a. The slope of the slow manifold in s X h is the positive root a+ of this quadratic, Le. a+ =
eo - Ks + {(eo- K s ) ~+ 4K~eol'/'
(3.6) 2KE Here Ks and a+ are experimentally measured, so that KE may then be calculated. The root a- corresponds to the (unphysical) fast manifold. Similarly, the slow manifold for (B) can be written as b*(s) = &s
+ O(s2)
(3.7) which leads to a quadratic in &. The slope of the slow manifold in s X 6 is the positive root of this quadratic, Le. &+
=
(1
+ ?)eo- R E + {[(I + :)eo 2kE
+ 4KEeo11/2
(3.8)
Here KEand &+ are experimentally measured, so that t may then be calculated. Again 6- corresponds to the fast manifold. The slow manifolds A Hand dBare arranged to coincide as s 0, by setting a+ = &+. Note that a, and &+ could be measured either by initial velocity or transient experiments. Together with eqs 3.1 and 3.3, this third relation between (kl and (k}makes the flows for (H) and (B) comparable in the sense described earlier. Figures 1 and 2 are drawn for such corresponding sets of rate constants. Comparison of Slow Manifolds and Reaction Velocities. Although we do not yet know whether (H) or (B) is correct from the experimental results that were used to fix the above relations (because the parametrization at very large and very small s was set up that way) accurate velocity measurements at intermediate s will now determine different functional forms for the slow manifolds of (H) and (B) and different steady-state velocity expressions. By superposing Figures 1 and 2, it is apparent that the slow manifolds dHand A Blie close together; these are redrawn in Figure 3, together with the common nullcline bE = hs. We can prove that, for s E [0, m), this nullcline and its first iterates hi and b,, under eqs 2.5 and 2.23, respectively, are ordered so that bl(S;[b~])2 h,(s;[hsl)2 ~ E ( S )= hs(s) (3.9) For s E [0, m), it can also be proved that the first-iterate velocities are ordered as follows: (3.10) ui(s;[hs]) 2 uS(s) = e,(#) 2 D,(s;[bE])
-.
8768 The Journal of Physical Chemistry. Vol. 95, No. 22, 1991
Roussel and Fraser
Thus uI for (H) lies aboue the standard approximate reaction velocity 0s = BE, and BI for (B)lies below this velocity. Theoretical arguments concerning flow topology and all our numerical evidence suggest that these orderings also hold for the exact manifolds and velocity functions. However, since time has been removed from the phase space picture, the slow manifolds do not immediately provide a “natural” way of comparing experimental results. But if we use the Eadie-Hofstee (or Woolf) graphical method, i.e., plot u vs. w = u / s or w vs. u, we restore a time scale into the representation by mapping the (exact) slow manifolds A H , B in r2into corresponding (exact) slow velocity manifolds Y H , B in the space Q = w X u. In Q the hyperbolic law is a straight line. The velocity manifolds ‘V contain more information than the phasespace manifolds & precisely because a time scale has been restored. Formally, we write this mapping as dH.B(Cr2)
-
‘VH,B(CQ)
(3.1 1)
We can carry out this mapping in two ways: Either (indirectly) we first find the slow manifolds &H,B in I’ and then calculate the reaction velocity expressions u*(w) for (H), and B*(w) for (B)in n, by change of variables; or (directly) we transform the corresponding trajectory differential equations and solve the resulting functional equation for the velocity manifolds ‘VH,B. We now describe this direct procedure. We explicitly define the variable changes for the mapping 3.1 1 for (H). With the relations u = k2h and w = u / s , i.e., s = u / w , the trajectory differential eq 2.4 can be rewritten as the velocity differential equation
N.B. From now on u and 0, their iterates and fixed points, will denote velocity functions in Q usually we will include the argument w for clarity, so u = u(w). Now as as ds(u,o) = - du - dw (3.13) a0
+ aw
so eq 3.1 3 in eq 3.12 can be explicitly evaluated and rearranged to give the velocity functional equation for (H): U(W)
1
= -(~G(k2 2k2
I
1.Oh
+ + k2(Vmax.- K ~ w +) W)
w
I4. EaditHofstee plots of reaction velocity u vs. w = u / s shown for mechanisms (H) and (B), at the rate constants given in Figures I and 2. The “diagonal” D is the (linear) hyperbolic law for both (H) and (B). All curves have the same intercept (V,, = p& and slope (-Ks= -Re) as D as w 0 (s m). Curve H I is the first iterate of D under eq 3.14 for (H), and H. is the (H)fixed-point function u*(w) of eq 3.14 corresponding to the velocity manifoid VH.All iterates Hjfor j = I, 2, ..., lie above D and so correspond to larger velocitiu than 0;they all intersect the waxis at k2a+ Curve BI is the first iterate of D under eq 3.17 for (B), and B, is the (B) fixed-point function C*(w) of eq 3.17 for the velocity manifold Yg.All iterates Bj for j = 1, 2, ..., lie below D and so correspond to smaller velocities than D, they all intersect the waxis at kze,-,. The B, curves have inflection points at small w and curve upward near the w axis. The Hjcurves have no inflaction point and always curve upward but arc rather straight near thew axis. While the slow manifolds are close together (scc Figure 3), the EadieHofstee velocity manifolds H. and Bo,Le., VH and VB,are quite distinct. F
- -
satisfied: First, B(0) = 8*(0) = PFX;second, provided B in the right-hand side of eq 3.17 is differentiable, i.e., C’exists, the iterate 0(w) on the left-hand side has a zero O(wo) = 0 at the correct position wo = k2eP Starting with the hyperbolic law for (B) as initial approximation, i.e. Bo(w) =
Pmax. - kEw
(3.18)
the first iterate of eq 3.17 can be simplified, by completing the square under the radical, to give
l[~’W(k2+ w ) + kAVmax - Ksw)12 4k2U’w[k2Vma,+ w(vmax- k2Ks) - w 2 K ~ ] ) ’ / ’(3.14) )
where u’ = du/do. The fixed point u*(o)of eq 3.14 is precisely the same as the velocity function u(w;h*) obtained from h* and represents the linelike (slow) velocity manifold YH C Q. Also, the boundary conditions for eq 3.14 are superstably satisfied, because at w = 0 (- s = m), the solution remains fixed at u(0) = u*(O) = Vmx: and at w+ = k2u+ (* s = 0), eq 3.14 can be proved to give the zero u*(w+) = 0 correctly for any differentiable function on the right-hand side. Starting with the hyperbolic law for (H) as initial approximation, i.e. UO(W) vmax - Ksw (3.15) the first iterate of eq 3.14 can be simplified, by completing the square under the radical, to give KSCJ = V,,, - K ~ w ++ 0 ( w 5 ) (3.16) k2 vmax so that, for mechanism (H), there is a positiue additive correction to the hyperbolic law. We can carry out a similar procedure for (B) to get the velocity manifold functional equation uI(o)
(B) there is a negatiue multiplicative correction to the hyperbolic law. Figure 4 illustrates these algebzaic results for the sets of corresponding rate constants (k)and {k]given earlier. The ‘diagonal” D represents the hyperbolic law in the Eadie-Hofstee linear form U(W) V, - KMw (3.20)
so that for mechanism
This provides the starting function for iteration for eqs 3.14 and 3.17, because the comparison of (H) and (B)earlier relates the maximum velocities by V,,, = kzeo =
p,,,,,,= L2kEe0
and the Michaelis constants by KM = Ks = k~
(3.21) (3.22)
(cf. eqs 3.15.3.18, and 3.20). which means uo(w) = Bo(w) = D(O) = D. If eqs 3.14 and 3.17 are iterated for the parameters in eqs 3.21 and 3.22, keeping all terms under the radical, for w E [O, k2eo]we can prove U I ( W ) 1 uo(0) = B,(o) L ;,(a) (3.23) To simplify the labeling in Figure 4, we use the notation BO B0(w), and Ho = uo(w), Le., Bo = Ho= D. For mechanism (H), using Ho = uo(w) as starting function, the first iterate of eq 3.14 is the curve HI = u I ( w ) . For mechanism (B), using BO = 00(o) as starting function, the first iterate of eq 3.17, is the curve B1
Accurate Steady-State Approximations
= uII(o).The characteristic shapes for these first iterates B , and HI also appear in all higher B and H iterates and in the exact answers shown as B, and H, in Figure 4. Thus B and H iterates always remain pinned at u = V , , at w = 0. For mechanism (H), HI and all higher iterates of eq 3.14 are pinned at w = k2a+for u = 0. Successive iterates H I , H2, etc., always lie above the hyperbolic law D and have no inflection point. This sequence converges very rapidly so that the third iterate H3and the exact result H, are indistinguishable on this scale. For mechapism (B) however, B I and all higher iterates are pinned a t w = k2eofor u = 0, which, because of the multiplicative correction, is also the intercept for the hyperbolic law at u = 0 on the w axis. Successive iterates, B , , B2, etc., all lie below the hyperbolic law D and have an inflection point. In fact the first iterate BI and the exact result B, are almost indistinguishable on the scale of this figure. We also remark that although there are many ways to write down a functional equation from a trajectory differential equation, the iteratively stable functional form for (planar) mechanisms (H) and (B) seems always to be that used above, with the derivative of the dependent function but not the function itself on the right-hand side. In suitable examples we believe the geometrical characteristics of the curves B I ,etc., and H I , etc., can be used to differentiate between mechanisms (H) and (B). Furthermore, since their analytic forms are known, using suitable fitting procedures, they allow all the rate constants to be extracted from steady-state experimental data. Although we do not discuss fitting methods here, we have found that simulated annealing3* finds the global minimum, whereas standard minimization methods39 often get trapped in local minima. Finally,we remark that in those cases where the rate constants { k )and {kl lead to very similar velocity manifolds YH,Blying close to D, accurate data would be required to decide mechanism, but even this proximity tells us about the domain in parameter space in which the rate constants must lie. 4. Discussion
Using two rival models of simple enzyme catalysis, we have shown that, in principle, accurate kinetic analysis leads to distinct steady-state behavior. In favorable circumstances this should allow mechanisms to be distinguished and all the rate constants found independently. To illustrate our method, we deliberately chose mechanisms that led to the same behavior (the hyperbolic law) in the usual steady-state analysis, but we resolved the precise steady-state description by choosing functional and graphical methods that emphasized deviations from the standard approximation. Enzyme kinetics results are commonly represented by linear functional forms derived from the hyperbolic law for mechanism (H).These standard forms are widely used even for more complicated mechanisms and are displayed as linear Of these standard representations, the Woolf or Eadie-Hofstee plots are most suitable for accurate steady-state’ analysis, because they emphasize the important crossover region where the order of reaction changes. Departures from linearity are necessary in order to extract rate constants and distinguish mechanism reliably where this can be done. Provided that reaction velocity is a monotonic increasing function of substrate concentration,”~~~ it is legitimate to use nonlinearity in the EadieHofstee plots to extract this additional rate constant information. The analysis given is global in the obvious sense that it is very accurate throughout the entire “slow-velocity” regime regardless of the particular values of the rate constants, because it precisely locates
(38) (a) Corana, A.; Marchesi, M.; Martini, C.; Ridella, S. ACM Trans. on M a l . Software lW, 13,263-280. Kirkpatrick. S.; Gelatt, C. D.; Vecchi, M. P. Science 1983,220,671-680. These methods use the basic Markov chain sampling described in: (b) Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A.; Teller, E. J . Chem. Phys. 1953, 21, 1087-1090. (39) See e&: (a) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes; Cambridge University Press: Cambridge, 1986. (b) Bard, Y. Nonlinear Parameter Estimation; Academic: New York, 1974.
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8169 the slow manifold Jt in phase space r2,or the corresponding velocity manifold Y in the plane s2 of the Eadie-Hofstee plot. The functional iteration method described here has been successfully applied to other two-step mechanisms: the reversible Michaelis-Menten mechanism where E and P may recombine, the Lindemann mechanism, and exciplex formationaM In these planar systems errors less than one-millionth of the SSA are obtained in three iteration^.'^^ The same approach has also been used to solve a vector functional equation for A, involving ordinary derivatives, in three-step mechanisms: the -extended” Michaelis-Menten mechanism where separate, reversibly coupled ES and EP complexes are identified,13band a competitive inhibition mechanism.41 We further note that for these three-dimensional flows, two-dimensional relaxation surfaces can be found by iteration of a scalar functional equation involving partial derivatives,lMso our method can be extended to higher dimensional sets. Even in these more complicated examples, the slow manifold can be located in a few steps with an error that is one-thousandth of the usual SSA. However, here the convergence rate does depend on the topology of the reaction network; for example, the convergence is fast for competitive enzyme inhibition where the network is highly connected but is slow in the extended, linear enzyme network with ES and EP complexes. Convergence of the functional equation can be very effectively accelerated in such cases by a generalized Aitken’s m e t h ~ d . l ~Indeed, * ~ ~ since only differentiation and simple algebraic rearrangement is required, the entire procedure can easily be programmed in an algebraic manipulative language such as Maple.43 We also note that the functional equation method finds the slow manifold 10-100 times faster than the best stiff differential equation solvers in the simple cases we have considered. It may sometimes be useful to have an explicit representation of the time evolution of chemical rate equations,26rather than the implicit phase-space representation of reaction velocity that is normally used in kinetic analysis. For the steady state, our analysis leads immediately to such a representation for progress along the slow manifold M : For any species X,motion on A can always be written in the single-variable form i = fA(x), where fA depends on X and the fixed-point function for M. This differential equation rearranges to give the separated form dt = dx/fA(x). Since fA is known as a series in k2, the solution of the separated form has a Laurent series33 in kz with leading term 0 ( k 2 - l ) . Such series could be treated by the methods of ref 10, for example; also, the pole at k2 = 0 implies that these series have a close connection with the singular perturbation series solutions in ref 9. For example, for mechanism (H) using h = h*(s) in eq 2.2 gives the closed center-manifold differential equation16aS = -kleg + (k,s + k-Jh*(s), which exactly expresses the slow time evolution of the system along AH.The nullcline functions h&) (eq 2.8), h&) (eq 2.10), or their iterates under the functional eq 2.5 may be used as approximations to h*(s) in this centermanifold differential equation for S; as the approximation to h* improves, the explicitly integrated form of this differential equation gives increasingly accurate relations between s and t for motion along AH. Many of the purely mathematical problems associated with steady-state kinetics seem now to be well understood to the extent (40) The Henri-Michaelis-Menten, Lindemann, and exciplex formation mechanisms all have different asymptotic behavior. In the Henri-Michaelis-Menten mechanism, dl is closer to SSA at large substrate concentration, as we discussed in ref 13b. In the Lindemann mechanism, as the excitable molecule concentration increases, dl,and the EA and S A nullclines become asymptoticallly parallel (in I’) with finite separation between them; see ref 13a. In exciplex formation the EA and SSA are parabolas, which separate as the relevant concentrations increase; in this case J(1 is closer to EA, but neither EA, nor SSA is an asymptote, Le., correct to O(conc-I), where conc denotes the independent variable. Fraser, S. J., unpublished. (41) Roussel, M. R.; Fraser, S.J., manuscript in preparation. (42) See e.&: Fox, L. Iterative Solution of Elliptic Finite-Difference Equations. In Fox, L. Ed. Numerical Solution of Ordinary and Partial Differential Equations; Pergamon: Oxford, 1962. (43) Char, B. W.; Geddes, K. 0.;Gonnet, G. H.; Monagan, M. 8.;Watt, S. M. Maple: Reference Manual; WATCOM: Waterloo, 1988.
8770 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991
that the slow manifold J# in r that characterizes such kinetics can be found accurately in a variety of cases. Correspondingly, the geometry of this manifold or the related velocity manifold in the (Woolf or EaditHofstee) space Q can used to extract detailed kinetic information.
Roussel and Fraser
Acknowledgment. This work was supported by the Natural Sciences and Engineering Rcsearch Council (NSERC) of Canada.
We also thank Paul Brumer, John Bunting, Jerry Kresge, Bob McClelland, Mitch Winnik, and Keith Yates for helpful conversations and criticism.
