[22] Diffusion-controlled enzymatic reactions

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

473

[22] D i f f u s i o n - C o n t r o l l e d E n z y m a t i c R e a c t i o n s By MALCOLM E. DAVIS, JEFFRY D. MADURA, JACQUELINE SINES, BROCK A. LUTY, STUART A. ALLISON, a n d J. ANDREW MCCAMMON Introduction The rate of diffusional encounter between reactant molecules in solution sets the ultimate limit on the speed of enzymatic and other reactions. If the reactant molecules are such that subsequent events develop very rapidly when the reactants come into contact, the net rate of the reaction will be equal to the rate of diffusional encounter. The reaction is then said to be diffusion-controlled. 1-5 Examples of such processes can be found among enzymatic and redox protein reactions, the binding of ligands to macromolecules and receptors, and the transport of ions or molecules by channels or other mechanisms. Because diffusion determines the maximum possible rate for such processes, the study of diffusional encounters will become increasingly important as protein engineering yields increasingly efficient systems. As discussed further below, a number of effects may operate to increase or decrease the rate of diffusion-controlled reactions. The magnitude of the rate constant, while suggestive, is therefore not a reliable indicator of diffusion control. Theoretical models of simple diffusioncontrolled reactions show that the rate constants for such reactions are inversely proportional to the viscosity of the solvent. Experimental studies of the viscosity dependency of the kinetics of enzymatic reactions, antigen-antibody binding, and other biological processes have therefore been made to probe for diffusion control. Such studies are of great value, but a variety of potential complications must be considered before the results can be unambiguously interpreted in terms of a diffusion-controlled mechanism. It should be established, for example, that the agents added to vary the solvent viscosity do not alter the rate of the reaction by binding to or altering the structure of the reactants or by changing the free energy of the unbound species. Even after these questions are resolved, a detailed interpretation may prove difficult to develop, since the viscosity depeni G. G. Hammes, "Principles of Chemical Kinetics." Academic Press, New York, 1978. 2 D. F. Calef and J. M. Deutch, Annu. ReD. Phys. Chem. 34, 493 (1983). 30. G. Berg and P. H. yon Hippel, Annu. Rev. Biophys. Biophys. Chem. 14, 131 (1985). 4 j. Keizer, Chem. Rev. 87, 167 (1987). 5 D. A. Case, Prog. Biophys. Mol. Biol. 52, 39 (1988).

METHODS IN ENZYMOLOGY, VOL. 202

Copyright © 1991 by Academic Press, Inc. All rights of reproduction in any form reserved.

474

ENZYMES

[22]

dence may arise from rearrangements within the reactants or the initial complex as well as from the diffusional approach of the reactants. 6 In the present chapter, we describe how computer simulations may be used to shed light on the nature of diffusion-controlled reactions. Such simulations can, in principle, aid in the detailed interpretation of experimental results or in the design o f molecules with prescribed kinetic properties. We also describe methods for calculating electrostatic interactions among solute molecules in solution. The long range of such interactions makes them particularly important in the consideration of diffusional encounters, but the methods we describe can also be used to study the role of electrostatics in determining equilibrium properties such as the stability of molecular folding and association. 7'8 Basic T h e o r y Diffusion Control

The kinetics of e n z y m e - s u b s t r a t e ( E - S ) encounters are typically explained using a model proposed by L e o n o r Michaelis and Maud Menten. 9 Their model is based on a simple e n z y m e - s u b s t r a t e reaction such as that found in Eq. (1): kcat

E + S *~-ES---> E + P

(I)

kr

The catalytic rate V for Eq. (1) can be expressed as V = kcat[ES]

(2)

Assuming the concentration of the encounter complex remains constant during the reaction, namely, d[ES]/dt = 0 (which is typically true for physiological conditions since [El 0 ~ [S]0), then the catalytic rate defined in Eq. (2) can be written in terms of the free e n z y m e concentration [E] and the substrate concentration [S]: V = ~

[E][S]

(3)

where Km, k n o w n as the Michaelis constant, is equal to (kcat + kr)/kD. If kcat >> kr, then Eq. (3) becomes 6j. A. McCammon and S. C. Harvey, "Dynamics of Proteins and Nucleic Acids." Cambridge Univ. Press, Cambridge, 1987. 7 M. Davis and J. A. McCammon, Chem. Rev. 90, 509 (1990). 8 K. A. Sharp and B. Honig, Annu. Rev. Biophys. Biophys. Chem. 19, 301 (1990). 9 L. Michaelis and M. L. Menten, Biochem. Z. 49, 333 (1913).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

V ~ kD[E][S]

475 (4)

and thus the catalytic rate is controlled by the rate of diffusional encounter of the enzyme and substrate. In other words, when the free energy barrier for the transformation of the encounter complex into products is negligible (so that the reaction effectively proceeds instantaneously once the encounter complex is formed) the reaction is diffusion controlled. One is usually concerned with more complicated enzyme-substrate encounters than that just described. In complicated enzymatic pathways, the maximal catalytic rate Vmaxwill depend on several rate constants and not on kcat alone. Therefore, when discussing catalytic rates the pertinent term is kcat/K m, which has dimensions of a second-order rate constant (M -~ sec-~). When [S] ~ Kin, it can be shown that k c a t / K m <~ kD; the reaction therefore cannot be faster than the diffusional encounter of the enzyme and substrate. ~°m The condition kcat/Km ~ 108-109 M -I sec-~ -~ kD has been used to indicate the enzyme has reached kinetic perfection, although it must be remembered that the rate of diffusioncontrolled reactions can be significantly increased or decreased by a variety of factors (see below). Determination of diffusion control in an enzymatic reaction is usually done using an experimental technique called the viscosity variation method. ~zThis method measures the catalytic rate as a function of viscosity. If the catalytic rate exhibits an inverse relationship to the viscosity of the solvent, the reaction is presumed to be diffusion-controlled because kD typically shows this same dependence on viscosity. Again, it must be remembered that other effects, such as reorganization of the initial enzyme-substrate complex, may lead to similar viscosity dependence in the net rate g . 6'13'14 The rate constant kD depends on viscosity because the frequency of diffusional encounter between the enzyme and substrate depends on how rapidly each of these molecules moves about in the solvent. The mobility of a molecule in a particular solvent is characterized by its diffusion constant. For large, neutral molecules, the diffusion constant is given by the Stokes-Einstein equation~4: D = kBT/(6~n)a)

(5)

l0 A. Fersht, " E n z y m e Structure and Mechanism," 2rid Ed. Freeman, San Francisco, California, 1983. H L. Stryer, "Biochemistry." Freeman, San Francisco, California, 1988. i2 S. C. Blacklow, R. T. Raines, W. A. Lira, P. D. Zamore, and J. R. Knowles, Biochemistry 27, 1158 (1988). 13 p. G. Debrunner and H. Frauenfelder, Annu. Rev. Phys. Chem. 33, 283 (1982). 14 j. T. Hynes, Annu. Rev. Phys. Chem. 36, 573 (1985).

476

ENZYMES

[22]

TABLE I TYPICAL DIFFUSION CONSTANTS FOR ENZYMES AND SUBSTRATES E n z y m e or substrate

D°

L y s o z y m e (egg white) b Hemoglobin c Superoxide dismutase d Triose-phosphate isomerase e Na + Superoxide (O~-) f Water g Sucrose h D-Glyceraldehyde i Glycine h

10.4 6.9 5.0 7.1 133 160 285 45.9 82.3 93.3

a D is the diffusion constant in units of 10 -7 cm 2 s e c - i. b j. R. Colvin, Can. J. Chem. 30, 831 (1952). c O. L a m m and A. Poison, Biochem. J. 30, 528 (1936). d W. W. Wilson and M. W. Salin, Mol. Cell. Biochem. 36, 157 (1981). e Estimated using the hydrodynamic radius a = 35.0 ~ which is 0.5 times the longest dimension. f Calculated from molecular dynamics simulations by J. Shen, C. F. Wong, and J. A. McCammon, J. Comput. Chem. 11, 1003 (1990). g R. Hauser, G. Maier, F. Noack, Z. Naturforsh. A: Phys. Phys. Chem. Kosmophys. 21, 1410 (1966); J. S. Murday and R. M. Cotts, J. Chem. Phys. 53, 4724 (1970). h C. R. Cantor and P. R. Schimmel, "Biophysical C h e m i s t r y . " Freeman, San Francisco, California, 1980. i Estimated using the hydrodynamic radius a = 3.0/~ which is 0.5 times the longest dimension.

where kBT is the Boltzmann constant multiplied by the absolute temperature, ~ is the solvent shear viscosity, and a is the radius of the molecule. Typical diffusion constants for various enzymes 15and substrates are listed in Table I. The observed diffusion constants for molecules that are comparable in size to the solvent molecules, or that interact strongly with the solvent, for example, ions in water, may differ from the Stokes-Einstein 15 L. J. Gosting, Adv. Protein Chem. 11, 429 (1956).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

477

value by a numerical factor that is usually in the range of 0.5-2.0.16,17 The effective diffusion constant for nonspherical molecules is well within an order of magnitude of that given by Eq. (5) if a is taken to be one-half the longest dimension of the molecule. Given the diffusion constant, the observed root mean squared displacement of a spherical molecule in a solvent in time t can be determined using the following relationship: (r2) I/2 =

(6Dt) 1/2

(6)

In enzyme-substrate encounters, the motion of a single molecule is not of interest, but rather the relative motion between the enzyme and substrate. To a first approximation, this relative motion is characterized by an effective diffusion constant: Dre 1 = Oenzyme + Osubstrat e

(7)

More accurate descriptions take account of the fact that the relative diffusion of the substrate to the enzyme is influenced by other dynamic effects of the solvent. As the substrate moves, it tends to move the surrounding solvent. The solute-induced solvent motion in turn tends to displace the enzyme. This coupling is known as the hydrodynamic interaction. Equation (7) neglects the effects of hydrodynamic interactions on the relative motion. To account for these interactions, one can rewrite Eq. (7) in the form Drel(r)

= (Oenzyme +

Dsubstrate)/h(r)

(8)

where h(r) increases from the limiting value h(~) = 1 as the separation r between the substrate and enzyme decreases. Approximate expressions for h(r) have been derived using continuum hydrodynamic theory. 16More refined expressions would have to recognize the molecular nature of the solvent, for example, the existence of solvation shells with a variety of relaxation times around each solute species.

Diffusion Equations In diffusion-controlled enzyme-substrate reactions, it is less likely to find a substrate molecule in a given volume element close to the enzyme than in a similar volume further away. This is because nearby substrate molecules will be consumed faster than they can be replaced by the diffusional motion of more distant substrate molecules. The distribution of the substrate molecules at any given time is described by an ensemble-averaged density function p(r,t). The average rate or flux $ at which substrate T6j. T. Hynes, Annu. Rev. Phys. Chem. 28, 301 (1977). i7 p. G. Wolynes, Annu. Rev. Phys. Chem. 31, 345 (1980).

478

ENZYMES

[22]

diffuses across a point, for example, a reaction site, is related to the density function p(r,t) by Fick's first law of diffusionl8: J

=

-DVp(r,t)

(9)

In the presence of external forces acting on the molecules, for example, electrostatic and solvent f o r c e s , 2 Eq. (9) becomes J(r,t) = - D V p ( r , t ) + (kaT)-lDFp(r,t)

(10)

The second term in Eq. (10) reflects the direct force F due to substrateenzyme interactions and the indirect counterforce due to the friction of moving through the solvent. The change in local substrate density with respect to time is equal to the change in flux with respect to position; this is the familiar Fick's second law of diffusionlS: Op(r,t____~) = _ V . $ = VDVp(r,t) Ot

(11)

Equation (1 l) merely reflects the law of conservation: in a constant volume system, the change in a density function, that is, the change in concentration of the substrate, with time must equal the divergence of the flux. That is, if the fluxes into and out of a volume element are different, the density function within that element must change as a function of time. Again, in the presence of external forces acting on the molecules, Eq. (11) becomes Op(r,t) = V. DVp - (kaT)-lV .DFp(r,t) Ot

(12)

These equations describe the relative motion of a substrate with respect to an enzyme ifr is the vector between them, D = Drel(r), and F is a mean force of interaction between the enzyme and substrate. The mean force of interaction in this case refers to an averaging of the solvent-solute interactions. A physical interpretation of the density function p(r,t) can be described based on the following one-dimensional example. An ensemble of systems representing the real system at time 0 is constructed. These systems are then allowed to evolve for time t according to the dynamic laws of the real system. A frequency distribution of the position of the molecule relative to its position at t = 0 for each system in the ensemble is constructed and normalized. This frequency distribution is the desired p(x,t), and the relative probability of finding the molecule in the interval x to x + dx at Is A. Fick, Philos. Mag. 10, 30 (1855).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

479

time t is just p(x,t)dx. An analytical form of p(x,t) for a noninteracting spherical molecule can be obtained by solving

Op(x,t) = D OZP(x't) Ot Ox 2

(13)

with the condition that all the particles are initially at x = 0. The solution to Eq. (13) is

p(x,At) = (4~-DAt)-l/2 e x p [ - (x - xo)2/4DAt]

(14)

which is the probability density of finding a particle at point x after a time At.

Analytic Theories o f Reaction Rates The experimental quantity of primary interest is the diffusion-controlled rate constant. The diffusion equations of the previous section are related to the reaction rate through the flux at the reactive surface. For a single reactant E (the enzyme), the integration of the flux over the reactive surface I = -fs l'dS

(15)

gives the rate of encounter of species S (the substrate) with a single molecule of E. [The negative sign in Eq. (15) is due to the standard definition of dS as the outward normal.] The rate for a given concentration [El is therefore given by R = I[E] as long as the system is sufficiently dilute that interactions between enzyme molecules can be neglected. Since the bimolecular rate constant is defined by R = k[S][E], the rate constant can be expressed as k = I/IS]. Therefore, rates can be determined for any system for which the diffusion equation can be solved. Analytic solutions, however, can only be obtained for systems of high symmetry. Also it is usually the steady state that is of interest, and so we concern ourselves only with that case in which Op/Ot = 0. We also consider only low concentrations, where only the relative motion of a pair of reactants is important and N-body correlations can be ignored. Reaction is assumed to occur on contact, so the concentration will vanish at the reactive surface while going to the bulk value at large separations. For the steady-state case, the integrated flux over any surface surrounding reactant E equals the rate of encounter, that is, I=

-q~ l . d S dS

= constant

This is just a restatement of the law of conservation of matter:

(16)

480

ENZYMES OplOt = - V . J

[22] (17)

In the case of a steady state, ~p/at = 0, so for a spherically symmetric system the integration of Eq. (16) over a sphere of radius r yields I = -47rr2j(r)

(18)

and the rate constant is given by k = -47rr2j(r)/[S]

(19)

where again k is a constant. For reactants with no interparticle forces the flux is related to the density by Eq. (9), which in spherical symmetry reduces to J(r) = - D do~dr

(20)

So the rate constant can be written as k = 4~rrZD[S] -I do~dr

(21)

The boundary conditions on O can be applied in integrating this expression from r = a (the contact separation) to r = ~: kr -2 dr =

f?

41rD[S]-l do

(22)

which, assuming D is independent of concentration and separation, yields k = 41rDa

(23)

or the Smoluchowski law. 4 Extending this to include interparticle forces, the flux is given by Eq. (10), which, again in the case of spherical symmetry, can be written J(r) = - D do~dr - (kBT)-~ D p d U / d r

(24)

k = 4~rrZD[S]-l[do/dr + (kBT)-IO dU/dr]

(25)

As before,

In this case, the integration is achieved by realizing that this equation can be rewritten as k = 47rr2D[S]-1 e- V/kBr[d( O eU/kaT)/dr]

(26)

So

fa~ [S]k eV/kardr 47r D r 2

or

=

[ [S] d(p e v/kBr )

Jo

(27)

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

481

= ~oo eU/kB T

k-1

3a ~

dr

(28)

which is also valid for a separation-dependent D ? For a Debye-H0ckel interaction potential U = qsqE e x p ( - r r ) 47re r

(29)

truncated to only the first two terms in the Taylor series expansion of the exponential U=

q_s4q7Er (e~

K)

(30)

the integral in Eq. (28) can be evaluated, assuming a constant diffusion coefficient and using the variable substitution u = r-~, to obtain k =

(qsqEekBT) 8 qsqEK/4~r~kaT e qsqE/4rr~akBT - 1

(31)

which is the Debye result 4 for ions, expressed in SI units. The simple cases considered above treat the reacting partners as uniformly reactive species with centrosymmetric direct forces and hydrodynamic interactions. Analytic theories have also been developed for more general models which account for orientation-dependent reactivity on one or both species, 19-22 nonspherical reactants, 23 active-site "gating", 24-~6 and noncentrosymmetric direct forces.E7 In addition, analytical-numerical formalisms developed by Wilemski, Fixman, and others ,28-33 Berg and yon Hippel, 3 and also Zientra, Freed, and co-workers 34 have been particularly 19 K. Solc and W. H. Stockmayer, Int. J. Chem. Kinet. 5, 733 (1973). 2o j. M. Schurr and K. Schmitz, J. Phys. Chem. 80, 1934 (1976). El R. Samson and J. M. Deutch, J. Chem. Phys. 68, 285 (1978). 22 D. Shoup, G. Lipari, and A. Szabo, Biophys. J. 36, 697 (1981). 23 p. H. Richter and M. Eigen, Biophys. Chem. 2, 255 (1974). 24 j. A. McCammon and S. H. Northrup, Nature (London) 293, 316 (1981). 25 A. Szabo, D. Shoup, S. H. Northrup, and J. A. McCammon, J. Chem. Phys. 77, 4484 (1982). 26 S. H. Northrup, F. Zarrin, and J. A. McCammon, J. Phys. Chem. 86, 2314 (1982). 27 j. W. van Leeuwen, FEBS Lett. 262, 156 (1983). 28 G. Wilemski and M. Fixman, J. Chem. Phys. 58, 4009 (1973). 29 G. Wilemski and M. Fixman, J. Chem. Phys. 60, 866 (1974). 3o G. Wilemski and M. Fixman, J. Chem. Phys. 60, 878 (1974). 31 M. Doi, Chem. Phys. 11, 107 (1975). 32 M. Doi, Chem. Phys. 11, 115 (1975). 33 C. James and G. T. Evans, J. Chem. Phys. 76, 2680 (1982). 34 G. P. Zientra, J. A. Nagy, and J. H. Freed, J. Chem. Phys. 73, 5092 (1980).

482

ENZYMES

[22]

useful in the study of intramolecular processes in chain molecules such as ring closure 33 and protein domain coalescence. 34 In order to obtain diffusion-controlled bimolecular rate constants using detailed models for actual e n z y m e - s u b s t r a t e reactions, one can employ a simulation p r o c e d u r e ) 5-41 This Brownian dynamics trajectory simulation method is the primary focus o f the present chapter. Electrostatics

Diffusional e n c o u n t e r involves the interaction of molecules at large separations. Electrostatic interactions are therefore particularly important because o f the long-range nature of the Coulombic potential. At small separations, the effects o f other interactions can often be represented by suitable boundary conditions (e.g., reflecting and absorbing). As a result, most attention has been paid to the calculation of electrostatic forces. One of the foremost considerations in simulations is the speed of calculating forces. Especially as the number o f charges increases, calculation of pairwise forces becomes prohibitive. To alleviate this, the force field can be p r e c o m p u t e d and tabulated on a grid. 4z An additional problem with pairwise potentials arises because Brownian dynamics simulations only treat the solute species explicitly. Such potentials cannot accurately model the many-body effects associated with, for example, solvent polarization. Finite difference methods 43-45 that determine the electrostatic potential in an inhomogeneous dielectric provide a better model o f the m a n y - b o d y effects, while providing a tabulated force field. The electrostatic potential in a region of spatially varying polarizability can be calculated from Gauss's law: -V.

eV4~ = p + p~

(32)

where e is the position-dependent permittivity, tb is the electrostatic potential, p is the explicit charge density associated with the enzyme, and PI is 35s. H. Northrup, S. A. Allison, and J. A. McCammon, J. Chem. Phys. 811,1517 (1984). 36S. A. Allison, S. H. Northrup, and J. A. McCammon, J. Chem. Phys. 83, 2894 (1985). 37K. Sharp, R. Fine, and B. Honig, Science 236, 1460(1987). 38T. Head-Gordon and C. L. Brooks, J. Phys. Chem. 91, 3342 (1987). 39S. A. Allison, R. J. Bacquet, and J. A. McCammon, Biopolymers 27, 251 (1988). 40S. H. Northrup, J. O. Boles, and J. C. L. Reynolds, Science 241, 67 (1988). 41j. D. Madura and J. A. McCammon, J. Phys. Chem. 93, 7285 (1989). 42j. A. Ganti, J. A. McCammon, and S. A. Allison, J. Phys. Chem. 89, 3899 (1985). 43j. Warwicker and H. C. Watson, J. Mol. Biol. 157, 671 (1982). I. Klapper, R. Hagstrom, R. Fine, K. Sharp, and B. Honig, Proteins 1, 47 0986). 45M. E. Davis and J. A. McCammon, J. Comput. Chem. 10, 386 (1989).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

483

the charge density of the mobile ions. At equilibrium the local density of the j t h ion species nj is given by a Boltzmann distribution: nj = cj exp( - GJkr~T)

(33)

where cj is the bulk number density of species j and Gj is the local value of its potential of mean force. The potential of mean force can be thought of as an effective potential energy; it describes the variation of the free energy of the system with the position of the ion of interest. The ionic charge density is then t91 = ~ qjcj exp( - Gj/kBT) J

(34)

Ignoring the finite size of the ions and correlation effects, the potential of mean force can be approximated as the charge times the electrostatic potential:

Pl = Z qjcj e x p ( - qjdp/kaT )

(35)

J

For a univalent salt, the ionic strength I is just the bulk concentration, so Eq. (35) becomes Pl

=

- -

2eI sinh(eqb/kB T)

(36)

If the electrostatic energy is small with respect to kBT, the exponentials can be Taylor series expanded and only the first terms kept, leaving PI =

(37)

-- eK2(~

where K (2IeZ/kBTe) 1/2 is the inverse of the Debye-Hi~ckel screening length. Using the full exponential form for PI in Eq. (32) results in the Poisson-Boltzmann (PB) equation, which for univalent salt has the form =

- V . eVck = p - 2 e l sinh(e4~/kBT)

(38)

Using Eq. (37) in Eq. (32) leads to the linearized Poisson-Boltzmann (LPB) equation, which has the form --V"

8Vt~

~" lO - -

8K2(~

(39)

The latter is more commonly utilized because of the difficulties inherent in solving the nonlinear form. Near highly charged species, for example, nucleic acids, however, the electrostatic energy e~ can readily exceed kBT so that the full PB equation should be used in accurate work. Solution of Eq. (39) for simple symmetric forms of s leads to Debye-H0ckel and Tanford-Kirkwood theory. 7 For arbitrary geometries,

484

ENZYMES

[22]

however, numerical techniques must be applied. The most popular methods in this field utilize finite differences, which reduce the partial differential equation to a large, sparse system of linear equations. The form of these equations is most readily seen by integrating Eq. (39) over a volume V with surface S and using the Divergence Theorem: - ~ s eV4~'dS +

fveK2t~dV = fvPdV

(40)

With space divided into a cubic grid of points with spacing h, Eq. (40) can be applied to the h × h × h cubes centered on each point. The surface integral over each face is approximated by the area of the face, times the permittivity assigned to the face, times the difference between the potential at the grid points on opposite sides of the face, divided by the separation. For the mobile ion volume integral, the value of the integrand at the center of the cube is multiplied by the volume of the cube. The explicit charge volume integral is taken to be the charge assigned to the grid point. Thus, the equation for grid point (i,j,k) is

q(i,j,k)/h = e i - ~,j,k [4~(i,j,k) - d~(i - 1,j,k)] + e(i + l,j,k)[4~(i,j,k) - 4~(i + + e(i,j-1,

1,j,k)]

k)[~(i,j,k)- dp(i,j-1,k)]

+ e(i,j+~,k)[d~(i,j,k)-d~(i,j+

1,k)]

+ e(i,j,k-1)[c~(i,j,k)-ch(i,j,k

-

I)]

+ e(i,j,k +l)[eb(i,j,k)-cb(i,j,k

+

1)]

(41)

+ e(i,j,k)rZ(i,j,k)cb(i,j,k)h 2 Various methods are available to solve this large, sparse, system of linear equations. Generally midpoints of grid lines (the lines between grid points) within an atomic radius of the position of an enzyme atom are set to the internal permittivity and the rest to that of the solvent. A range of radii definitions

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

485

have been used, from one to two times the van der Waals radii to a fixed value of 4/~. Born solvation energies suggest that ionic radii would be most appropriate for anions and covalent radii for cations. 46 Excluded volumes defined this way, however, leave "solvent" cavities and clefts too small for a solvent molecule, thus leading to the use of expanded radii. Alternately, the permittivity assignment can be made using solventaccessible surface algorithms. Techniques are under development to provide for a smoother, and so more physically appropriate, transition for the dielectric constant at the surface. Charges are assigned to grid points using smoothing criteria as well. The charge of an atom is distributed over the eight grid points forming the box around the actual charge location. The trilinear weighting function ql = (1 - a)(1 - b)(1 - c)q

(42)

has been shown by Edmonds et al. to minimize the errors in the lower multipole moments of the charge. 47 This relation takes the atomic charge q and assigns the partial charge qi to the ith grid point, with a, b, and c being the fractional distances along the grid axes that the actual charge location is from the grid point. Once the potential is determined for the enzyme, the force on a singly charged substrate at each time step of the dynamics simulation is calculated as the charge multiplied by the negative gradient of that potential. The process of interpolating the gradient from the grid is most easily visualized as placing the charge of the substrate, distributing it onto the grid as above, calculating the force at each grid point, and summing to obtain the net force. A number of shortcomings should be pointed out. The current status of computer hardware limits the grids to at best around 1003 cubic elements, which limits the accuracy of the molecular surface description. Grid spacings must be of the order of 1 ,~. Inherent in discretizing space is the limitation that charges must be at least a grid point away from the surface to avoid "errors" arising from the charge "spreading" into the solvent. The method of force calculation is also an approximation. Treating the substrate as a test charge in the potential of the enzyme does not include polarization charges that may be induced in the enzyme, or at its surface, by the substrate charge. 48 Also neglected are any effects arising from the solvent-excluded volume of an extended substrate. To overcome these approximations in the force calculation requires that the potential 46 A. A. Rashin and B. Honig, J. Phys. Chem. 89, 5588 (1985). 47 D. T. Edmonds, N. K. Rogers, and M. J. E. Sternberg, Mol. Phys. 52, 1487 (1984). 48 M. E. Davis and J. A. McCammon, J. Comput. Chem. 11, 401 (1990).

486

ENZYMES

[22]

be recalculated for each configuration and the stresses at the dielectric boundaries accounted for. Calculation speed and accuracy currently preclude the more accurate approach. In spite of these difficulties, however, reasonable results can be obtained.

Brownian Dynamics Simulation of the diffusional motion of an assembly of interacting solute particles can be accomplished by Brownian dynamics. The algorithm developed by Ermak and McCammon 49 for computing Brownian dynamic trajectories is best introduced by considering the following one-dimensional example. A particle is placed at an initial position x0. From the probability distribution function Eq. (14) the new position at At is randomly chosen. Therefore the first step of a diffusional trajectory is x = x0 + R

(43)

where R is a random number with the statistical properties (R) = 0 and (R 2) = 2DAt. The trajectory is extended to t = 2At, t = 3At, etc., by sampling distributions that are centered at the position found during the previous time step. This procedure yields one representative diffusional trajectory. By computing a large number of such trajectories with different sets of random numbers, one generates a description of how an ensemble of diffusing molecules behaves. The generalization of Eq. (43) for the case of a particle diffusing in three dimensions relative to a second particle and subject to an arbitrary force is r = r0 + (kBT)-lDF(ro)At + R

(44)

where r 0 is the position before a step is taken, ka is Boltzmann's constant, T is the temperature, D is the relative diffusion constant, F(r 0) is the force between the two particles at r 0, At is the time step, and R is a random vector with statistical properties (R) = 0 and (RiRj) : 2DSijAt, (i,j = x,y,z). Equation (44) can be further generalized to treat structured entities, hydrodynamic interactions, e t c . 49-51 The time step At in any region should be small enough that the force on the particle changes by only about 10% or less during any step. Longer time steps can be achieved using a secondorder algorithm 52 instead of the first-order algorithm in Eq. (44). 49 D. L. Ermak and J. A. McCammon, J. Chem. Phys. 69, 1352 (1978). 50 S. A. Allison and J. A. McCammon, J. Phys. Chem. 89, 1072 (1985). 51 S. A. Allison, G. Ganti, and J. A. McCammon, Biopolymers 24, 1323 (1985). 52 A. Iniesta and J. G. de la Torre, J. Chem. Phys. 92, 2015 (1990).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

487

Simulations o f Diffusional Encounter

The Brownian dynamics method described in the previous section can be used to generate diffusional trajectories of a substrate in the field of an e n z y m e target. 36-42'5°'51 To calculate a rate constant, the probability/3 of a substrate diffusing from some given initial separation and reacting with the rigid enzyme is calculated. It is this probability that is used to calculate the diffusion-controlled rate constant. (If the enzyme is small, it may be necessary to allow for rotational diffusion of the target e n z y m e ) 5,4° The space around the enzyme is divided by a spherical surface of radius b into an outer region (r > b) and an inner region (r < b). Trajectories of the substrate are initiated on the surface of this sphere. The value of b is chosen to be sufficiently large that the interparticle forces between the pair are approximately centrosymmetric for r > b. Thus, any effect of charge asymmetry of the species has a negligible effect on diffusion at distances r > b. This condition, which allows for the simplest simulation procedure, can be relaxed to increase computational efficiency. An outer spherical surface of radius q is defined as a truncation surface, at which any trajectories that wander too far from the target particle are terminated. A correction must be applied to the combination probability to allow for the possibility of trajectories recrossing the q surface and reacting with the target particle. As shown by Zhou, q should be chosen to be sufficiently large that the inward reactive flux at this distance would be isotropic for the system being modeled. 53 Calculation of the reaction probability for this system is accomplished using the following algorithm. Starting with the diffusing substrate on the b surface, a step is made using Eq. (44). The trajectory is continued until the substrate reacts with the enzyme or is truncated at the q surface. Several thousand of these trajectories are used to compute the bimolecular diffusion-controlled rate constant: k

=

kD(b)/3[1 - (1 - /3)11]-1

(45)

Where ko(b) is the steady-state rate at which mobile reactants with r > b would first strike the spherical surface at r = b and/3 is the computed reaction probability. With the assumptions that the interparticle potential U is centrosymmetric for r > b, kD(b) is given by [cf. Eq. (28)]

53 H. X. Zhou, J. Chem. Phys.

92,

3092 (1990).

488

ENZYMES

kD(b) = 4zr

(etV(r)/~sr~/r2D)dr

[22]

(46)

The quantity f~, given by l~ = ko(b)/kD(q)

(47)

is the probability that a particle at r = q will eventually return to r = b. The most straightforward way of determining the reaction probability /3 is simply to carry out a large number of trajectories (say, N) and equate /3 to the fraction of trajectories which reach the active site (free diffusion approach). Implicit in Eq. (44) and the statistical properties of the random displacement vector, R, is the assumption that diffusional displacements in At are small, on average, compared to distances from absorbing and/or reflecting boundaries. This condition can be satisfied by choosing At small enough. An alternate approach developed by Lamm and Schulten 54 and Northrup et al. 55 uses the exact one-dimensional solutions of the diffusion equation [Eq. (13)] near boundaries, which makes longer time steps permissible. Northrup et al. showed that the reaction probability/3 can be related to the survival probability for each trajectory i, wi, by N

/3 = 1 - ~ w i / N i=1

(48)

The survival probability is computed from Wi = H Wik k

(49)

where Wik is the survival probability of the kth dynamics step of the ith trajectory. The wik, in turn, are given by Wik = Preact/Preflect

(50)

Where, Preactis the probability of the given step using the reactive boundary condition and Pretlect is the probability of the step using the reflective boundary condition. Analytic expressions for these probabilities and methods for generating steps from these distributions can be found elsewhere. 54-56Northrup and co-workers have applied this method, with per-

54 G. Lamm and K. Schulten, J. Chem. Phys. 78, 2713 (1983). 55 S. H. Northrup, M. S. Curvin, S. A. Allison, and J. A. McCammon, J. Chem. Phys. 84, 2196 (1986). 56 S. A. Allison, J. A. McCammon, and J. J. Sines, J. Phys. Chem. 94, 7133 (1990).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

489

fectly absorbing active sites, to carbonate dehydratase 57 and heme proteins. 58 Recently the partially absorbing case has been studied for simple model systems. 56 Gated reactions can also be studied by simulation. The dynamics of the gate can be described most simply by the scheme where kc and ko are rate constants for closing and opening. Restricting the discussion to stochastic gating, suppose the gate was open/closed during the previous dynamics step and that a dynamics step of duration At is now to be taken. The probability pg that the gate remains open/closed is ~{ko + kcexp[-(ko + kc)At]}/(ko + kc) Pg = ~{kc + k o e x p [ - ( k o + kc)At]}/(k o + k c)

open closed

(51)

where " o p e n " and "closed" refer to the initial states. If pg exceeds a uniform random number selected from the interval (0,1), the gate remains unchanged. Otherwise the gate is switched. In the free diffusion approach, the active site is simply turned offwhen the gate is closed. In the probability distribution approach, Wik is set to 1 when the gate is closed and calculated according to Eq. (50) when it is open. Gated reaction simulations for simple cases have been reported elsewhere. 56

Applications of Brownian Dynamics In this section, we briefly describe two recent studies in order to illustrate the use of the methods described in this chapter. Superoxide D i s m u t a s e

Simulations of the encounter of superoxide (. O2-) with bovine erythrocyte Cu/Zn-superoxide dismutase (SOD) were instrumental in the development of Brownian dynamics techniques for enzyme kinetics.37-39'42'5°'51'59-62 The catalytic reaction involves the rate-limiting diffusion of the negatively charged superoxide substrate into either of two active-site channels on the

57 j. C. Reynolds, K. F. Cooke, and S. H. Northrup, J. Phys. Chem. 94, 985 (1990). 58 S. H. Northrup, J. O. Boles, and J. C. Reynolds, J. Phys. Chem. 91, 5991 (1987). 59 S. A. Allison, S. H. Northrup, and J. A. McCammon, Biophys. J. 49, 167 (1986). 60 K. Sharp, R. Fine, K. Schulten, and B. Honig, J. Phys. Chem. 91, 3624 (1987). 61 S. A. Allison, J. J. Sines, and A. Wierzbicki, J. Phys. Chem. 93, 5819 (1989). 62 j. Sines, S. A. Allison, and J. A. McCammon, Biochemistry 29, 9403 (1990).

490

ENZYMES

[22]

surface of the SOD dimer. (See Refs. 63-65 for comprehensive review.) Although the enzyme carries a net negative charge, positive residues in the active site regions appear to steer the superoxide toward the reactive copper ions .44,66,67 Current detailed models account for electrostatic forces using full 6L68o r l i n e a r i z e d 37'39'6°'62'69 Poisson-Boltzmann potentials as disc u s s e d a b o v e . 7 Brownian dynamics simulations have been successful in reproducing the negative salt dependence of the rate constant 37'39'6°which is observed experimentally 7°-73 and suggests electrostatic control of the reaction. Simulated rates for the current models are about twice the experimentally observed rates. 39 Recent work from our laboratories addressed the roles of individual amino acid charges in the SOD electrostatic steering mechanism. 62 The enzyme representation was based on the 2-,~ resolution crystal structure as solved by Tainer and co-workers. TMWe assigned charges, in whole units of protonic charge, to charged amino acid residues. This simplified charge model yields the same results as a model with complete Amber partial charges 59 and is particularly apt for alteration of single point charges associated with amino acids. For purposes of calculating linearized Poisson-Boltzmann potentials, SOD was centered in a 100 × 100 × 104 ,~3 grid with charges and dielectric constants assigned to grid points at 1-,~ resolution. Dielectric constants for enzyme interior and solvent were 2 and 78, respectively. Superoxide was treated as a sphere containing - le point charge. Ionic strength ranged from 0.005 to 0.437 M. We calculated boundary potentials for the grid by using Debye-Hiickel theory based on a finite ion radius of 30 ,~ and the net charge of the enzyme model under study. A finite difference solution of the linearized Poisson-Boltzmann equation yielded potentials at the remaining grid points. 63 B. Halliwell and J. M. C. Guneridge, Biochem. J. 219, 1 (1984). 64 I. Fridovich, Ado. Enzymol. 58, 61 (1986). 65 j. V. Bannister, W. H. Bannister, and G. Rotilio, Crit. Rev. Biochem. 22, 111 (1987). 66 W. H. Koppenol, in "Oxygen and Oxy-Radicals in Chemistry and Biology" (M. A. Rodgers and E. L. Powers, eds.), p. 671. Academic Press, New York, 1981. 67 E. D. Getzoff, J. A. Tainer, P. K. Weiner, P. A. Kollman, J. S. Richardson, and D. C. Richardson, Nature (London) 306, 287 (1983). 68 B. Jayaram, K. Sharp, and B. H. Honig, Biopolymers 28, 975 (1989). 69 j. Sines, S. A. Allison, A. Wierzbicki, and J. A. McCammon, J. Phys. Chem. 94, 959 (1990). 70 A. Rigo, P. Viglino, G. Rotilio, and R. Tomat, FEBS Lett. 50, 86 (1975). 71 A. Cudd and I. Fridovich, J. Biol. Chem. 257, 11443 (1982). 72 E. Argese, P. Viglino, G. Rotilio, M. Scarpa, and A. Rigo, Biochemistry 26, 3224 (1987). 73 p. O'Neill, S. Davies, E. M. Fielden, L. Calabrese, C. Capo, F. Marmocchi, G. Natoli, and G. Rotilio, Biochem. J. 251, 41 (1988). 74 j. A. Tainer, E. D. Getzoff, K. M. Beem, J. S. Richardson, and D. C. Richardson, J. Mol. Biol. 160, 187 (1982).

[22l

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

491

Dynamics simulations were carried out as discussed above. We usually ran 12,000 trajectories per charge/ionic strength model. The relative diffusion constant Dr was defined using hydrodynamic radii of 28.5 and 2.05 ,~ for SOD and superoxide, respectively, with stick boundary conditions. Trajectories began at 41.5 A from the SOD center of charge and terminated at 500 A. We defined the excluded volume of the enzyme to disallow overlap between nonhydrogen SOD atoms and superoxide, which was assumed to be spherical. Further than 9 ,~ from the reactive coppers, the sum of the radii of SOD atoms and superoxide was assumed to be 4 A. For atoms closer to the copper ions than this, 4 ,~ was too large to allow diffusion down the active-site channel. A value of 2.75 A was used for these active-site atoms instead, accounting for channel flexibility in an approximate way. We assumed that reaction occurred when the distance between the centers of superoxide and one of the two reactive coppers was less than 4.5 ~.. A distance-dependent time step was small when superoxide was close to absorbing or reflecting boundaries and large in the extensive region of free diffusion. When superoxide was within 30 ~, of the SOD center of charge, At (in psec) depended on the distance d (in ~) to the nearest copper: At = 0.02; At = 0.15; At = 0.50;

d < 6.0 6.0 --< d < 7.5 7.5 --< d

(52)

When the distance r (in A) from superoxide to the SOD center of charge was greater than or equal to 30 A, At (in psec) was given by At = 0.0325(r - 25)2; At = 0.130r2; At = 0.130(500 - r) 2 + 32.48;

30 --< r < 100 100 -- r --< 270.75 270.75 < r--< 500.0

(53)

This scheme allows very large steps to be taken (time steps up to about 10 nsec, and displacements of the order of 100 A) when the substrate is far from any boundaries. These values of At, as well as values for b, q, and atomic radii, were chosen to be within ranges where simulation results were independent of small changes in their values. We ran simulations on numerous SOD models which differed only in the enzyme charge distribution. For each model we derived not only reaction rates, but also rates for collision of superoxide anywhere on the SOD surface. The latter, nonspecific collision rates proved to be determined only by the net charge on the enzyme molecule, independent of the details of the point charge distribution. At physiological salt, this long-range monopole-monopole effect was negligible; only 6% variation in collision rates was observed in the + 2e to 10e net charge range. Rates -

492

ENZYMES

[22]

-0.5

v

o

-1.0

J

o.o

&2

i

J

i

o.6

4T FIG. 1. Simulation results, k/k c is the ratio of the rate of superoxide encounter at copper ions to the rate of superoxide encounter for the entire SOD surface. Standard deviations are about +0.02 log units. I is in molar units. The dashed line shows results from unmodified SOD simulations. Symbols are results from simulations of SOD with charge modifications: + I or - 1 charge introduced at valine-5 ( x ); neutralization of glutamate- 131 (A); neutralization of lysine-134 (11); neutralization of arginine-141 (V).

for the specific encounter of superoxide at the active site, however, were indicative of electrostatic steering roles for the residues under study. For a large region around the active site channels, a + le or - le modification of the point charges associated with one residue of the primary sequence (one amino acid in each subunit) resulted in only a 10-15% increase or decrease in reaction rates; the only exceptions are residues that are located in or very close to the active-site channels. This effect simply disappeared as the modifications were made more distant from the channel. Figure 1 illustrates this, as well as the more dramatic effects seen at three residues. We have plotted the collision efficiency factor, that is, the ratio of reaction rate to nonspecific collision rate. This factor is more indicative of electrostatic steering effects than is the reaction rate, since long-range monopole effects (as reflected in the nonspecific collision rate) have been factored out. The three residues exhibiting greater effects on the steering mechanism are arginine- 14 I, lysine- 134, and glutamate- 131. Arginine- 141 lies deep in the SOD active-site channel. It has been experimentally modified, resulting in large decreases in SOD a c t i v i t y . 71'75-79 This residue and lysine-134, 75 D. P. Malinowski and I. Fridovich, Biochemistry 18, 5909 (1979). 76 D. COCCO, t . Rossi, D. Barra, and G. Rotilio, FEBS Lett. 150, 303 (1982).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

493

which is at the lip of the channel, appeared to play major roles in the electrostatic steering mechanism in previous theoretical w o r k . 37'66'67 The third residue, glutamate-131, is situated near lysine-134 at the lip of the channel. Although it is negatively charged, it has also been suggested as a contributor to the electrostatic steering mechanism. 67 We see large decreases in the reaction rate in the cases of neutralization or reversal of the two positive charges, consistent with roles in the steering mechanism. A pronounced increase in rate is seen in the case of glutamate-131, in agreement with a previous simulation at zero ionic strength 37 and not consistent with any electrostatic facilitating role.

Triose-Phosphate Isomerase Another recent application of the Brownian dynamics method was to the diffusional encounter of o-glyceraldehyde phosphate (GAP) with triose-phosphate isomerase (TIM).41 These simulations made use of the new, integrated computer program U H B D for electrostatic and diffusional calculations. 8° The diffusion-controlled rate constant for the reaction has been carefully analyzed experimentally by Knowles and co-workers. ~2The three-dimensional Cartesian coordinates for chicken muscle TIM were taken from the 1988 Brookhaven National Laboratory Protein Data Bank tape. 81 Hydrogens attached to polar atoms (N and O) were added using CHARMm from Polygen. 82Charges were taken from the OPLS parameter set. 83 The overall charge of TIM was zero. The substrate was treated in this preliminary work as a sphere with isotropic reactivity, hydrodynamic radius 2.5 A, and charge - 2 e . The volume excluded to the center of the diffusing GAP molecule was determined using twice the standard van der Waals radii for the protein atoms. 84The hydrodynamic radius for TIM was taken to be 35.0 ,~, which is approximately one-half the largest dimension 77 F. Marmocchi, I. Mavelli, A. Rigo, and G. Rotilio, Biochemistry 21, 2853 (1982). 78 C. L. J. Borders, J. E. Saunders, D. M. Blech, and I. Fridovich, Biochem. J. 230, 771 (t985). 79 W. F. Beyer, Jr., I. Fridovich, G. T. Mullenbach, and R. Hallewell, J. Biol. Chem. 262, 11182 (1987). 8o M. E. Davis, J. D. Madura, B. A. Luty, and J. A. McCammon, Comput. Phys. Commun. 62, 187 (1991). 81 F. C. Bernstein, T. F. Koetzle, G. J. B. Williams, E. F. Meyer, Jr., M. D. Brice, J. R. Rogers, O. Kennard, T. Simanouchi, and M. Tasumi, J. Mol. Biol. 112, 535 (1977). 81 B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus, J. Comput. Chem. 4, 187 (1983). 83 W. L. Jorgensen and J. Triado-Rives, J. Am. Chem. Soc. 110, 1657 (1988). 84 L. Jarvis, C. Huang, T. Ferrin, and R. Langridge, " U C S F MIDAS User's Manual." University of California, San Francisco, 1986.

494

ENZYMES

[22]

of the enzyme. TIM was placed at the center of a 58 x 58 × 58 grid of points which had a grid spacing of 1.6 .~. The electrostatic potential was calculated for this grid by solving the linearized Poisson-Boltzmann equation via a method developed in this laboratory. 45 This equation was solved using dielectric constants of 2 and 78 for the interior and exterior regions of the enzyme, respectively, and zero ionic strength. The b and q surfaces had radii of 80 and 300 A, respectively. Reaction was assumed to occur whenever the substrate diffused to within 6.0/~ of HG of residue Ser-2 I0 or H of residue Gly-232 in the active site of either monomer of the dimeric enzyme. A flexible loop (residues 168-177) that may intermittently cover the active sites was held fixed in an open conformation in the preliminary study. Finally, a variable time step algorithm was used in this study with the time step determined by At = 0.02psec, At = 2.00 psec, At = 0.10 psec,

r < 50 or r > 2 9 5 85 < r < 290 otherwise

(54)

where r is the center-to-center distance between the target and substrate. These values were chosen so that (1) the forces do not change significantly on displacement and (2) the displacements are small compared to distances from the absorbing and reflecting boundaries. It would probably be possible to use larger steps both in the middle region and near the q surface. The diffusion-controlled rate constant was determined using the Brownian dynamics method described in the previous section. Five independent runs of 400 trajectories each were used to determine the rate constant. One set of 2000 trajectories had the interparticle forces set to 0 (i.e., no electrostatics), while a second set of 2000 trajectories included full electrostatics. A summary of the reaction probabilities and the diffusion-controlled rate constants calculated from the Brownian dynamics simulations is found in Table II. The rate constant determined from the five runs for the system with no electrostatic interactions is 0.306 -- 0.114 × 101°M-J sec -j, and that for the system with such interactions is 1.48 -- 0.13 × l0 ~°M - l ; the error ranges correspond to the standard deviations of the five runs for each system. Both of these rate constants are higher than the experimental rate constant of 4.8 × 108 M - 1sec- 1.85The higher rates in the calculations are likely due in part to the neglect of the loop motion, which should " g a t e " the reaction; the generous conditions assumed for reaction, in terms of the substrate location and orientation; and the neglect of hydrody85 D. C. Phillips, P. S. Rivers, M. J. E. Sternberg, J. M. Thorton, and I. A. Wilson, Biochem. Soc. Trans. 5, 642 (1977).

[22]

DIFFUSION-CONTROLLED ENZYMATIC REACTIONS

495

T A B L E II CALCULATED RATE CONSTANTS FOR TRIOSEPHOSPHATE ISOMERASE--D-GLYCERALDEHYDE PHOSPHATE REACTION /3

run a

1 2 3 4 5 1 2 3 4 5

N o electrostatics 0.0475 0.0200 0.0525 0.0475 0.0700 Electrostatics 0.275 0.242 0.215 0.235 0.260

kb

0.307 0.130 0.338 0.307 0.448 1.64 1.46 1.31 1.42 1.56

A run is defined as 400 trajectories. h k is in units o f 10 I° M -j s -I.

namic interactions. More detailed simulations are planned to explore these effects. The preliminary results make clear, however, that the electrostatic field of the enzyme tends to steer the diffusion of the GAP substrate toward the active sites. This is reminiscent of the steering effects observed in Brownian dynamics simulations of superoxide diffusion to the active sites of superoxide dismutase) 7-39,5°'86 As in the case of superoxide dismutase, it is likely that the effectiveness of the steering will decrease with increasing ionic strength, resulting in decreases in the computed rate constants. These effects, and the possibility that the electrostatic fields might lead to orientational as well as translational steering of the substrate, will also be explored in future simulations. Current and Future Directions As with all computer simulations, accuracy and efficiency are two primary concerns that are often at odds. Efficiency favors larger time steps, but the accuracy of the constant force assumption is reduced for electrostatic, intramolecular, and exclusion forces. The variation of the electrostatic forces with position is the smallest of these and therefore is the only limiting factor in volumes away from exclusion or reactive sur86 j. A. M c C a m m o n , Science 238, 486 (1987).

496

ENZYMES

[22]

faces if diffusing entities are rigid. Exclusion forces are not typically included explicitly, but are handled by rejecting steps that would violate the excluded volume conditions. In either case, however, large steps might allow the particles to jump by each other. Similarly, the truncation sphere needs to be represented reasonably accurately, and so large steps may be precluded near this surface. A number of different approaches are being explored to overcome these obstacles. For the problem of surfaces (exclusion, reaction, and truncation) the simplest approach is to provide for variable time steps. In this manner time steps can be made larger in regions away from critical surfaces. A more interesting approach is to handle the diffusion near a surface by an approximate analytic method, namely, by sampling displacements from non-Gaussian distributions derived for simple moUels of the region near the surface. This also has the advantage of permitting calculations of the probability of reaction from each step, providing better statistics for low-probability events. For the intramolecular forces, implicit algorithms are being developed that permit more accurate integration of vibrational motion at larger time steps than standard methods. As details of the enzyme are more accurately represented through the electrostatics and exclusion grid, it also becomes desirable to add increased structure to the diffusing unit. As a first step, the structure of the substrate can be represented by a collection of spherical, hydrodynamically interacting subunits. These subunits are not necessarily selected to represent each atom of the substrate but should be chosen to represent most accurately the domains or functional groups of the diffusing unit. With this increase in resolution of the substrate comes the ability to add orientational constraints to the reaction criteria. This will allow examination of the possibility that the enzyme may "steer" the substrate into a productive orientation for reaction. This increase in detail also brings in some added complexities. As mentioned, the subunits will interact hydrodynamically. Ermak and McCammon 49 have shown that this implies that the stochastic terms will be correlated. Generation of correlated Gaussians is accomplished by Choleski factorization of a matrix which grows as the number of particles squared. Therefore, the generation of the correlated Gaussians can quickly become a bottleneck. A second problem which arises is in maintaining the structure of the diffusing unit. If allowed to diffuse according to Eq. (44) with no interactions, the subunits would eventually diffuse apart. Therefore, some method must be employed to constrain the subunits to the proper geometry. One method is to use a harmonic potential between neighboring subunits and incorporate this added force into the equations. A second method

[23]

CALCULATIONS

OF RELATIVE BINDING AFFINITIES

497

is to attempt to return all the constrained subunits to prescribed separations after an unconstrained step. This can be accomplished with a shaketype algorithm appropriate to Brownian dynamics, s7 Finally, the area of error estimates, both systematic (model dependent) and statistical is being investigated. Acknowledgments Work at the University of Houston was supported by grants from the National Institutes of Health, the Robert A. Welch Foundation, the Texas Advanced Research Program, and the U.S. Office of Naval Research. Supercomputer access was provided by the National Center for Supercomputing Applications, the San Diego Supercomputer Center, and HNSX Supercomputers. Stuart A. Allison would like to acknowledge the National Science Foundation for support through a Presidential Young Investigator Award. J. A. McCammon is the recipient of the 1987 George H. Hitchings Award from the Burroughs Wellcome Fund. Brock A. Luty is supported by the NIH Molecular Biophysics Training Program. 87 S. A. Allison and J. A. McCammon, Biopolymers 23, 167 (1984).

[23] T h e o r e t i c a l C a l c u l a t i o n s o f R e l a t i v e Affinities o f B i n d i n g

By T. P. STRAATSMAand J. A. MCCAMMON Introduction The analysis and prediction of enzyme activity by means of computer simulation have become possible in recent years as a result of advances in theoretical and computational chemistry. The new computational tools allow the calculation of fundamental thermodynamic and kinetic quantities such as those displayed in Schemes I and II. Here, E, S, and P represent enzyme, substrate, and product, respectively, X represents a reaction intermediate, I represents a competitive inhibitor, K m is the Michaelis constant (commonly approximated by the equilibrium constant for dissociation of the enzyme-substrate complex), and K I is the equilibrium constant for the dissociation of inhibitor. Calculations of rate constants for initial binding, for example, kl and k4, can be accomplished in some cases by simulations of the corresponding diffusional encounters. This is discussed elsewhere in this volume. 1 The present chapter treats the calculation of thermodynamic quantities such as the equilibrium constants in Schemes 1 M. E. Davis, J. D. Madura, J. Sines, B. A. Luty, S. A. Allison, and J. A. McCammon, this volume [22].

METHODS IN ENZYMOLOGY, VOL. 202

Copyright © 1991 by Academic Press, lnc. All rights of reproduction in any form reserved.