Energetics of Infinite Homopolypeptide Chains: A New Look at

6872

J. Phys. Chem. B 2008, 112, 6872–6877

Energetics of Infinite Homopolypeptide Chains: A New Look at Commonly Used Force Fields Evgeni Penev,† Joel Ireta,*,‡ and Joan-Emma Shea*,† Department of Chemistry and Biochemistry, and Department of Physics, UniVersity of California, Santa Barbara, California 93106, and Departamento de Quı´mica, DiVisio´n de Ciencias Ba´sicas e Ingenierı´a, UniVersidad Auto´noma Metropolitana-Iztapalapa, A.P. 55-534, Me´xico D.F. 09340 ReceiVed: January 4, 2008; ReVised Manuscript ReceiVed: March 4, 2008

We present a novel method for comparing the long-range part of force fields in the presence of a maximally cooperative network of nonbonded interactions. The method is based on mapping the potential energy surface of an infinite polypeptide chain in the gas phase by using cylindrical coordinates (the twist and pitch) as geometry descriptors. We apply our method to an infinite polyalanine chain and consider the AMBER99, AMBER99SB, CHARMM27, and OPLS-AA/L fixed partial-charge force fields and the protein-specific version of the AMOEBA polarizable force field. Results from our analysis are compared to those obtained from high-level density-functional theory (DFT) calculations. We find that all force fields produce stronger stabilization of the helical conformations as compared to DFT, with only AMBER99/AMBER99SB satisfactorily reproducing all three helical conformations (π, R, and 310). I. Introduction Proteins play a critical role in most cellular processes. To perform their biological function, proteins must first fold to a specific, stable three-dimensional structure. The large sizes of proteins, coupled with the relatively long time scales (∼µs, ms or greater)1 associated with the folding process, render their study extremely challenging from a computational standpoint. While the subtle interplay between strong (covalent) and weak (hydrogen bond, van der Waals, etc.) interactions that govern folding mechanism and protein stability calls for a first principles approach, a quantum mechanical treatment is in practice only possible for very small peptides. The current method of choice for studying protein folding and stability is through the use of molecular dynamics simulations,2 which employ empirical force fields3,4 parametrized from gas-phase experiments and quantum mechanical calculations. These force fields come in a variety of flavors,5–8 each with slightly different functional forms and parametrization. The assessment of the quality of these force fields is typically gauged by their ability to reproduce known experimental features for very small peptide systems. A standard test involves the calculation of the (φ, ψ) Ramachandran map for dipeptides. Although such studies are computationally tractable, the singleresidue energetics approximated by dipeptides may not provide meaningful information for larger systems where the conformational degrees of freedom of a given residue are coupled to the conformations of its nearest neighbors and/or residues far apart along the polypeptide chain. Even if a force field can reproduce the energetics of dipeptides, this does not imply that it is adequate for larger polypeptides in which long-range interactions are present. In particular, the “gold standard” benchmark of the alanine dipeptide is inherently unable to describe the network of hydrogen bonds present in helical motifs, the most common structural element in proteins.9 * Correspondingauthors.E-mail:[email protected]; [email protected] † University of California. ‡ Universidad Auto ´ noma Metropolitana.

The aim of the present work is to provide a new means of evaluating the ability of force fields to accurately describe the energetics of polypeptide chains. Our approach is based on the use of infinite homopolypeptide chains which, unlike the dipeptide system, incorporate the effects of a maximally cooperative network of nonbonded interactions. We focus here on the infinite right-handed Ala polypeptide chain (abbreviated hereafter Ala∞) in the gas phase as a model system. The energetics of this system have been previously studied using density-functional theory (DFT) by one of the authors of this paper, as well as by other researchers.10–13 In particular, Ireta et al.11,14,15 showed that the use of an infinite chain allows the characterization of a complete pathway (a “continuum” of conformations) on its energy landscape connecting the unfolded polypeptide chain to the folded, helical states, i.e., the process of hydrogen bond network development. They11,14 identified the helical twist θ and pitch L, shown in Figure 1a, as suitable descriptors of the polypeptide chain conformation in cylindrical coordinates, and mapped its complete potential energy surface (PES) in the (L, θ) space. The resulting PES accounts for the allowed conformationssnamely, π-helix, R-helix, 310-helix, 27 structure, and the fully extended structure (FES), Figure 1sof a residue embedded in a long chain under the influence of both longrange and nearest-neighbor interactions. It can thus be considered as an effective Ramachandran map in cylindrical coordinates. Building on the work of Ireta et al., we aim at gaining new insight into the ability of commonly used force fields to describe polypeptide energetics. We wish to (i) compare the description of van der Waals and long-range electrostatic interactions by different force fields and (ii) compare the available DFT PES of Ala∞ to that obtained from both fixed-charge and polarizable force fields. We note that protein force fields are generally intended to mimic behavior in aqueous solution. In this work, we do not include a model for the solvent environment as this would entail further uncertainty to the force field comparison, given that force fields and solvent models tend to be parametrized independently. Indeed, subtle changes to the solvent parametrization may strongly affect the relative stabilities of

10.1021/jp800058f CCC: $40.75  2008 American Chemical Society Published on Web 05/14/2008

Energetics of Infinite Homopolypeptide Chains

J. Phys. Chem. B, Vol. 112, No. 22, 2008 6873 assuming the helical axis oriented along eˆz (cf. Figure 1a), the coordinates ri, n of the atom i in the nth residue can be written as

ri, n ) ui + r cos(nθ)eˆx + r sin(nθ)eˆy + nLeˆz

(1)

where r is the radial position with respect to the chain axis of an internal (fixed) reference point in the residue (e.g., its center of mass or specific atom), and ui is the displacement in a local frame with origin at that point. Due to the periodicity along eˆz, the number of residues N, the number of helical turns per supercell m, and the helical twist θ are related via θ ) 2πm/N (to decouple periodic images of the system in the orthogonal directions, a sufficiently large vacuum space is used). Hence, the PES is calculated on an irregular grid of points (Li, θi), relaxing all other internal degrees of freedom, i.e.,

∆E(Li, θi) ) min E(R) - E0 R

Figure 1. (a) Geometry of the infinite Ala chain in the right-handed R-helical conformation: L, helix length per residue along the helix axis (pitch); θ, helical twist; r, helix radius. Hydrogen bonds are indicated by dashed lines (only amide hydrogens are shown). Helix axis is shown by a long green arrow, and the red arrows indicate two consecutive backbone CR atoms. Other characteristic conformations of the polypeptide chain: (b) π-helix; (c) 310-helix; (d) 27 structure (repeated γ′-turn); (e) fully extended structure (FES).

helical motifs predicted by a given force field.16 Calculating the PES in a solvent-free environment is hence an appropriate computational scheme to characterize differences due solely to the force fields. In addition, we should stress that to ensure comparability to the ground-state, gas-phase DFT calculations, the force field results discussed in the following sections do not take into account thermal fluctuations. While vibrational entropy contribution to the free energy of the isolated infinite homopolypeptide chains are within the reach of DFT,17,18 an explicit account of the solvent still poses a major challenge. The fixed-charge protein force fields, on the other hand, are mainly targeted at liquid-phase simulations, but part of their parameter sets3,4 (e.g., partial atomic charges, torsional parameters) is fitted against high-level gas-phase quantum mechanical calculations, such as Hartree-Fock or second-order Møller-Plesett perturbation theory (MP2). Comparison of the force-field peptide energetics in the gas phase is thus not only legitimate, but also offers a means for addressing the crucial question of how well a given force field can describe properties that are intrinsic to the peptide structure. For some short polypeptides gas-phase experimental data are available19–22 and have been used to address the ability of a protein force field to accurately capture intramolecular interactions.23 Our paper is organized as follows. The description of the computational methods is presented in section II. In section III, we provide a detailed analysis of the potential energy landscape of Ala∞, the nonbonded energy components along the minimum energy pathway and the energetic convergence from finite AlaN fragments to the infinite Ala∞ chain. A discussion and concluding remarks are given in section IV. II. Methods The computational setup for mapping the PES of the infinite polypeptide chains with empirical force fields follows closely the one used in the first-principles calculations of Ireta et al.11,14 To represent the infinite polypeptides we use an orthorhombic box (“supercell”) and apply periodic boundary conditions in all three dimensions. In the Cartesian frame with basis {eˆx, eˆy, eˆz},

{

L(R) ) Li θ(R) ) θi

(2)

where E(R) is the total potential energy per residue for particular atomic configuration R ≡ {ri, n}, and E0 is a reference energy. Following refs 10 and 11, E0 is chosen to be the potential energy per residue in the FES of the polypeptide chain, Figure 1e. The final map of the PES is obtained by interpolating the calculated points with bicubic splines.24 A more detailed characterization of the PES is then carried out along the minimum energy pathway s ) s(L, θ) connecting all local minima. All calculations are carried out with the Tinker 4.2 molecular modeling package25 that allows for simulations of infinite polymer systems. The package offers access to a rich spectrum of force fields from which we have selected AMBER99,5 CHARMM27,6 OPLS-AA/L,26 and the (preliminary) proteinspecific versions of the AMOEBA polarizable force field, denoted hereafter as AMOEBApro.3,27 We have also implemented the very recent modification of the AMBER99 set, AMBER99SB,28 that was shown to improve the description of Ala and Gly backbone conformations. Electrostatic interactions were treated with the Ewald summation method, while the default cutoff of 9 Å was used for the van der Waals potential energy interactions (having the default Lennard-Jones functional form). For given (Li, θi), the initial coordinates of the polypeptide chain were constructed from the corresponding DFT-optimized geometries.11,14 Geometry optimization was performed with the limited-memory Broyden-Fletcher-Goldfarb-Shanno nonlinear optimization routine, constraining (Li, θi) and relaxing all other internal degrees of freedom (cf. eq 2), using a rms gradient tolerance of 0.005 kcal mol-1/Å. The functional form of the total potential energy is common for all fixed partial charge force fields considered here:3,4

E ) Ebonds + Eangles + Etors + EvdW + Eel

(3)

where Ebonds and Eangles represent the bond stretching and angle bending contributions, and EvdW, Eel, and Etors are the van der Waals, electrostatic, and torsional terms, respectively. The AMOEBApro force field extends the above form by adding an explicit polarization term.3,27,29 III. Results A. Potential Energy Landscape of Ala∞. The PES’s of the infinite Ala homopolypeptide chain, calculated with the AMBER99SB, CHARMM27, OPLS-AA/L, and AMOEBApro force fields, are shown in Figure 2. The topology of the PES’s is qualitatively similar to that obtained from DFT calculations,11,14 and can be described as consisting of two broad regions: a helical basin (L j 2.5 Å) where hydrogen bonds form between non-nearest neighbor peptide units, and an extended region for

6874 J. Phys. Chem. B, Vol. 112, No. 22, 2008

Penev et al.

Figure 2. Potential energy surface, eq 2, of the right-handed Ala∞ chain in vacuum. Energies are given in kcal/mol per residue; all panels show an energy range of 15 kcal/mol measured from the global PES minimum.

Figure 3. Total potential energy along the minimum energy pathway connecting the local minima on the PES of the Ala∞ chain, Figure 2. For reference, the DFT-calculated path according to refs 11 and 14 is shown as well, and the shaded regions indicate the key conformations of the polypeptide chain.

L > 2.5 Å where only nearest-neighbor interactions occur. The fine structure of the helical basin for the right-handed Ala∞ chain consists of two or three local minima of increasing L, depending on the force field. The PES profile along the minimum energy pathway s is presented in Figure 3 and geometry parameters of the helical conformations are compared in Table 1. We find that only the AMBER99 force field can satisfactorily reproduce all the π, R, and 310 helical minima, preserving also the relative stability calculated with DFT: ∆ER < ∆E310 j ∆Eπ. In contrast, no energy landscape feature can be associated with the 310 conformation on the PES calculated with CHARMM27 and AMOEBApro, and only a shallow metastable minimum is observed for OPLS-AA/L. Figure 3 also clearly reveals the differences between the “ff99” and “ff99SB” modifications28 of the AMBER force field. The optimization of the dihedral part of the potential corresponding to the φ and ψ backbone torsional angles in AMBER99SB considerably destabilizes the 310 minimum, somewhat less the R-helix, and hardly affects the π minimum.

For all force fields considered, the helical domain is significantly more stable relative to the FES than observed in the DFT PES. With the exception of CHARMM27, the force fields and DFT all agree in finding the R-helix to be the most stable conformation of the infinite Ala chain. The global minimum on the CHARMM27 PES corresponds to the π-helix, which also represents the strongest stabilized conformation with respect to the FES with ∆Eπ = -8 kcal/mol per residue. Figures 2 and 3 also show that although ∆E may differ substantially, the structural characteristics of the stable conformations in the helical region, given in Table 1, are overall similar for all force fields and the DFT. For the π- and R-helices, however, AMOEBApro gives longer O · · · H distance d and 10-20° smaller N-H · · · O angle ϑ. The helical region is separated from the extended conformations (unfolded region) by a barrier at L = 2.5 Å and twist θ = 130°. In this transition region the minimum energy pathways obtained with AMOEBApro and CHARMM27 agree better with DFT. Further increase in the twist and pitch transforms the polypeptide chain into the 27 conformation, Figure 1d, which is only marginally represented by AMBER99. This structure on the force-field PES appears more stable than the FES by =1 kcal/mol per residue. It is interesting to note that all PES’s of Ala∞, Figure 2, in agreement with the DFT results,14 reveal a slight left-handed twist of the FES, i.e., θFES ) 180° + δθ, with some small δθ > 0, a feature that is actually a fingerprint of sheet conformations.30 B. Nonbonded Components Along the Minimum Energy Pathway. The conformational preferences of the polypeptide chain are ultimately influenced by the nonbonded and torsional energy terms in the empirical potential, and their parametrization is crucial for the force-field performance.29 All force fields considered here include two- and three-body terms (stretching of covalent bonds, and bending of angles between nearestneighbor covalent bonds, respectively), apart from the nonbonded (EvdW, Eel) and torsional (Etors) contributions to the total potential energy E, eq 3. In Figure 4, we present these

Energetics of Infinite Homopolypeptide Chains

J. Phys. Chem. B, Vol. 112, No. 22, 2008 6875

TABLE 1: Geometry Parameters of the Helical Conformations of Ala∞, Figures 2 and 3a R

π AMBER99SB CHARMM27 OPLS-AA/L AMOEBApro DFT (ref 11)

310

L/θ

φ/ψ

d/ϑ

L/θ

φ/ψ

d/ϑ

L/θ

φ/ψ

d/ϑ

1.15/80 1.15/80 1.15/80 1.17/83 1.17/80

-76/-54 -72/-57 -77/-54 -80/-47 -76/-53

1.92/162 1.91/165 1.94/158 2.13/146 1.95/164

1.47/98 1.47/98 1.46/98 1.52/103 1.50/98

-65/-42 -62/-45 -65/-42 -64/-38 -63/-42

1.91/164 1.89/166 1.94/159 2.08/154 1.95/164

1.95/120

-60/-17

1.97/168

1.95/120

-68/-21

2.32/157

1.95/120

-60/-20

1.92/169

a All angles (θ, φ, ψ, ϑ) are measured in degrees, and distances (L, d) are given in Å. φ and ψ denote the standard backbone dihedral angles C-N-CR-C and N-CR-C-N, respectively. Hydrogen bonds are characterized by the O · · · H distance d, and the N-H · · · O angle ϑ.

Figure 5. Potential energy ∆E˜ along the minimum energy pathway s, excluding the EvdW contribution.

Figure 4. Nonbonded and torsional energy contributions (per residue relative to the FES) to the total potential energy of the Ala∞ chain. Full curves are just those for the corresponding force fields shown in Figure 3. Eel for AMOEBApro refers to the atomic multipoles energy component. Principal conformations of the polypeptide chain are schematically indicated as shaded regions on one of the panels, similar to Figure 3.

contributions to the total potential energy along the same pathways shown in Figure 3. This figure clearly demonstrates that although the total potential energies may show an overall similar profile, the contribution of individual energy components may differ considerably among the force fields. We find that the characteristic conformations of the Ala∞ chain, shown in Figure 1, are all local minima only for the Eel component in the force fields considered here. The 310-helix does not emerge as a stable structure on the CHARMM27 and AMOEBApro PES (cf. Figure 2). The lack of propensity toward the 310-helical conformation for the later two force fields may thus be attributed mostly to unfavorable two-body and three-body bonding terms (not shown in Figure 4). Furthermore, for all force fields, the π- and R-helices are local minima for both Eel and EvdW. It is also worth pointing out that the Etors contribution of AMBER99SB differs from the other force fields in that it has a destabilizing (>0) contribution along the whole energy pathway. Another qualitative difference arises between the Etors term of CHARMM27 and the other force fields: while for the former it has a maximum in the transition region, for the latter it is a minimal value. We note that DFT may encounter problems in correctly accounting for van der Waals/dispersion interactions,31–33 whereas all force fields considered here explicitly include a Lennard-Jones potential function. As a simple test for the importance of EvdW for Ala∞ we have compared in Figure 5 the DFT results to the force field calculated potential energy with the EvdW contribution omitted,34 ∆E˜ ) ∆E - ∆EvdW. Exclusion of the van der Waals term does indeed result in closer agreement

between the force fields and DFT in the helical region of the Ala∞ PES, although the helical conformations are still less stable in the DFT calculations. The agreement between the DFT and force field results, excluding CHARMM27, is particularly good around the π-helix minimum. This agreement lessens when moving along the minimum energy pathway toward more extended (larger L) conformations. The trend observed in Figure 5 may point to differences in the description of hydrogen bonds by DFT and force fields. According to ref 11, DFT predicts that hydrogen bonds are stronger in π-helices than in R-helices, and the later is stronger than in 310-helices, despite the fact that all three helices have hydrogen bonds with similar bond distances. The orientation of the groups forming such hydrogen bonds, however, varies between the helical conformations. As the force fields studied here do not include an explicit term to represent hydrogen bonding, they may not capture the hydrogen bonding orientational dependence.35 Hence, one would expect the attractive energy between the atoms forming the hydrogen bonds to be similar for the three helical conformations in force field calculations, as opposed to the DFT results presented here. Figure 5 also points to the importance of including van der Waals interactions in the electronic structure studies, which limits the DFT-force field comparison presented here. We do note, however, that it is still not entirely clear whether force fields are capable themselves of correctly capturing van der Waals interactions.36 C. Energetic Convergence: From Finite AlaN Fragments to Ala∞. In the previous section we have investigated the conformational space accessible to a residue in an infinite chain. In contrast to the study of dipeptides, an infinite polypeptide chain allows us to study the effect of long-range interactions on the conformations of a residue. For the case of finite chains it is important to known how large a finite helix should be in order to fully include the long-range contributions, i.e., what is the minimum length (in terms of number of residues) of the chain so that a “bulk-like” behavior is exhibited. To evaluate whether such a critical length is force-field dependent, we

6876 J. Phys. Chem. B, Vol. 112, No. 22, 2008

Penev et al. 6 also indicates that, for both R- and π-helices, chains of N g 8 residues are required to fully include the long-range effects on, e.g., a residue in the middle of the chain, and this critical length appears independent of the force field. IV. Discussion and Concluding Remarks

Figure 6. Cooperativity of the finite Ace-(Ala)N-Nme fragment as a function of N in (a) R-helical conformation, and (b) π-helical conformation. Cartoons indicate the formation of the first and second turn for the corresponding helix [note that the capping groups can contribute hydrogen bonds (hb)].

consider the energetics of a sequence of finite fragments of increasing length N, truncated from the Ala∞ chain in a given conformation, and capped with Ace and Nme at the N- and C-terminus, respectively. The study of helices of finite length has been shown to be a useful tool to separate the long- and short-range contributions to the cooperativity of R-helices37 as well as for assessing the accuracy of the “gas phase part” of empirical force fields.38 To assess the convergence of the contributions from nonbonded interactions, we calculate changes in chain strength as it grows, ∆εN ) εN - εN-1, where εN is the strength of the nonbonded interactions in a chain of N residues. Since we consider fragments of fixed geometry, the strength of nonbonded interactions is defined simply as εN ) EN - EN-1 - E0, where EN is the potential energy of Ace-(Ala)N-Nme and E0 is a reference energy. This leads to the simple difference expression,

∆εN ≡ ∆∆EN ) EN + EN-2 - 2EN-1

(4)

The quantity ∆∆EN measures the cooperative effect of the nonbonded interactions. In DFT calculations, ε measures primarily the hydrogen bond strength10 as current exchangecorrelation functionals do not describe van der Waals interactions. Hence, in DFT calculations ∆∆EN is a measure of the hydrogen bond cooperative effect. Since none of the force fields studied here contain an explicit term to account for hydrogen bonding, the cooperativity is merely due to the electrostatic interactions and, in the case of AMOEBApro, also to polarization effects. In Figure 6, ∆εN is plotted for capped finite fragments of the Ala∞ chain in R- and π-helical conformations since they are local minima for all force fields considered. Clearly, the formation of the first hydrogen bond (helix “nucleation”) leads to a substantial energy gain. Upon further increase of N the cooperativity ∆εN within the nonpolarizable force fields is a monotonically increasing function of N for the R-helix, and displays a narrow step for the π-helix in the range N ) 6-8. The explicit inclusion of multipole moments and polarization in the AMOEBApro force field gives rise to two distinct features of ∆εN: (i) the second hydrogen bond formed for both helices has practically the same strength as the first one, ∆εN = 0, and (ii) ∆εN decreases again for the next few hydrogen bonds formed. Thus, although all force fields show a cooperativity effect on ε, only AMOEBApro can make a distinction between the helix nucleation and at least the initial growth process. Figure

AMBER, CHARMM, and OPLS are among the most popular fixed atomic point-charge force fields for biomolecular simulations, and AMOEBA is one of the more elaborate polarizable force fields. The performance of these force fields is generally considered in the context of their ability to reproduce the conformational preferences of single residue dihedral angles (φ, ψ) from the study of dipeptide systems. In the present work, we introduce a new method to evaluate force fields using infinite chains. This approach allows a conformational characterization of a residue that incorporates the effects of a maximally cooperative network of nonbonded interactions. Our calculations of the Ala∞ system show that by using infinite chains, welldefined minima associated with helical conformations can be easily identified, and differences between the force fields can be readily illustrated. Ireta et al. have shown using DFT that such infinite systems reproduce the allowed local conformations of a residue in the bulk of a protein.14 The use of infinite chains hence also allows a comparison of the performance of ab initio versus force-field methods. Since polarization may be critical for the correct description of nonbonded interactions, e.g., hydrogen bonding,11,35 it is particularly interesting to examine the differences between the nonpolarizable force fields and AMOEBApro. From our calculations it appears that differences mainly show up in the context of forming finite helices. As DFT does include polarization effects, one would expect a better qualitative agreement between the DFT PES11,14 and the AMOEBApro PES (cf. Figures 2 and 3). Interestingly, this does not appear to be the case, with the later lacking the minimum associated with the 310-helix. In terms of the ability of the force fields to capture the characteristic conformations of infinite chains (cf. Figure 1), we find AMBER to be in closest qualitative agreement with DFT, reproducing the π-, R-, and 310-helices. It is interesting to note that these three conformations are all local minima only for the electrostatic energy term, while the van der Waals contribution displays local minima only for the π- and R-helices for this force field. In the case of the CHARMM27 force field, the π-helix is notably overstabilized. This bias can be traced back to the underlying CHARMM22 force field and has been discussed in detail, e.g., in refs 39-41 where a map-based spline interpolation correction was introduced (at present, this correction scheme is not available in the Tinker package and was not considered in this study). In conclusion, the use of infinite helices provides a reference limiting case that can be used to critically assess, compare, and improve the parametrization of empirical force fields commonly used in biomolecular simulations. By comparing DFT and force field result we also illustrated in a simple manner the importance of improving the description of the van der Waals interactions in both methodologies. Acknowledgment. This work was supported by the David and Lucile Packard Foundation and the NSF (MCB no. 0642086). J.I. acknowledges support by “programa de repatriación 2007” of the CONACYT.

Energetics of Infinite Homopolypeptide Chains References and Notes (1) Jackson, S. E. Folding Des. 1998, 3, R81-R91. (2) Frenkel, D.; Smith, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 2001; Vol. 1. (3) Ponder, J. W.; Case, D. AdV. Protein Chem. 2003, 66, 27–85. (4) MacKerell, A. D., Jr. J. Comput. Chem. 2004, 25, 1584–1604. (5) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179–5197. (6) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wio´rkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586–3616. (7) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (8) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. J. Comput. Chem. 2001, 22, 1205–1218. (9) Stickle, D. F.; Presta, L. G.; Dill, K. A.; Rose, G. D. J. Mol. Biol. 1992, 226, 1143–1159. (10) Ireta, J.; Neugebauer, J.; Scheffler, M.; Rojo, A.; Galva´n, M. J. Phys. Chem. B 2003, 107, 1432–1437. (11) Ireta, J.; Neugebauer, J.; Scheffler, M.; Rojo, A.; Galva´n, M. J. Am. Chem. Soc. 2005, 127, 17241–17244. (12) Improta, R.; Barone, V.; Kudin, K.; Scuseria, G. E. J. Am. Chem. Soc. 2001, 123, 3311–3322. (13) Kaschner, R.; Hohl, D. J. Phys. Chem. A 1998, 102, 5111–5116. (14) Ireta, J.; Scheffler, M. In preparation. (15) The DFT calculations of refs 11 and 14 employed ab initio pseudopotentials in conjunction with a plane-wave basis set of 70 Ry energy cutoff. The Perdew-Burke-Ernzerhof generalized gradient approximation was used to describe the electronic exchange and correlation interaction and the Brilloin-zone integration was replaced by the Γ-point. Greater details can be found in ref 11. (16) Sorin, E. J.; Rhee, Y. M.; Shirts, M. R.; Pande, V. S. J. Mol. Biol. 2006, 356, 248–256. (17) Ismer, L.; Ireta, J.; Boeck, S.; Neugebauer, J. Phys. ReV. E 2005, 71, 031911. (18) Ismer, L.; Ireta, J.; Neugebauer, J. J. Phys. Chem. B. 2008, 112, 4109–4112. (19) Hartings, M. R.; Kinnear, B. S.; Jarrold, M. F. J. Am. Chem. Soc. 2003, 125, 3941–3947.

J. Phys. Chem. B, Vol. 112, No. 22, 2008 6877 (20) Kohtani, M.; Jarrold, M. F. J. Am. Chem. Soc. 2004, 126, 8454– 8458. (21) Chin, W.; Piuzzi, F.; Dognon, J.-P.; Dimicoli, I.; Tardivel, B.; Mons, M. J. Am. Chem. Soc. 2005, 127, 11900–11901. (22) Abo-Riziq, A.; Bushnell, J. E.; Crews, B.; Callahan, M.; Grace, L.; de Vries, M. S. Chem. Phys. Lett. 2006, 431, 227–230. (23) Wei, Y.; Nadler, W.; Hansmann, U. H. E. J. Chem. Phys. 2007, 126, 204307. (24) Preusser, A. ACM Trans. Math. Software 1989, 15, 79–89. http:// www.fhi-berlin.mpg.de/grz/pub/xfarbe. (25) Ren, P.; Ponder, J. W J. Phys. Chem. B 2003, 107, 5933–5947. http://dasher.wustl.edu/tinker/ (26) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. J. Phys. Chem. B 2001, 105, 6474–6487. (27) Ren, P.; Ponder, J. W. J. Comput. Chem. 2002, 23, 1497–1506. (28) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, A. R. C. Proteins 2006, 65, 712–725. (29) Rasmussen, T. D.; Ren, P.; Ponder, J. W.; Jensen, F. Int. J. Quantum Chem. 2007, 107, 1390–1395. (30) Chothia, C. J. Mol. Biol. 1973, 75, 295–302. (31) Jurecˇka, P.; Cˇerny´, J.; Hobza, P.; Salahub, D. S. J. Comput. Chem. 2007, 28, 555–569. (32) Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2007, 3, 289– 300. (33) Improta, R.; Barone, V. J. Comput. Chem. 2004, 25, 1333–1341. (34) Exclusion of the EvdW contribution implies neglect of both the repulsive, ∼ 1/r12, and attractive, ∼-1/r6, terms, r being the interatomic distance. Note, however, that the range of r relevant to our discussion is large enough to make the repulsive term negligible. (35) Morozov, A. V.; Kortemme, T.; Tsemekhman, K.; Baker, D. Proc. Natl. Acad. Sci. USA 2004, 101, 6946–6951. (36) Sorin, E. J.; Pande, V. S. J. Comput. Chem. 2005, 26, 682–690. (37) Morozov, A. V.; Tsemekhman, K.; Baker, D. J. Phys. Chem. B 2006, 110, 4503–4505. (38) Jagielska, A.; Skolnick, J. J. Comput. Chem. 2007, 28, 1648–1657. (39) Feig, M.; MacKerell, A. D., Jr.; Brooks, C. L., III J. Phys. Chem. B 2003, 107, 2831–2836. (40) MacKerell, A. D., Jr.; Feig, M.; Brooks, C. L., III J. Comput. Chem. 2004, 25, 1400–1415. (41) The energy grid correction map approach of ref 40 was used to modify CHARMM22 so that it can reproduce the gas-phase (φ, ψ)-PES of the Ala dipeptide calculated at the MP2 level of theory.

JP800058F