The role of electrostatics in the binding of small ... - ACS Publications

3342

J . Phys. Chem. 1987, 91, 3342-3349

The Role of Electrostatics in the Binding of Small Ligands to Enzymes Teresa Head-Gordon and Charles L. Brooks III* Department of Chemistry, Carnegie- Mellon University, Pittsburgh, Pennsylvania 1521 3 (Received: November 14, 1986)

In this paper we present a study of the effects of electrostatic models on the rate of ligand-protein association. The rates of association were calculated by using Brownian dynamics on a system representing superoxide dismutase (SOD)/superoxide, although the general form of these models will be applicable to other systems. Four analytical dielectric models were investigated, including a continuous dielectric model and a dielectric discontinuity model. Salt effects on the rate of binding were also examined through the introduction of extended Debye-Hiickel theory. Our results suggest that these simple models are sufficient in elucidating the primary mechanisms of electrostatic “steering” in diffusion-controlledreactions. The calculated ionic strength dependence is consistent with the experimentally observed trend. In addition, the more realistic dielectric discontinuity model provides an enhancement in the association rate of 45%-60% over models which do not account for the protein-solvent interface.

I. Introduction The rate of catalysis of many enzyme-substrate systems appears anomalously large when one considers that the reactive site is buried in the bulk of the protein, and is small relative to the size of the protein. One example of such a system is superoxide dismutase (SOD), which oxidizes the superoxide radical to hydrogen peroxide and oxygen.’ The small activation energy of these enzymatic reactions suggests that the reaction mechanism is diffusion-controlled and that proteinsubstrate recognition may actually begin in the initial stages of diffusion. Simulation studies have attempted to provide an understanding of the complex interaction of macromolecules diffusing through solution, starting with the Smoluchowski model of two uniformly reactive spheres with no forces,2 and continuing with the work of Debye on diffusion in the field of centrosymmetric force^.^ Subsequent model studies have incorporated more realistic molecular detail, which includes internal motions of the p r ~ t e i n , ~ J details of the active site,6 and complex charge distributions emEnzyme-substrate interactions which bedded in the include hydrodynamic effectsgJOand orientational dependence on one or both species”-13 have also been investigated. Relevant examples of model studies which include a complex charge distribution embedded in the enzyme are those by Allison and co-workers’ and Northrup et a].*. These studies have focused on the rate of binary encounter in models representative of the SOD/superoxide system and cytochrome c, respectively. Their results suggest that the inclusion of long-range electrostatic forces (monopole), with shorter range orientating forces (quadrupole) in a spherically modelled enzyme with reactive patches, may enhance the rate two- to twentyfold over the Smoluchowski model. The inclusion of an extended electrostatic dipole in the latter study of cytochrome c gave a fivefold enhancement when the dipole vector was directed along the axis of the reactive patch and a tenfold retarding effect in the rate when the dipole was directed perpendicular to the site, relative to the model of Solc and Stockmeyer, which treats the dipole as a point source.14 As a whole, these studies relay the importance of the following features for understanding the resultant activity, namely the disposition of charges in the protein and the need for an accurate representation of the electrostatic forces. In this paper we have taken this problem a step further to distinguish the protein interior from that of the solvent for a spherical enzyme. The simplification of spherical symmetry allows us to analytically solve for the electrostatic forces which include a dielectric discontinuity at the protein-solvent interface for any given charge distribution, while maintaining an accurate representation of the electric field distant from the protein surface. Previous work which has incorporated a dielectric discontinuity into model chemical studies include that of K i r k w ~ o d , ’ ~whose J~ *Author to whom correspondence should be addressed.

0022-3654/87/2091-3342$01.50/0

potential functions we utilize; Gilson et ai.,” who presented a general analysis of the electrostatic problem applied to proteins, as well as numerical calculations for proteins with an arbitrary geometry;18 and very recently Honig and c o - w o r k e r ~ , ’ who ~,~~ examined specific effects of the dielectric discontinuity and ionic strength on the electrostatic fields around superoxide dismutase. Our results for the SOD/superoxide system show that the inclusion of the discontinuity enhances the rate 45% to 60% over the continuous dielectric model and about 40% over the diffusion-controlled rate for the native protein. This enhancement is consistent with conclusions drawn by Honig et a1.,20 based solely on the electrostatic potential they calculated for the enzyme. Another anomalous feature peculiar to the SOD/superoxide system is the relationship between kinetic activity and added salt. A naive argument based on the overall negative charge of the protein (4e7 and the superoxide molecule (le-), would predict an increase in the rate with increasing ionic strength due to the screening of the repulsive forces. However, experimental studies2’ of activity of the native enzyme with addition of 0.1 and 0.5 M NaCl have shown the opposite to be true; the rate decreases as the ionic strength increases in this range of salt concentrations. In this study, we have modified the basic electrostatic models to include the effect of added salt using extended Debye-Huckel theory. Our results are in qualitative agreement with the decrease in rate with increasing ionic strength found experimentally2’ for concentrations of salt greater than 0.005 M NaC1. In the remainder of this paper we present the methods and (1) Fridovich, I. Adu. Enzymol. Relat. Areas Mol. Biol. 1974, 41, 33. (2) Smoluchowski, M. Phys. Z. 1916, 17, 557. (3) Debye, P. Trans. Electrochem. SOC.1942, 82, 265.

(4) Szabo, A.; Shoup, D.; Northrup, S.; McCammon, J. J . Chem. Phys. 1982, 77, 4484.

(5) Northrup, S.; Zarrin, F.; McCammon, J. J . Phys. Chem. 1982, 86, 2314. (6) Samson, R.; Deutch, J. J . Chem. Phys. 1978, 68, 285. (7) Allison, S.;Ganti, G.; McCammon, J. J . Phys. Chem. 1985,89, 3899. ( 8 ) Northrup, S.;Smith, J.; Boyles, J.; Reynolds, J. J . Chem. Phys. 1986, 84, 5536. (9) Friedman, H. J. Phys. Chem. 1966, 70, 3931. (10) Deutch, J.; Felderhof, B. J . Chem. Phys. 1973, 59, 1669. (1 1) Solc, K.; Stockmayer, W. J. Chem. Phys. 1971, 54, 2981. (12) Schmitz, K.; Schurr, J. J . Phys. Chem. 1972, 76, 534. (13) Schurr, J.; Schmitz, K. J. Phys. Chem. 1976, 80, 1934. (14) Solc, K.; Stockmayer, W. Int. J . Chem. Kinet. 1973, 5, 733. (15) Scatchard, G.; Kirkwood, J. Phys. Z. 1932, 33, 297. (16) Kirkwood, J. J. Chem. Phys. 1934, 2, 351. (17) Gilson, M.; Rashin, A.; Fine, R.; Honig, B. J . Mol. Biol. 1985, 183, 503. (18) Warwicker, J.; Watson, H. J . Mol. Biol. 1982, 157, 671. (19) Sharp, K.; Fine, R.; Schulten, K.; Honig, B., submitted for publication in J. Phys. Chem. (20) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; Honig, B. Proteins 1986, I, 47. (21) Cudd, A.; Fridovich, I. J . B i d . Chem. 1982, 257, 11443.

0 1987 American Chemical Society

The Journal of Physical Chemistry, Vol. 91, No. 12, 1987 3343

Binding of Small Ligands to Enzymes

that of position relaxation, one may approximate eq 1 by the Brownian equation of motion, which in turn leads to a Brownian, or diffusion-limited, algorithm for numerical integration

1+z

aDmno

Rm(t

+ At) = Rmo(t)+ E-At n aRn + En-

Dmno Fn0'k T At + Q m o ( f ) (2)

I

where Dm,Ois the relative diffusion tensor for particles m and n, Tis the temperature, Qmo(t) is a random displacement vector, and the superscript o denotes the evaluation of the quantity at the previous time step. Hydrodynamic interactions were ignored in this study, in which case eq 2 breaks down to an algorithm for the motion of two Brownian particles

\

where Dd = (al + a2)/4~qala2 is the relative diffusion coefficient, 7 is the solvent viscosity, a , and a2 are the radii of the spherical particles, and Q is a vector of Gaussian random numbers where

( Q , )= 0 (QiQj)

Figure 1. Superoxide dismutase model. The enzyme is represented as a sphere of radius 28.5 8, with two reactive patches defined by the +/z axis and a polar angle, 4, of loo. The fixed charges denoted by X reproduce the monopole, vanishing dipole, and quadrupole of the truncated 76-charge SOD model7 The vector R is from the origin to the superoxide ligand; r is the vector from the origin to one of the fixed charges in the enzyme and B the angle between R and r.

results of our studies of the model SOD/OT system (Figure 1). The Brownian dynamics algorithm of Ermak and McCammon22 and the trajectory analysis method of Northrup et al.23were used in this work and are briefly reviewed in section 11. In section 111, the results for the rate of association for the SOD/superoxide system for both a continuous dielectric and dielectric discontinuity, and as a function of ionic strength, are presented. We also include in this section some discussion of the role of electrostatic potentials in determining the calculated rates. In section IV we conclude with a discussion of our model results and the applicability of such models in the diffusive regime. In the Appendix A, we derive the hierarchy of electrostatic models used in the present study. Explicit results for the electrostatic potential in the presence of added salt and a dielectric discontinuity at the protein surface are given. 11. Methods and Models In this section we describe the Brownian dynamics and trajectory analysis methods used in the calculation of the association rate. Both of these methodologies have been presented in detail in earlier work, so the present discussion is brief. The interested reader is referred to ref 2 2 and 23. In addition to the methods, we also provide a description of the model systems we have studied. Brownian Dynamics. To describe the trajectories of particles in a condensed phase, we begin with the Langevin equation of motion (LE) for an interacting pair of particles

miv,(t) = -mipvi(t) + S , ( t ) + F,(t)

(1)

where i labels one of the pair of Brownian particles, m denotes its mass, and F , ( t ) represents the systematic force, such as an electrostatic lorce. The host medium degrees of freedom are represented in an average way through a frictional force, mi/3v,(t), and a random force, Si(t). The frictional force represents the drag exerted on the Brownian particle by the solvent molecules, where p is the friction coefficient. The random force is a nonsystematic contribution which represents random displacements of the Brownian particle due to collisions with many host particles. If the time scale of the momentum relaxation is very different from (22) Ermak, D.; McCammon, J. J . Chem. Phys. 1978, 69, 1352. (23) Northrup, S.; Allison, S.;McCammon, J. J . Chem. Phys. 1984,80, 1517.

= 2Dre1At

(4)

This algorithm was developed by Ermak and McCammon and is discussed more thoroughly in ref 2 2 . We have also employed the variable time step algorithm described by Allison et al.,7 with time steps ranging from 0.01 ps near the reactive patch to 1.O ns away from the reactive patch. In addition, a viscosity of 1 CPwas used to represent liquid water at room temperature in all calculations. Trajectory Analysis. Extracting a diffusion-controlled, bimolecular rate constant from eq 3 would require the evaluation of a large number of trajectories at initially large relative separation distances and is not feasible in practice. We therefore resort to the trajectory analysis method of Northrup et aL2j In this method the configuration space for a reactive pair of particles is divided into an inner sphere of radius b such that the forces between the pair are isotropic for R > b (for species with reactive patches the angular dependence of the reactive flux through the surface defined by R = b must vanish). The rate constant is then given as a product of a conditional flux, k(b),that the pair achieve a separation b for the first time, and a probability P, that the pair will ultimately react and not escape to infinite separation.

K = k(b)P

(5)

Because the forces are isotropic in the region of space R > b, k(b) can be evaluated analytically from the Smoluchowski relationship.2

The probability, P,is defined as a product of a purely diffusional quantity B, and a reactive quantity A . The former is the probability that particles diffusing in an infinite domain will reach separation R = b and subsequently reach the reactive surface at least once. The latter quantity is the probability that, once the reactive surfaces of the particles are touching, the required activation energy barrier for reaction is overcome. For diffusioncontrolled reactions the activation energy is negligible and A = 1. Equation 5 now reads K = k(b)B, (7) This is still not a tractable entity to simulate for the reason that a particle diffusing in an infinite domain may do so indefinitely. Therefore a truncation sphere is defined at some separation distance q, where q > b. This offers a means of truncating a trajectory after the separation distance reaches R = q. The probability B , is now considered to be proportional to (1 - a), the probability that the particles achieve a separation distance R = q, and (1 - y ) , the probability that once the separation distance of the particles reaches R = q the particles will escape

Head-Gordon and Brooks

3344 The Journal of Physical Chemistry, Vol. 91, No. 12, 1987 to infinite separation. Since the particle may diffuse back and forth across the surface R = q, B , is not just a simple product of the two probabilities. Instead, all recrossings must be counted in computing B,, the probability of reaching the surface R = b. When these recrossings are taken into account, the final expression of the rate becomes23 K = k(6)6/(1 - (1 - 6 ) ~ )

TABLE I: Rates of Ligand Association for the SOD/Superoxide Model with 10-deg Reactive Patches,at Zero Ionic Strength" charge model cont dielec (KIK,,) dielec discont ( K / K d ) -4 MQ 0 MQ 4 MQ

0.888 1.735 3.850

1.404 2.720 5.529

(8)

The rates are normalized by the diffusion-limitedresult (no forces and IO-deg reactive patches), Kd = 1.15 X IO') A3/s.

(9)

TABLE 11: Rates of Ligand Association for Electrostatic Models with Uniform Reactivity at Zero Ionic Strength' charge model cont dielec ( K / K o ) dielec discont ( K / K o )

The quantity 6 is defined as the ratio of successfully reactive trajectories to the total number of trajectories, all of which begin at an initial separation distance R = b. It is 6 which is calculated from the Brownian dynamics simulation. Model System. Superoxide dismutase is modelled as a sphere of radius 28.5 A (Figure 1) with two reactive patches whose centers define the z axis. The extent of each reactive site is defined by a 10-deg polar angle. A five-charge electrostatic model, which represents reasonably well the monopole and quadruple moments of SOD,7 was employed. This charge distribution is analogous to that used by Allison and co-workers.' The reduced charge model is parametrized so that the monopole, dipole, and quadrupole moments for the full enzyme are reproduced; the dipole is essentially zero for the full e n ~ y m eleaving ,~ only the monopole and quadrupole. The parameters used to describe the quadrupole were the polar and azimuthal angles, O and 4, and the product of the magnitude of the charge with the square of the moment arm. It was found that a superposition of two quadrupoles was necessary to accurately reproduce the quadrupole moment for the full e n ~ y m e .The ~ parameters describing the two quadrupoles are O1 = Oo, 41 = Oo, and q , r I 2= 1205 A2 for the first, and O2 = 18.26O, 42 = 149.09O, and q2r2 = 1197 A2 for the second where r , = r2 = 10 A. The monopole was then reproduced by placing a charge at the origin so that the overall charge of the enzyme (-4e in the case of the SOD enzyme, for example) was realized. By referring to the potential energy maps (see section 111, Figures 2 and 3), the values of b and q used in the trajectory analysis and dynamics were taken to be 200 and 400 A, respectively. The initial positions for the superoxide radical anion were uniformly sampled over the surface of the sphere defined by the radius b. The electrostatic models considered in this paper include a continuous dielectric (CD) model where the protein dielectric constant is assumed to be the same as water (e = 78); a dielectric discontinuity (DD) model where the protein is distinguished from the solvent by a dielectric constant of 2; the above continuous dielectric model with the addition of extended Debye-Huckel theory (CD-DH) to account for salt effects exterior to the protein; and the dielectric discontinuity model considered above with inclusion of extended Debye-Huckel theory (DD-DH) to account for salt effects exterior to the protein. For the CD and DD models a total of 25 000 trajectories were analyzed for all cases. For the CD-DH and DD-DH models the number of trajectories was increased to 60 000. The potential energy maps displayed in Figures 2 and 3 and 5-7 were generated by calculating the potential on a grid with 2.0-A spacing, using the exact potentials outlined in Appendix A. It was found that the 2.0-A maps used in displaying the potentials provided sufficient detail. However, we note that the continuous, analytical, potentials were used in the calculations. Features of the electrostatic fields in the protein interior were eliminated when the figures were made. 111. Results

In this section we present the results for calculations on the rate of the diffusional encounter in the SOD/superoxide system. In Tables I and 11, the rate constants calculated from our analysis of Brownian dynamics simulations of the SOD/superoxide system using the C D and DD electrostatic models in the absence of salt are presented. Table I presents the rate constants for the case of reactivity confined to two 10-deg reactive patches on the protein surface. Table I1 presents rate constants for a uniformly reactive

-4 MO 0 MQ 4 MQ

0.652 1.065 1.638

0.660 1.086 1.604

"The rates are normalized by KO = 4~D,,,a,where a is the radius of the enzyme. surface. The level of precision was found to be i5% for these calculations. There are four broad classes of calculations represented here: (1) a protein with reactive patches and inclusion of a dielectric discontinuity (DD); (2) a protein with reactive patches and a continuous dielectric (CD), (3) a protein with uniform reactivity and inclusion of a dielectric discontinuity, and (4) a uniformly reactive protein and a continuous dielectric. Each of these classes contain the following subcategories: an electrostatic charge distribution representing a monopole and quadrupole (MQ) of overall charge -4,0, and 4. Note that the MQ -4 charge model is representative of the native protein. The association rates for the reactive patch cases in Table I are normalized by the diffusion-limited (no forces and reactivity confined to the 10-deg reactive patches) rate constant. This rate was evaluated by using the Brownian dynamics algorithm and trajectory analysis method for our model enzyme with no embedded charges and was found to be Kd = 1.15 X 1013 A3/s. The uniform reactivity rates of association listed in Table I1 are normalized by the Smoluchowski result2 for a uniformly reactive sphere with no forces KO = 4aD,, a, (where a is the radius of the spherical enzyme). As illustrated in Table I for the reactive patch case, the inclusion of a dielectric discontinuity (DD) enhances the rate 60% for the M Q model with overall charge of -4 and zero, and 45% for the MQ model with overall charge of +4, over that of the CD model. The dramatic enhancement in the more realistic DD model may be interpreted by examining the potential energy surfaces for the DD and CD models with a MQ and overall charge of -4. These are displayed in Figures 2 and 3, respectively. In Figure 2, a and b, we plot a contour map and the corresponding energy surface for the DD M Q -4 case: a cut of the potential through the Y-Z plane, (Oy,z), is taken. In Figure 3, a and b, analogous plots for the CD MQ -4 model are displayed. Both electrostatic surfaces show similar characteristics, with repulsive lobes perpendicular to the reactive patch regions, which define a channel into the active site; each channel exhibits a small energy barrier and finally a region of attractive forces near the surface of the protein. However, in contrasting these two models, a number of features are outstanding. First, we note the extended range of the repulsive energy contours away from the reactive patch in the DD case (Figure 2a,b). In addition to their spatial extent, these positive contours are energetically steeper for the DD model, where the maximum contour has a value of 2.04kT compared to the C D model, for which the corresponding contour is 1.64kT (Figures 2b and 3b). This suggests that the equilibrium concentration of 02-perpendicular to the reactive patch will be less for the DD case than for the C D case, for distances less than 100 A from the protein center. In addition, the attractive channel into the reactive site is deeper for the DD model (the minimum energy contours for the DD and CD models are -0.8kT and -0.1 kT, respectively) and is also longer ranged. Finally, the barrier leading into the attractive region is lower for the DD case (O.1kT compared with 0.2kT for the C D case). In fact, the C D model shows a small repulsive barrier surrounding the entire protein, a reason why this

The Journal of Physical Chemistry, Vol. 91, No. 12, 1987 3345

Binding of Small Ligands to Enzymes

I '

I

I

I

I

t

loo 50

N

0

-50

t-

-100

-100

-50

0

50

100

Y

Figure 2. Potential energy for the SOD dielectric discontinuity (DD) model. (a, top) Potential energy contours for the case of the monopole/quadrupole with overall charge -4e are plotted vs. y and z, in units of angstroms. Each contour energy is normalized by k T and each contour level represents 0.2kT. (b, bottom) The potential energy surface for the SOD/superoxide model described in (a).

100

50

N

o

-50

-100 -100

-50

0

50

100

Y

Figure 3. Potential energy for the SOD continuous dielectric (CD) model. (a, top) Potential energy contours for the case of monopole/ quadrupole of overall charge of -4. All other details are given in the caption of Figure 2a. (b, bottom) The same as Figure 2b except for the CD model.

model does not attain, let alone surpass, the diffusion-controlled result (Figure 3a). This repulsive barrier is absent in the DD model and offers further explanation as to why the DD rate surpasses the diffusion-controlled (no forces) rate. Thus, the influence of a differing dielectric constant between the interior of the protein and the surrounding solvent is to enhance the rate of ligand association as a result of more directed channeling of the ligand into the active site region (also, see below) and a larger effective target area, which is ultimately due to the greater extent of the attractive channel in the DD m ~ d e l . ' ~A*similar ~~ analysis is valid for the other MQ cases, since the same features of the

TABLE III: Rate of Ligand Association for the DD-DH and CD-DH 10-deg Reactive Patch SOD/Superoxide Model Systems at Finite Ionic Strength" charge modelb

CD-DH(K/Kd)

DD-DH(K/Kd)

0.000 0.000001 0.000005 0.000075 0.00001 0.00005 0.00075 0.000 1 0.0005 0.00 1 0.005 0.0 1 0.05 0.1 0.2 0.3 0.4 0.5

0.888

1.404 1.235 1.291 1.235 1.227 1.275 1.227 1.25 1 1.324 1.403 1.609 1.586 1.548 1.357 1.299 1.303 1.184 1.035

C

0.856 C

0.888 0.944 C

0.856 0.912 0.877 1.046 1.192 1.292 1.184 1.173 1.234 1.08 1 1.150

.

=The rates are normalized by the diffusion-limited rate (no forces with 10-deg reactive patches), Kd = 1.15 X 1013A3/s. bConcentration in units of moles/liter. 'Indicates no calculations were made at this salt concentration.

potential energy surface that determine the overall rate enhancement for the -4 charge model appear in the other two charge cases as well. The results for uniform reactivity presented in Table I1 are essentially the same for both the CD and DD models, within statistical error. This is attributed to a cancellation effect, in which favorable energy regions near the reactive site are balanced by unfavorable energy regions perpendicular to the protein reaction center. Thus the enhanced geometric restrictions provided by the extended repulsive electrostatic field in the DD model makes the probability of reaction perpendicular to the active site region small, despite the greater intrinsic probability of reaction provided at the surface of the molecule. For the CD model, after the ligand has surmounted the small barrier surrounding the protein, the repulsive forces away from the reactive sites are approximately 0.4kT lower than in the DD model. In the simplest argument, one would thus predict that the ratio of encounter rate away from the active site would be -K(CD)/K(DD) = exp(0.4), or a 60% enhancement of the rate of association in regions roughly perpendicular to the active site for the CD model. Therefore the CD model participates in a larger number of productive encounters with the protein surface away from the active site region, and these productive encounters counteract the smaller encounter probability at the active site, equalizing the overall rate between the DD and CD models. In Table I11 the Debye-Hiickel model results for varying concentrations of monovalent salt, such as NaC1, are presented. The precision in these rates range from &5% to &lo%. Experimental studies21of the SOD/superoxide system assayed the effect of increasing ionic strength on the activity of the native enzyme for both the trivalent salt sodium phosphate and the monovalent salt sodium chloride. We have focused on the results for the monovalent ions because of the simplicity of our ionic strength treatment (Deb~e-Hiickel)~~ and potential problems with multivalent salts.21 The rate for the native enzyme was found to decrease on addition of NaCl from 0.1 to 0.5 M. We have carried out calculations on the -4 MQ case with reactive patches for the DD and CD models, in which the former is most representative of the native enzyme, at concentrations of the monovalent salt between and 0.5 M. We observe a trend in the rate of association which is consistent with experimental studies in NaCl solutions at concentrations greater than 0.1 M when we employ the DD-DH model. (24) Lewis, G.; Randall, M.; Thermodynamics, 2nd ed.; McGraw-Hill: New York, 1956; see Chapter 23.

3346 The Journal of Physical Chemistry, Vol. 91, No. 12. 1987 0.50

0.25

Head-Gordon and Brooks

I -/-

3 \ M

N

v

3

0.00

100.

1

50

.

I

10-~ 1 Salt (M) Figure 4. log-log plot of the normalized rate of ligand association vs. ionic strength. The curved denoted by A represents results for the DDDH model, and the curve denoted by 0 is that for the CD-DH model. The filled symbols, A and are the zero ionic strength results from Table I.

-loo

The inclusion of extended Debye-Huckel theory into the formulation of the potential was investigated in order to explain the puzzling fact of decreasing rate with increasing ionic strength. The observed behavior emphasizes the importance of positive charges near the active site of SOD. From eq A9 in Appendix A, it is clear that the potential is dominated by the Debye-Hiickel screening exponential, which screens out the influence of the electric moments at distances large relative to the protein radius, a. As the ionic strength increases (note that K , the inverse of the Debye screening length, is proportional to the square root of the ionic strength), the screening is even more pronounced. Figure 4 presents the results for the rate of ligand binding vs. the concentration of monovalent salt in the DD-DH and CD-DH models; log-log plots of the normalized rate vs. ionic strength are displayed. As shown here for the DD-DH model, the rate is relatively flat at very low salt concentrations and then begins to rise at concentrations greater than 0.001 M. On further addition of salt, beginning at 0.005 M NaCI, the rate begins to fall, and continues to do so until a concentration of 0.5 M is reached, where the rate is virtually diffusion-controlled due to the effective screening of all forces. The flatness at very low salt concentrations is not surprising since the argument of the screening exponential is nearly zero, and consequently the potential is essentially that of the DD case. This is evident from the potential energy contour plot displayed in Figure 5 for the DD-DH model at 0.001 M salt concentration. By comparing with Figure 2a, one sees that the general features are the same between the two surfaces. For concentrations of salt greater than 0.001 M, the rate begins to increase because the monopole (the repulsive contribution) is screened more effectively than the quadrupole, with the result that the electrostatic barrier now disappears and is replaced by an extended attractive channel. On addition of salt greater than or equal to 0.1 M in concentration, the rate decreases below that of the no-salt case. The diminished rate is due to the full screening of the monopole, quadrupole, and all higher electric moments, with the result that the electrostatic forces are nearly zero, and the trajectory of the superoxide ligand is virtually diffusion-controlled. The screening of the higher moments reduces the extent of the repulsive contours, as seen in Figure 6 (DD-DH with 0.1 M salt), which leads to loss of the steering feature important to enhancement. We note that the potential energy contours radiating from the protein surface at this salt concentration are, at the lowest contour level, -0.2kT. The fact that the repulsive lobes may be important for steering the ligand into the active site is now clearly illustrated by noting three points of comparison between the DD-DH (0.1 M) salt results just discussed (Figure 6 and Table 11) and the earlier results for the DD no-salt case (Figure 2 and Table 11). First, the

.

, 1 -100

10-8

.,

/

0 .

-50

-0.25

1

-50

50

0

100

Y Figure 5. Potential energy contours for the SOD DD-DH model with 0.001 M salt. Potential energy contours for the case of monopole/ quadrupole with overall charge of -4e. Other details are given in Figure

2a. I\

100

-

I

I

1

\ \

-

/ /

\

N

I

/ \

50

,

/

\

\

-

/ \

cz24b \

-

/ /

-

0 e?'

-50

/

-

-

I

\ \

-

\

/

\

\

/

\

/

-100

1 1 -100

\

/

-

/ I

/

I

-50

I

I

0

50

L

I

100

Y Figure 6. Potential energy contours for the SOD DD-DH model with 0.1 M salt. Potential energy contours for the case of the monopole/ quadrupole with overall charge of -4e. Other details are given in Figure

2a. repulsive lobes are much less extended in the 0.1 M salt case. Second, the effective target area, the surface area defined by the outermost attractive contour,19is essentially the same in both cases. Finally, the 0.1 kT barrier blocking the entrance to the attractive channel in the DD case. is absent a t 0.1 M salt concentration. To reconcile the anticipated increase in rate due to the lack of any barrier at 0.1 M salt, with the observed decrease in the rate for this case, we conclude that the repulsive lobes must play a role in steering the ligand into the channel. The CD-DH model is similar in behavior to the DD-DH model at concentrations of salt less than 0.1 M (although the actual numbers differ by 60% between the two models). At salt concentrations greater than or equal to 0.005 M, the CD-DH model has screened out the repulsive barrier surrounding the protein in the no-salt case. This is evident when Figures 3a and 7 (CD-DH at 0.005 M salt) are compared. Between salt concentrations of 0.05 and 0.5 M the rate falls with increasing ionic strength, which is consistent with the experimental trend.21 However, at low salt concentration the rate of ligand association is lower than the diffusion-controlled rate which is nor consistent with experiment2' and points to a fatal flaw in the C D model^.^^^^ Another C D study on SOD/superoxide2s used a simple Debye-Hiickel screening potential of the (25) Allison, S.;McCammon, J. J . Phys. Chem. 1985, 89, 1072

The Journal of Physical Chemistry, Vol. 91, No. 12, 1987 3347

Binding of Small Ligands to Enzymes

loo

-100

t

i -100

-50

0

50

100

Y Figure 7. Potential energy contours for the SOD CD-DH model with 0.005 M salt. Potential energy contours for the case of the monopole/ quadrupole with overall charge of -4e. Other details are given in Figure 2a.

type exp(-KR?, where R’ is the distance from the embedded charge to the ligand. In the results of this work the rate initially increases with addition of salt due to the fact that the quadrupole is screened less than the monopole, and only at high salt concentrations, when K dominates, does the rate begin to decrease. However, the rate remains above the no-salt case. Thus, one is unable to correctly predict salt effects with a continuous dielectric model, and the use of a dielectric discontinuity model seems warranted for a more reliable study of the effect of ionic strength on the rate of association.

IV. Discussion and Conclusion Addition of the more realistic feature of a dielectric discontinuity to the superoxide dismutase/superoxide model fits another piece into the puzzle aimed at understanding the role of electrostatics in enzyme reactions. For the SOD/superoxide system, the change brought about by the dielectric discontinuity is to remove the repulsive barrier surrounding the protein and enhance the definition of the channel leading into the active site by extending the region of repulsive forces perpendicular to this site. In addition, these lobes are steeper energetically than the corresponding lobes in the CD model, making the approach of the ligand perpendicular to the reactive site less likely. The extended repulsive lobes are a result of a “lensing effect”, whereby the lines of force are focused due to a change of dielectric constant in going from the protein interior to the solvent, via the curved boundary at the protein surface. Once the ligand finds itself positioned in the vicinity of the reactive site, it is drawn into the site by a strongly attractive energy channel which is -0.8kT more attractive than the C D model due to the dielectric discontinuity. Again, the smaller dielectric constant in the interior of the protein provides less screening of the fixed charges, so that the positive charges near the active site exert a larger Coulomb attraction when the ligand is a t the mouth of the reactive channel. It is these synergistic features which are responsible for the rate enhancement for the M Q reactive patch case. The calculated results for uniform reactivity are the same for both models. The previous synergistic effects of the high-energy lobes with the deep energy channel into the reactive site now cancel one another for the M Q uniform reactivity case. The DD model still has a larger number of productive reactions in the region of the deep energy channel which leads to the surface of the protein. However, other regions of the surface will remain for the most part unreactive due to the increase of repulsive forces at the surface. The C D model, which has about 60% fewer successful collisions in the active site region than does the DD model, makes up this deficiency with successful encounters at the surface away from the channel. The final result is that these two effects cancel

to leave the overall reaction rate in the two models the same. The inclusion of extended Debye-Huckel theory into the formulation of the potential provided a means of exploring the puzzling observation that the rate of association decreases with increasing ionic strength, even though the overall charge of the protein is negative. It was found that small additions of salt resulted in an increased rate, due to the fact that the repulsive monopole is screened out more effectively than the attractive quadrupole, and thus the ligand moves in a more attractive electrostatic field. At higher salt concentrations, the quardupole is screened out, and the rate is reduced due to the diminishing of the steering feature, which was found to guide the ligand into the enzyme for a productive encounter at the active site at zero ionic strength. Only the DD-DH model was able to correctly account for the trend of diminished rate with increasing salt concentration. The CD-DH model is inherently unable to predict this trend due to the fact that the no-salt rate constant is less than the reactive patch diffusion-controlled rate; the latter rate is that which must be found for the limit of infinite salt addition. This model was similar in behavior to the CD model of Allison et aLZ5 which used a simpler screening function; their results showed an initial increase of rate with addition of small concentrations of salt, followed by a decrease as the ionic strength passed 0.05 M salt. Again, the rate did not dip below the no salt case. In the study of ionic strength effects on the rate of ligand association, the employment of dielectric discontinuity models over continuous dielectric appears to be important. Finally, we comment on the validity and utility of the simplified models we have employed. There are a great many advantages in maintaining a simplistic description of the enzyme in numerical studies of dynamical ligand association, and these advantages must be weighed against the possible misinterpretation of fundamental physical characteristics for the reacting system. In cases where the diffusional time scale is truly the dominant time scale in the process of interest, and the initial stages of encounter dynamics appear to be such a case, the framework of rigorous Brownian motion theory dictates that all atomic level detail is inherently absent in an explicit form and is included in an average way only. Therefore, the use of simple electrostatic models such as those employed here are wholly consistent with this description. The question then becomes, do we lose any of the basic phenomenology when such models are employed? The main conclusion reached in this study on SOD/superoxide is that the dielectric discontinuity enhances the steering of the ligand into the active site region. It is interesting to contrast the potential energy maps of SOD, generated by the simple electrostatic models of this study, with the potential energy maps generated numerically by Honig and co-workers.20 Their study evaluated several dielectric models of the full SOD enzyme, where each atomic position was placed on a grid and assigned a full or partial charge as required and the potential at each grid point was then evaluated from the linearized Poisson-Boltzmann equation. The potential energy maps in the present study compare quite favorably in general features to the full map within the same dielectric model. One sees that the phenomenology predicted by the simple models, Le., a favorable energy channel and repulsive lobes, is evident in the potential maps of the full enzyme. Honig et al. have argued that the extent of the repulsive lobes does not contribute to the rate enhancement; they contend that the enhancement is a result of the increased extent of the favorable energy channel for the DD model which provides a larger “target” with which the ligand can react. This conclusion is reasonable when one considers that the enzyme is elongated in shape along the direction perpendicular to the active site. Our simple model, which assumes a spherical shape, will probably overemphasize the repulsive lobes which run perpendicular to the reactive site and will underestimate the extent of the energy channel when compared to the elongated protein. However, it should be pointed out that some artifacts may also be present in the calculations of Honig and co-workers.20 In particular, the truncation of the potential at 95 A in their study will alter the extent of the high-energy lobes. This can be shown if we place a new boundary

3348 The Journal of Physical Chemistry, Vol. 91, No. 12, 1987

Head-Gordon and Brooks and that due to ions outside the protein surface. We take the origin to be at the center of the sphere and the charge to lie along the z axis. The potential may be written as

.(

\kin = m

-100

I -100

-50

0

50

100

Y Figure 8. Potential energy contours for the SOD DD model with an additional boundary at R = 100.0A. By comparing with Figure 2a one sees that the effect of the second boundary is to truncate the spatial extent of the repulsive lobes and enhance the attractive channel.20

in our model at R = 100 A and require the potential to vanish a t this surface. The potential energy surface with this new boundary is shown in Figure 8. By comparing the potential energy surface of Figure 2a with Figure 8, it is clear that the additional boundary condition manifests itself by truncating the extent of the high-energy lobes, while extending the range of the attractive energy channel and hence the “target” area. Thus the contention that there is no steering involved, and that the extent of the energy channel around the active site is solely responsible for the rate enhancement, may be incomplete. Instead, a more complete description of the association dynamics should include influences from the repulsive lobes as well as the effect of the channel barrier and the effective target area. In this scenario, the ligand first encounters the repulsive lobes away from the channel region and more rapidly diffuses into the channel due to the reduction in dimensionality (from three to something less than three) provided by motion along the surface of the lobes. Once positioned in the channel above the reactive patch, the barrier and effective target area (as defined by the extent of the attractive forces) play their anticipated roles in determining the overall rate.I9q2O Thus the inclusion of the dielectric discontinuity when evaluating the electrostatic forces, which leads to more extended repulsive lobes, a lower barrier, and a larger effective target size, provides a more realistic picture of enzymesubstrate encounter dynamics and offers further evidence that enzyme-substrate recognition begins in the early stages of diffusional encounter. It is hoped that the results and general conclusions reached in this present study will find application to systems other than SOD/superoxide. Acknowledgment. Acknowledgment is due to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for partial support of this work. We thank the Cottrell Research Corporation and the Samuel and Emma Winters Foundation for partial support of this work. Also, we acknowledge Professor B. Honig for providing ref 19 prior to publication. Appendix A Presented here is the derivation of the potential for an arbitrary distribution of point charges in a sphere of radius a with dielectric constant el, placed in a medium of dielectric constant e2. The effect of ionic strength is incorporated by using extended Debye-Huckel theory. Most of the derivation is due to Kirkwood,I6 although the present derivation employs a multisite angular expansion, which provides azimuthal symmetry around each charge, and focuses on the potential external to the protein. We will begin by deriving the DD-DH model, which yields the CD-DH, DD, and CD models when the appropriate limits are taken. The electrostatic potential inside the protein is a sum of contributions from the charge distribution inside the protein itself

+ BmRm)P,,,(cos

Rm+l

0)

R

< a (Al)

Azimuthal symmetry around each charge embedded in the sphere is broken with respect to the center of the sphere when considering the charge distribution as a whole. However, this form of the potential is still applicable to an arbitrary charge distribution by rotating the z axis so that it lies along the line joining the origin and the embedded charge. The forces are calculated by using ( A l ) and the Cartesian forces recovered by rotating back to the lab frame. A superposition of the forces for each charge results in the forces due to the charge distribution as a whole. Referring to Figure 1, R is the distance between the origin and the test charge, r is the distance between the origin and a charge embedded in the sphere, and 6 is the angle made by the vectors R and r. The coefficients B, will be determined from the boundary conditions as discussed below, and the functions Pm(cos 6) are the Legendre polynomials, which satisfy the following recursion relation

( m + l ) P , ( x ) = (2m

+ l)xP,(x)

- mP,-,(x)

where Pl(X) = 1.0 P2(x) = x

In addition, the following derivative recursion relation will be needed when evaluating the forces. dPm(x) dx

(2- 1) -- mxP,(x)

- mP,_,(x)

(A3)

The form of these potential functions is based on the model originally proposed by Kirkwood16 for zwitterions. To find the potential outside the protein surface we assume the extended theory of Debye-Hiickel is adequate. The charge distribution surrounding an ion in the solvent is assumed to be of the form of a Boltzmann distribution. In order to solve the Poisson equation with this form of the charge distribution, one must expand the exponential with respect to the mean force potential. If the salt concentration is dilute, higher order terms in the expansion may be ignored, and we are left with a linear equation in \k, where K is the inverse of the Debye screening length and is related to the salt concentration through the following relation

where e is the electronic charge, N is Avogadro’s number, I is the ionic strength in units of molarity, k is Boltzmann’s constant, T is the temperature, and t2 is the dielectric constant of the solvent. The general form of the potential is found to be A , eXp(-KR)K, (KR)P,( COS 8) R >a (A6) *,,ut = ,E Rm+1 where the K,(KR) are the modified spherical Bessel functions which obey the following recursion relationships.

+

- (2m 1 + t)K,(t) - (2m + l)Km+l(t)(A71 dt t The quantities A , and B, are unknown coefficients to be determined by imposing the appropriate boundary conditions at the interface of the two dielectric media. The appropriate conditions dK,(t) --

3349

J . Phys. Chem. 1987, 91, 3349-3354 to satisfy at the boundary are that the tangential component of the electric field and the normal component of the displacement join smoothly at the interface.

After solving for A , and B,, one arrives at an expansion for the potential for R > a qP(2m

*out

=

+ ~)K,(KR)P,(cos

c Rm+'[e2(2m + l)K,+,(Ka) m

potential precludes the penetration of ions into the protein interior as is evident by the arguments of the exponential term. This form of the Debye-Huckel treatment is different than that employed in the earlier work of McCammon et al.25 The potential found for the DD-DH model, (A9), can be used to find the potential for the other three models considered. If the dielectric inside the protein, el, is allowed to be that of the solvent, e2, the DD-DH model breaks down to the continuous dielectric model with inclusion of salt effects exterior to the protein, namely the CD-DH model.

0) exp(-K(R - a ) )

+ mKm(Ka)(e, - e2)]

('49) The expansion is appropriate for a single charge located on the z axis in the interior of the protein; the superposition principle allows one to evaluate (A9) for an arbitrary charge distribution embedded in the spherical protein as discussed above (note that if a charge is positioned at the origin, the potential is the first term (m = 0) of the sum in (A9)). When the quantity in (A9) is differentiated with respect to Cartesian coordinates, we arrive at the systematic electrostatic force required in the Brownian dynamics algorithm described in section 11. The series was found to converge for m 2 8, and the evaluation of the forces was truncated at m = 10. We also note that the

Again, the superposition principle allows one to assume azimuthal symmetry around each charge. If we now assume that e l and c2 are different, but let the Debye screening length go to infinity, we arrive at the DD model with no-salt effects *out

=

q P ( 2 m + l)Pm(cos 6) cRm+'[e2(2m + 1) + m(cl m

el)]

(A1 1)

Finally, it is assumed that el = t2, and ( A l l ) reduces to the continuous dielectric model of Allison et a].'

Simple Intramolecular Model Potentials for Water Liem X. Dang and B. Montgomery Pettitt* Department of Chemistry, University of Houston, Houston, Texas 77004 (Received: December 12, 1986)

An effective intramolecular potential is presented for use in conjunction with existing three-site models of water. Two commonly used internal geometries were fit to the same form yielding slightly different parametrizations. By including a Urey-Bradley-like term in an otherwise standard molecular mechanics form it was found that the experimental transition frequencies of water monomer can be reproduced accurately. Good qualitative agreements for spectral shifts were subsequently found for the models in condensed-phase applications. Harmonic analysis of clusters indicates good qualitative agreement with experimental environmental shifts in frequencies at low temperatures for these models. This model should be useful for a wide variety of applications including simulations of biopolymers and ionic solutions.

I. Introduction In recent years there has been an increased interest in the vibrational motions of water from the point of view of both ex~ r i m e n t l -and ~ theory.- In principle, given the potential energy surface as a function of the nuclear coordinates, the observable properties are determined from quantum mechanical and statistical mechanical methods. For the study of condensed-phase systems (1) Page, R. H.; Frey, J. G.;Shen, Y. R.; Lee, Y. T. Chem. Phys. Lett. 1984, 108, 1373. Vernon, M. F.; Lisy, J. M.; Krajnovich, D. J.; Tramer, A,; Kwok. H.-S.: Sen. Y. R.: Lee. Y. T. Faradav Discuss. Chem. SOC.1982. 73.

387. Vernon, M.'F.; Krajnovich, D. J.; Kwok, H A ; Lisy, J. M.; Shen, Y. R.; Lee, Y. T. J. Chem. Phys. 1982, 77,47. (2) Chen, S. H.; Toukan, K.; Loong, C. K.; Price, D. L.; Teixeira, J. Phys. Rev. Lerf. 1984, 53, 1360. (3) Walfren, G . E. J. Chem. Phys. 1967, 47, 116. Walfren, G. E. In Water: A Comprehensive Treatise; Franks, F., Ed.; Plenum: New York, 1972; Vol. 1; p 151. (4) Bentwood, R. M.; Barnes, A. J.; Oroille-Thomas, W. J. J. Mol. Spectrosc. 1980, 84, 391. (5) Coker, D. F.; Miller, R. E.; Watts, R. 0. J . Chem. Phys. 1985, 82, 3554. Reimers, J. R.; Watts, R. 0.Chem. Phys. 1984,91, 201 and references therein. (6) Toukan, K.; Rahman, A. Phys. Rev. E 1985, 31, 2643. (7) Postma, J. P. M. Ph.D. Thesis, Rijksuniversiteit te Groningen, 1985. (8) Berens, P. H.; MacKay, D. H.; White, G . M.; Wilson, K. R. J. Chem. Phys. 1983, 79, 2375. (9) Bopp, P.; Jansco, G.;Heinzinger, K. Chem. Phys. h i t . 1983,98, 129.

0022-3654/87/2091-3349$01.50/0

two techniques of choice have been the Monte Carlolo and molecular dynamics" methods. In view of the importance of water and aqueous solutions and the computational resources devoted to their study, it is significant to have effective classical models of water which give reasonable properties with a minimum computational effort. This is most imperative in large-molecule studies where the majority of the time is spent simulating the environment rather than the solute. There now exist a number of relatively simple intermolecular potential surfaces for water, none of which are ideal for all the properties determined by the intermolecular interactions but many of which reproduce selected properties. In this effort to construct a simple intramolecular potential we restrict our attention to the simplest intermolecular models of water, namely, the three-site models,'2J3 as these widely used models are most computationally efficient at providing an aqueous environment in molecular solvation studies. The intramolecular oscillatory motions of an (10) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. A,; Teller, E. J . Chem. Phys. 1951, 55, 4291. (11) Alder, B.; Wainwright, T. E. J. Chem. Phys. 1957, 27, 1208. (1 2) Berendsen, H. J. C.; Postma, J. P. M.; van Gunstren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (13) Jorgensen, W. L. J . Am. Chem. SOC.1981, 103, 335.

0 1987 American Chemical Society