3886
J. Phys. Chem. B 1997, 101, 3886-3890
On the Presence of Rotational Ewald Artifacts in the Equilibrium and Dynamical Properties of a Zwitterionic Tetrapeptide in Aqueous Solution Paul E. Smith* Department of Biochemistry, Willard Hall, Kansas State UniVersity, Manhattan, Kansas 66506-3702
Herb D. Blatt and B. Montgomery Pettitt Department of Chemistry, UniVersity of Houston, 4800 Calhoun Road, Houston, Texas 77204-5641 ReceiVed: NoVember 12, 1996X
The presence of Ewald artifacts in a molecular dynamics simulation of a zwitterionic tetrapeptide in aqueous solution has been investigated. Both equilibrium and dynamical aspects of the rotational behavior of the peptide were examined. The equilibrium distribution of rotational states obtained from a 10 ns simulation were consistent, within statistical errors, with a random distribution free from Ewald artifacts. The rotational dynamics of the peptide was observed to obey simple Debye diffusion with a rotational diffusion rate in agreement with that predicted by a simple rotational diffusion model assuming no substantial long-lived water hydration shell. These results suggest that rotational Ewald artifacts will be negligible for small biomolecular simulations in solvents with high relative permittivities. The data support the results obtained from a previous study of simpler model systems.
I. Introduction The Ewald potential has fast become an accepted treatment of electrostatic interactions in biomolecular systems.1-3 This is largely due to the presence of faster implementations of this traditionally computationally intensive technique.4-6 However, the possibility of artificial periodicity effects in these systems is still the subject of much discussion.7-12 Artificial periodicity effects are not confined to the Ewald potential. The use of periodic boundary conditions per se produces subtle system size dependent orientational effects even in nonpolar systems such as liquid argon.13 Fortunately, these effects are less important for biological systems because of the larger numbers of solvent molecules involved. Ewald artifacts are also system size dependent. However, the exact system size dependence of these effects in the context of biomolecular simulations is still unclear. Artifacts are not expected to be significant for simple dilute solutions because of the free mobility of the molecules. In contrast, peptides and protein simulations involve systems containing very polar, often charged, molecules in which large permanent electric moments exist. Here, many-body correlation effects might be expected to increase the possibility of artifacts when using the Ewald potential. Recently, Smith and Pettitt investigated the dependence of the Ewald potential on the orientation of a simple dipole, a linear quadrupole, and a rigid protein within a central cubic cell in order to obtain more information concerning the presence of Ewald artifacts in biomolecular systems.12 An artifact of the Ewald potential is that a dipole prefers to orientate along one of the axes of the lattice (θ ) 0°, θ/φ ) 90°/0° or 90°/90°, in spherical polar coordinates), with a maximum in the rotational potential surface at θ/φ ) 54.74°/45° corresponding to the dipole pointing toward the corner of the cell. The situation is reversed for a linear quadrupole. To a good approximation, the Ewald rotational potential energy surface scales with q2r4L-5-1, where q is the charge, r is the charge separation, L is the box length, and the relative permittivity of the surrounding system.12 X
Abstract published in AdVance ACS Abstracts, April 15, 1997.
S1089-5647(96)03764-9 CCC: $14.00
For the simple model used in the previous work the value of was assumed to be that of the solvent.12 Since the relative permittivity of water is high, it was argued that Ewald artifacts on the rotational properties of most molecules in aqueous solution would be negligible. However, since the dipoles in the Ewald lattice are, by definition, perfectly correlated and no molecular solvent effects were considered, a more rigorous test under typical simulation conditions is required. To determine the magnitude of Ewald artifacts in a realistic simulation, we have analyzed the molecular dynamics trajectory of a zwitterionic alanine tetrapeptide in an explicit solvent performed using the Ewald potential.14 The peptide may be considered as a flexible dipole. Rotational Ewald artifacts are expected to produce at least two effects. First, the equilibrium distribution describing the orientation of the peptide dipole in the simulation cell will exhibit features corresponding to the underlying Ewald rotational potential energy surface. Second, a reduced rotational diffusion rate would be expected due to the presence of underlying potential energy barriers to rotation. Therefore, the rotational dynamics of the peptide would not correspond to that of simple diffusion. Here, both the distribution of rotational states and the rate of rotational diffusion of the peptide are investigated and compared with the distributions expected for a simple dipole in solvents of varying relative permittivity and the predicted rate of rotational diffusion assuming a simple Stokes-Einstein-Debye model. The results present a picture of the extent to which Ewald artifacts affect the equilibrium and dynamical properties of peptides in typical liquid state simulations of charge-bearing biomolecules. II. Method A. Simulation Details. A molecular dynamics simulation of the zwitterionic form of Ala4 in explicit water was performed. The system contained 1 peptide and 496 waters in a cubic box of length 2.478 nm. The CHARMM23 all-atom force field was used to describe the solute interactions,15 with the TIP3P water model for the solvent.16 The simulation was performed in the microcanonical (NVE) ensemble resulting in an average temperature and pressure of 298 K and -216 atm, respectively. © 1997 American Chemical Society
Rotational Ewald Artifacts
J. Phys. Chem. B, Vol. 101, No. 19, 1997 3887
All bonds were constrained to their equilibrium values using SHAKE17 and a tolerance of 10-5 nm. A time step of 2 fs was used to integrate the equations of motion. Electrostatic interactions were treated using the Ewald technique18 with convergence parameter R ) 2.5 nm-1, including all lattice vectors with n2 e 36, and using a real space cutoff of Rc ) 1.2 nm. No nonbonded pair list was used. The initial configuration of the peptide was a fully extended conformation (φ ) ψ ) 180°). The peptide was solvated using a pre-equilibrated box of 512 waters, removing any waters for which the water oxygen to any peptide atom (excluding hydrogens) was less than 0.23 nm. The system was then minimized with 10 steps of steepest descent. An initial 40 ps of equilibration was performed with the peptide fixed and intermittent reassignment of the water velocities corresponding to a temperature of 300 K. A further 10 ps of simulation with the peptide free to move completed the equilibration. A production run of 10 ns was then performed, saving configurations every 0.1 ps for analysis. All 100 000 configurations have been used for the results presented here. Error estimates were obtained using 10 000 configuration subaverages. B. Equilibrium Model. By use of the molecular dynamics simulation data, the orientation of the peptide dipole in spherical polar coordinates was determined by defining cos θ ) µz and tan φ ) µy/µx, where µR are the peptide dipole unit vectors. Symmetry was then used to restrict the range of interest of both θ and φ to 0-90°, representing one octant of the simulation box. By use of methods described previously,12 the rotational potential energy surface for a rigid linear dipole can be mapped as a function of the same set of polar angles. The potential energy maps were obtained using the full Ewald potential and the same parameters as used for the simulation. Maps were calculated using ranges of θ ) 0-90° and φ ) 0-90° in increments of 0.25°. The probability distribution for θ and φ corresponding to a dipole of fixed length r at a temperature T can be written
Pr(θ,φ,T,) dθ dφ ) sin θe-Er(θ,φ,)/(kBT) dθ dφ
(1)
where is the relative permittivity, Er(θ,φ,) is the Ewald energy of the system, and kB is the Boltzmann constant. Probability distributions for θ and φ as a function of dielectric permittivity were then obtained by Boltzmann averaging over the potential surface maps assuming a temperature of 298 K. C. Dynamical Model. In order to model the rotational dynamics of the Ala4 zwitterionic peptide in solution, we have assumed that the peptide is essentially spherical and that rotational diffusion is isotropic and hence characterized by a single rotational diffusion constant. This is the Debye theory of rotational diffusion.19 Application of Debye theory to the rotational diffusion of peptides and proteins in solution is a common procedure and good results have been obtained.20,21 Although the peptide is clearly a nonspherical, flexible, asymmetric rotor, we will restrict our attention to the simplest model and leave the extension to asymmetry and flexibility for another study. For our purposes the orientation of the peptide has been defined by the vector between the terminal N and C atoms. If θ(t) is the angle between this vector at time 0 and at time t, then we have
〈Pl [cos θ(t)]〉 ) e-l(l+1)DRt ) e-t/τl
(2)
where Pl is a Legendre polynomial of rank l, DR is the rotational diffusion constant, and τl is a rotational relaxation time associated with each of the Legendre polynomial correlation
Figure 1. Radial distribution function between the terminal nitrogen and carbon atoms.
functions. The angular brackets denote an averaging over time origins. For molecules undergoing Debye diffusional rotation a plot of 1/τl against l(l + 1) should be linear with a slope equal to DR. There are many models to estimate rotational correlation times of peptides and proteins in solution.19,20,22-25 The simplest approximation assumes that the peptide is a sphere of volume V rotating in a solvent of viscosity η at a temperature T. This is the Stokes-Einstein-Debye model of rotational friction from which one obtains the rotational correlation time τ2 given by26
τ2 )
ηV kBT
(3)
Using this relationship, we have compared the calculated and observed values of τ2 to determine if the rate of rotational diffusion of the peptide has been affected. III. Results The radial distribution function between the terminal N and C atoms is displayed in Figure 1. One observes that the peptide adopted two distinct conformational states characterized by end to end distances of 0.45 and 0.72 nm corresponding to dipole moments of 22 and 35 D, respectively. A barrier at 0.6 nm separated the two states. Henceforth, we shall refer to the two states as the low dipole state and the high dipole state, respectively. Our analysis has been performed on both states separately. The high dipole state was the most highly populated at 70%. It is unclear if this distribution represents the long time equilibrium distribution for the system.14 However, for the present analysis it is sufficient that one can separate the two states. The Ala4 peptide contains three peptide groups in addition to the charged termini. In the CHARMM all-atom model the dipole moment of a peptide group is ∼4.7 D. In order to determine if the contribution of the peptide groups to the total dipole moment can be ignored, the correlation between the magnitude of the end to end distance and the magnitude of the total dipole moment was determined. In addition, the average cosine of the angle between these two vectors was also determined. The results for each state are shown in Table 1. It is apparent that the high dipole state was totally dominated by the charge separation of the two termini. However, the low dipole state contained a sizeable contribution from the peptide groups. Therefore, comparison of the results with model linear dipole calculations will be less reliable for the low dipole state.
3888 J. Phys. Chem. B, Vol. 101, No. 19, 1997
Smith et al.
TABLE 1: Simulation and Modeling Results dipole state
low
high
MD Simulation population (%) 30 N-C distance, r (nm) 0.45 correlation, Ca 0.979 correlation angle, Rb (deg) 40 radius of gyration, Rg (nm) 0.32 ( 0.01 τ1 (ps) 84 ( 4 τ2 (ps) 29 ( 2 τ3 (ps) 15 ( 1 τ4 (ps) 10 ( 1 τ5 (ps) 6(1 DR (ns-1) 5.1 ( 0.2
70 0.72 0.998 8 0.40 ( 0.01 236 ( 13 61 ( 4 31 ( 2 19 ( 1 13 ( 1 2.6 ( 0.3
Modeling 0.41 ( 0.01 0.30 ( 0.03 28 ( 3
0.52 ( 0.01 0.58 ( 0.04 54 ( 4
radius, Rc (nm) volume, V (nm3) τ2d (ps)
a C ) σ(|r||µ|)/(σ(r2)σ(µ2)), where σ(ab) ) x〈ab〉 - 〈a〉〈b〉 and r and µ are unit vectors describing the orientation of the end to end distance and dipole vector, respectively. The angular brackets denote a time average. b R ) a cos(〈r‚µ〉). c Assuming a sphere where R ) x5/3Rg. d Calculated using eq 3 at 298 K assuming the viscosity of TIP3P water to be 0.39 cP.
Figure 3. Normalized probability distributions Pr for the polar angles θ (top) and φ (bottom) corresponding to the high dipole state at 298 K. The solid line is the distribution obtained from the molecular dynamics simulation. The corresponding distributions calculated for a model dipole (r ) 0.72 nm) under the same conditions as the simulation are also displayed assuming dielectric permittivities of 1 (dotted), 10 (dashed), and infinity (thin solid).
Figure 2. Normalized probability distributions Pr for the polar angles θ (top) and φ (bottom) corresponding to the low dipole state at 298 K. The solid line is the distribution obtained from the molecular dynamics simulation. The corresponding distributions calculated for a model dipole (r ) 0.45 nm) under the same conditions as the simulation are also displayed assuming dielectric permittivities of 1 (dotted), and infinity (thin solid).
In contrast, the high dipole state should afford a good comparison between theory and simulation. The probability distributions of the polar angles θ and φ for the low and high dipole states are displayed in Figures 2 and 3, respectively. Distributions were restricted to one dimension to improve both the statistical error and the clarity of presentation. The predicted probability distributions obtained from the full Ewald rotational potential energy surface assuming different
values of the relative permittivity are also shown. Using the equilibrium model the barrier to rotation of a diple in a dielectric continuum is calculated to be 0.3/ and 2.0/ kJ/mol for the low and high dipole states, respectively. Unfortunately, to our knowledge the relative permittivity of TIP3P water has not been determined and hence a direct comparison with the distribution corresponding to this water model is not possible. However, Ewald artifacts are predicted to be small for both dipole states. It is clear that for both conformational states of the peptide the angle distributions obtained from the simulation were in agreement with a calculated distribution corresponding to a relative permittivity between 10 and infinity (or random). It is apparent that, within the accuracy of this 10 ns simulation, no significant Ewald artifacts on the equilibrium distribution of rotational states were observed. This agrees with our previous conclusions assuming a simple nonmolecular dielectric continuum representation of the solvent.12 The first five Legendre polynomials corresponding to changes in orientation of the peptide have been calculated using the dipole moment vector time series. In generating the polynomials for each state, data were only included for time periods in which the peptide remained in that specific state. The correlation times τl have been obtained from the resulting log-time plots. Figure 4 displays the inverse of the correlation time as a function of l (l + 1) for both states. The plots are linear, indicating that the rotational dynamics of both states were consistent with pure Debye rotational diffusion. The average radius of gyration for each state has been determined from the simulation (see Table 1). Assuming the peptide is essentially spherical in each state, one can estimate the volume of the peptide in each state. However, to obtain
Rotational Ewald Artifacts
Figure 4. Inverse of the decay time characterizing the Legendre polynomials of rank l obtained from the orientational dynamics of the dipole moment vector as a function of l(l + 1). Results are for the low dipole (solid) and high dipole (dashed) states.
the rotational correlation time τ2 one also needs to know the viscosity of the solvent. To our knowledge the viscosity of TIP3P water is unknown. Therefore, we have estimated the value by reference to the self-diffusion constant of TIP3P water. From the simulation the self-diffusion constant of the water was determined to be (4.9 ( 0.1) × 10-9 m2 s-1. This is high compared to the experimental value but in good agreement with other literature values for this model.6,27 Assuming a StokesEinstein relation (D ) kBT/(6πηa)) between the diffusion constant D and the viscosity η, and knowing the experimental diffusion constant (2.2 × 10-9 m2 s-1)28 and viscosity (0.851 cP)29 of water at 300 K, the estimated viscosity of TIP3P water is 0.39 ( 0.01 cP. The low viscosity correlates with the high diffusion rate and is a property of the water model rather than a consequence of the Ewald summation.27 With our approximation for the viscosity of TIP3P one finds that the predicted values (eq 3) for τ2 of 28 ( 3 and 54 ( 4 ps, for the low and high dipole states, respectively, closely match the values determined from the simulation of 29 ( 2 and 61 ( 4 ps, respectively (see Table 1). Interestingly, good agreement is obtained without requiring the presence of a long-lived hydration shell surrounding the peptide. This is in agreement with other simulation results30 but contradicts the view of protein hydrodynamics that often invokes a 0.3 nm water hydration shell.31,32 The presence of a solvent shell around the peptide would result in a larger predicted value of τ2, suggesting the simulation results display an increased rate of rotational diffusion. This is opposite that expected from the presence of any Ewald artifacts. Therefore, within the approximations of the present calculations, there appears to be no evidence that an underlying Ewald induced rotational potential results in a significantly decreased rotational diffusion rate. IV. Conclusions It has been shown that the equilibrium and dynamical rotational properties of a small zwitterionic tetrapeptide in aqueous solution were not significantly affected by the use of the Ewald technique for the treatment of electrostatic interactions. The equilibrium orientational distribution of the peptide within the central box was consistent with that of an equivalent dipole in a medium of relative permittivity equal to or greater than 10. Both the simulation results shown here and the previous model calculations12 indicate that rotational Ewald artifacts are negligible for relative permittivities of this mag-
J. Phys. Chem. B, Vol. 101, No. 19, 1997 3889 nitude. The rotational dynamics of the peptide were observed to be consistent with that of a simple spherical particle, with no significant hydration layer, in a medium of viscosity 0.39 cP. This also suggests that the rotational diffusion of the peptide was not significantly hindered by artificial potential barriers arising from application of the Ewald technique. Recently, Luty and van Gunsteren presented simulation evidence for the presence of Ewald artifacts.11 Their system involved two charges placed at a separation of L/2 on the surface of a sphere that excluded water. For this system, the relative permittivity of the system is not so easy to assign. Since there is no dielectric screening across the sphere, the effective relative permittivity is considerably lower than that in bulk solution, and therefore, Ewald artifacts are observed. This is still consistent with the model presented by Smith and Pettitt12 but requires a less intuitive guess for . For charged residues interacting across a globular protein the protein itself will provide some degree of dielectric screening. Unfortunately, exactly how much screening is the subject of much debate.33-36 However, the alternative approaches for the treatment of electrostatic interactions in biomolecular simulations usually involve cutoffs or switching functions and would probably neglect this interaction altogether. It is not clear that this is more satisfactory than the slight system size effect introduced by the Ewald potential. The current results, together with previous work,37,38 suggest that the simulation of biomolecular systems in solvents with high relative permittivities using the Ewald technique provides a reasonable description of the structure and dynamics of these systems and is free from significant rotational artifacts. Acknowledgment. This project was supported by the Kansas Agricultural Experimental Station (Publication 97-190-J). We would also like to thank the Robert A. Welch Foundation and the National Science Foundation for partial support. H.D.B. thanks the Keck Center for a graduate fellowship under Grant No. BIR-92-56580 from the National Science Foundation. References and Notes (1) Smith, P. E.; Dang, L. X.; Pettitt, B. M. J. Am. Chem. Soc. 1991, 113, 67-73. (2) Forester, T. R.; McDonald, I. R. Mol. Phys. 1991, 72, 643-660. (3) York, D. M.; Wlodawer, A.; Pedersen, L. G.; Darden, T. A. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 8715-8718. (4) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 1008910092. (5) Fincham, D. Mol. Simul. 1994, 13, 1-9. (6) Smith, P. E.; Pettitt, B. M. Comput. Phys. Commun. 1995, 91, 339344. (7) Card, D. N.; Valleau, J. P. J. Chem. Phys. 1970, 52, 6232-6240. (8) Valleau, J. P.; Whittington, S. G. A Guide to Monte Carlo for Statistical Mechanics: 1. Highways. In Statistical Mechanics A. Modern Theoretical Chemistry; Berne, B. J., Ed.; Plenum Press: New York, 1977; pp 137-168. (9) van Gunsteren, W. F.; Berendsen, H. J. C.; Rullmann, J. A. C. Faraday Discuss. Chem. Soc. 1978, 66, 58-70. (10) Figueirido, F.; Del Buono, G. S.; Levy, R. M. J. Chem. Phys. 1995, 103, 6133-6142. (11) Luty, B. A.; van Gunsteren, W. F. J. Phys. Chem. 1996, 100, 25812587. (12) Smith, P. E.; Pettitt, B. M. J. Chem. Phys. 1996, 105, 4289-4293. (13) Pratt, L. R.; Haan, S. W. J. Chem. Phys. 1981, 74, 1873-1876. (14) Blatt, H. D.; Smith, P. E.; Pettitt, B. M. To be submitted. (15) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Wiorkiewicz-Kuczera, J.; Karplus, M. FASEB J. 1992, 6, A143. (16) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (17) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327-341.
3890 J. Phys. Chem. B, Vol. 101, No. 19, 1997 (18) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, Ser. A 1980, 373, 27-56. (19) Debye, P. Polar Molecules; Dover: New York, 1945. (20) Mu¨ller, J. J. Biopolymers 1991, 31, 149-160. (21) Kakalis, L. T.; Kumosinski, T. F. Biophys. Chem. 1992, 43, 3949. (22) Perrin, F. J. Phys. Radium 1934, 5, 497-511. (23) Gierer, A.; Wirtz, K. Z. Naturforsch., A 1953, 8, 532-538. (24) Hu, C.; Zwanzig, R. J. Chem. Phys. 1974, 60, 4354-4357. (25) Brune, D.; Kim, S. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 38353839. (26) Boere, R. T.; Kidd, R. G. In Rotational Correlation Times in Nuclear Magnetic Spectroscopy; Webb, G. A., Ed.; Annual Reports on NMR Spectroscopy 13; Academic Press: New York, 1982. (27) Tasaki, T.; McDonald, S.; Brady, J. W. J. Comput. Chem. 1993, 14, 278-284. (28) Gillen, K. T.; Douglas, D. C.; Hoch, M. J. R. J. Chem. Phys. 1972, 57, 5117-5119.
Smith et al. (29) Weast, R. C. Handbook of Chemistry and Physics; CRC Press: Cleveland, 1975. (30) Smith, P. E.; van Gunsteren, W. F. J. Mol. Biol. 1994, 236, 629636. (31) Squire, P. G.; Himmel, M. E. Arch. Biochem. Biophys. 1979, 196, 165-177. (32) Kakalis, L. T.; Baianu, I. C. Arch. Biochem. Biophys. 1988, 267, 829-841. (33) Gilson, M. K.; Honig, B. H. Biopolymers 1986, 25, 2097-2119. (34) Harvey, S. C. Proteins: Struct. Funct. Genet. 1989, 5, 78-92. (35) Smith, P. E.; Brunne, R. M.; Mark, A. E.; van Gunsteren, W. F. J. Phys. Chem. 1993, 97, 2009-2014. (36) Simonson, T.; Brooks, C. L., III. J. Am. Chem. Soc. 1996, 118, 8452-8458. (37) Smith, P. E.; Pettitt, B. M. J. Chem. Phys. 1991, 95, 8430-8441. (38) Schreiber, H.; Steinhauser, O. Biochemistry 1992, 31, 5856-5860.
