10059
2009, 113, 10059–10062 Published on Web 07/02/2009
Unravelling the Conformational Dynamics of the Aqueous Alanine Dipeptide with First-Principle Molecular Dynamics M.-P. Gaigeot* UniVersite´ d’EVry Val d’Essonne, LAMBE UMR8587 Laboratoire Analyse et Mode´lisation pour la Biologie et l’EnVironnement, BlVd F. Mitterrand, Baˆt. Maupertuis, 91025 EVry, France ReceiVed: April 23, 2009; ReVised Manuscript ReceiVed: June 15, 2009
First principle DFT-based molecular dynamics simulations are performed in order to bring new insights into the conformational dynamics of the alanine dipeptide analogue Ac-Ala-NHMe (with methyl group caps at the extremities) immersed in liquid water at ambient temperature. Two simulations have been run for a total of 100 ps, which allows for a relevant statistical sampling of the phase space, at the ab initio level. PII-β equilibrium and (PII-β)-RR conformational interconversions are obtained, without a free-energy barrier for the Φ angle and with a rather low barrier of 2-3 kcal/mol for the Ψ angle, easily overcome from solute-solvent energy transfers. We furthermore give first insights into the rather strong zwitterionic character of the peptide bonds of the dipeptide when immersed in the liquid. The structural and zwitterionic properties extracted from first-principle dynamics in the liquid phase will be useful as benchmarks for force field developments. Introduction Though the alanine dipeptide is one of the simplest dipeptide motifs (the alanine dipeptide analogue Ac-Ala-NHMe with methyl group caps at the extremities is considered throughout the present work, and it will be referred to as 2-Ala in the text meaning that it is composed of two peptide bonds and one alanine residue), its structure in liquid water is still a source of debate. The most probable conformations that have been suggested in the literature are the polyproline II (PII), the β-sheet (β), and the right-handed R-helix (RR), but experiments and theory fail at providing a definitive answer on the actual conformation in liquid water. Some experiments provide a view where the PII and the RR structures coexist in the liquid,1,2 others show that only the PII structure can be seen.3,4 Raman spectroscopy further provides indications that the C7eq conformation, well-known to be the conformation of least energy adopted by the peptide in the gas phase, also coexists with the PII and RR conformations in liquid water.5 The calculations performed up to now have consisted of either quantum chemistry geometry optimizations of the alanine dipeptide hydrogen bonded with a few water molecules6 or classical molecular dynamics simulations of the alanine dipeptide immersed in liquid water.7-12 While the static calculations cannot provide a proper view of the bulk solvent surrounding the alanine dipeptide, this is not the case with molecular dynamics simulations in the liquid phase, but the results from classical simulations strongly depend on the force fields used, as elegantly summarized by Cho et al.7 Note that inclusion of the Charmm22/CMAP correction leads to a better samplings of conformations.9 We provide here, for the first time, DFT-based Car-Parrinello molecular dynamics simulations of 2-Ala immersed in liquid water in order to characterize the conformational dynamics and zwiterrionic properties of the aqueous dipeptide from first* E-mail:
[email protected].
10.1021/jp903745r CCC: $40.75
principle calculations. In this work, no force field is required, and the bulk liquid is explicitely taken into account at the same level of theory as the solute peptide. We employ the CPMD method,13,14 using the DFT/BLYP electronic representation,15,16 70 Ryd Cut-off expansion of the plane wave basis set, Martins-Trouillier pseudopotentials,17 a 500 au fictitious electron mass, and a time step of 5 au. The system is composed of the alanine dipeptide analogue 2-Ala immersed in 118 water molecules, in a cubic box length of 15.52 Å. Two CPMD dynamics, respectively of 44 (run1) and 54 ps (run2) durations have been performed, originating from two different initially thermalized conformations (nuclei positions and velocities). Average temperatures are, respectively, 312 ( 18 and 301 ( 23 K. Details on the general setup applied can be found in our previous publications.18-20 Briefly, the thermalization step was achieved with a 100 ps classical MD using the CHARMM22 force field for the solute and the TIP3P water model, followed by a 5 ps CPMD relaxation in the microcanonical ensemble. The CPMD simulations over which data are collected are subsequently performed in the microcanonical ensemble, and their time lengths give access to relevant statistical averages at this level of representation. Two configurations have been extracted from the classical MD that serve as starting points for CPMD, whose geometries are, respectively, of RR (run1, Φ ) -52°, Ψ ) -54°) and PII (run2, Φ ) -89°, Ψ ) 164°) types. According to the literature, they are expected to be the main types of conformations adopted by the alanine dipeptide immersed in liquid water and have therefore been purposely chosen here as starting conformations in order to get insights on the stability of these conformations in the liquid, on the organization of the water solvent around both types of conformations in relation to their relative stability, and on the interconversion between these conformations during the dynamics. Running two separate dynamics from two distinct conformations enhances sampling and statistics. 2009 American Chemical Society
10060
J. Phys. Chem. B, Vol. 113, No. 30, 2009
Figure 1. Evolution with time of the Φ (top) and Ψ (bottom) dihedral angles of the aqueous alanine dipeptide from the two Car-Parrinello dynamics performed here (run1, left, run2, right). Angles (vertical axes) are in degrees, and time (horizontal axes) is in picoseconds. The angle discontinuities arise from definitions of the angles at (180°.
Letters
Figure 2. 2D representation of the (Φ, Ψ) dihedral angles of the aqueous alanine dipeptide from Car-Parrinello dynamics. Φ is on the horizontal axis, and Ψ is on the vertical axis. The angle discontinuities arise from definitions of the angles at (180°. The top spot corresponds to PII/β and the bottom spot to RR conformations.
Results and Discussion In the following, we adopt the nomenclature of Cho et al.7 in order to define the dipeptide structures. The β conformation corresponds to -180 e Φ e -120° and 60 e Ψ e 180°, the right-handed RR helix corresponds to -180 e Φ e 0° and -150 e Ψ e 60°, and the polyproline PII conformation corresponds to -120 e Φ e -30° and 60 e Ψ e 180°. The evolution with time of Φ and Ψ angles of the alanine dipeptide are presented in Figure 1 for the two CPMD dynamics. The dynamics originating from the alanine RR conformation (run1; Figure 1, left) grossly remains in that state, with Φ/Ψ average values of -106/-7°, which correspond to an RR helix arrangement of the dipeptide, on average. The second dynamics originates from a PII conformation (run2; Figure 1, right), and one can see that Φ oscillates between two average values along all of the dynamics, that is, -120 and -150°, while Ψ evolves from ∼150 (first half of the dynamics) to ∼20° (second half of the dynamics). Therefore, the alanine dipeptide alternates between the polyproline PII and the more extended β conformations during the first half of the dynamics. The alanine dipeptide subsequently undergoes a drastic conformational change with the Φ angle still adopting the -120/-150° alternation, but the Ψ angle is now at 20°. This corresponds to a change for an RR helix arrangement of the dipeptide, nicely spontaneously sampled during the dynamics. The Φ/Ψ 2D plot is presented in Figure 2, which summarizes the two distinct RR/(PII-β) domains sampled during the dynamics and the transition regions in between. Free-energy curves calculated as G(Φ) ) -kT ln P(Φ) and G(Ψ) ) -kT ln P(Ψ), where P(Φ) (P(Ψ)), is the histogram probability of a certain angle Φ(Ψ) in the course of the dynamics,21 are presented in Figure 3. The free-energy curve G(Φ) presents a primary minimum corresponding to Φ ) -106° and a second minimum for Φ close to -150°. The potential is going uphill rather smoothly between these two minima, with no evident energy barrier, and yielding an overall free-energy difference of 0.2 kcal/mol between these two Φ regions. Going further uphill, the free-energy curve toward Φ angles in the range of -50/0°, as would be required for a full exploration of the angle in RR conformations of the dipeptide, necessitates an increase of free energy of 2-3 kcal/mol. From the point of view of the Φ angle coordinate, this rather broad free-energy curve explains the easy conversions between β, PII, and RR conforma-
Figure 3. Free-energy profiles along the Φ (full line) and Ψ (dashed line) dihedral angles obtained along the two CPMD trajectories recorded in the present work. Angles are reported in degrees, and the energy is in kcal/mol.
tions along the time. On the contrary, the free-energy curve G(Ψ) presents two separate wells, one centered on -7° (and grossly related to RR helix types of conformations) and one centered on 148° (grossly related to a mixing of β and PII types of conformations), the first one being the broadest. The freeenergy barrier to overcome for the conversion between β/PII and RR conformations is ∼2 kcal/mol. Beware that this energy barrier is not well-defined in our calculations as there is only one spontaneous passage of the barrier during the dynamics. Longer MD and subsequent recrossings of the barrier would be required in order to refine the energy barrier, but this would undoubtely lead to a decrease of the value obtained here. The gain of energy needed to cross the RR/(PII-β) barrier is very likely obtained via solute-solvent energy transfers. The broader basin associated with RR conformations indicates a slightly higher energy preference, presumably favored by entropic effects. Radial distribution functions (RDF) are presented in Figure 4. Though the 100 ps time length of the DFT-based MD reported here is long in comparison to the usual time lengths of these simulations in the literature, it remains short in comparison to what is affordable in classical MD. The following conclusions would therefore be improved if longer time scales could be
Letters
J. Phys. Chem. B, Vol. 113, No. 30, 2009 10061
Figure 4. Radial distribution functions of alanine dipeptide solute-water solvent. Full line: run1-RR; dashed line: run2-(PII-β); dashed-dotted line: run2-RR (see text for the signification). (A) N-H(Nter) · · · OW, (B) N-H(Cter) · · · OW, (C) CdO(Cter) · · · HW, (D) CdO(Nter) · · · HW.
TABLE 1: Hydrogen Bonds Formed between the Alanine Dipeptide and the Surrounding Water Solvent Moleculesa run1-RR
run2-PII-β
run2-RR
Positions of RDF Peaks: first peak/second peak (Å) N-H(Cter) N-H(Nter) CdO(Cter) CdO(Nter)
1.91/3.78 2.00/3.86 1.80/3.35 1.71/3.61
2.00/no peak 1.90/no peak 1.73/3.28 1.74/3.43
1.92/no peak 2.15/3.71 1.74/3.28 1.71/no peak
Strengths of H Bonds N-H(Cter) · · · OW N-H(Nter) · · · OW CdO(Cter) · · · HW CdO(Nter) · · · HW
5.7 4.6 17.4 29.5
1.5 5.0 16.4 13.0
3.5 3.0 10.7 50.3
Number of H-Bonded Water N-H(Cter) N-H(Nter) CdO(Cter) CdO(Nter) a
0.9 ( 0.3 0.8 ( 0.4 2.3 ( 0.6 2.1 ( 0.4
0.5 ( 0.5 1.0 ( 0.3 3.0 ( 0.5 2.9 ( 0.4
0.7 ( 0.5 0.8 ( 0.4 2.8 ( 0.5 2.8 ( 0.4
See text for calculation details.
achieved. Note that we have calculated two sets of RDFs over the 54 ps CPMD (run2), one associated with the PII-β conformational sampling and one associated with the RR conformational sampling over the time. This therefore provides two sets of RDFs for the RR type of conformations, from run1 and run2, respectively. The RDFs of the water solvent are not displayed here, but they are identical to the ones published in our previous investigations,18,19 corresponding to the expected liquid water. Positions of the RDF peaks are reported in Table 1. Strengths of the hydrogen bonds that are formed between the N-H or CdO groups of 2-Ala and the surrounding water molecules can also be found in Table 1, calculated as gmax/gmin, where gmax and gmin are, respectively, the amplitudes of the first maximum and first minimum of the RDF. All N-H and CdO groups of 2-Ala are involved in hydrogen bonds with the surrounding solvent, as indicated by the first maxima in the radial distribution functions that are located at distances in the 1.70-2.00 Å range. All RDF of N-H(Nter) display a rather organized second shell of surrounding water molecules, as revealed by the presence of a second peak in the radial distribution functions located at around 3.7-3.8 Å. This
second peak disappears in the RDF of N-H(Cter). The same picture can be drawn for the CdO groups, with a second shell of water molecules rather well organized around the CdO(Nter) and not organized around the CdO(Cter). Overall, the water liquid is thus organized into H-bonded (first shell) and wellorganized second shell molecules around the N-terminus of the alanine dipeptide, while it only participates in the first hydrogenbonded shell at the C-terminus. Intramolecular CdO · · · N-H H bonds are never formed along the dynamics. We found that the two RDFs related to the RR conformation obtained from the two separate dynamics performed here are not strictly equivalent (run1-RR versus run2-RR). In particular, the amplitudes of the first peak of the radial distribution functions obtained in run2-RR are systematically lower than thoseobtainedinrun1-RR,whilethesecondpeakofCdO(Nter) · · · HW is not present in run2-RR. Though caution should be taken when analyzing such “short dynamics” (in comparison to the length scales affordable in classical MD), such differences are expected considering that the average RR geometries sampled during the two runs are not strictly equivalent (see discussion above about Φ/Ψ domains), therefore giving rise to slightly different details of hydration. We have checked that the RDFs do not evolve when they are calculated over different periods of time along the trajectories, which does indicate that the observed geometrical differences arise from the sampling of the large RR geometrical domain. The biggest difference arising in the solute/solvent hydrogen bonds between the RR (run1 and run2) and PII-β (run2) alanine conformations can be seen in the N-H(Cter) group as the first peak of the N-H(Cter) · · · OW RDF of run2-(PII-β) has a much lower amplitude than that in the two other dynamics. Moreover, the second shell of water solvent around the (PII-β) conformations can be seen as less well organized than that around the RR conformations, as can be inferred from less well defined second peaks in all radial distribution functions associated with run2-(PII-β). The number of water molecules hydrogen bonded to the N-H and CdO groups are reported in Table 1. The N-H groups of the dipeptide are involved in one (or close to one) H bond on average, apart from the PII-β conformation for which the N-H(Cter) does not form such a stable hydrogen bond along the dynamics (leading to the smaller value of 0.5 H-bonded water molecules on average). Both CdO groups accommodate two hydrogen bonds on average along the run1-RR trajectory, while they accommodate up to three H bonds on average in the run2-PII-β and run2-RR trajectories. Note that the ratio of 5 between CdO(Cter) · · · HW and CdO(Nter) · · · HW hydrogen bond strengths (run2) is not reflected in the number of H bonds formed as both carbonyls are involved in an average of 2.8 H bonds. The strengths of the solute-solvent H bonds that are formed depend on the solute site (N-H versus CdO) but also vary for a similar site depending on the underlying average conformation of the peptide. For comparisons, the OW · · · HW strength in liquid water is on the order of 7,19,22 though the liquid is slightly more structured in the present simulations, with values actually close to 10. Carbonyl groups of 2-Ala form the strongest H bonds, with the CdO(Nter) of the run2-RR dynamics giving a strength up to 50. H-bond strengths arising from the N-H groups are on the same order as the solvent-solvent H bonds, with the exception of the N-H(Cter) in the PII-β dynamics, which displays a significantly lower H-bond strength. The very high values of CdO · · · HW H-bond strengths very likely arise from the CO-dNH+ zwitterionic form of the two
10062
J. Phys. Chem. B, Vol. 113, No. 30, 2009
Letters
TABLE 2: Zwitterionic Character of the Alanine Dipeptide Immersed in Liquid Water Investigated through Bond Length Values of CdO, N-H, and N-C Groups of the Two Peptide Bonds As Obtained in the Present Calculations (column 2), Compared to the Values of the Gas-Phase Alanine Dipeptide (values from ref 23, columns 3 and 4), and Aqueous trans-NMA (N-methyl-acetamide) (from ref 18, column 5)a bond type
2-Ala liquid phase
2-Ala23 gas phase (C7eq)
2-Ala23 gas phase (C5)
CdO(Nter) CdO(Cter) N-H(Nter) N-H(Cter) N-C(Nter) N-C(Cter)
1.267 ( 0.023 1.265 ( 0.024 1.030 ( 0.028 1.034 ( 0.026 1.340 ( 0.027 1.359 ( 0.029
1.239 ( 0.005 1.246 ( 0.006 1.027 ( 0.009 1.020 ( 0.011 1.364 ( 0.006 1.369 ( 0.008
1.241 ( 0.005 1.241 ( 0.005 1.019 ( 0.004 1.023 ( 0.005 1.366 ( 0.008 1.372 ( 0.008
trans-NMA18 liquid phase 1.268 ( 0.022 1.031 ( 0.032 1.351 ( 0.028
a Nter refers to the N-terminal of the dipeptide, and Cter refers to the C-terminal of the dipeptide; NMA has only one peptide bond. Bond-lengths are reported in Å.
peptide bonds of the peptide when immersed in the aqueous solution. In a previous publication on the prototype N-methylacetamide peptide bond model, we have already shown the appearance of such a zwitterionic state in liquid water.18 The zwitterionic character of the two peptide bonds of 2-Ala can be inferred from the behavior of the CdO, N-H, and N-C bond lengths as obtained during the dynamics. The average values obtained over the two trajectories are reported in Table 2 (there are no significant differences between the two runs), together with the values from the 20 K trajectories of the isolated alanine peptide (at the same level of theory) extracted from our previous investigation23 and with the values from the 300 K aqueous trajectory of N-methyl-acetamide (again at the same level of theory) extracted from our previous investigation.18 In the case of isolated alanine, we report the bond lengths from the C7eq trajectory, where an internal H bond between the Cter CdO and Nter N-H is observed on average, and from the C5 conformation, where no internal H bond is formed over the time. As can be immediately observed, all bond lengths obtained for the aqueous alanine peptide here are very similar to the ones observed for the aqueous N-methyl-acetamide for which we have further demonstrated the zwitterionic character of the peptide bond using the electronic Wannier Centers (see ref 18). Clearly, we observe that the CdO and N-H bonds of 2-Ala are lengthened and weakened by the formation of hydrogen bonds with the liquid, while the N-C peptide bonds are strongly shortened. It is remarkable that the length of the N-C peptide bond of the dipeptide at the N-terminal is even more strongly shortened than that of the C-terminal one. This correlates very well with the surrounding solvent, which is more structurally organized around the N-terminal of the peptide, as shown above. As a result, the N-terminal peptide bond displays a more pronounced (C-O)-dNH+ zwitterionic character than the C-terminal peptide bond. Note that, although an intramolecular H bond is formed in the gas phase C7eq peptide, this does not induce such tremendous CdO/N-H lengthenings and N-C shortenings. The zwitterionic character can only be obtained once the peptide is immersed in the liquid phase. Conclusion In conclusion, first-principle DFT-based molecular dynamics simulations show that the three PII, β, and RR conformations of the alanine dipeptide (alanine dipeptide analogue Ac-Ala-NHMe with methyl group caps at the extremities) can coexist in the liquid phase, with no (PII/β interconversion) or low (PII/β to RR interconversion) free-energy barriers. The N-terminal part of the peptide was found to have a surrounding solvent that is well organized into first- and second-shell H-bonded water molecules, while the C-terminal part display only first-shell H-bonded water
molecules. Evidence of the strong zwitterionic character of the peptide bonds of the aqueous alanine dipeptide was given, with an enhancement of the zwitterionic character of the N-terminal peptide bond. As a consequence, an average of three water molecules are hydrogen bonded to the CO- entities of the peptide chain, and strong CO--water H-bond strength enhancements have been observed. We offer here 100 ps first-principle dynamics of the alanine dipeptide immersed in explicit liquid water (118 water molecules) that can be used as benchmarks for structure, free-energy, and zwitterionic properties. Longer dynamics would be required in order to get more definite conclusions, and we hope that the present simulations will stimulate further developments and investigations on these important topics. Acknowledgment. IDRIS (Orsay, France) is acknowledged for generous access to the BlueGene platform. Support from Genopole-France through the program ‘ATIGE’ Action The´matique Incitative de Ge´nopole is acknowledged. References and Notes (1) Madison, V.; Kopple, K. D. J. Am. Chem. Soc. 1980, 102, 4855. (2) Mehta, M. A.; Fry, E. A.; Eddy, M. T.; Dedeo, M. T.; Anagnost, A. E.; Long, J. R. J. Phys. Chem. B 2004, 108, 2777. (3) Weise, C. F.; Weisshaar, J. C. J. Phys. Chem. B 2003, 107, 3265. (4) Kim, Y. S.; Wang, J. P.; Hochstrasser, R. M. J. Phys. Chem. B 2005, 109, 7511. (5) Takekiyo, T.; Imai, T.; Kato, M.; Taniguchi, Y. Biopolymers 2004, 73, 283. (6) Han, W. G.; Jalkanen, K. J.; Elstner, M.; Suhai, S. J. Phys. Chem. B 1998, 102, 2587. (7) Kwac, K.; Lee, K.-K.; Han, J. B.; Oh, K.-I.; Cho, M. J. Chem. Phys. 2008, 128, 105106. (8) Seabra, G. D.; Walker, R. C.; Elstner, M.; Case, D. A.; Roitberg, A. E. J. Phys. Chem. A 2007, 111, 5655. (9) MacKerrel, A.; Feig, M.; Brooks, C. J. Comput. Chem. 2004, 25, 1400. (10) Drozdov, A. N.; Grossfield, A.; Pappu, R. V. J. Am. Chem. Soc. 2004, 126, 2574. (11) Hu, H.; Elstner, M.; Hermans, J. Proteins 2003, 50, 451. (12) Brooks, C. L.; Case, D. A. Chem. ReV. 1993, 93, 2487. (13) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (14) CPMD; International Business Machines Corporation: Zurich, Switzerland, 1990-2008 and Max Planck Institute fuer Festkoerperforschung: Stuttgart, Germany, 1995-2001. (15) Becke, A Phys. ReV. A 1988, 38, 3098. (16) Lee, C.; Yang, W.; Parr, R. Phys. ReV. B 1988, 37, 785. (17) Trouillier, N.; Martins, J. Phys. ReV. B 1991, 43, 1993. (18) Gaigeot, M.-P.; Vuilleumier, R.; Sprik, M.; Borgis, D. J. Chem. Theory Comput. 2005, 1, 772. (19) Gaigeot, M.-P.; Sprik, M. J. Phys. Chem. B 2004, 108, 7458. (20) Gaigeot, M.-P.; Sprik, M. J. Phys. Chem. B 2003, 107, 10344. (21) Marinica, C.; Gre´goire, G.; Desfrancois, C.; Schermann, J.; Borgis, D.; Gaigeot, M.-P. J. Phys. Chem. A 2006, 110, 8802. (22) Sprik, M.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1996, 105, 1142. (23) Gaigeot, M.-P. J. Phys. Chem. A 2008, 112, 13507.
JP903745R
