J . Phys. Chem. 1986,90, 3901-3905
3901
FEATURE ARTICLE Diffusionai Dynamics of Ligand-Receptor Association J. Andrew McCammon,* Department of Chemistry, University of Houston-University Park, Houston, Texas 77004
Scott H. Northrup, Department of Chemistry, Tennessee Technological University, Cookeville, Tennessee 38505
and Stuart A. Allison Department of Chemistry, Georgia State University, Atlanta, Georgia 30303 (Received: February 7, 1986)
Biological molecules generally have complicated shapes and charge distributions. Their binding sites for ligands or substrates are often geometrically restrictive and may display fluctuating steric and electrostatic properties. Detailed studies of biomolecular associations that involve such complications can be carried out with the aid of new Brownian dynamics simulation methods. These methods should also prove useful in studies of diffusion-influencedprocessq in other areas of chemistry. This article provides a simple and pedagogic treatment of diffusional bimolecular association and describes some of the initial applications of Brownian dynamics to biochemical systems.
Introduction The first step in many biological processes is the diffusional encounter of ligand and receptor molecules.’ Among the many important types of ligand-receptor pairs involved in biological action, familiar examples include substrate-enzyme, hormonereceptor, and antigenantibody pairs. The rate of initial diffusional encounter determines or controls the overall rate of a number of biological processes. For example, a number of enzymes catalyze reactions of their substrates at diffusion-controlled rates. Such enzymes are said to have been “perfected” by evolution. Because the net catalytic rate is already limited by the diffusional transport of substrates to the active site of such an enzyme, there is no selective advantage associated with increased reactivity within the enzyme-substrate complex.* In most biochemistry textbooks, the rate of diffusional encounter between ligands and receptors is discussed within the framework of the traditional Smoluchowski t h e ~ r y . ~In this theory, the reaction partners are assumed to be spherical, with uniformly reactive surfaces. Apart from their instantaneous reactivity upon contact, interactions between the reaction partners are ignored. For the enzyme-substrate association
Structural fluctuations in the ligand or receptor may result in a time-dependent reactivity upon contact. Effects such as these influence biological processes in a variety of interesting ways. For example, ligands and receptors that have complementary distributions of electric charge may be “steered” into productive orientations during diffusional encounters, resulting in increased reaction rates. In such cases, it can be said that molecular recognition begins before the ligand and receptor actually come into ~ontact.~ Because of the complications outlined above, reactions between biological molecules may be diffusion-controlled and yet have rate constants that differ by several orders of magnitude from what is expected on the basis of the simple Smoluchowski theory. Analytic theories that take some of these complications into account have been developed and are described in several recent Howeyer, it appears likely that the most ’detailed descriptions of diffusional biomolecular encounters will be provided by new computer simulation methods.*-15 In such a simulation, one computes representative diffusional trajectories of the molecules of interest and then analyzes these trajectories to determine
k
E+S-ES (1) and typical molecular sizes, the Smoluchowski theory yields a bimolecular rate constant k of about 1O’O M-’s-l. It is frequently supposed that bimolecular reactions with rate constants more than about 2 orders of magnitude below this value are not diffusioncontrolled. Of course, biomolecular encounters are in fact much more complicated than what is pictured in the Smoluchowski theory. Typically, only a small part of the molecular surface is active. Electrostatic and other interactions between the reaction partners will generally produce attractive or repulsive forces that can increase or decrease the rate of diffusional encounter.
1984,80, 1517. (9) Allison, S. A.; Srinivasan, N.; McCammon, J. A.; Northrup, S. H. J . Phys. Chem. 1984,88,6152. (10) Allison, S. A.; McCammon, J. A. J . Phys. Chem. 1985, 89, 1072. (1 1) Allison, S. A.; Ganti, G.; McCammon, J. A. Biopolymers 1985, 24, 1323. (12) Ganti, G.; McCammon, J. A.; Allison, S. A. J . Phys. Chem. 1985,
(1) McCammon, J. A.; Harvey, S. C. Dynamics ofProreins and Nucleic
89, 3899. (13) Allison, S. A.; Northrup, S.H.; McCammon, J. A. J. Chem. Phys.
Acids; Cambridge University: London, in press. (2) Knowles, J. R.; Albery, W . J. Acc. Chem. Res. 1977, 10, 105. (3) Cantor, C. R.; Schimmel, P. R. Biophysical Chemistry, Parr IIR Freeman: San Francisco, 1980.
1985,83,2894. (14) Dickinson, E. Chem. SOC.Rev. 1985, 14, 421. (15) Northrup, S. H.; Curvin, M. S.;Allison, S.A,; McCammon, J. A. J. Chem. Phys. 1986,84,2196.
(4) Neumann, E. In Structural and Functional Aspects of Enzyme Catalysis: Eggerer, H., Huber, R., Springer-Verlag: Berlin, 1981. ( 5 ) Calef, D. F.; Deutch, J. M. Annu. Rev. Phys. Chem. 1983, 34, 493. (6) Berg, 0. G.; von Hippel, P.H. Annu. Reu. Biophys. Chem. 1985, 14, 131. ( 7 ) Keizer, J. Acc. Chem. Res. 1985, 18, 235. (8) Northrup, S.H.; Allison, S.A.; McCammon, J. A. J. Chem. Phys.
ms.;
0022-3654/86/2090-3901$01.50/00 1986 American Chemical Society
3902 The Journal of Physical Chemistry, Vol. 90, No. 17, 1986 the mechanism and rate constant of molecular association. This approach may be viewed as a dijfwional analogue of the molecular dynamics simulation of inertial motions in proteins and nucleic acids.] This article provides an introductory overview of the theory of biomolecular diffusion-controlled reactions in dilute solution. After a review of the basic concepts of such reactions, the new trajectory approach to such reactions is described and illustrated. The present treatment is intended to be simple and pedagogic. Diffusion Control A quantitative measure of the extent of diffusion control in a simple bimolecular reaction of the form A
+ B ckD; l AB k,
kf
McCammon et al. molecule. For small molecules (comparable to the size of the solvent m ~ l e c u l e s or ) ~ molecules ~ that interact strongly with the solvent (e.g., ions in water):, the observed diffusion constant may differ from the Stokes-Einstein value by a numerical factor that is usually in the range 0.5-2.0. For nonspherical molecules, the effective diffusion constant is usually well within an order of magnitude of that given by eq 6 if a is taken to be half the longest dimension of the molecule. Accurate values of the diffusion constant can be calculated for such molecules by numerical techniques.21-22 In bimolecular reactions, one is not interested in the motion of a single molecule, but rather in the relative motion of two reactant molecules. To a first approximation, this relative motion is characterized by an effective diffusion constant
products
Drel = DI+ D2
can be obtained by consideration of the rate constants kD, k,, and kP Here, AB indicates the diffusional encounter complex with the reactants A and B properly positioned for subsequent transformation into products. Under typical reaction conditions (e.g., when A and B represent enzyme and substrate, respectively, and the initial concentrations satisfy [A], > k,, ken = kD;the rate of the reaction then equals the rate of diffusional encounter, and the reaction is said to be diffusion-controlled. Physically, diffusion-controlled reactions are obtained when there is a negligible free energy barrier for the transformation of the encounter complex into products. The reaction then effectively proceeds instantaneously once the encounter complex is formed. In the opposite case where k, >> kf, the effective rate constant becomes k , ~ =([AB]/[A][B])kf. = These two cases can usually be distinguished experimentally by measuring the dependence of keffupon the solvent viscosity, since kD typically depends inversely on solvent viscosity (see below) and kf is typically independent of solvent viscosity. It should be noted, however, that kf may exhibit a viscosity dependence in some processes. 1*16,17 Diffusion “Constants” The frequency of diffusional encounter between reactants A and B clearly depends on how rapidly each of these molecules moves about in the solvent. Considering for the moment a single spherical molecule, this mobility is characterized by the molecular diffusion constant D. For such a molecule, the root mean square (rms) displacement observed during a time t is (6Df)1/2.18 The diffusion constants of glycine and the small enzyme lysozyme are respectively 9.3 X 10“ cm2 s-l and 1.0 X 10” cm2 s-l a t 20 O C in water. If t = lo-” s, the rms displacements of these molecules are respectively 2.4 and 0.77 A. For large, neutral molecules, the diffusion constant is given by the Stokes-Einstein equation” D =k,T/(6~7~)
(6)
Drei(4 = (Dl + D 2 ) / h ( 4
(8)
where h(r) increases from the limiting value h(m) = 1 as the separation r of the reactants decreases. The exact magnitude of the effect depends on the detailed nature of the solvent-reactant interactions and the concentration of the solution. For proteins with no electrostatic interactions in dilute aqueous solution, hydrodynamic interactions are expected to reduce the diffusional encounter frequency by perhaps 30%.23 Diffusion Equations In a diffusion-controlled bimolecular reaction, one is less likely to find a reactant molecule in a given volume element close to the reaction partner than in a similar volume farther away. This is because nearby reactants will be consumed more rapidly than they can be replaced by the diffusional motion of fresh rea~tants.~’ The distribution of reactant molecules a t any given time is described by an ensemble-averaged density function p. Consider for simplicity the diffusion of a single spherical molecule in one dimension. Imagine that one prepares a large number (ensemble) of systems that are representative of the real system a t time 0. These systems are allowed to evolve for a time t according to the dynamic laws of the real system. The position of the molecule in each system is recorded, and the frequency with which molecules are found in discrete intervals along the x axis is calculated. This frequency, normalized in any convenient manner, gives the desired density function; the relative probability of finding the molecule in the interval, x, x dx, at time t is just p(x,t)dx. The average rate at which molecules diffuse across a point x in the above ensemble is given by a flux functionj(x,t). The flux depends on the gradient of the density and any external force F acting on the molec~les:~
+
j(x,t) = -Dap/dx
+ (kBT)-IDFp
(9)
According to this equation, molecules will exhibit a net tendency to diffuse to the left if they are more concentrated on the right or if there is a force pulling the particles to the left. The rate of change of density a t a given point is just the difference between the flux into and out of the point:
where kBT is the Boltzmann constant multiplied by absolute temperature, 7 is the solvent viscosity, and a is the radius of the (16) Debrunner, P. G.; Frauenfelder, H. Annu. Rev. Phys. Chem. 1982, 33, 283. (17) Hynes, J. T. Annu. Rev. Phys. Chem. 1985, 36, 573. (18) Friedman, H. L. A Course in Statistical Mechanics; Prentice Hall: Englewood Cliffs, NJ, 1985.
(7)
where D, and D2 are the individual diffusion constants of the two reactants. In addition to the purely kinematic effect included in eq 7 , more accurate descriptions consider the effects of hydrodynamic interactions on the relative m ~ t i o n . ~As a reactant molecule moves, it tends to move the surrounding solvent. The solvent motion in turn tends to displace the reaction partner. This hydrodynamic interaction slows the relative motion of the reactants. Thus, the relative diffusion “constant” takes the general form
ap(x,t)/at = -aj/ax (19) (20) (21) 81. (22) (23)
(10)
Hynes, J. T. Annu. Reu. Phys. Chem. 1977, 28, 301. Wolynes, P. G. Annu. Reu. Phys. Chem. 1980, 31, 345. Garcia de la Torre, J.; Bloomfield, V. A. Q. Reu. Biophys. 1981, 14, Allison, S.A,; McCamrnon, J. A. Biopolymers 1984, 23, 167. Wolynes, P. G.; McCammon, J. A. Macromolecules 1977, 10, 86.
The Journal of Physical Chemistry, Vol. 90, No. 17, 1986 3903
Feature Article
Equations 9 and 10 can be combined to give the diffusion equation in its familiar form for one-dimensional systems: For a molecule diffusing in three dimensions, eq 9-1 1 become j(r,t) = -DVp
+ (kBT)-'DFp
(12)
ap(r,t)/dt = -V.j
(13)
ap(r,t)/at = V*DVp- (kBT)-'V.DFp
(14)
The above equations also describe the relative motion of two spherical reactant molecules if r is the vector from the "target" molecule to the "mobile" reactant and D = Dm,(r). In this case, j(r,t).dA gives the ensemble-average rate of diffusion of mobile reactants through any small element dA in the diffusion domain. Also in this case, F is usually a mean force of interaction between the reactant molecules in which solvent-molecule effects appear in averaged form (e.g., through a dielectric constant if the reactant molecules are electrically ~ h a r g e d ) . ~ Analytic Theories of Reaction Rates Simple models of diffusion-controlled reactions can be studied by analytic r n e t h o d ~ . ~ - ' , For ~ ~ the case of dilute solutions, it is sufficient to consider the relative motion of a single reactant pair. Then the encounter frequency and reaction rate of spherical molecules can be obtained by solving eq 14 supplemented by appropriate boundary conditions. Under the usual steady-state approximation, the density and flux become time-independent quantities, p(r) and j(r), and the reaction rate is just the normal component of the flux integrated over any surface that encloses the target molecule. In the Smoluchowski theory, it is assumed that the reactant spheres have uniformly reactive surfaces and that reaction occurs instantaneously upon collision. Forces and hydrodynamic interactions between the reactants are neglecttd. Thus, the equation to be solved is the familiar Laplace equation DV2p(r) = 0 subject to the boundary conditions
(15)
Arc) = 0 (16a) corresponding to total depletion of reactants at the collision distance r,, and A m ) = Po (16b) where po is the density of reaction partners that would be found everywhere if the target were unreactive. Integration of the corresponding flux j(r) = -DVp(r) over a surface around the target yields eq 4 with
kerf= kD = 4?rr$ (Smoluchowski) (17) In this equation, the usual units of p [B]and [ A ]are molecules per cubic centimeter, and those of r, and D are centimeters and centimeters squared per second, respectively. Inclusion of hydrodynamic interaction and a centrosymmetric force
F = -VU(r) where U is a potential of mean force, yields5
(18)
Approximate analytic expressions have also been obtained for models that depart from the above in one respect or another. These include models for spherical reactants with reactive "patches", in some cases including simple intermolecular interactions, models that include the general effect of binding site fluctuations, etc. References to this work may be found in other recent reviews and paper^.^'"^ Such work allows for useful quick estimates of the (24) Wilemski, G.; Fixman, M. J . Chem. Phys. 1973, 58, 4009.
rates of diffusion-controlled reactions. However, for more detailed analyses that include realistic descriptions of the irregular surface topography, multiple interaction sites, etc. that characterize biological molecules, one must rely upon computer simulations of the diffusional motion of correspondingly detailed models. Brownian Dynamics Simulation Brownian dynamics is a method that allows one to simulate the diffusional motion of an assembly of interacting solute molecules.26 The method is most easily illustrated by considering its application to a single spherical molecule moving along the x axis in the absence of any applied force. Suppose that the molecule is initially at the point x o . Then, from the analytic solution of eq 11, the probability density for finding the molecule at a point x after a time At is just the Gaussian function :,At) = (4?rDAt)-'/*exp[-(x - ~ ' ) ~ / 4 D A t l (20) The first step of the simulation is completed by choosing the position of the molecule at random from this probability distribution. In other words, the algorithm for the first step of a diffusional trajectory is x=xo+R (21) where R is a random number with the statistical properties p(
(R)= 0 ( R 2 )= 2DAt (22) The trajectory is extended to times t = 2At, t = 3At, etc. by repeated application of the above algorithm, taking the position at the start of each step to be that chosen in the preceding 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. For example, the ensemble-average position distribution at time At approaches the result given in eq 20. The above method has been extended to systems that are not amenable to exact analytic analysis. For example, for the relative motion of a pair of spherical reactant molecules that interact by a noncentrosymmetric force F(r) but not by hydrodynamic interactions (cf. eq 14), one has the algorithm'0J1-26 r = ro (kBT)-'DF(ro)At R (23) where
+
+
(R) = 0 (RiRj)= 2D6,,At (24) and i and j label the Cartesian components of R. In the generation of trajectories with eq 23, the time step At must be short enough that F(r) = F(ro). Similar algorithms are available for the cases of many spherical particles with arbitrary forces and hydrodynamic interactions,26nonspherical particles modeled as flexible or rigid assemblies of spherical subunits,22and spheres with translational and rotational motion explicitly coupled by hydrodynamic intera~tion.~' Because the Brownian dynamics algorithm given in eq 23 applies to motions of particles in the absence of boundaries, it will give inaccurate results when particles are near reactive or reflecting surfaces unless the time step is kept exceedingly small. In the past, this restriction hampered efficient generation of the trajectories required for calculation of reactive rate constants, but improvements on the free diffusion algorithm for motion near boundary surfaces have been developed. l5 In particular, Lamm and SchultenZ8 proposed some very useful one-dimensional Brownian algorithms for generating the displacements of particles near reactive and reflecting boundaries. These are based on the exact distribution functions for diffusion in one dimension analogous to eq 20 except that they describe diffusion in the presence of constant forces near reactive and relecting boundary surfaces. These were extended by Lamm to treat three-dimen(25) Temkin, S. I.; Yakobson, E. I. J . Phys. Chem. 1984,88, 2679. (26) Ermak, D.L.; McCammon, J. A. J . Chem. Phys. 1978, 69, 1352. (27) Dickinson, E.; Allison, S. A.; McCammon, J. A. J. Chem. Soc., Faraday Trans. 2 1985,81, 591. (28) Lamm, G.; Schulten, K. J . Chem. Phys. 1983, 78, 2713.
3904 The Journal of Physical Chemistry, Vol. 90, No. 17. 1986
sional spherically symmetric diffusion ~ a s e s . 2 ~In a recent paper,I5 it has been shown how the use of such algorithms increases the efficiency of rate constant calculations by allowing much larger step sizes to be taken near boundaries. Despite the progress outlined above, there remain a number of unsolved problems related to the implementation of Brownian trajectory methods for calculating rate constants. Boundary curvature places limitations on time steps, even when improved algorithms are being used that treat diffusion near absorbing or reflecting surfaces, since these algorithms are essentially onedimensional while the surfaces generally considered in bimolecular reactions are curved. Angular constraints on reactions (reactive patch type surfaces) introduce discontinuities that are not handled by algorithms based on one-dimensional diffusion distributions. Irregular surface topography (such as buried active sites) require special provisions which must be treated on a case by case basis. Also, it would be desirable to incorporate hydrodynamic interactions between the reactant and surface. Simulated Reactions In principle, the Brownian dynamics approach offers an attractive route to the detailed analysis of diffusion-controlled reactions. For example, for the case of an enzymatic reaction in dilute solution, one can imagine generating diffusional trajectories of the substrate in the field of a fixed enzyme target. From the frequency of collisions of properly oriented substrates with the active site of the enzyme, one would then calculate the rate constant. In practice, this approach encounters an obvious difficulty. Many trajectories may wander quite far from the enzyme. To determine the ultimate fate of such trajectories (Le., whether they return and lead to reaction or escape reaction altogether) would require infinitely long simulations. Northrup et al. have recently devised a method to circumvent this difficulty.s In this method, one essentially divides the diffusion space around the enzyme into two regions. The %mer" region is finite and comprises that volume adjacent to the enzyme in which the interactions are complicated and best dealt with by numerical techniques. The "outer" region is of infinite volume but is everywhere far enough from the enzyme that the diffusional behavior can be described analytically. Trajectories then need be computed only in the inner region. Any trajectory that reaches the outer region is truncated; the contribution of such trajectories to the rate is given by an analytic correction factor. The factors needed in the rate constant calculation can be introduced successively as follows. Let the target molecule (e.g., an enzyme) be surrounded by a spherical surface of radius b, where b is large enough that the surface lies just outside the inner region. Then the rate constant can be written as
Here, k,(b) is the steady-state rate at which mobile reactants with r > b would first strike the surface and is given by an expression such as eq 17 or 19 with r, = b. The factor p, is the probability that a reactant starting at r = b, and free to diffuse in the inner and outer regions, will react rather than escape. The probability of escape from r = b is then 1 - P,. Now let the b surface be enclosed by a larger spherical surface of radius q. The probability of escape from r = b can be written as 1 - P, = (1 - P)Pm (26) where 1 - ,8 is the probability that a reactant starting at r = b will reach r = q before reacting and P, is the probability that a reactant starting at r = q will escape rather than react. Although evaluation of P, formally still requires consideration of reactant motion in the inner and outer regions, it can be shown that the contribution of the inner region dynamics can be expressed in terms of @.* Here, we simply quote the result P , = (1 - Q)[l - (1 - p)n]-' (27) where (29) Lamm, G. J . G e m . Phys. 1984, 80, 2845.
McCammon et al.
Combining eq 25-28, one obtains k = kD(b)P[1 - (1 - P)Q]-'
(29)
Thus, to calculate the rate constant for a bimolecular diffusion-controlled reaction in dilute solution, one need only compute a large number of trajectories of one reactant diffusing in the vicinity of the other fixed reactant. Trajectories are initiated at r = b and terminated upon reaction or upon reaching r = q. The fraction of trajectories that react is 6, and this quantity yields the rate constant through eq 29. Analysis of the trajectories also provides information on the mechanistic details of the reaction, e.g., whether reactants tend to be steered into productive collision geometries during the diffusional encounter. The general approach described here can be extended in many ways. Some of these, e.g., consideration of partial reactivity at the active sites and the use of intermediate boundary surfaces to simplify the consideration of active site models that differ in their details,13 have been described elsewhere. A particularly powerful extension of the rate constant procedure has been possible by the employment of exact reflecting and absorbing boundary algorithms near collision surfaces.'* This technique is based on a statistical weighting method which allows the simultaneous solution of several different reactive boundary problems in one simulation. Trajectories are generated in a reference unreactive system, in which every encounter of particles leads to reflection, and trajectories are eventually truncated at large separation r = q. For a given boundary model, a statistical weight is calculated for each trajectory. This weight is the product of ratios of the value of the true (reactive) distribution function for each step divided by the reference (reflecting) distribution function. The value of the statistical weight of the trajectory at truncation is its escape probability (1 - @),which is subsequently related to a rate constant. Escape probabilities may be tallied for an entire series of reactive boundary conditions during one simulation. This technique is valuable for performing angle averaging over trajectory starting angles or for determining rate constants as a function of the size or intrinsic reactivity of the reactive surfaces on the particles. Example: Superoxide Dismutase The diffusion-controlled reaction of superoxide (OF) catalyzed by the enzyme superoxide dismutase (SOD) is anomalous in a t least two respect^.^^*^' First, the reaction rate ( k = 2 X lo9 M-' s-' at T = 25 OC and ionic strength of 0.02 M) is larger than what one might expect, given that (a) the reactants have net electrostatic charges of the same sign (-4 for SOD and -1 for 02-, in units of proton charge) and hence suffer some electrostatic repulsion and (b) the two specific binding pockets for substrate account for only 0.1%; of the surface area of the dimeric enzyme. Second, the reaction rate decreases with increasing ionic strength of the solvent at moderate salt concentrations, whereas the opposite behavior would be expected for electrostatically repelling reactants. Chemical and structural studies of SOD suggest that these anomalies arise from a highly asymmetric charge distribution on SOD that acts to steer 02-into the active site^.^^,^] Allison et al. have initiated a series of Brownian dynamics simulation studies to examine the diffusional encounter dynamics in increasingly realistic models of the SOD-02- system.'0,'' In the initial studies, SOD has been represented by a sphere of 30-8, radius. Two small circular reactive patches (each subtending an azimuthal angle of loD)were defined on opposite sides of the enzyme. A cluster of five charges was placed in the enzyme to produce an electrostatic field that had monopole, dipole, and quadrupole components outside the enzyme similar to those produced by all 76 ionic groups in the native structure. The 02molecule was represented as a sphere of 1-5-8, radius with a single central charge. The dielectric constant was taken to be 78 (30) Cudd, A.; Fridovich, I. J . B i d . Chem. 1982, 257, 11 443. (31) Getzoff, E. D.; Tainer, J. A.; Weiner, P. K.; Kollman, P. A.; Richardson, J. S.; Richardson, D. C. Nature (London) 1983, 306, 287.
The Journal of Physical Chemistry, Vol. 90, No. 17, 1986 3905
Feature Article
I Figure 1. Schematic illustration of the S O H 2 - reaction model used in the initial simulation studies.IoJ1Crosses indicate positions of charges. Active sites are indicated by the dark caps on the SOD sphere; Bo = loo.
throughout the system. Figure 1 depicts the model used in the initial studies. The diffusional encounters were simulated by computing trajectories of 0,in the field of the fmed enzyme. Trajectories were typically initiated at b = 300 A and terminated upon collision with an active site or with a truncation surface at q = 500 A. This simple model system successfully reproduced the qualitative features of the S O W 2 - system. In particular, the rate constant was found to increase by nearly half when the noncentral components were added to the monopole field of the enzyme. Also, when salt effects were represented by a crude Debye-HUckel type model, it was observed that the reaction rate first increased and then decreased to a plateau as the solvent ionic strength was increased. The initial increase corresponds to the screening of the long-range monopole repulsion, while the subsequent decrease corresponds to the screening of the shorter range noncentral “steering fields”. The initial calculations are being extended in a variety of ways. Ganti, et al. have developed a simulation method in which the net intermolecular forces are computed once and saved at lattice points throughout the diffusion space; these forces then need not be recomputed during the subsequent trajectory calculations.12 This method has been used to consider the diffusion of 02-in the field of partial charges associated with each of the 2196 nonhydrogen atoms of SOD.The results in this case are similar to those obtained by using the quadrupole model. Other studies that incorporate dielectric inhomogeneity, more sophisticated treatments of solvent ions, and topographic details of the enzyme surface are also in progress. Example: Cytochrome c Quite recently, Northrup and co-workers have begun a series of studies on the association of the electron transport protein cytochrome c with its redox partnerseg2 In the initial work, cytochrome c is modeled as a spherical dipole with orientationdependent reactivity; the protein reacts with a spherical monopolar partner that has isotropic reactivity. Simulations were performed at various ionic strengths and were treated according to extended Debye-Hiickel theory and with varying reactive patch size on the protein. A series of these model dipolar proteins was considered in which the reactive surface is variously disposed relative to the dipole moment vector. This initial system includes some of the features thought to be important in the kinetics of cytochrome electron transfer reactions;33in particular, the dipole moment of cytochrome c is thought to play a significant role in the association rates of the electron-transfer complexes. Two models of the protein dipole were compared, a surface-charge-dipole model, in which opposite charges are placed on the surface of the protein, and a (32) Northrup, S. H.; Smith, J. D.; Boles, J. 0.; Reynolds, J. C. L. J. Chem. Phys., in press. (33) Margoliash, E.;Bosshard, H. R. Trends Biochem. Sci. (Pers. Ed.) 1983, 8, 316.
point-dipole model, in which the charge separation distance is infinitesimal but the dipole moment magnitude is held constant. This study also represents the first application of the extended methodology now available in Brownian dynamics simulations. The calculation (1) takes advantage of the increased efficiency of trajectory-initiating surfaces in close proximity to the reaction target, but with the added feature of performing the angle averaging and truncation correction exactly; (2) uses the correct absorbing and reflecting step algorithms near boundaries; and (3) utilizes the statistical weighting method to enable a treatment of a large set of boundary problems in a single simulation (thus enabling the comprehensive study of varying patch size and dipole orientation). The most significant findings of this work may be summarized as follows. For the special case of an isotropically reactive dipolar protein, it is found that the electrostatics have virtually no effect on the diffusional rate constant, despite the fact that the equilibrium-averaged interaction potential between the species is attractive and would yield a factor of 20 increase in the equilibrium population of the reactive complex. On the other hand, when the dipolar protein has a small axially symmetric reactive surface (of 1Oo extent), the surface-charge-dipole produces a 5-fold increase in the rate when the dipole vector is along the reactive patch axis and retards the rate by a factor of 10 when the dipole is oppositely disposed. A significant electrostatic effect persists through the physiological range of ionic strength. The positioning of dipolar charges relative to the protein surface has an important effect on the electrostatic influence of the dipole. Thus, in this system the electric field in’the neighborhood of the reactive surface of the protein appears to play a more important role than the long-range steering effect of the dipole moment. The surface-charge-dipole enhancement of the rate is comparable to that calculated for two small isotropic monopolar charged species, one of which has a charge of -le like the ligand while the other has a charge of +2e, equal to the charge at the active site of the protein. This again indicates the greater importance of the field strength near the active surface compared to the global electrostatic properties. Variations of the cytochrome c dipole relative to its reactive surface have been achieved experimentally by chemical modific a t i o n ~and ~ ~could also be obtained by site-directed mutagenesis. As with superoxide dismutase, work is continuing to develop more detailed theoretical models for comparison with such experimental systems. Concluding Remarks The introduction of simulation methods is expected to open the way for detailed study of a wide range of diffusional phenomena in cellular and molecular biology. For example, the calculations on enzymes can easily be extended to predict the effects of amino acid sequence changes. Internal dynamics of the enzyme (e.g., binding-site fluctuations that modulate the reactivity) and substrate can be incorporated in the calculations in a natural and straightforward way. Other simple association phenomena (e.g., antigen-antibody, hormone-receptor) can be handled in the same manner. Moreover, the specific methods outlined here are based on quite general principles and could be adapted for application to nonbiological systems. The influence of diffusion on fluorescence energy transfer is one such problem currently being studied. Acknowledgment. This work has been supported by the N I H (at UH and TTU), the Robert A. Welch Foundation (at UH), and the Texas Advanced Technology Research Program (at UH). Acknowledgment is made to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for partial support (through TTU) of this work. J.A.M. is the recipient of a Camille and Henry Dreyfus Teacher-Scholar Award. S.H.N. is the recipient of an NIH Research Career Development Award. S.A.A. is the recipient of a Camille and Henry Dreyfus Grant for Newly Appointed Faculty in Chemistry and an N S F Presidential Young Investigator Award. (34) Koppenol, W. H.; Margoliash, E.J . Biol. Chem. 1982, 257, 4426.
