Baixe A Second Generation Force Field for the Simulation of e outras Notas de estudo em PDF para Engenharia Elétrica, somente na Docsity!
J. Am. Chem. SOC. 1995, 117, 5179-5197 5 179
A Second Generation Force Field for the Simulation of
Proteins, Nucleic Acids, and Organic Molecules
Wendy D. Cornell? Piotr Cieplak,’ Christopher I. Bayly,s Ian R. Gould,l
Kenneth M. Merz, Jr.,” David M. Ferguson,&David C. Spellmeyer: Thomas Fox,
James W. Caldwell, and Peter A. Kollman*
Contribution from the Department of Pharmaceutical Chemistry, University of California, San Francisco, California 94143
Received November IO, 1994@
Abstract: We present the derivation of a new molecular mechanical force field for simulating the structures, conformational energies, and interaction energies of proteins, nucleic acids, and many related organic molecules in condensed phases. This effective two-body force field is the successor to the Weiner et al. force field and was developed with some of the same philosophies, such as the use of a simple diagonal potential function and electrostatic potential fit atom centered charges. The need for a 10-12 function for representing hydrogen bonds is no longer necessary due to the improved performance of the new charge model and new van der Waals parameters. These new charges are determined using a 6-31G* basis set and restrained electrostatic potential (RESP) fitting and have been shown to reproduce interaction energies, free energies of solvation, and conformational energies of simple small molecules to a good degree of accuracy. Furthermore, the new RESP charges exhibit less variability as a function of the molecular conformation used in the charge determination. The new van der Waals parameters have been derived from liquid simulations and include hydrogen parameters which take into account the effects of any geminal electronegative atoms. The bonded parameters developed by Weiner et al. were modified as necessary to reproduce experimental vibrational frequencies and structures. Most of the simple dihedral parameters have been retained from Weiner et al., but a complex set of 4 and yj parameters which do a good job of reproducing the energies of the low-energy conformations of glycyl and alanyl dipeptides has been developed for the peptide backbone.
Introduction The application of computer-based models using analytical potential energy functions within the framework of classical mechanics has proven to be an increasingly powerful tool for studying molecules of biochemical and organic chemical interest. These applications of molecular mechanics have employed energy minimization, molecular dynamics, and Monte Carlo methods to move on the analytical potential energy surfaces. Such methods have been used to study a wide variety of phenomena, including intrinsic strain of organic molecules, structure and dynamics of simple and complex liquids, ther- modynamics of ligand binding to proteins, and conformational transitions in nucleic acids. In principle, they are capable of giving insight into the entire spectrum of non-covalent interac- tions between molecules, and, when combined with quantum mechanical electronic structure calculations, modeling covalent
bonding changes, essentially all molecular reactions and interac-
tions. Given their importance, much effort has gone into consideration of both the functional form and the parameters that must be established in order to apply such analytical potential energy functions (or “force fields”).
t Graduate Group in Biophysics, University of California, San Francisco.
* Permanent address: Department of Chemistry, University of Warsaw,
8 Current address: Merck Frosst Canada, Inc., C.P. 1005 Pointe Claire-
Current address: Deuartment of Chemism. Universitv of Manchester.
Pasteura 1, 02-093, Warsaw, Poland.
Domal, Quebec H9R 4P8, Canada.
Lancs M13 9PL, U.K. Current address: Department of Chemistry, The Pennsylvania State University, State College; PA 16802.
Minnesota, Minneapolis, MN 55455.
Bi Current address: Department of Medicinal Chemistry, University of
Current address: Chiron Corporation, Emeryville, CA 94608.
*Author to whom correspondence and reprint requests should be @ Abstract published in Advance ACS Abstracts, April 15, 1995.
addressed.
In the area of organic molecules, the book by Allinger and
Burkert’ provides a thorough review pre-1982 and the subse- quent further development of the MM2* and MM33 force fields by Allinger and co-workers has dominated the landscape in this area. The number of force fields developed for application to biologically interesting molecules is considerably greater, prob- ably because of the greater complexity of the interactions which involve ionic and polar groups in aqueous solution and the difficulty of finding an unequivocal test set to evaluate such force fields. Many of these force fields developed prior to 1987 are described briefly by McCammon and Harvey. Given the complexities and subjective decisions inherent in such biological force fields, we have attempted to put the development of the force field parameters on a more explicitly stated algorithmic basis than done previously, so that the force field could be extended by ourselves and others to molecules and functional groups not considered in the initial development. This is important, because, if the assumptions, approximations, and inevitable imperfections in a force field are at least known, one can strive for some cancellation of errors. Approximately a decade ago, Weiner et al.596 developed a force field for proteins and nucleic acids which has been widely
(1) Burke& U.; Allinger, N. J. Molecular Mechanics; American Chemical Society: Washington, DC, 1982. (2)Allinger, N. L. J. Am. Chem. SOC. 1977, 99, 8127-8134 and subsequent versions, e.g. MM2-87, MM2-89, MM2-91. (3) Allinger, N. L.; Yuh, Y. H.; Lii, J.-H. J. Am. Chem. SOC. 1989,ll I ,
(4) McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, 1987. ( 5 ) Weiner, S. J.; Kollman, P. A,; Case, D. A,; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. J. Am. Chem. SOC. 1984,106, 765-
(6) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comp. Chem. 1986, 7, 230-252.
8551-8566, 8566-8576, 8576-8582.
0002-7863/95/1517-5179$09.00/0 0 1995 American Chemical Society
5180 J. Am. Chem. SOC., Vol. 117, No. 19, 1995 Comell^ et^ al.
used. Important independent tests of this force field were performed by Pavitt and Hall for peptides’ and Nilsson and Karplus8for nucleic acids and it was found to be quite effective. Nonetheless, it was developed in the era before one could routinely study complex molecules in explicit solvent. Weiner et al. attempted to deal with this issue by showing that the same force field parameters could be effectively used both without explicit solvent (using a distance-dependent dielectric constant ( E = Ru)) and with explicit solvent ( E = 1) on model systems. Further support for this approach was provided by molecular dynamics simulations of proteins9-” and DNAl2.l3 which compared the implicit and explicit solvent representations. As computer power has grown, it has become possible to carry out more realistic simulations which employ explicit solvent representations. It is therefore appropriate that any new force field for biomolecules focus on systems modeled in the presence of an explicit solvent representation. This approach has been pioneered by Jorgensen and co-workers in their OPLS (Optimized Potentials for Liquid Simulations) m0de1.I~ In particular, the development of parameters which reproduce the
enthalpy and density of neat organic liquids as an essential
element ensures the appropriate condensed phase behavior. The OPLS non-bonded parameters have been combined with the Weiner et al. bond, angle, and dihedral parameters to create the OPLS/Amber force field for peptides and proteins,I5 which has also been effectively used in many systems.I We have been influenced by the OPLS philosophy of balanced solvent-solvent and solute-solvent interactions in our thoughts about a second-generation force field to follow that of Weiner et aL5v6 The Weiner et al. force field used quantum mechanical calculations to derive electrostatic potential (ESP) fit atomic centered charges, whereas the OPLS charges were derived empirically, using mainly the liquid properties as a guide. For computational expediency, Weiner et al. relied principally on the STO-3G basis set for their charge derivation. This basis set leads to dipole moments that are approximately equal to or smaller than the gas-phase moment but tends to underestimate quadrupole moments. Thus, it is not well balanced with the commonly used water models (SPC/E,” TIP3P,I8 TIP4PI8) which have dipole moments that are about 20% higher than the gas-phase value for water. These water models, which have empirically derived charges, include condensed-phase electronic polarization implicitly. Kuyper et aZ.l9 suggested that the logical choice of a basis set for deriving ESP-fit partial charges for use in condensed phases is the 6-3 lG* basis set, which uniformly overestimates molecular polarity. Standard ESP charges derived with that basis set were shown ~ (7) Pavitt, N.; Hall, D. J. Compur. Chem. 1984, 5, 441- (8) Nilsson, L.; Karplus, M. J. Comput. Chem. 1986, 7, 591-616. (9) Tilton, R. F.; Singh, U. C.; Weiner, S. J.; Connolly, M. L.; Kuntz, I. D., Jr.; Kollman, P. A.; Max, N.; Case, D. J. Mol. Eiol. 1986, 192, 443-
(10) Guenot, J. M.; Kollman, P. A. Prorein Sci. 1992, 1, 1185-1205. ( 1 1) York, D. M.; Wlodawer, A.; Redersen, L.; Darden, T. A. Proc. Narl. (12) Sinsh. U. C.: Weiner. S. J.: Kollman. P. A. Proc. Natl. Acad. Sci.
Acad. Sci. U.S.A. 1994, 91, 8715-8718.
U.S.A’l983, 82, 755-759. ’ (13) Seibel, G. L.; Singh, U. C.; Kollman, P. A. Proc. Natl. Acad. Sci. U.S.A. 1985, 82, 6537-6340. (14) Jorgensen, W. L.; Pranata, J. J. Am. Chem. SOC. 1990, 112,2008-
(15) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. SOC.1988, 110, 1657 - 1666. (16) (a) Tirado-Rives, J.; Jorgensen, W. L. J. Am. Chem. SOC. 1990,112, 2773-2781. (b) Orozco, M.; Tirado-Rives, J.; Jorgensen, W. L. Eiochem- istry 1993, 32, 12864-12874. (17) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6271. (18) Jorgensen, W. L.; Chandreskhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1982, 79, 926-935. (19) Kuyper, L.; Ashton, D.; Men, K. M., Jr.; Kollman, P. A. J. Phys. Chem. 1991, 95, 6661-6666.
to lead to excellent relative free energies of solvation for benzene, anisole, and trimethoxyani~ole.’~ A 6-31G* based ESP-fit charge model, like the OPLS model, is capable of giving an excellent reproduction of condensed- phase inter molecular properties such as liquid enthalpies and densities and free energies of solvation.20 A major difference between such a model and most others is the magnitude of the charges on hydrocarbons. For example, 6-3 lG* standard ESP charges derived from the trans conformation of butane have values of -0.344 for the methyl carbon and 0.078 for the methyl hydrogen. In both cases, however, the carbon and hydrogen charges offset each other, resulting in small net charges on the methyl groups of -0.1 10 and -0.059 for the trans and gauche charges, respectively. Furthermore, free energy perturbation calculations involving the perturbation of methane with standard ESP charges (qc = -0.464 and q H = 0.116) to methane with charges of 0.0 in solution yield essentially no change in free energy.21 The standard ESP charges also result in conforma- tional energies for butane which are in reasonable agreement with experiment, when used with a 1-4 electrostatic scale factor of m.2. Nevertheless, the 6-3 lG* standard ESP charges are less than ideal for two reasons. First, when charges generated using different conformations of a molecule are compared, there is often considerable variation seen. This was demonstrated by Williams, who studied the conformational variation of ESP-fit charges in alanyl dipeptide for 12 different conformations. Butane is another example, where charges from the gauche conformation have values of -0,197 and 0.046 for the methyl carbon and hydrogen, respectively. Another example is pro- pylamine, which was studied at length by Comell et aL2O Five low-energy conformations can be identified for propylamine, and the 6-31G* standard ESP charges calculated for each conformation show significant variation. The average and standard deviation for the charge on a given atom over the five conformations are as follows: a-carbon qav = 0.339 and IJ =
0.059, /3-carbon qav = 0.033 and u = 0.060, and y-carbon qav
= -0.205 and u = 0.146. This inconsistency is potentially problematic in terms of deriving other force field parameters which may be sensitive to the variation. Furthermore, it reduces the reproducibility of a particular calculation, which is not a problem in other force fields where the charges are assigned empirically. The second reason that the 6-3 lG* standard ESP charges are less than ideal is that the charges on “buried” atoms (such as the sp3 carbons described above for butane and propylamine) are statistically underdetermined and often assume unexpectedly large values for nonpolar atoms. Bayly et aLZ3 found that the electrostatic potential of methanol could be fit almost equally well using either the standard ESP charges determined by the linear least-squares fit or an altemative set of charges derived with the methyl carbon constrained to have a much smaller value. Considering the problems associated with the standard ESP charge model, it might seem tempting to adopt the OPLS approach of empirically derived charges. However, any empiri- cally derived charge model cannot easily describe transition states and excited states, as can an electrostatic potential fit
(20) Cornell, W.; Cieplak, P.; Bayly, C.; Kollman, P. A. J. Am. Chem. SOC. 1993, 115, 9620-9631. (21) (a) Sun, Y.X.; Spellmeyer, D.; Pearlman, D. A.; Kollman, P. A. J. Am. Chem. SOC. 1992, 114, 6798-6801. (b) Sun, Y. X.; Kollman, P. A. Hydrophobic Solvation of Methane and Nonbond Parameters of the TIP3P Water Model. J. Cornput. Chem., in press. Pang, Y. P.; Kollman, P. A., unpublished. (22) Williams, D. E. Biopolymers 1990, 29, 1367-1386. (23) Bayly, C.; Cieplak, P.; Comell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269-10280.
5182 J. Am. Chem. SOC., Vol. 117, No. 19, 1995
Table 1. List of Atom Types“
atom type description
Come11 et al.
carbon CT C CA CM cc
cv
cw
CR
CB
C*
CN
CK
CQ
nitrogen N NA
NB
NC
N*
N
N oxygen OW OH
os 0 0 2 sulfur S SH phosphorus P hydrogen H Hw HO HS HA HC
H
H
H
HP
H
H
any sp3carbon any carbonyl sp2carbon any aromatic sp2carbon and ( C E of Arg) any sp2 carbon, double bonded sp2aromatic in 5-membered ring with one substituent + next to nitrogen ( C y in His) sp2aromatic in 5-membered ring next to carbon and lone pair nitrogen (e.g. C6 in His (6)) sp2aromatic in 5-membered ring next to carbon and NH (e.g. C6 in His ( E ) and in Trp) sp2 aromatic in 5-membered ring next to two nitrogens ( C y and C E in His) sp2aromatic at junction of 5- and 6-membered rings (C6 in Trp) and both junction atoms in Ade and Gua sp2aromatic in 5-membered ring next to two carbons (e.g. C y in Trp) sp2junction between 5- and 6-membered rings and bonded to CH and NH (Ce in Trp) sp2carbon in 5-membered aromatic between N and N-R (C8 in purines) sp2carbon in 6-membered ring between lone pair nitrogens (e.g. C 2 in purines) sp2nitrogen in amides sp2nitrogen in aromatic rings with hydrogen attached (e.g. protonated His, Gua, Trp) sp2 nitrogen in 5-membered ring with lone pair (e.g. N7 in purines) sp2 nitrogen in 6-membered ring with lone pair (e.g. N3 in purines) sp2 nitrogen in 5-membered ring with carbon substituent (in purine nucleosides) sp2 nitrogen of aromatic amines and guanidinium ions sp3nitrogen sp3oxygen in TIP3P water sp3oxygen in alcohols, tyrosine, and
sp3 oxygen in ethers sp2oxygen in amides sp2oxygen in anionic acids sulfur in methionine and cysteine sulfur in cysteine phosphorus in phosphates H attached to N H in TIP3P water H in alcohols and acids H attached to sulfur H attached to aromatic carbon H attached to aliphatic carbon with
H attached to aliphatic carbon with
H attached to aliphatic carbon with
H attached to aliphatic carbon with
H attached to carbon directly bonded to
protonated carboxylic acids
no electron-withdrawing substituents
one electron-withdrawing substituent
two electron-withdrawing substituents
three electron-withdrawing substituents
formally positive atoms (e.g. C next to NH3+ of lysine)
electronegative neighbor (e.g. hydrogemon C5 of Trp, C6 of Thy)
electronegative neighbors (e.g. H8 of Ade and Gua and H2 of Ade)
H attached to aromatic carbon with one
H attached to aromatic carbon with two
a See refs 5 and 6.
determined by linear interpolation between pure single and double bond values using the observed bond distances and the
KO value taken from vibrational analysis of a simple sp2 atom
containing fragments such as benzene and NMA. That this approach was reasonably successful is supported by the reason- able agreement found in nucleic acid base vibrational analysis
and suggested by the critical analysis of Halgren of the diagonal force constants used in different force fields. One “difficulty” arose in the development of this new force field compared to that of Weiner et al. which was related to the switch to the 6-31G* basis set for charge derivation. With 6-31G* standard ESP charges and a 1-4 electrostatic scale
factor of M. 2 rather than U2.0 (see below), we found that the
exocyclic -NH2 groups of the bases moved considerably away
from their res and 8 , values upon energy minimization. This
problem was considerably reduced with RESP charges and a 1-4 electrostatic scale factor of M.2, so we chose not to
selectively increase the KO values around the - N H 2 group to
force it to more “canonical” geometries. In general, however, one might have resorted to a more
complex optimization of re,.,, Oeq, K,, and KO to ensure that the
geometries of simple fragments were as close as possible to experiment after energy minimization, rather than taking r, and 8 , from experiment and assuming little distortion would occur (which is generally the case, with the slight exception of the case of the - N H 2 groups noted above). We chose not to undertake a more time-consuming iterative self-consistent derivation of geometrical parameters, because of our assumption that any such errors which we were making were of much smaller consequence for accurately representing conformations and intermolecular interactions than the inaccuracies remaining in the dihedral and non-bonded (charge and VDW) parameters.
C. Dihedral Parameters. Weiner et ~ 1. ~ developed a 9 ~
limited set of general and specific dihedral parameters which were appropriate for the functionalities found in proteins and DNA and calibrated to adjust the energies of small model compounds. In this strategy, a dihedral parameter is optimized on the simplest molecule possible and then applied to larger and more complex molecules. This approach is in contrast to one employed by many other force field developers where the parameters are optimized to best reproduce the conformational energies of a large number of molecules. An advantage of our approach is the lack of dependence of the resulting parameters on the particular molecules chosen for the test set. For the most part, a minimalist approach has been retained
with regards to dihedral parameters. For example, we have only
a 3-fold Fourier component (V3) for dihedrals around -C-C- bonds, with the exception of cases such as E-C-C-E’ where E and E’ are electronegative atoms like 0 or F. In these cases, there is a “gauche” effect which stabilizes the gauche conforma- tion over the trans and this can be modeled with a 2-fold Fourier component (V2). The rotation around phosphorus-ester bonds
(CT-OS-P-OS) also requires a 2-fold component. In these
cases, we have been able to go beyond the Weiner et al. force field by making use of reasonably high level ab initio models (MW6-31G*) to fit the values of such V2 Fourier components. Two exceptions were made to the principle of adding extra Fourier terms to the dihedral energies only in the presence of a compelling physical basis. These exceptions are the dipeptide
v and 4 and the nucleoside x dihedrals. Here we used additional
Fourier components to try to reproduce as well as possible the relative energies of the alanyl and glycyl dipeptides and a model nucleoside fragment calculated at a high level of theory without the requirement of “a physical picture’’. An altemative approach would be to empirically adjust the atomic partial charges to
achieve the same aim. Given the power of the RESP methodol-
ogy for deriving atomic partial charges which lead to good representations of intermolecular interactions and the importance of maintaining an accurate balance between intra- and inter- molecular interactions, we chose to empirically adjust the terms
in the Fourier series for v and 4 as well as x.
(31)Halgren, T.A. J. Am. Chem. SOC. 1990, 112, 4710-4723.
Simulation of Proteins and Nucleic Acids
Table 2. Standardized Parameters for Scaling Algorithms
J. Am. Chem. SOC., Vol. 117, No. 19, 1995 5183
unavailable for the OPLS force field. These Monte Carlo simulations were the f i s t calculations carried out as part of the development of this new force field, and as such employed 6-31G* standard ESP charges. The electrostatic contribution for the n-alkanes was very small regardless of the charge model-at most a few tenths of a kcdmol. We note that the standard ESP charges for benzene (qc= -0.145 and qH= 0.145) accurately reproduce the quadrupole moment of that molecule. We have taken most of the remaining VDW parameters from the OPLS modelI5-sp2 and sp3 N; sp2 0, ether ester (OS), hydroxyl (OH) and TIP3P water (OW) sp3 oxygens; and sulfur (SH and S)-since it has been optimized for reproducing liquid properties. The Weiner et ~ 1. phosphorus (P) parameters were ~ 9 ~ not re-optimized since that atom is most frequently found buried inside of four other heavy atoms. The VDW model is minimalist as well, with some exceptions. A standard VDW parameter is used for a given atom and hybridization, e.g. all sp2 carbons have the same VDW parameters. The only heavy atom exceptions are sp3 0, where oxygens in water (OW), alcohol (OH), and ether (OS) have slightly different parameters, as found in OPLS. We suspect that this is due to the use of a zero VDW radius on hydrogens bound to oxygen, so that an effectively larger _R_* is required for a water oxygen than alcohol than ether. A significant departure has been made from the previous model in the treatment of hydrogens. The current model does not employ 10-12 hydrogen bonding Ha .X parameters, although these are still supported within the AMBER software. The original Hagler et and OPLS a p p r o a ~ h ' ~? ' ~suggested a zero _R_* and 6 for hydrogen binding hydrogens. Thus the TIP3P water model has _R_* and 6 equal to 0.0 for its hydrogen
(HW). We opted not to develop a new water model, but to use
the T P 3 P one. Hydrogen and helium are unique in the periodic table in not having an inner shell of electrons. Consequently, it makes physical sense for the hydrogen VDW radius, unlike other atoms, to be very sensitive to its bonding environment. This has been extensively analyzed for the hydrogen _R_* in X-C-H systems by Gough et al. and Veenstra et al.,26,27 who demon-
strated the sensitivity of R* to the electron-withdrawing proper-
ties of X. For example, a ''normal" C-H has VDW R* = 1.
8,; whereas in CF3-H it is -0.3 A shorter and in CH3NH3+ it
is -0.4 8, shorter still.
We have employed the following approach here. A C-H
has R* = 1.487 8, and, based on nucleic and base pairing energy
minimization, an N-H has R* = 0.6 8,. This qualitative
dependence on electronegativity makes physical sense. Based on the Veenstra et aL2I studies we have chosen to reduce the
R* on sp3 C-H atoms by 0.1 8, for each electronegative (0, N,
F, S) substituent. The hydrogen atom types are then defined as H1, H2, and H3 for 1, 2, and 3 electronegative groups,
respectively. The hydrogen R* is reduced by 0.4 8, for each
neighboring positively charged group (atom type HP). For sp
C-H, R* has been reduced by 0.05 8, for each electronegative
neighbor (atom types H4 and H5). Given our retention of the simplicity of a 6-12 rather than a 6-exponential VDW representation, we have continued to reduce 1-4 VDW interactions since the 6- 12 approximation and the lack of polarization in the model both will lead to exaggerated short-range repulsion. It is difficult to determine the scale factor unambiguously so we have retained the value of 112.0 used by Weiner et a1.s. E. Electrostatic Energies. In Cornel1 et aL20 and Cieplak et we have extensively analyzed the development of our
~~ bond r," KP
pure C=C 1.336e 57of pure C-N 1.4499 337h pure C=N 1.273' 570,
pure C-C 1.507' 317d
torsion re( v2k pure X-C-C-X 1.507c 0.0' partial X-C=C-X 1.397" 14.5" pure X-C=C-X 1.336' 30.0" pure X-C-N-X 1.4499 0.w partial X-C=N-X 1.3354 10.0' pure X-C=N-X 1.273' 30.0" a In A. In kcaV(mo1 A2). Microwave data from acetone (ref 32). Value taken from MM2, ref 2. e Microwave data from propene (ref 32). f Default from NMA normal mode analysis for carbonyl force constant. 9 Benedetti structural^ data^ (ref^ 33).^ *^ Value derived from normal mode analysis on NMA. ' Microwave data from methylenimine (ref 32). J Default value, see footnote f. I; In kcaVmo1. Assumed free rotation about pure C-C single bond. Structural data^ from benzene (ref 32). From normal modes analysis of^ benzene.^ Approximate rotational barrier of ethylene is -60 kcaVmol (see ref 34). p Assumed free rotation about a^ pure single^ C-N^ bond.^^4 Benedetti structural data (ref 33). 'Reference 35. Calculated rotational barrier in methylenimine is 57.5 kcaUmol (see ref 36).
In our previous force field, the bond length and V2 parameters for X-C-N-X and X-C-C-X fragments involving sp hybridized atoms were determined by a linear interpolation approach (according to the experimental bond length) between the known barriers of pure single, pure double, and partial double bonded systems (benzene for X-C-C-X and NMA for X-C-N-X). We have used the same approach here, but have adjusted the V2 term of benzene to more accurately describe its out-of-plane frequencies (Weiner et ~ 1. had used the ~ 3 ~ V derived for a united atom model of benzene, which was significantly different). Table 2 presents the parameters used. For example, given a C(sp2)-C(sp2) bond length, its bond stretching force constant is linearly interpolated between the values for pure single bond and double bond given in Table 2. Its V2 torsional potential is interpolated between the values for pure double and partial double or between partial double and single, dependin on whether the bond length is greater or less than the 1.397 f of benzene. This is exactly the procedure used by Weiner et ~ 1. ~ 9 ~ D. VDW Parameters. Given the success of the OPLS approach in modeling liquids, we have developed all-atom sp carbon and aliphatic hydrogen VDW parameters by canying out Monte Carlo simulations on C b , C2H6, C3H8, and C4H
liquids and empirically adjusting R* and E for the C and H to
reproduce the densities and enthalpies of vaporization of these Such parameters have also been employed in calcula- tions of relative free energies of solvation of CQ, C2H6, and C3H8.21,38We also derived VDW parameters for sp2 C and aromatic H employing Monte Carlo simulations on benzene liquid and adjusting the _R_* and 6 of these atoms to reproduce the density and enthalpy of liquid benzene.37 At the time these parameters were developed, such all-atom parameters were
(32) Harmonv. M.: Laurie. V.: Kuezkowski. R.: Schwendeman. R.: Ramsay, D.; Livas, F.; Lafferty, W.; Maki, A. J. Phys. Chem. Re$ Data 1979, 8, 619-721. (33) Benedetti. E. In Peutides-Proceedings o f the 5th American Peutide Symposium; Goodman, MI, Meienhofer, J.,-Edi.; J. Wiley and Co.: 'New York, 1977; pp 257-273. (34) Douglas, J.; Rabinovich, B. S.; Looney, F. J. Chem. Phys. 1955,
23. 315-323.
Symposium; Goodman, MI, Meienhofer, J.,-Edi.; J. Wiley and Co.: 'New York, 1977; pp 257-273. (34) Doudas, J.; Rabinovich, B. S.; Looney, F. J. Chem. Phys. 1955,
23. 315-323. (35) Momany, F.; McGuire, R.; Burgess, A,; Scheraga, H. J. Phys. Chem. 1975, 79, 2361-2381. (36) Lehn, J.; Munsch, B.; Millie, P. H. Theor. Chim. Acta 1970, 16, (37) Spellmeyer, D., unpublished.
35 1-372.
(38) Sun, Y.; Kollman, P. A. Hydrophobic Solvation of Methane and Nonbond Parameters of the TIP3P Water Model. J. Comput. Chem., accepted for publication.
Simulation of Proteins and Nucleic Acids
Table 4. Results for Alcohols and Ethers (Energies in kcal/mol, Angles in deg) parameter this work MM3" experiment
J. Am. Chem. Soc., Vol. 117, No. 19, 1995 5185
Dimethyl Ether AE(ec1ipsed-staggered) 2.74 2. O(C-0-C)(staggered) 112.3 111. B(C-0-C)(eclipsed) 113.
AE(Cs-Cz) 0.12 0. AE(c2v-c~) 3.98 4. structural parameters C2 conformation
Tetrahydrofuran
qe 0. e(c-o-c) 108.8 108. e(c -0 -0) 106.8 106. e(c-c-c) 100.4 101.
qe 0. e(c-o-c) 105.4 104. e(c-0-0) 105.1 105 .O e(c-c-c) 103.6 103.
C , conformation
Methyl Ethyl Ether AE(gauche-trans) 1.46 1.
4(gauche) 76.0 74. B(C-0-C)(trans) 112.3 112.
AE(cis-trans) 6.46 6. structural parameters
O(C-C-O)(trans) 108.3 108. Methanol AE(ec1ipsed-staggered) 1.03 0.
2.72b 111.8b
0 f 0.3" 3.5d
0.364,' 0.38f 106.2' 105.0' 104.1'
1.5 f 0. 7.01h
84 f 6' 111.7j 108.9k
1.06' a Reference 56. Reference 57. Reference 58. Reference 59. e Reference 60. fReference 61. g Reference 62. hAb initio MF'2/6-31G/ /HF/6-31G calculations. Reference 63. Reference 64. Reference
experiment. The parameters in MM3 were derived by fitting to a wide variety of data for hydrocarbons, whereas our approach is to start with ethane as the simplest model and add additional dihedral parameters in a conservative way. As one can see, the barriers and geometry of n-butane are well described with such a model, as is the energy to eclipse the first and second methyl group of propane with the methylene.
We next tum to the alcohols and ethers (Table 4 ). Here we
begin with only two general V3 dihedrals, as in Weiner et al.,5X for X-CT-OH-X and X-CT-OS-X. This leads to es- sentially exact reproductions of the dihedral barriers in methanol and dimethyl ether. The cis-trans energy difference is about 0.5 k c d m o l greater than that calculated by the Weiner et al. force field; however, the Weiner et al. value matched the ex- perimental data originally used. When these dihedral parameters are applied to methyl ethyl ether (MEE) and tetrahydrofuran (THF), one finds that a small Vz(CT-CT-OS-CT) dihedral of 0.1 k c d m o l (Weiner et al. had such a parameter with mag- nitude 0.2 kcal/mol) leads to an excellent reproduction of the g/t energy difference in MEE and a slight preference for C
(56) Allinger, N. L.; Rahman, M.; Lii, J.-H. J. Am. Chem. SOC. 1990,
(57) Blukis, U.; Kasei, P. H.; Myers, R. J. J. Chem. Phys. 1963, 38, ( 5 8 ) Almenningen, A,; Seip, H. M.; Willadsen, T. Acta Chem. Scand. (59) Engerholm, G. G.; Luntz, A. C.; Gwinn, W. D.; Harris, D. 0. J.
(60) Cremer, D.; Pople, J. J. Am. Chem. SOC. 1975, 97, 1354-1358. (61) Geise, H.; Adams, W.; Bartell, L. Tetrahedron 1969, 25, 3045-
(62) Kitagawa, T.; Miyazawa, T. Bull. Chem. SOC. Jpn. 1968, 41(8), 1976-1976; MP2/6-31G//HF/6-31G ab initio calculations lead to AE = 1.4 kcal/mol. (63) Oyanagi, K.; Kutchitsu, K. Bull. Chem. SOC. Jpn. 1978,51,2237-
(64) Hayashi, M.; Adachi, M. J. Mol. Struct. 1982, 78, 53-62. (65) Lees, R. M.; Baker, J. G. J. Chem. Phys. 1968, 48, 5299-5318.
Chem. Phys. 1969, 50, 2446-2457.
Table 5. Low-Frequency Vibrational Modes
Dimethyl Phosphate Energies, Structures, and
Relative Energies" (kcaVmo1) conformationb E(MM) E(QM) 0.00 0. 1.42 1. 2.83 3. Geometrical Parameters (Angles in dedb.' (^) ~~ MM QM X-rayd $%42z(g,g) 67.7, 67.7 75.2, 75.2 73. 4Jlr42(grt)
e(opo')(g,g) 108.
e(cop)(g,g) 122.1 118.5 121. e(opo)(g,g) 103.8 99.3 104. 107.5 110. e(o'po')(g,g) 119.3 124.9 119. e(cop)(g,t) 120.5 118. e(oPO)(g,t) 102.5 96. e(o'po)(g,t) 108.2 108. e(o'Po')(g,t) 120.1 122. e(cop)(t,t) 120.2 116. 103.0 94. 108.2 109. 119.9 120. Vibrational Frequencies < 500 cm-' (cm-I)
e(opo')(t,t) e(o'Po)(t,t) e(o'po')(t,t)
MM exp' MM expe MM expe 78 262 210 359 357 109 295 321 383 393 196 195 302 345 42 1 503 a Absolute energies for g,g conformations are -40.77 kcaVmol (MM) and -720.606019 au (QM). The quantum mechanical calculations used the model MP2//6-31G//HF/6-31G. Dihedral angles around C-O- P-0. Bond angles, 0 is ester oxygen and 0' is anionic oxygen. See Table 4 in ref 5. e Reference 66.
THF over C,, as inferred from experiments. The calculations
overestimate the barrier to planarity of THF, but not by as much
as MM3. We next turn to dimethyl phosphate, the model for the backbone of nucleic acids. We have carried out ab initio calculations (MP2/6-31G//HF/6-31G) on dimethyl phosphate in its g,g; g,t; and t,t conformations and adjusted the Vz(0S- P-OS-CT) parameter to reproduce the (g,g)/(g,t) energy difference of 1.41 kcdmol. These results are reported in Table 5. The reoptimized VZparameter has a value of 1.20 as opposed to the value of 0.75 determined by Weiner et al. with the V parameter of 0.25 left unchanged. Reasonable agreement with ab initio calculations and consensus structural values from X-ray data has been achieved. The normal mode frequencies calcu- lated with such a model are also compared with those developed based on experimental frequencies of diethyl phosphate.& Given the difference in molecules, the agreement between calculation and experiment for the low-frequency modes reported in Table 5 is acceptable. The low-frequency modes for the simple hydrocarbons, alcohols, ethers, and thio compounds are presented in Table 6. The average error between the calculated and experimental frequencies is 31 cm-' for the 36 low-frequency examples where experimental data are available, compared to an error of 21 cm-' with MM3. Again, it should be noted that our parameters have been optimized using this limited set of simple molecules whereas the test set of molecules used to derive the MM parameters is much larger. Next to consider in the development of a force field for nucleic acids are the bases. Elsewhere, we have reported the (66) Brown, E.; Peticolas, W. Biopolymers 1975, 14, 1259-1271. (67) Allinger, N. L.; Quinn, M.; Rahman, M.; Chen, K. J. Phys. Org.
(68) Allinger, N. L.; Quinn, M.; Rahman, M.; Chen, K. J. Phys. Org.
Chem. 1991, 4 , 647-658.
Chem. 1991, 4, 659-666.
5186 J. Am. Chem. SOC., Vol. 117, No. 19, 1995
Table 6. Low-Frequency (< 1000 cm-I) Vibrational Modes for Small Hydrocarbons, Ethers, Alcohols, and Sulfur Compounds
symm @(thiswork) O(MM3). @(exp)b,c moded trans-NMA
Come11 et al.
Table 7. Normal Modes of trans-NMA and Benzene (cm-I) nmode no. symm this work experimenta mode
312 811 811 898
23 1 275 356 733 809
866
877
127 236 27 1 272
364
297 867
212 279 416 798
123 225
27 1
283
404
755
806
707 80 1
279 69 1 720
105 236 275 509 710 713
Ethane 279 283 CH3-CH3 torsion 908 822 CH3 asym rocking 908 822 CH3 asym rocking 962 995 C-C stretch
208 217 CH3-CH2 torsion 255 265 CH3-CH2 torsion 375 379 C-C-C bend
803 748 CH2 rock + CH3 def
850 868 CH3 rock + sym C-C str/str
938 921 CH3 rock +
asym C-C str/str
961 899 CH2 twist + CH3 def
122 121 CHI-CH~ torsion 216 CH3-CH2 torsion 245 266 CH3-CH2 torsion
287 asym C-C-C bend +
C-C -C bend
394 427 sym C-C-C bend +
C-C-C bend
Propane
Butane
Methanol 263 270 CH3-0 torsion 1052 1034 C-0 stretch Dimethyl Ether 188 198 CH3-O sym torsion 273 242 CH3-0 asym torsion 400 424 C - 0 - C bend 924 918 C-0 sym stretch Methyl Ethyl Ether 114 C2H5-0 torsion
216 CH3-C torsion +
CH3-O torsion
257 238 CH3-0 torsion +
CH3-C torsion
296 308 C-0-C bend +
C-C-0 bend
420 472 C-C-0 bend +
C-0-C bend
870 820 CH3 rock + CH2 rock +
CH2 twist 897 855 C - 0 str + CH3 wag + C-C str Methanethiol 695 704 c-s 823 803 C-S-H Dimethyl Sulfide 285 282(285) C-S-C 683 691 (683) S-C sym 702 741 (704) S-C asym Dimethyl Disulfide 116 102 (106) C-S-S-C torsion 241 239 (242) S-S-C bend 279 272 S-S-C bend
1 A” 44 2 A” 97 3 A“ 184 192 4 A’ 286 289 5 A’ 440 439 6 A” 587 600 7 A’ 591 628 8 A” 696 725 9 A’ 80 1 883 10 A’ 963 99 1 11 A” 1037 1044 12 A” 1046 13 A’ 1075 1114 14 A‘ 1082 1161 15 A‘ 1209 1300 16 A‘ 1395 1374 17 A’ 1398 1414 18 A” 1402 1441 19 A” 1407 1451 20 A’ 1428 1458 21 A’ 1516 1471 22 A’ 1614 1569 23 A’ 1693 1660 24 A’ 2868 2935 25 A’ 2869 2935 26 A” 2980 298 1 27 A’ 2982 298 1 28 A’ 2982 2994 29 A’ 2983 2994 30 A’ 3304 3307
1 e2u 410 410 2 e2g 609 606 3 a2u 66 1 673 4 b2g 704 703 5 e l g 900 849 6 e2u 979 975 7 alg 941 992
9 blu 1167 1010 10 elu 1124 1038 11 b2u 1194 1150 12 e2g 1129 1178 13 b2u 1331 1310 14 a2g 1729 1326 15 elu 1493 1486 16 e2g 1706 1596
18 alg 3062 3062 19 elu 3064 3063 20 blu 3068 3068
Benzene
8 b2g 947 995
17 e2g 3064 3047
ring def ring def CH bend ring def CH bend CH bend ring stretch (breathing) CH bend ring def CH bend CH bend CH bend ring stretch (kekule) CH bend
ring stretch + def
ring stretch CH stretch CH stretch CH stretch CH stretch a Reference 70 for trans-NMA. Reference 71 for benzene.
the bases is the dihedral potential for out-of-plane motion, as discussed by Weiner et al. As in the development of our previ- ous force field, normal mode analyses of benzene and NMA are imuortant. The results for the normal mode analyses auulied 514 509 (514) S-S stretch (^) to thesk molecules are presented in Table 7. We ha& readjisted 701 689 (694) S-C stretch 703 (694) S-C stretch a References 2, 56,67, and 68. See references 2,5,56,67, and 68 for experimental frequencies. Experimental frequencies given in parentheses refer to those used as reference for MM3 values. See references 2, 5, 56, 67, and 68 for’the mode assignments.
hydrogen bond energies and structures of A:T and G:C pairs and these appear to be in good agreement with the highest level of ab initio data currently available.69 However, a critical element in the development of planar functionalities such as
(69) Gould, I. R.; Kollman, P. A. J. Am. Chem. SOC. 1994, 116, 2493-
the X-CA-CA-X V 2 value and the improper out-ofiplane dihedral X-X-CA-HA to ensure correct representation of the lowest frequency modes of benzene, with the four lowest modes (1700 cm-’) in good agreement with experiment.” We next turn to NMA, the model for the peptide backbone. With a few adjustments to the Weiner et al.53b bonded parameters, the agreement with experiment70for the six lowest frequency modes is again excellent. In NMA, a key adjustment (70) Rey-Lafon, M.; Ford, M. T.; Garrigen-Lagrange, C. Specrrochim. Acta, Part A 1973, 29A, 471-486. (7 1) Shimananchi, T. Tables of Molecular Vibrational Frequencies; National Stand. Ref. Data Ser.; National Bureau of Standards: Washington, DC, 1967; Parts^ 1-3.
5188 J. Am. Chem. Soc., Vol. 117, No. 19, 1995
Table 9. Conformational Energies for Deoxyadenosine (Angles in deg, Energies in kcal/mol)”
pucker qa wd y b zb 3‘OHC SOHd Ee A d Sugar Pucker Profile E = 18 C2’endo 0.40 146.1 51.3 -158.4 176.5 171.4 -52.40 0 C3’endo 0.37 5.7 56.9 -162.3 -178.3 -179.3 -51.87 0. 04’endo 0.38 65.3 54.1 -156.7 -175.5 -175.4 -49.53 2. O4’exo 0.29 276.6 42.6 -178.7 175.8 178.5 -46.54 5.
E = 48 C2’endo 0.39 144.9 55.3 -153.1 176.2 -179.2 -1.65 0 C3’endo 0.38 14.1 56.3 -156.1 178.1 179.9 -0.61 1. 04’endo 0.39 56.6 55.6 -153.1 179.1 179.5 0.21 1. 04’exo 0.30 285.0 42.7 176.9 177.3 179.8 4.03 5.
Gamma Dependence E = 18 C2’endo 0.40 146.1 51.3 -158.4 176.5 171.4 -52.40 0 C2’endo 0.42 141.9 -168.6 -168.6 179.6 -179.9 -50.39 2. C2’endo 0.42 151.7 -62.3 -169.9 -179.7 -179.7 -50.92 1.
E = 48 C2’endo 0.39 144.9 55.5 -153.1 176.2 179.2 -1.65 0 C2’endo 0.40 148.0 -172.9 -162.4 180.0 -179.9 -1.31 0. C2’endo 0.40 150.0 -66.6 -169.9 -179.9 -179.9 -0.13 1.
x Dependence E = 18 C2’endo 0.40 146.1 51.3 -158.4 176.5 171.4 -47.46 0 C2’endo 0.40 166.7 63.3 60.8 -179.1 179.8 -41.52 5.
E = 48 C2’endo 0.39 145.4 55.3 -141.6 176.5 179.2 3.24 0 C2’endo 0.40 144.0 52.6 37.6 180.0 179.6 1.84 -1.
a q and W defined in ref 76. y and x defined in ref 77. Above, the
first entry corresponds to an “anti” conformation for x, the second to
“syn”. 3’OH refers to C4’-C3’-03’-H03’ dihedral. 5’OH refers to H05’-05’-C5’-C4’ dihedral. e Absolute molecular mechanical energy. f Relative conformational energy. g E is the dielectric screening factor used in eq 1.
Cornel1 et al.
*H61* N H 6 2 N
HE-CE
H04’-04’-Cl ‘-HI ’
HZ’l-CZt-H2’
H2’ Figure 1. Model of deoxyadenosine employed in the quantum mech- anical and molecular mechanical conformational studies reported in Table 10. In the quantum mechanical calculations, the H04’-04’- Cl’-N9 and H2’3-C2’-Cl’-N9 dihedrals were held fixed at values characteristic of a C2’-endo sugar, in order to mimic the conformation of the sugar ring. In the molecular mechanical calculations, the dihedrals were restrained to those values with dihedral restraints of 500 k c d mol.
agreement with high-level ab initio data is very good for all but alanyl dipeptide C7ax and glycyl dipeptide a R. The ala C,, conformation is rarely found in proteins and gly occurs relatively infrequently in a-helices, due to the loss of conformational entropy, so these conformations were reasonable ones in which to tolerate any error.
(80) Ben Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016-2027. (81) Wolfenden, R. Biochemisrry 1978, 17, 201-204.
Table 10. x Angle Profile for Base with Sugar Fragment (kcdmol) AMBER (E = 1) ab initio no specific with specific pb MP2/6-31G//HF/6-31G dihedral dihedral‘
60 min 120 180 min 210 240 300 360
60 min 120 180 min 210 240 300 360
Model of Deoxyadenosine 0.94 4. 0.63 (74.7°)d 4.62 (61.1°)d 3.37 5. 0.06 0. 0.00 (198.2°)d 0.00 (196.7°)d 0.22 0. 1.45 1. 5.33 6.
2.27 6. 2.02 (72.3°)d 5.99 (61.0°)d 7.02 8. 0.74 1. 0.00 (210.00)d 0.00 (205.2°)d 0.00 0. 1.29 1. 8.15 8.
Model of Deoxythymidine
0.00 (197.2°)d
1.45 (54.50)d
2.83 (55.4°)d
0.00 (205.5°)d
Reference 77. Degrees. Specific VI and V2 dihedral terms were added for OS-CT-N-CK (purines) and OS-CT-N-CM (pyri- midines) dihedral angles. Minimized value of x.
Table 11. Conformational Energies of Glycyl and Alanyl Dipeptides (kcdmol) glycyl dipeptide alanyl dipeptide E(MM) UQMY E(MM) E(QM)” C I 0.0 0.0 c 7 q 0.0 0. c 5 1.9 2.0 c7, 1.5 2. aR 6.0 4.0 c 5 1.5 1. aR 3.9 3. Quantum mechanical energies calculated at the MP2/TZP//HF/6- 31G* level on methyl-blocked versions of the dipeptides. See ref 28 for further details.
One of the important features in our force field is the attempt to reproduce the solvation free energies of a representative set of molecules. In Table 12, we present such a representative set. As one can see, the absolute solvation free energy of methane is somewhat (0.5 kcdmol) too large with our model, but the relative solvation free energies of methane, ethane, and propane are within 0.3 kcdmol of experiment. For our protypal polar molecules, methanol and NMA, the agreement with experimental solvation free energies is within -0.5 kcdmol. We wished also to assess the solvation free energies for sulfur compounds and the relative solvation free energies of those are in reasonable agreement with experiment (again within 0.5 kcaV mol). The calculated free energy of 9-methyladenine is a
prediction, because there are no precise experiment^,^^ but the
relative free energies of this force field and that of Weiner et al.536 suggest that the experimental determination of this quantity would be of great interest. Turning to the ionic molecules, our results make clear that a typical two-body additive force field will tend to overestimate ion solvation (when corrected for long- range cut-off) unless its parameters are significantly modified, but fully non-additive calculations with exactly the same parameters reproduce experiment very well. (82) Meng, E.; Cieplak, P.; Caldwell, J.; Kollman, P. Accurate Solvation Free Energies of Acetate and MethylammoniumCalculated with a Polarized Water Model. J. Am. Chem. Soc., in press. (83) Analyzed by Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1991, 113, 8305-8311. (84) Hine, J.; Mookejee, P. K. J. Org. Chem. 1975, 40, 292-298. (85) Ferguson, D. M.; Pearlman, D. A.; Swope, W. C.; Kollman, P. A. J. Comput. Chem. 1992, 13, 362-370.
Simulation of Proteins and Nucleic Acids J. Am. Chem. SOC., Vol. 117, No. 19, 1995 5189
0 PSI
PHI
0 PSI
PHI
Figure 2. (a) The molecular mechanical (4.v) map for methyl-blocked
glycyl dipeptide generated using the force field presented here. Contours are drawn every 2 kcdmol. (b) The molecular mechanical (4,q) map for methyl-blocked alanyl dipeptide generated using the force field presented here. Contours are drawn every 2 kcal/mol.
Table 12. (kcal/mol’,
Solvation Free Energies for Model Compounds
molecule AAG(ca1c) AAG(exp) C& - nothing
C3Hs - C2H
NMA - CHq CH3NH3+ - nothing
CH3C02- - nothing
C2H6 - c&
CH30H - CH3CH
CH3SCH3 - CH3OCH CH30H A CHsSH 9-CH3 adenine - CHq
-2.5 f 0.1‘ -0.1 f 0.1‘ -0.2 f 0.1c 6.9 f O.ld 11.6 f 0.2d 87.6 f 2. (75.4 f 1.7)’ 87.1 f 1. (71.6 f 1.0)’ 0.9 f 0.1” 3.5 f 0.1” 18.3 f 2.6, 13.9 f 0.4’
**-2.P
-0. 6.9* 12.1e 77-**
70-71h
0.4‘ 3.7‘ 15.6 2r 1.1) a Reference 21b. Reference 80. Reference 21a. Because of the uncertainty in the electrostatic potential derived charges for ethane, the average of^ the free energies for the electrostatic potential derived and Mulliken charges for these free energy calculations are presented. Reference 20. Reference 81. f Reference 82, additive potential, values in parentheses are for nonadditive potential. 8 Reference 83. This paper. Reference 84. j See Note Added in Proof.
The results described above were obtained on model systems that were relatively very simple3’ (neat liquids) andor small (dipeptides and nucleosides). In order to test the performance of the new force field on a more complex system, we carried
1
0. 0 60.0 120.0 180. time [ps] Figure 3. RMS deviation (A) between the crystal structure of ubiquitin and structures along an MD trajectory as modeled by the Weiner et aL6 and Cornell et al. (this work) force fields. The lower lines correspond to the RMS deviation of the heavy backbone atoms only and the upper lines to the RMS deviation for all heavy atoms.
out an MD simulation of ubiquitin in water with periodic boundary conditions. The RMS difference was calculated for structures along the trajectory relative to the crystal structures for (1) the backbone atoms and (2) all of the heavy atoms. These results were then compared to those obtained with the Weiner et ~ 1. force field (Figure ~ 3 ~ 3). The RMS values are reported for the first 72 residues only, since the four residues of the carboxy terminus were mobile. The behavior of the new force field presented here is better in two ways. First, the protein structure seems to have stabilized after 50 ps of simulation with the new force field, while the RMS deviation continues to increase throughout the trajectory with the Weiner et al.536 force field. Second, the RMS deviation for all of the heavy atoms
after 180 ps of simulation is about 2.0 8, with the force field
presented here and about 2.5 8, with the Weiner et al.5.6 force
field. Alonso and Daggett have also reported the results of a long MD simulation of ubiquitin, and they found a backbone
RMS deviation of 1.4 8, from the crystal structure, comparable
to the deviation found here.87 A referee has pointed out that smaller deviation from a crystal structure could simply be a consequence of an “unrealistically stiff’ force field. We cannot rule this out, but stress that we did not, in our force field derivation on the fragments described above attempt to add “stiffness”. Even closer agreement with a protein crystal structure has been obtained by York et al.,” who carried out a 1000 ps MD simulation of BPTI with the long-range electrostatic forces of the crystal environment treated using the particle mesh Ewald method and the Weiner et ~ 1. force field. ~ 5 ~ With this model they obtained an RMS deviation from the crystal structure of
0.33 8, for backbone atoms. These results serve to illustrate
the difference between errors arising from the force field itself and those arising from its implementation in a given calculation.
Currently, most MD simulations employ an 8 or 9 8, cutoff for
nonbonded interactions in order to reduce this rate-determining part of the calculation. In systems where long-range electrostat- ics play an important role, this approximation is clearly inadequate. Although the Ewald method is only fully appropri- ate for periodic crystal systems, other methods also exist which allow for the more accurate treatment of long-range electrostat- (86) Vijay-Kumar, S.; Bugg, C. E.; Cook, W. J. J. Mol. B i d. 1987,194, 531-544. (87) Alonso, D. 0. V.; Daggett, V. Molecular Dynamics Studies of Partially Unfolded Conformations of Ubiquitin in Methanol and Their Refolding in Water, submitted for publication.
Simulation of Proteins and Nucleic Acids
0.0978 0.
J. Am. Chem. SOC., Vol. 117, No. 19, 1995 5191
0.2719 H H1 0-0.
- -0.4157^ N-CT-C-I^ I^ I^ 0. I -0'
0.2719 H Hl 0-0.
- -0.41 57^ N-CT-C-I^1 I^ 0.
I -o'w
0.1 123 0.2719 H h1 0-0.
HP 0 -0. 0.1984 M \ 0.1592 1 10, 0.1984 H- N3 - CT - C- 0.1984 H (^) ' I o'
-^ -0.4157^ N-CT-C-I^ I^ i^ 0. 0.0295 HC - CT - HC 0.
0.1330HA
0.0295 HC - CT - HC 0.
0.1699 HA
0.0125 HC-CT-HC 0, I 0'
I 0'
I
I -0'
0.0292 H1- CT - H1 0.
S-0.
0.0597 H1-CT-Hl 0.
-0.1704CA CA -0.1 704 0.1430HA '\CAY \ H A 0.
I -0'
-0.2341CA CA -0 2341 0.1656 HA / \ c / \ H A 0 1 6 5 6
I 0'3226 HA 0.1572I H 0.
HA 0.1 297 -0.5579 .OH HO 0.3992 Hl 0. PHE TYR TRP MET-nt
0.2681 H H -0,382, I I /02 -0'
- N-CT- I - 0. 2 ~ 9 : ~ ' ~ 02 -0. H1 0.
GLY-ct
Figure 4. Charges for the peptide fragments. All the charges for the non-terminal amino acids are presented. For histidine, the three protonation states are presented (HID, HIE,HIP). For cysteine, both the disulfide bonded forms (CYX) and the reduced form (CYS) are presented. The N terminal and C terminal blocking groups are presented (ACE and NME, respectively). The N and C terminal amino acids for ubiquitin (MET-nt and GLY-ct) are shown; the remaining N and C terminal residues are available by anonymous ftp as are the protonated forms of GLU and ASP and the deprotonated form of LYS. See ref 39 for a description of how these charges were derived.
A B \ / H 0.
0.4167 H N2-0.
0.1877 H5 \ -0.6175 / C IA Y 0.1607 CK 0.0725 CB NC -0. I I I
0.4314 H N2-0.
0.1 963HA CA 0.8439^ I d , 5 k C r ' 'NC-0. I^ I -0.0268 N-* ' CB^ I CQ 0. 0'3800 \NC/ \H5 0. -0.
ADE
-0.0183& 0. 0.2293H4 ' \y-i39'0 -0.
0.0754 H 1. -0.0069 Base 0.0754 *HINCT* 7; / 0.1176 H1^ ,CTO.l629 H ,CT-CT^ \0.0713^ (0, 0.0985 H
C M 0 -0. -0.5725 C I 0,4918 H 0. 0.1997 H5 \ / \ / 0.0736 c,K 0.1 991 c r NA-0.
I
I HC 0. 0.0577 pi' - 0.1814 CB (^) LNC/CA 0.7432\ /H 0. N2-0.
7 - H 2 0.1746 H^ I 0.
(or\ 1 o,0358 -0. OH -0. \ HO 0.
D
-0'5232 O r C
I GLIA
THY
Figure 5. Charges for DNA. The four bases and the C1' and H 1 ' charges are shown separately. These are combined with any of four combinations of sugarhackbone charges. A nucleoside corresponds to fragments A and D with the sugar. A 5' terminal residue corresponds to fragments A and C and the sugar; a 3' terminal residue corresponds to combining B and D with the sugar and a central residue corresponds to combining B and C with the sugar. See ref 39 for how these charges were derived.
ics.88 Thus, it appears that the way electrostatic interactions are handled is significantly more important than the detailed force field parameters in ensuring that a molecular dynamics trajectory stays near an experimental (X-ray or NMR) structure. We suggest, however, that comparing two force fields with the same cutoff protocol can be illustrative and we conclude, on that basis, that the new force field performs at least as well as, or slightly better than, that of Weiner et ~ 1. for full solution ~ 3 ~ simulations.
Discussion
We have presented the development and the description of a new force field for proteins, nucleic acids, and organic molecules. Previously, we have attempted to give a coherent description of the underlying basis for the Weiner et al. force field,5s6 in order that it could be extended by others as well as ourselves for studies of molecular interactions and conforma- tions. We should emphasize again that our goal is to describe molecular conformational energies and structures as accurately (88) Saito, M. J. Chem. Phys. 1994, 101, 4055-4061. as possible in condensed phases with a simple, transferable, and
5192 J. Am. Chem. SOC., Vol. 117, No. 19, 1995 Comell et al.
A B
-0.7760 02 - P^ I^ l.l -^662 02 -0.
I
0.4295 HO
-0.6223^ \ OH^ OS^ -0.
0.0679 H1.CT0.0558 Base
0.0679 H 1 ’ \ 7; /
TC , **0.
- 1 1 7 4 H l** \
\cTo.z02z 10.
00615 H 1 ’ I - T L H l 0.
(^1) -0.61. 9 OH 3 HO 0.41 86
-0.0251 N- CB* CQ 0. 0’3053 \NC’ \H50.
(J‘ \cT0,0394l H 2 0.
ADE 0 -0.
4.1709 I
0.1640 H 5 , ,C.y77o,H 0. o,l 374 c,K 0.1 744c,@ NA -0.
I
0.0492 N’ - 0.1222CB \ N C / CA 0.7657 \ (^) /H 0. N24.
(M \ c ~ 0. 0 1 9 1 ).^ H2 0.2006 -0’6323 H^ I 0.
I GUA
0.4234 H N2-0.
0.1928 HA C-^ I^ 0.
-o,5kC,M’ \ N C 4. 7 S 8 4
0.0053 CM C^ I 0.
0.1958 H4^ /^ \ * / \
CYT
0 -0.
C D
Figure 6. Charges for RNA. Notation is the same as in Figure 5.
Table 13. Comparison of Cornell et aL, Weiner et aL, CHARMm, OPLSIAMBER, and GROMOS Force Fields
C^ I 0.5952 H 0.31 54
NA-0.
\ / \ /
0.1811 HA
-0’3635 -0.1 126 ‘I*CM cI 0.
0.2188 H 4 ’ \yA18‘o -0.
URA
force field electrostatics van der Waals VDW combining rules“ torsions CHARMm93(1983) empirical fit to QM dimers empirical (x-tals) _R_* arithmetic mean; E geometric mean single bond pathb GROMOS92 empirical empirical (x-tals) A and B “non-standard” geometric meanC user specified OPLS/AMBERIS(1990) empirical (MC on liqs) empirical (liquids) A and B geometric means equal division among
Weiner et a1.s,6 ESP fit (STO-3G) empirical (x-tals) R* arithmetic mean; E geometric mean equal division among
this work E S P fit (6-31G) empirical (liquids) _R_ arithmetic mean: E geometric mean equal division among
equiv bond paths
equiv bond paths
equiv bond paths A = E R * ’ ~ and B = ~ E R * ~. In CHARMm22, the torsion representation was changed to the more commonly used equal division of the energy along equivalent bond paths. GROMOS employs the geometric mean method for calculating VDW interactions, but for water-methyl interactions, for example, a smaller VDW radius is assumed for the water since it is no longer in a hydrogen bonding interaction. This has been shown to result in a “too hydrophilic” methyl g r o ~ p. ~ ~ , ~ ~
general model. This goal has framed our approach, which has been to focus mainly on the electrostatic, VDW, and dihedral energies and use both ab initio calculations, empirical liquid and solvation data, and experiment to calibrate the model. However, our approach differs significantly from that of many in building from the ground up with the simplest model and defining relatively few general principles, which are elucidated in the section General Description of the Model above. We will attempt to summarize the salient features of some of the more commonly used force fields here, in order to compare and contrast our approach with theirs. They can be roughly grouped into four different categories, depending on the nature and complexity of the force field equation: (1) those with rigid or partially rigid geometries, (2) those without electrostatics, (3) simple diagonal force fields, and (4) more complex force fields. The ECEPP force field of Scheragas9 employs rigid internal geometries which allow a more efficient exploration of con-
formational space. This approach has the disadvantage that it
can cause certain conformations and conformational barriers to be too high in energy. A second force field which uses only partially rigid geometries is developed by Lavery
(89) (a) Roterman, I. K.; Lambert, M. H.; Gibson, K. D.; Scheraga, H. A. J. Biomol. Srruct. Dyn. 1989, 7, 421-453. (b) See also: Kollman, P. A.: Dill, K. A. 1991, 8, 1103- 1107. Gibson, K. D.; Scheraga, J. Biomol. Struct. Dyn. 1991, 8, 1109-1111. (90) Lavery, R.; Hartmann, B. Biophys. Chem. 1994, 50, 33-45.
and co-workers. This force field has been developed for nucleic acids and allows flexibility in the sugar ring but uses mainly internal geometries and keeps the bases rigid. The SYBYL force field9’ has been developed for the calculation of internal geometries and conformational energies. Because it contains no electrostatic term, it is inappropriate for studying detailed condensed-phase properties. The YETI force field?2 developed by Vedani and Huhta, is a modification of the Weiner et al. force field with highly damped electrostatics and an angular dependent hydrogen bond (and metal ligation) potential added. This approach could be valuable in some modeling situations, where large and difficult to handle elec- trostatic energies are present, but it is also unlikely to be general and extendable to condensed-phase phenomena. The category of simple diagonal force fields includes the Weiner et ~ l. , ~ , ~ GROMOS,93 CHARMm?4 and OPLS/AM- BERI5 force fields. All of these force fields employ a simple harmonic diagonal representation for the bond and angle terms. Descriptions of the nonbonded and dihedral energies are given in Table 13. The Weiner et al. force field derived charges from
(91) Clark, M.; Cramer, R. D.; van Oppenbosch, N. J. Comput. Chem. (92)Vedani, A,; Huhta, D. A. J. Am. Chem. SOC. 1990, 112, 4759-
(93) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen Molecular (94) Brooks, B. R.; Brucoleri, R. E.; Olafson, B. D.; Slater, D. J.;
1989, 10, 982-1012.
Simulations (GROMOS) Library Manual; Biomos: Groningen, 1987.
Swaminathan, S. ; Karplus, M. J. Compur. Chem. 1983, 4 , 187-217.
5194 J. Am. Chem. SOC., Vol. 117, No. 19, 1995
Table 14. Molecular Mechanical Parameters"
Bond Parameters
Come11 et al.
C-CA C-CB C-CM C-CT C-N C-N* C-NA C-NC c- c - 0 2 C-OH C-CB C-CT c * - c w C*-HC CA-CA CA-CB CA-CM CA-CN CA-CT CA-H
CA-HA CA-N CA-NA CA-NC CB-CB CB-CN CB-N* CB-NB CB-NC CC-CT cc-cv cc-cw CC-NA CC-NB CK-H CK-N* CK-NB CM-CM CM-CT CM-H CM-H
1.080 CM-HA 1.340 CM-N* 1.381 CN-NA 1.339 CQ-H 1.370 CQ-NC 1.419 CR-H 1.374 CR-NA 1.391 CR-NB 1.354 CT-CT 1.504 CT-F 1.375 CT-H 1.371 CT-H 1.385 CT-H 1.394 CT-HC 1.080 CT-HP 1.371 CT-N 1.304 CT-N* 1.350 CT-N 1.510 CT-N 1.080 CT-OH 1.080 CT-OS Angle Parameters
1.47 1
CT-S CT-SH CV-H CV-NB CW-H CW-NA H-N H-N* H-N H-N H-NA HO-OH HO-OS HS-SH 02-P OH-P o s - P o w - H w s-s
227.0 1. 237.0 1. 367.0 1. 410.0 1. 367.0 1. 427.0 1. 434.0 1. 434.0 1. 434.0 1. 434.0 1. 434.0 1. 553.0 0. 553.0 0. 274.0 1. 525.0 1. 230.0 1. 230.0 1. 553.0 0. 166.0 2.
" angle Ked @eo' angle Ksd @eae angle Ked angle Ked OeQe C-CA-CA C-CA-HA C-CB-CB C-CB-NB C-CM-CM C-CM-CT C-CM-H C-CM-HA C-CT-CT C-CT-H C-CT-HC C-CT-HP C-CT-N C-CT-N C-N-CT C-N-H C-N-CM C-N-CT C-N-H C-NA-C C-NA-CA C-NA-H C-NC-CA C-OH-HO C-CB-CA C-CB-CN C -CT-CT C * -CT-HC C * -CW -H C* -CW -NA CA-C-CA CA-C-OH CA-CA-CA CA-CA-CB CA-CA-CN CA-CA-CT CA-CA-H CA-CA-HA CA-CB-CB CA-CB-CN CA-CB-NB CA-CM-CM CA-CM-H CA-CM-HA CA-CN-CB CA-CN-NA
CA-CT-HC CA-N2-CT CA-N2-H CA-NA-H CA-NC-CB CA-NC-CQ CB-C-NA CB-C- CB-C-CT CB-C-CW CB-CA-H CB-CA-HA CB-CA-N CB-CA-NC CB-CB-N* CB-CB-NB CB-CB-NC CB-CN-NA CB-N-CK CB-N-CT CB-N-H CB -NB -CK CB -NC-CQ CC-CT-CT CC-CT-HC CC-CV-H CC-CV-NB CC -CW -H CC -CW -NA CC -NA-CR CC-NA-H CC-NB -CR CK-N-CT CK-N-H CM-C-NA CM-C- CM-CA-N CM-CA-NC CM-CM-CT CM-CM-H CM-CM-HA CM-CM-N CM-CT-HC CM-N* -CT CM-N*-H CN-CA-HA
CN-NA-H CR-NA-CW CR-NA-H CR-NB-CV CT-C-N CT-C- CT-C- CT-C-CW CT-CC-CV CT-CC-CW CT-CC-NA CT-CC-NB CT-CT-CT CT-CT-H 1 CT-CT-H CT-CT-HC CT-CT-HP CT-CT-N CT-CT-N CT-CT-N CT-CT-N CT-CT-OH CT-CT-OS CT-CT-S CT-CT-SH CT-N-CT CT-N-H CT-N2-H CT-N3-H CT-OH-HO CT-OS-CT CT-OS-P CT-S-CT CT-S-S CT-SH-HS CV -CC -NA CW -CC-NA CW-CC-NB CW-NA-H F-CT-F F-CT-H H-N-H H-N2-H H-N3-H HI-CT-H H1 -CT-N
CA-CT-CT 63.0 114.00 CN-NA-CW 70.0 111.60 HI-CT-N* 50.0 109.
H1 -CT-N H1 -CT-OH H1-CT-OS H1-CT-S H1-CT-SH H2-CT-H H2-CT-N* H2-CT-OS H4 - CM - N* H4-CV-NB H4-CW-NA H5-CK-N* H5 - CK - NB H5 -CQ-NC H5-CR-NA H5 -CR-NB HC - CT - HC HO-OH-P HP-CT-HP HP-CT-N HS-SH-HS HW-OW-HW N-C- N-C-NA N-C-NC N-C- N -CB-NC N-CK-NB N-CT-OS N2-CA-N N2-CA-NA N2-CA-NC NA-C- NA-CA-NC NA-CR-NA NA-CR-NB NC-C- NC-CQ-NC 0-c- 0 2 - c - 0 2 02-P- 02-P-OH 02-P-os OH-P-OS o s - P - o s P- o s - P
50.0 109. 50.0 109. 50.0 109. 50.0 109. 50.0 109. 35.0 109. 50.0 109. 50.0 109. 35.0 119. 35.0 120. 35.0 120. 35.0 123. 35.0 123. 35.0 115. 35.0 120. 35.0 120. 35.0 109. 45.0 108. 35.0 109. 50.0 109. 35.0 92. 100.0 104. 80.0 122. 70.0 115. 70.0 118. 80.0 120. 70.0 126. 70.0 113. 50.0 109. 70.0 120. 70.0 116. 70.0 119. 80.0 120. 70.0 123. 70.0 120. 70.0 120. 80.0 122. 70.0 129. 80.0 126. 80.0 126. 140.0 119. 45.0 108. 100.0 108. 45.0 102. 45.0 102. 100.0 120.
Simulation of Proteins and Nucleic Acids
Table 14 (Continued)
J. Am. Chem. Soc., Vol. 117, No. 19, 1995 5195
Torsional Parameters torsion no. of pathd VJ2g Y h n' torsion no. of pathd' V,,l2g Y h n' X-C-CA-X X-C-CB-X X-C-CM-X X-C-CT-X X-C-N-X X-C-N-X X-C-NA-X X-C-NC-X X-C-OH-X X-C -CB -X X-C-CT-X x-c -cw-x X-CA-CA-X X-CA-CB-X X-CA-CM-X X-CA-CN-X X-CA-CT-X X-CA-N2-X X-CA-NA-X X-CA-NC-X X-CB-CB-X X-CB-CN-X X-CB-N-X X-CB-NB-X X-CB-NC-X X-CC-CT-X x-cc-cv-x x-cc-cw-x X-CC-NA-X X-CC-NB-X X-CK-N-X X-CK-NB-X X-CM-CM-X X-CM-CT-X X-CM-N-X X-CN-NA-X X-CQ-NC-X X-CR-NA-X X-CR-NB-X X-CT-CT-X X-CT-N-X X-CT-N-X X-CT-N2-X X-CT-N3-X
X-CT-OH-X
X-CT-OS-X
X-CT-S-X
X-CT-SH-X
X-CV-NB-X
X-CW-NA-X
X-OH-P-X
x-os-P-x C-N-CT-C C-N-CT-C C-N-CT-C C-N-CT-C CT-CT-C-N CT-CT-C-N CT-CT-C-N CT-CT- C-N CT-CT-N-C CT-CT-N-C CT-CT-N-C CT-CT-N-C CT-CT-OS -CT CT-CT-OS-CT CT- S -S -CT CT-S -S -CT H-N-C- H-N-C- N-CT-C-N N-CT-C-N N-CT-C-N N-CT-C-N OH-CT-CT-OH OH-CT-CT-OH OH-P-OS-CT OH-P-OS-CT OS-CT-CT-OH OS-CT-CT-OH OS-CT-CT-OS OS-CT-CT-OS OS-CT-N-CK OS-CT-N-CK OS-CT-N-CM OS-P-OS-CT OS-P-OS-CT S -CT-N-CM
Improper Torsions
1 .OO
1 .oo
1 .oo
1 .oo
-4. -3. -2.
-4. -3. -2. 1 .o -4. -3. -2. 1 .o -3.
-2. -2. 1 .o -4. -3. -2. 1 .o -3.
-3.
-3.
-3.
-2. 1 .o -2. -3.
1 .o
torsion VJ29 n' torsion X-CT-N-CT X-N2-CA-N x-02-c- x-x-c- X-X-CA-H X-X-CA-H X-X-CA-HA X-X-CK-H X-X-CM-H X-X-CM-HA
1 .o
X-X-CQ-H
X-X-CR-H
X-X-CV-H
X-X-CW-H
X-X-N-H
X-X-N2-H
X-X-NA-H
CA-CA-C-OH
CA-CA-CA-CT
CB-NC-CA-N
VJ
1 .o 1 .o 1 .o
~~ CK-CB-N-CT CM-C-CM-CT CM-C-N-CT CT-CM-CM-C CW-CB -C*-CT NC-CM-CA-N NA-CV-CC-CT NA-CW-CC-CT NA-NC-CA-N NB-CW-CC-CT
VJ
1 .o
1 .o
Y h
ni
Van der Waals Parameters atom type _RJ_* €k atom type R*I ek atom type R*I €k atom type R*J €k C' CA CM c s CT F H H
0.0860 H
0.0860 H
0.0860 H
0.0000806 H
0.1094 HA
0.061 HC
0.0157 HO
0.0157 HP
HS
Hw IP K Li N" N3" 0
O.oo
OH
os ow P Rb S SH
- m 2.oo
Simulation of Proteins and Nucleic Acids
involve simulations up to 800 ps in length. On the other hand,
the second leg of the simulation, which involves the disappear- ance of the van der Waals interactions of the base atoms and changing the N-C bond to the H-C bond of methane, gives very different free energies for the two protocols. The calculated value is 6.16 kcdmol with SPASMS and 1.94 kcdmol with GIBBS, for simulations as long as 500 ps (with SPASMS) and 800 ps (with GIBBS). This leads to the values of 18.3 and 13.9 reported in Table 12. Using the calculated free energy of solvation of methane (2.5 kcdmol), this leads to an absolute
J. Am. Chem. SOC., Vol. 117, No. 19, I995 5197
solvation free energy for 9-methyladenineof 15.8 kcdmol with SPASMS and 11.4 kcdmol with GIBBS. Interestingly, these two values bracket the extrapolated experimental value of 13.
f 1.1 kcdmol (see Ferguson et al. (Ferguson, D. M.; Pearlman,
D. A.; Swope, W. C.; Kollman P. A. J. Comp. Chem. 1992,
13,362-372) for a discussion on how this "experimental" value was determined). Further calculations are required to sort out this issue.
JA943664F
