J. Phys. Chem. 1996, 100, 5149-5154
5149
Time-Dependent Rate Coefficients from Brownian Dynamics Simulations Michael J. Potter* Department of Chemistry and Biochemistry and Department of Pharmacology, UniVersity of California at San Diego, La Jolla, California 92093-0365
Brock Luty† Physical Chemistry, ETH Zentrum, UniVersitatstrasse 6 CH 8092 Zurich, Switzerland
Huan-Xiang Zhou‡ Department of Biochemistry, The Hong Kong UniVersity of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
J. Andrew McCammon§ Department of Chemistry and Biochemistry and Department of Pharmacology, UniVersity of California at San Diego, La Jolla, California 92093-0365 ReceiVed: October 31, 1995; In Final Form: January 10, 1996X
The kinetics of a diffusion-influenced reaction is described by a time-dependent rate coefficient k(t). In this paper, we implement an algorithm for calculating k(t), introduced by one of us, in the UHBD software package. In the implementation, the electrostatic force exerted on a diffusing substrate by an enzyme is obtained from a finite-difference solution of the Poisson equation. The rate coefficient is obtained by starting the substrate in the active site and monitoring its survival probability using Brownian dynamics simulations. The technique is applied to the binding of superoxide to Cu/Zn superoxide dismutase. The long-time limit of k(t) is found to be in agreement with both the experimental value and that calculated using an algorithm designed specifically for finding k(∞).
Introduction The description of the kinetics of diffusion-influenced reactions in terms of a time-dependent rate coefficient, k(t), is a problem known to be analytically solvable only for cases of idealized geometry. There are, however, many important systems too complicated for such representations. Examples include ligand binding to biological macromolecules, enzymecatalyzed reactions, and the self assembly of multi-subunit proteins. It is thus desirable to develop a numerical approach for the study of such systems. Recently such an algorithm was introduced for the calculation of the time-dependent rate coefficient for systems of arbitrary geometry using Brownian dynamics.1 The use of Brownian dynamics simulations in the calculation of steady-state rate coefficients is a well-known technique.2-5 Its ability to deal with complex geometry and potentials, particularly when combined with a Poisson-Boltzmann treatment of electrostatics, has led to its application to large, diffusion-controlled, enzymatic systems, such as Cu/Zn superoxide dismutase (Cu/Zn SOD).6-12 Cu/Zn SOD is a metalloprotein that plays an important role in the control of oxygen toxicity through its conversion of the superoxide radical to less dangerous species according to the overall reaction
2O2- + 2H+ f H2O2 + O2
(1)
Because of this important role, there has been much interest in * Corresponding author. E-mail:
[email protected]. † E-mail:
[email protected]. ‡ E-mail:
[email protected]. § E-mail:
[email protected]. X Abstract published in AdVance ACS Abstracts, March 1, 1996.
0022-3654/96/20100-5149$12.00/0
this class of enzymes.13,14 Cu/Zn SOD’s steady-state kinetics have been successfully modeled using Brownian dynamics simulation, but no attempt to model non-steady-state effects has been reported. In this paper we report the implementation and testing of the algorithm of Zhou1 for the calculation of k(t), and its application to the Cu/Zn SOD system. Theory The algorithm has been described in detail previously.1,15,16 Here we include a summary for the paper to be self-contained. We are interested in calculating the time-dependent rate coefficient describing the kinetics of the reaction between two particles, A and B. We define the pair distribution function p(r,t) as the concentration of B at the position r, normalized by its bulk concentration. It is assumed that p(r,t) satisfies the Smoluchowski equation17,18
∂p(r,t)/∂t ) -∇‚J(r,t) ≡ ∇‚De-βU(r)∇eβU(r)p(r,t)
(2)
where β ) (kBT)-1, D is the relative diffusion coefficient of the species (i.e., D ) DA + DB), and U(r) is the potential of mean force between A and B. Until t ) 0, when the reaction is initiated, it is assumed that B is present in a Boltzmann distribution about A
p(r,0) ) e-βU(r)
(3)
After t ) 0, reaction between the particles is described by a radiation boundary condition19
-n‚J(r,t) ≡ n‚De-βU(r)∇eβU(r)p(r,t) ) κ(r)p(r,t), r ∈ Γ (4) © 1996 American Chemical Society
5150 J. Phys. Chem., Vol. 100, No. 12, 1996
Potter et al.
Here Γ is the surface around A at which B reacts, n is the normal of Γ, and κ(r) is the space-dependent intrinsic bimolecular rate constant. The outer boundary condition is
p(∞,t) ) 1
(5)
Our goal is the calculation of the time-dependent rate coefficient defined by the equation
k(t) ) -∫ΓdA n‚J(r,t)
(6)
The most direct method for obtaining k(t) would be to solve eq 2, subject to the given initial and boundary conditions, and then integrate to calculate k(t). This is very difficult for cases in which the geometry or potential is complicated. Brownian dynamics simulations are capable of treating such systems, and the relationship between the pair distribution function p(r,t) and the survival probability at time t of a particle started at r20,21 -βU(r)
p(r,t) ) e
S(t|r)
(8)
Thus, k(t) can be found by starting particles on the boundary Γ with a distribution of κ(r)e-βU(r) and recording the fraction of surviving particles at time t. Away from any boundary, if the diffusion coefficient and force are constant over a small time step ∆t, the motion of a particle can be simulated using the Brownian dynamics algorithm22
r ) r0 + D∆tβF(r0) + (2D∆t)1/2R
(9)
where R is a vector of Gaussian random numbers with the properties 〈Ri〉 ) 0, 〈RiRj〉 ) δij for each of its three components. However, as a particle nears a boundary, care must be taken to ensure that boundary conditions are satisified. This is accomplished by converting the radiation boundary condition to a reflecting boundary condition and adding a sink term to the equation governing the conditional probability.23-26 The reflecting boundary condition is modeled by making a new attempt whenever a particle crosses Γ. Reactions are allowed within the volume Λ between Γ and a secondary surface Γ′ (see Figure 1). At each step in which the particle is within Λ, a random number, uniformly distributed between 0 and 1, is generated. If this number is greater than or equal to exp{-[∆tκ(r0)/ + ∆tκ(r)/]/2}, the particle is considered to have reacted and the trajectory is terminated. Here is the distance separating the Γ and Γ′ surfaces. The modifications given above lead to a slightly different expression for k(t)
kF(t) ) ∫ΛdV
κ(r) -βU(r) SF(t|r) e
approximated as
(7)
forms the basis of the algorithm used here. The survival probability is easily obtained from a Brownian dynamics simulation by calculating the fraction of particles started at r that have not reacted by time t. We can write k(t), as defined by eq 6, in terms of S(t|r) using eqs 4 and 7:
k(t) ) ∫ΓdA κ(r)e-βU(r)S(t|r)
Figure 1. Illustration of the boundaries Γ and Γ′, their separation , and the reactive volume Λ.
(10)
Here the subscript F emphasizes that a different conditional probability is being used. Note that the surface integration in eq 8 has been replaced by an integral over the reactive volume Λ. It can be shown that as f 0, kF(t) f k(t). In practice, κ/D , 1 yields adequate convergence. Equation 10 may be
κ kF(t) ) V e-βUh 〈SF(t|r)〉
(11)
Here V is the volume of the reactive region Λ, U h is the average of the potential of mean force between the particles over Λ, and 〈SF(t|r)〉 is the average survival probability of particles started with distribution κ(r) exp[-βU(r)] within the reactive volume. To summarize, particles are started within the reactive volume Λ according to the distribution κ(r) exp[-βU(r)]. The trajectory is propagated using eq 9. If the particle crosses a reflective barrier, a new attempt is made. Trajectories are terminated if the reaction conditions are satisfied or if a preset lifetime cutoff is exceeded. From the particle lifetimes, 〈SF(t|r)〉 is calculated and eq 11 is used to obtain k(t). Technical Details The k(t) algorithm was implemented within the UHBD software package.27 The implementation was tested through its application to cases for which analytic or near analytic solutions are known. Specifically, the reaction of two uniformly reactive spheres with radiation boundary conditions, both with and without Coulomb forces, was examined. The radii of the spheres were chosen to be 28.55 and 2.05 Å to approximate the characteristics of SOD and superoxide. This gives a relative diffusion constant of 0.1290 Å2/ps, assuming no hydrodynamic interactions between the spheres and stick boundary conditions. Testing was conducted under conditions of high and low reactivity, κ ) 0.042 and 0.0042 Å/ps. High-reactivity cases are more challenging numerically due to the small size of and the small time steps that are thus required. For the simulations with a detailed model of SOD, the coordinates of the protein were taken from the 2 Å crystal structure of Tainer et al.28 The structures of the subunits of the SOD homodimer are slightly different in the crystal structure, particularly in the area of the active sites. In order to facilitate comparisons of the rates of reaction at each active site (to judge the statistical quality of the results), a “symmetric” dimer composed of identical monomers was constructed. This was accomplished by overlaying a copy of one monomer on the other, such that the root mean square deviation between corresponding atoms was at a minimum. The original monomer was then deleted, leaving a dimer composed of structurally identical monomers in the correct relative orientation. Polar hydrogens were added to the structure using the molecular modeling package CHARMm 22.0.29 Titratable sites were
Rate Coefficients from Brownian Dynamics Simulations
J. Phys. Chem., Vol. 100, No. 12, 1996 5151
considered to be in the same protonation state as their equivalent amino acids at pH 7.0. The generation of correctly distributed starting points within the reactive volume Λ was accomplished analytically in the test cases. For SOD, a 65 × 65 × 65 grid of the finest spacing possible (∼0.1 Å), such that the entire reactive volume was still included, was overlaid on Λ. Grid points within Λ but not within the excluded volume of the protein were saved. The excluded volume of the protein was defined as the area within the surface generated by the sum of the Van der Waals radii of the atoms and the radius of the diffusing substrate. A Monte Carlo procedure was then used to select a set of points distributed according to κ(r) exp[-βU(r)]. It was necessary to handle the test cases analytically, as the large size of the reactive volume forced the use of a very coarse grid. This led to a biased distribution of points. This was not a problem for SOD, as the small size of the reactive volume allowed a very fine grid to be used. Electrostatic forces were calculated using finite difference solutions of the Poisson equation (i.e., with zero ionic strength). Details of the procedure have been described elsewhere.30 Briefly, the test sphere, or detailed model of SOD, was placed at the center of a 110 × 110 × 110 grid with a spacing of 1.0 Å. Uniform dielectrics of 40 and 78 were used for the tests with Coulombic interactions. For SOD, the interior of the protein was set to a dielectric of 2 and the solvent to 78. Charges of -4.0 and -1.0 were assigned to the target and diffusing particle, respectively, for the test cases. Charges and radii for SOD were taken from the GROMOS parameter set31 except for the metals and their ligands. Charges for these atoms were assigned on the basis of earlier quantum mechanical calculations.32 The electrostatic potential at the boundary of the grid and beyond was approximated using a dielectrically screened ( ) 78) Coulomb potential for an ion of radius 30 Å. Trajectories were propagated using the Brownian dynamics algorithm (eq 9). To increase the efficiency of the algorithm, the time step was increased as the separation between particles increased:
{
10-22 reb 2D ∆t ) (r - b)2 10-2 2 r>b +c 2D D
(12)
This allowed large time steps to be taken when the diffusing particle was far away from the target. In all cases presented here, b ) 41.5 Å and c ) 0.01. Variation of either parameter was determined to have little effect on the results. Trajectories were terminated either when the reaction criteria given above were satisfied or if they exceeded a preset lifetime cutoff. The lifetime cutoffs were 100 ns for the low-reactivity test cases, 1 ns for the high-reactivity test cases, and 50 ns for SOD. The recorded lifetimes for the trajectories were then used to calculate average survival probabilities. The time-dependent rate constant, k(t), was then calculated using eq 11. It was, however, necessary to modify eq 11 for use with the test cases. This was the result of the geometry and large size of the reaction volume. As was mentioned above, kF(t) f k(t) as f 0 and κ/D , 1 generally yields adequate convergence. This was not the case for the test cases. The survival probability was converged but the V(κ/) portion of eq 11 was not. This lack of convergence was due to the spherical shell nature of the reactive volume and the large radius of the target particle. This difficulty can be circumvented by modifying eq 11. Equation 8 can be written in the following way:1
k(t) ) k(0)〈S(t|r)〉
(13)
k(0) ) ∫ΓdA κ(r)e-βU(r)
(14)
∫ΓdA κ(r)e-βU(r)S(t|r) 〈S(t|r)〉 ) ∫ΓdA κ(r)e-βU(r)
(15)
where
This suggests replacing V(κ/)e-βUh with k(0) for cases in which the integral can be easily evaluated. This was straightfoward for the test cases, so eq 13 was used to calculate k(t). It was not feasible to do the same for SOD, so eq 11 was used. It is not likely that this problem of convergence occurs with the much smaller, nonspherical reactive volume of SOD in any case. For the SOD simulations, it was necessary to define the reaction volume and choose a value for the intrinsic reactivity, κ. The reaction volume was defined as the volume enclosed by 7 Å spheres about each active site Cu and outside the excluded volume of the atoms lining the active site. This roughly corresponds to a minimum in the O2- potential of mean force, as determined by earlier MD studies, from which the reaction is thought to take place.33 With regard to the value of κ, the algorithm itself places a practical upper bound on κ once the value of has been defined. In the case of SOD, by defining the reactive volume as described above, the value of has been set to about 2 Å. This is approximately the largest distance between the boundary of a 7 Å sphere centered on Cu and the Van der Waals surface of the protein. In order to be assured of convergence, κ can be no larger than 0.1D/ ) 0.0032 Å/ps. Given the nearly diffusion-controlled nature of SOD, it is reasonable to use the largest value of κ allowed by the algorithm and our choice of reaction volume. Values of 0.0022 and 0.0012 Å/ps were also considered for κ. The steady-state rate constant obtained from the above simulations can be compared with the results obtained in previous studies of SOD.9,10 To make a fair comparison, the procedure in these studies was used to find k(∞) under the conditions described above. The electrostatic potential was treated in the same way as for the k(t) calculations. Trajectories were randomly initiated on a sphere of radius b ) 41.5 Å centered on SOD. Termination occurred either when the particle diffused to within 6 Å of an active site Cu or when it reached a distance of q ) 300 Å. A total of 15 000 trajectories were used to calculate β, the probability that a particle started on the b surface will react rather than escape by crossing the q surface. Given β, the bimolecular rate constant can be obtained from the following equations2
k(∞) )
(
kD(b)β
)
(16)
kD(b)
1 - (1 - β) kD(q)
kD(a) ) 4πD(∫a r-2eβU(r) dr)-1 ∞
(17)
Results and Discussion For spheres with radiation boundary conditions, in the absence of forces, it is known that18
{
}
κa D 1 + Ψ[(1 + κa/D)(Dt)1/2/a] k(t) ) k(0) D + κa D k(0) ) 4πκa2
(18) (19)
5152 J. Phys. Chem., Vol. 100, No. 12, 1996
Potter et al.
Figure 2. Results for uniformly reactive spheres with radiation boundary conditions and without Coulomb forces. The dashed line was obtained from simulation (κ ) 0.0042 Å/ps). The solid line is the analytical result calculated from eq 18.
where Ψ(x) ) exp(x2) erfc(x) and a is the radius of the reaction surface. To calculate k(t) for the test cases by Brownian dynamics, 40 000 trajectories were run. This required approximately 2.2 h on a SGI 200 MHz R4400 workstation. In Figure 2 the calculated k(t) for low reactivity (κ ) 0.0042 Å/ps) is compared to the analytical solution. It can been seen that the simulated and analytical values for k(t) match very well throughout the entire range of the simulation. The results were similar for the high-reactivity case (not shown). There is no exact solution for k(t) for spheres with radiation boundary conditions that are subject to Coulomb forces. There is, however, an approximate solution that is accurate at long times18,24
aeffD k(t) ) k(0) 2 -r /a[1 + (πDt)-1/2aeff] κa e c
(20)
aeff ) rc[(1 + Drc/κa2)erc/a - 1]-1
(21)
k(0) ) 4πκa2e-rc/a
(22)
Figure 3. Results for uniformly reactive spheres with radiation boundary conditions including Coulomb forces. The dashed lines were obtained by simulation (κ ) 0.0042 Å/ps). The solid lines are the near analytical results calculated from eq 20. Note that the labels on the graph refer to the dielectric used not to the separation between Γ and Γ′.
Here
(z1z2e2/4π0kBT).
Figure 3 and rc is the Onsager distance compares the results of simulations using dielectric constants of 40 and 78 with the curves calculated from eq 20. For each case, approximately 40 000 trajectories were run. A dielectric constant of 40 was included in the testing to provide a more stringent test of the Poisson representation of the potential. As can be seen from the graph, close agreement between the simulations and eq 20 was obtained throughout the time range considered. Figure 4 shows the k(t) curve for SOD and its long-time limit. A total of 40 000 trajectories were run. Also shown is the rate constant obtained from steady-state Brownian dynamics simulations with κ ) ∞. The long-time limit of k(t) was obtained by fitting the curve to a function of the form1
k(t) ) a +
b t
1/2
(23)
Figure 4. Simulation results for Cu/Zn SOD. The “perfect absorber” rate constant obtained from steady-state Brownian dynamics simulations and the long-time limit of the k(t) curve are given for comparison.
The long-time limit of k(t) is thus a. In this case, the curve was fitted over the range 20-50 ns to yield a value of 2.97 Å3/ps. The fit was very good, with a standard error of 8.4 × 10-4. This long-time limit corresponds to a steady-state rate constant of 1.79 × 109 M-1 s-1, which is in good agreement with the reported experimental rate of ∼4 × 109 M-1 s-1.34-36 The experimental result is, however, for an ionic strength greater than zero. As can be seen from the graph, k(t) is very large during the initial stages of the reaction, as substrate molecules very near or within the active site react quickly. However, once this supply has been depleted, convergence to steady-state occurs quite rapidly and is essentially complete by 20 ns. Despite the speed of convergence to steady state, it is clear that there are significant deviatons from steady-state behavior at short times. Figure 5 illustrates the effect this short-time behavior has on
Rate Coefficients from Brownian Dynamics Simulations
J. Phys. Chem., Vol. 100, No. 12, 1996 5153
Figure 5. Ratio of product formed for the non-steady-state reaction to that formed if steady-state is assumed to hold from t ) 0. This was calculated using eq 24.
Figure 7. Plot illustrating the linearity of the dependence of 1/k(∞) on 1/κ.
reactivity in various cases where results are known analytically, we propose the following relation
1 B )A+ κ k(∞)
Figure 6. Cu/Zn SOD simulation results for a variety of intrinsic reactivities.
(25)
If this relation holds, a graph of 1/k(∞) vs 1/κ would be linear and the intercept A should give the inverse of the steady-state rate constant for infinite κ (the “perfectly” absorbing case). Figure 7 illustrates that 1/k(∞), obtained from the long-time limits of the k(t) curves, is indeed linearly dependent on 1/κ. From the intercept in Figure 7, k(∞) in the perfectly absorbing case was found to be 5.16 Å3/ps. In comparison, the result obtained by directly simulating k(∞) in this case using eq 16 was 7.68 Å3/ps. If the steady-state rate constants for κ ) 0.0022 and 0.0012 Å/ps were obtained by fitting k(t) at longer times (as opposed to the 5-10 ns range), the resulting k(∞) for κ ) ∞ would likely be even closer to the value of 7.68 Å3/ps. In summary, the technique developed in this study yields both the full time dependence of the rate coefficient for arbitrary intrinsic reactivity and the steady-state rate constant for the perfectly absorbing case. Conclusion
product formation. The curve is the ratio of product formed when the time-dependent rate constant is used to that formed if one assumes steady state from the beginning of the reaction. The curve was obtained using the following formula
ratio )
∫0tk(t) dt k(∞)t
(24)
It can be seen that the large size of k(t) early in the reaction results in significant deviations from steady state with respect to product formation. To study the dependence of the rate coefficient k(t) on the intrinsic reactivity κ, two other values of κ, 0.0022 and 0.0012 Å/ps, were also considered. Figure 6 presents the results of these calculations. The steady-state rate constants, obtained by fitting the simulated k(t), in the 5-10 ns range, to eq 23, were found to be 2.68 and 1.85 Å3/ps, respectively. By examination of the dependence of the steady-state rate constant on intrinsic
The Brownian dynamics algorithm for the calculation of the time-dependent rate coefficent, k(t), was implemented and tested in a computer program that allows for detailed descriptions of biomolecular structures and interactions. The method accurately reproduces analytical and near analytical results for simple test cases. Most importantly, the algorithm is not limited by the need for idealized geometries or simplified interaction potentials and thus is applicable to a wide range of important biological systems. This is illustrated by the successful application of the method to the Cu/Zn SOD system. The long-time limit of k(t) is in good agreement with the experimentally measured steadystate rate constant. The steady-state rate constant for a “perfect” absorber, calculated from the k(t) data for three values of κ, also agrees well with the value calculated using an established steady-state Brownian dynamics procedure. The steady-state constant for κ ) ∞ is approximately twice that calculated for κ ) 0.0032 Å/ps. An earlier study, which coupled Brownian dynamics to molecular dynamics in order to
5154 J. Phys. Chem., Vol. 100, No. 12, 1996 allow for an explicit treatment of interactions within the active site, found the same relationship between the “perfectly absorbing” rate constant and one that included the influence of events within the active site.37 In our case these events are described by the intrinsic reactivity, κ. This agreement is encouraging and supports the assumptions made in the definition of κ for SOD. The success of the method as applied to Cu/Zn SOD, combined with its flexibility with regard to electrostatics, opens the possibility of extending theoretical studies of the effects of enzymatic electrostatic potentials on their reaction rates beyond previous steady-state studies. For example, it has been suggested that the steep nature of the electrostatic potential gradient that results from the presence of anionic residues near the active site might act to speed convergence to steady state and thus play an important role in the function of SOD in vivo, where superoxide concentrations can fluctuate very rapidly.11,14 Acknowledgment. This work has been supported in part by the National Institutes of Health and the Metacenter Program of the National Science Foundation’s supercomputer centers. M.J.P. is a Howard Hughes Medical Institute Predoctoral Fellow. References and Notes (1) Zhou, H.-X. J. Phys. Chem. 1990, 94, 8794. (2) Northrup, S. H.; Allison, S. A.; McCammon, J. A. J. Chem. Phys. 1984, 80, 1517. (3) McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, 1987. (4) Zhou, H.-X. J. Chem. Phys. 1990, 92, 3092. (5) McCammon, J. A. Curr. Biol. 1992, 2, 585. (6) Allison, S. A.; McCammon, J. A. J. Phys. Chem. 1985, 89, 1072. (7) Allison, S. A.; Ganti, G.; McCammon, J. A. Biopolymers 1985, 24, 1323. (8) Sharp, K.; Fine, R.; Honig, B. Science 1987, 236, 1460. (9) Allison, S. A.; Bacquet, R. J.; McCammon, J. A. Biopolymers 1988, 27, 251.
Potter et al. (10) Sines, J. J.; Allison, S. A.; McCammon, J. A. Biochemistry 1990, 29, 9403. (11) Getzoff, E. D.; Cabelli, D. E.; Fisher, C. L.; Parge, H. E. Nature 1992, 358, 347. (12) Sergi, A.; Ferrario, M.; Polticelli, F.; O’Neill, P.; Desideri, A. J. Phys. Chem. 1994, 98, 10554. (13) McCord, J. M.; Fridovich, I. Free Radical Biol. Med. 1988, 5, 363. (14) Bannister, J. V.; Bannister, W. H.; Rotilio, G. CRC Crit. ReV. Biochem. 1988, 22, 111. (15) Zhou, H.-X. Biophys. J. 1993, 64, 1711. (16) Zhou, H.-X.; Szabo, A. J. Phys. Chem., in press. (17) Smoluchowski, M. Z. Phys. Chem. 1917, 92, 129. (18) Rice, S. A. In ComprehensiVe Chemical Kinetics; Bamford, C. H., Tipper, C. F. H., Compton, R. G., Eds.; Elsevier: Amsterdam, 1985; Vol. 25. (19) Collins, F. C.; Kimball, G. E. J. Colloid Sci. 1949, 4, 425. (20) Onsager, L. Phys. ReV. 1938, 54, 554. (21) Tachiya, M. J. Chem. Phys. 1978, 69, 2375. (22) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352. (23) Szabo, A. J. Phys. Chem. 1989, 93, 6929. (24) Zhou, H.-X.; Szabo, A. J. Chem. Phys. 1990, 92, 3874. (25) Wilemski, G.; Fixman, M. J. Chem. Phys. 1973, 58, 4009. (26) Northrup, S. H.; Hynes, J. T. J. Stat. Phys. 1978, 18, 91. (27) Davis, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A. Comput. Phys. Commun. 1991, 62, 187. (28) Tainer, J. A.; Getzoff, E. D.; Beem, K. M.; Richardson, J. S.; Richardson, D. C. J. Mol. Biol. 1982, 160, 187. (29) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminatham, S.; Karplus, M. J. Comput. Chem. 1982, 4, 187. (30) Madura, J. D.; Davis, M. E.; Gilson, M. K.; Wade, R. C.; Luty, B. A.; McCammon, J. A. ReV. Comput. Chem. 1994, 5, 229. (31) van Gunsteren, W. F. GROMOS Molecular Dynamics Package. (32) Shen, J.; Wong, C. F.; Subramaniam, S.; Albright, T. A.; McCammon, J. A. J. Comput. Chem. 1990, 11, 346. (33) Shen, J.; McCammon, J. A. Chem. Phys. 1991, 15, 191. (34) Bray, R. C.; Cockle, S. A.; Fielden, E. M.; Roberts, P. B.; Rotilio, G.; Calabrese, L. Biochem. J. 1974, 139, 43. (35) Fielden, E. M.; Roberts, P. B.; Bray, R. C.; Lowe, D. J.; Mautner, G. N.; Rotilio, G.; Calabrese, L. Biochem. J. 1974, 139, 49. (36) McAdam, M. E. Biochem. J. 1977, 161, 697. (37) Luty, B. A.; Elamrani, S.; McCammon, J. A. J. Am. Chem. Soc. 1993, 115, 11874-11877.
JP953229N
