Peptide folding driven by Van der Waals interactions
Shen-Shu Sung* Department of Pharmacology, College of Medicine, Penn State University, Hershey, Pennsylvania 17033 Received 20 December 2014; Accepted 13 May 2015 DOI: 10.1002/pro.2710 Published online 26 May 2015 proteinscience.org
Abstract: Contrary to the widespread view that hydrogen bonding and its entropy effect play a dominant role in protein folding, folding into helical and hairpin-like structures is observed in molecular dynamics (MD) simulations without hydrogen bonding in the peptide-solvent system. In the widely used point charge model, hydrogen bonding is calculated as part of the interaction between atomic partial charges. It is removed from these simulations by setting atomic charges of the peptide and water to zero. Because of the structural difference between the peptide and water, van der Waals (VDW) interactions favor peptide intramolecular interactions and are a major contributing factor to the structural compactness. These compact structures are amino acid sequence dependent and closely resemble standard secondary structures, as a consequence of VDW interactions and covalent bonding constraints. Hydrogen bonding is a short range interaction and it locks the approximate structure into the specific secondary structure when it is included in the simulation. In contrast to standard molecular simulations where the total energy is dominated by chargecharge interactions, these simulation results will give us a new view of the folding mechanism. Keywords: intramolecular interaction; compact structure; secondary structure; covalent bonding constraint; hydrogen bonding; folding mechanism Introduction Historically, protein intramolecular hydrogen bonding was a determining factor in the original prediction of protein secondary structures.1,2 Currently, the hydrophobic effect is considered to be the most important single factor in protein folding in our textbooks,3,4 and the hydrophobic effect is considered to arise from entropy changes in water hydrogen bonding around water-exposed protein nonpolar groups. However, the helical structure could be predicted based on the energy minimum of peptide intramolecular VDW interactions.5 The entropy and enthalpy effect of VDW interactions contribute importantly to protein folding in the VDW dynamic hydration shell proposal.6,7 Studies also showed that hydrophobic interactions are an important driving force for forming the transition state8 and for molecular recogniGrant sponsor: College of Medicine, Penn State University. *Correspondence to: Shen-Shu Sung, Department of Pharmacology, R130, Penn State College of Medicine, 500 University Drive, Hershey, PA 17033-0850. E-mail: [email protected] or [email protected]
C 2015 The Protein Society Published by Wiley-Blackwell. V
tion.9 It may be difficult to separate these contributing factors in laboratory tests, but it is feasible to change or remove hydrogen bonding effects in computations. In the widely applied atom-based point charge model, each atom is modeled by a point in space with a partial charge representing the molecular charge distribution and with a 6–12 VDW potential representing dispersion forces and electron shell repulsions. VDW interactions and covalent constraints of bond lengths, bond angles, and dihedral angles are necessary for keeping atomic distances realistic and their parameters are largely derived from experimental structural data. Atomic partial charges are mainly based on calculated values.10 Hydrogen bonding is calculated as part of Coulomb interactions between atomic charges. Setting the atomic charges to zero is a simple way to hypothetically remove hydrogen bonding and other charge– charge interactions from a computation, and it is tested in this study focusing on the qualitative nature of the peptide folding mechanism.
PROTEIN SCIENCE 2015 VOL 24:1383—1388
1383
Figure 1. From the left, the helix turn formation observed at 0.12 ns, a helix-turn-helix at 0.94 ns, an approximate helix at 6 ns, and a partial helix at 10 ns during the first simulation with AMBER force field parameters without atomic charges, and a helical structure at 17.75 ns during the simulation with CHARMM parameters without atomic charges. The carbon atoms are colored green, oxygen red, nitrogen blue, and polar hydrogen white. Water molecules and hydrogen atoms bound to carbon atoms are not shown. The molecular graphics was generated using Maestro display software from Schrodinger LLC.
Results Helix folding Early peptide folding simulations11,12 were carried out successfully on synthetic peptide sequences. Among them, a well characterized alanine-based ahelical peptide13 Ac(AAQAA)3Y-NH2 is studied here in helix folding simulations without atomic partial charges of the peptide and water molecules. The first simulation was carried out using the AMBER force field14 covalent bonding constraints and VDW parameters. The initial structure was a fully extended peptide in a solvent box of 4695 water molecules, built using the AMBER software package.15 This simulation was carried out at 273 K for 20 nanoseconds (ns). Helical turns formed at both ends of the peptide at 0.12 ns, and tertiary interactions between helical segments were observed between 0.6 and 1.1 ns, as shown in Figure 1. Helical structures with partially unfolded N- or C-termini were observed repeatedly during several periods of time, including 2.1-2.3 ns, 3.8-4.3 ns, and 4.9-5.6 ns. More complete helical structures were observed between 5.8 and 6.2 ns. The structure at 6.0 ns in Figure 1 is an example of an approximate helical structure. Partial unfolding and refolding were observed repeatedly. The structure at 10 ns in Figure 1 is a partially unfolded helical structure. Approximate helical structures were observed again at 11.4–11.6 ns, 16.3–16.7 ns, and 17.5–17.9 ns. Helix-turn-helixlike structures were observed in several periods of time, including 10.3–10.6 ns and 19.7–19.9 ns.
1384
PROTEINSCIENCE.ORG
Two views of a helical structure observed at 6.05 ns are shown in Figure 2, one perpendicular to the helix axis and the other along the axis from the C-terminus. Its backbone amide oxygen and hydro˚ distance are connected by gen atoms within 2.6A dotted lines. There are 10 (i, i14) oxygen-hydrogen pairs corresponding to the a-helix structure. There
Figure 2. Two views of the helical structure observed at 6.05 ns. One or two a-carbon atoms are labeled in each turn for the visualization of the number of residues per turn. The view perpendicular to the helix axis shows backbone oxygen and hydrogen atom pairs within 2.6 A˚, including ten (i, i14) pairs corresponding to a-turns and four (i, i13) pairs corresponding to b-turns, including 310 helix turns. For clarity, sidechains are not shown in the view along the helical axis from the Cterminus (lower right). Atom colors and other descriptions are the same as in Figure 1.
Driving Forces of Peptide Folding
are 4 (i, i13) pairs, 2 of them are b-turns near both ends of the helix, and the other two share atoms with (i, i14) pairs, indicating possible 310 helix turns. The oxygen-hydrogen atom pairs not connected by the dotted lines have their distances ˚ . One or two a-carbon atoms per larger than 2.6 A turn are labeled for visualizing the number of residues per turn, which is close to 3.6 residues per turn of the a-helix. The pitch of the helix is between ˚ , depending on the atom pairs used in 5.3 and 5.5 A the distance measurement. Without hydrogen bonds connecting the backbone hydrogen and oxygen atoms, amide groups tilt with oxygen atoms outwards, as shown in both views in Figure 2. This amide tilt makes the oxygen-hydrogen distance larger than that in a standard a-helix. During this simulation, atomic coordinates were recorded and the backbone /, w angles calculated every 10 ps. Several backbone dihedral angle criteria have been used in literatures16,17 for calculating the helical content or the helical ratio, such as / 5 2578 6 308 and w 5 2478 6 308, / 5 2658 6 358 and w 5 242.58 6 37.58, / 5 2658 6 358, and w 5 2378 6 308, etc. In the current study, when the /, w angles of two or more consecutive amino acid residues are within 308 of the standard a-helical angles (/ 5 2578 and w 5 2478), these residues are assumed to be in the helical conformation, which includes the majority of a-helix and 310 helix residues. The ratio of the number of helical residues to the total number of the residues (16 residues of the alanine based peptide) as a function of time in the simulation is shown in Figure 3. The average helical ratio during the 20 ns is 18.66%. In comparison, a standard MD simulation was carried out with atomic charges, which shows a higher helical ratio of 35.03%, indicating that hydrogen bonds stabilize helical structures. Amide tilting changes the backbone dihedral angles from those of the standard helix, which is another reason for the lower helical ratio without hydrogen bonding. To confirm these qualitative observations, multiple simulations with different initial structures were carried out. The initial structures were unfolded coil structures, generated from a high temperature (900 K) simulation at 10 ps time intervals. AMBER force field parameters were used without atomic charges of the peptide and water molecules. A group of 16 simulations was carried out for 20 ns each at 273 K with 1529 solvent water molecules. Two helical segments with tertiary interactions forming helix-turn-helix-like structures were observed in all these simulations, and the peptide formed single helices in 11 of the 16 simulations. The results of multiple simulations are qualitatively in agreement with those of the previous simulation with the extended initial structure. The helical ratio of each simulation averaged over the 20 ns ranges from
Sung
Figure 3. The helical ratio (the number of residues in helical conformations to the total number of residues of the peptide) as a function of time in simulations with AMBER force field parameters. The solid line shows the helical ratio in the simulation with the extended initial structure and without atomic charges, the dashed line shows the average helical ratio of 16 simulations with various unfolded initial structures and without atomic charges, and the dotted line shows the helical ratio in a standard MD simulation with the extended initial structure and with atomic charges.
10.15% to 20.21%. The helical ratio of the previous simulation (18.66%) is within this range. The average helical ratio of all 16 simulations as a function of time is shown in Figure 3. To see the force field dependence of these observations, CHARMM force field18 covalent bonding constraints and VDW parameters have been applied to a simulation of the same alanine-based peptide without atomic charges of the peptide and water molecules. This simulation was carried out at 273K for 20 ns using the NAMD program.19 The initial structure was a fully extended peptide in a solvent box of 4,427 water molecules, built using the VMD program.20 During the first 10 ns, the peptide was in random coil conformations with isolated helical turns. After the first 10 ns, various helical structures were observed, including partial helices and helix-turn-helix-like structures. Between 17.20 and 17.75 ns, the peptide formed a single helix. One of these structures observed at 17.75 ns is shown in Figure 1. Compared with the structure with AMBER parameters, the backbone amide tilt is more pronounced with CHARMM parameters. The periodicity of the turns is between 3.5 and 4 residues. The pitch ˚ , because of the larger of the helix is shorter at 5.2 A amide tilt. Helix folding occurred with different force field parameters without peptide intramolecular and water hydrogen bonding.
Hairpin folding Without hydrogen bonding, is the secondary structure amino acid sequence dependent? Besides the ahelix, the b-sheet is another common secondary structure. With limited computing resources, a small peptide is the first choice for a folding simulation. Blanco et al.21 have successfully designed a small b-
PROTEIN SCIENCE VOL 24:1383—1388
1385
Figure 4. A hairpin-like structure observed in the simulation using CHARMM force field parameters without atomic charges (left) and that using AMBER parameters without atomic charges (right), and a b–hairpin structure using CHARMM parameters with a fraction of atomic partial charges (center). Dotted lines connect the potential hydrogen bonding atom pairs with their distances shown in A˚. Atom colors and other descriptions are the same as in Figure 1.
hairpin peptide YQNPDGSQA, and Wu et al.22 have carried out a folding simulation of this b-hairpin peptide using the SGMD method23 with explicit solvent. Can this peptide fold into a hairpin structure without hydrogen bonds? A MD simulation was carried out with CHARMM force field parameters18 without atomic charges of the peptide and water molecules. The initial structure was a fully extended peptide in a solvent box of 1002 water molecules, built using the VMD program.20 The simulation was carried out at 300 K for 10 ns using the NAMD program.19 During the simulation, structures were recorded every 5 ps. Various transient hairpin-like structures were observed during the time period of 4,320 to 4,360 ps. A structure at 4340 ps is shown in Figure 4. Hydrogen bonding atom pairs indicated by NOE data21 are connected by dotted lines. Three out of the four distances between the hydrogen and the ˚ , about 0.5 A ˚ oxygen atoms are 2.93, 2.47, and 2.58 A larger than normal hydrogen bond distances. The backbone NH group of the amino acid A9 is pointing away from the Y1 carbonyl group. To test the effect of hydrogen bonding on the hairpin structure, a short simulation is carried out starting from the hairpin-like structure with a fraction (1/2) of atomic partial charges of the CHARMM18 force field. Within 50 ps, all four hydrogen bonds formed with the hydrogen-oxygen distan˚ , consistent with the NOE ces smaller than 2.6 A 21 data, as shown in Figure 4. These results show hairpin-like structures can form rapidly as a consequence of VDW interactions and covalent constraints, and hydrogen bonding locks them into the specific secondary structure. This observation may suggest a computational approach using mainly VDW interactions to direct folding into compact structures and including certain charge-charge interactions to stabilize the folded structure.
1386
PROTEINSCIENCE.ORG
AMBER force field14 parameters have also been tested in hairpin folding simulations. 20 unfolded structures were taken from a simulation at 900 K in 10 ps time intervals as initial structures of the following hairpin folding simulations. Each simulation was carried out for 20 ns at 300 K with 1,273 water molecules, using AMBER force field14 parameters without atomic charges of the peptide and water molecules. During these simulations b-turns were observed repeatedly, but they are not always flanked by antiparallel b-strands forming the hairpin structure. For example, in the 12th simulation hairpinlike structures were observed at 0.80–1.25 ns, but the two strands moved apart after 1.25 ns while the b-turn at the P4-D5 remained. A hairpin-like structure observed at 1.14 ns in this simulation is included in Figure 4. In this structure, three potential hydrogen bonding atom pairs are at distances ˚ , but the N-terminal NH3 between 2.5 and 3.5 A group does not point towards the oxygen atom in the C-terminus. Using the structures saved every 5 ps, the / and w angles were calculated. Standard type I bturn angles are /i11 5 2608, wi11 5 2308, / i12 5 2908, and w i12 5 08. The difference of these /, w angles between i11 and i12 positions is 308. The largest dihedral angle difference between the type I b-turn and the a-turn is 478 in wi12. To clearly distinguish between the a-turn and the b-turn and between the i11, i12 positions of the type I b-turn, the type I, II, I’, and II’ b-turns and the a-turn are assigned using a smaller variation range of 6158 from their standard /, w angles. The percentage of a type of turns at a (i11, i12) site in the total number of structures recorded during the simulation is the ratio of that type of turns at the given site. The type I b-turn ratios at i 1 1 5 2 to 7 (i 1 2 5 3 to 8) sites (the 20 simulation average) are 0.000, 0.117, 3.625,
Driving Forces of Peptide Folding
2.320, 0.870, and 1.420%, and helical ratios at the same sites are 0.052, 0.085, 0.695, 0.273, 0.140, and 0.442%. The first and the last residues do not have sufficient numbers of atoms for their dihedral angle calculations. The type II, I’, II’ b-turn ratios are all lower than 0.1%. The highest turn ratio is 3.625% of the type I b-turn at the P4-D5 site, consistent with the model structure based on NOE data.21 During the simulation, the b-hairpin conformation lasted shorter time than the b-turns, which is consistent with findings in previous studies that the hairpin conformation may not be the only structure of the peptide21 and that cooperative folding remains a rare event in the simulation.22 The helix and hairpin folding simulations show that without hydrogen bonding peptide structures are amino acid sequence dependent based on VDW interactions and covalent bonding constraints.
Discussion The electric polarity of water molecules and peptide groups has been the main focus of the protein folding theory. Without hydrogen bonding and other charge–charge interactions, what forces drive secondary structure folding in simulations? The covalently bound groups of the peptide chain have very different thermal motions from those of individual water molecules, including the translational motion, the rotational motion, and vibrational frequencies and amplitudes. Water molecules near the peptide chain are restricted by the covalently connected chain structure to move differently from those in bulk water, and consequently push the peptide chain toward intramolecular interactions forming compact structures. These intermolecular collisions in the thermal motion are repulsive VDW interactions. Furthermore, peptide covalent bonds are shorter than intermolecular distances. Within a given cutoff distance, the attractive VDW forces between two peptide segments include more atom pairs contributing to stronger attractions than those of each segment with water, which leads to peptide compact structure formation. Therefore, both repulsive and attractive VDW interactions contribute to peptide intramolecular interactions leading to compact structures. This rationalization agrees with the VDW dynamic hydration shell proposal.6,7 The compact structures are amino acid sequence dependent and closely resemble the a-helix and the b-hairpin, as a consequence of VDW interactions and covalent bond constraints, including torsional constraints. Both hydrogen bonding and VDW interactions have been considered important driving forces for compact packing.24 This MD study demonstrates that VDW interactions with covalent constraints are sufficient to fold these peptides into the proximity of their known structures. Hydrogen bonding is a short range interaction and it locks the
Sung
approximate structure into the specific secondary structure. In contrast to standard molecular simulations where the total interaction energy is dominated by atomic charge interactions, these results without atomic charges and hydrogen bonding provide the evidence for a new view of the protein folding mechanism.
Materials and Methods MD simulations were carried out without atomic charges and without any hydrogen bonding potential functions for all atoms of the peptide and water molecules. VDW parameters and covalent bonding constraints of AMBER14 and CHARMM18 force fields with the TIP3P water model were tested. Periodic boundary conditions were included in constant vol˚, ume MD simulations with the cutoff distance of 8 A the time step of 0.002 picoseconds (ps), and the bond length connecting to hydrogen atoms constrained using SHAKE algorithm.25
Acknowledgments The author is grateful to Dr. R. L. Baldwin and Dr. G. Matters for helpful discussions and valuable suggestions.
References 1. Pauling L, Corey RB, Branson HR (1951) The structure of proteins: two hydrogen-bonded helical configurations of the polypeptide chain. Proc Natl Acad Sci USA 37: 205–211. 2. Pauling L, Corey RB, Branson HR (1951) Configurations of polypeptide chains with favored orientations around single bonds: two new pleated sheets. Proc Natl Acad Sci USA 37:729–740. 3. Nelson DL, Cox MM (2008) Lehninger principles of biochemistry, 5th ed. New York: W. H. Freeman, pp 114– 115. 4. Berg JM. Tymoczko JL, Stryer L (2002) Biochemistry, 5th ed. New York: W. H. Freeman, pp 14. 5. De Santis P, Giglio E, Liquori AM, Ripamonti A (1965) Van der waals interaction and the stability of helical polypeptide chains. Nature 206:456–458. 6. Baldwin RL (2014) Dynamic hydration shell restores Kauzmann’s 1959 explanation of how the hydrophobic factor drives protein folding. Proc Natl Acad Sci USA 111:13052–13056. 7. Jorgensen WL, Gao J, Ravimohan C (1985) Monte Carlo simulations of alkanes in water: hydration numbers and the hydrophobic effect. J Phys Chem 89:3470– 3473. 8. Anil B, Sato S, Cho JH, Raleigh DP (2005) Fine structure analysis of a protein folding transition state: distinguishing between hydrophobic stabilization and specific packing. J Mol Biol 354:693–705. 9. Bachmann A, Wildemann D, Oraetorius F, Fischer G, Kiefhaber T (2011) Mapping backbone and side-chain interactions in the transition state of a coupled protein folding and binding reaction. Proc Natl Acad Sci USA 108:3952–3957. 10. Schulz GE, Schirmer RH (1990) Principles of protein structure. New York: Springer-Verlag, pp 30–33.
PROTEIN SCIENCE VOL 24:1383—1388
1387
11. Sung SS (1994) Helix folding simulations with various initial conformations. Biophys J 66:1796–1803. 12. Wang H, Sung SS (2000) Molecular dynamics simulations of three-strand b-sheet folding. J Am Chem Soc 122:1999–2009. 13. Scholtz JM, York EJ, Steward JM, Baldwin RL (1991) A natural water-soluble, a-helical peptide: the effect of ionic strength on the helix-coil equilibrium. J Am Chem Soc 113:5102–5104. 14. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Jr., Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA (1995) A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc 117:5179–5197. 15. Pearlman DA, Case DA, Caldwell JW, Ross WS, Cheatham TE, III, DeBolt S, Ferguson D, Seibel G, Kollman PA (1995) AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comp Phys Commun 91:1–41. 16. Daggett V, Levitt M (1992) Molecular dynamics simulations of helix denaturation. J Mol Biol 223:1121–1138. 17. Best RB, Hummer G (2009) Optimized molecular dynamics force fields applied to the helix-coil transition of polypeptides. J Phys Chem 113:9004–9015.
1388
PROTEINSCIENCE.ORG
18. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M (1983) CHARMM: a program for macromolecular energy, minimization, and dynamics calculations. J Comp Chem 4:187–217. 19. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kale L, Schulten K (2005) Scalable molecular dynamics with NAMD. J Comput Chem 26:1781–1802. 20. Humphrey W, Dalke A, Schulten K (1996) VMD—visual molecular dynamics. J Mol Graph 14:33–38. 21. Blanco FJ, Jimenez MA, Herranz J, Rico M, Santoro J, Nieto JL (1993) NMR evidence of a short linear peptide that folds into a b-hairpin in aqueous solution. J Am Chem Soc 115:5887–5888. 22. Wu X, Wang S, Brooks BR (2002) Direct observation of folding and unfolding of a b-hairpin in explicit water through computer simulation. J Am Chem Soc 124:5282–5283. 23. Wu X, Wang S (1998) Self-guided molecular dynamics simulation for efficient conformational search. J Phys Chem B 102:7238–7250. 24. Dill KA, Ozkan SB, Shell MS, Weikl TR (2008) The protein folding problem. Annu Rev Biophys 37:289–316. 25. Ryckaert JP, Ciccotti G, Berendsen HJC (1997) Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of nalkanes. J Comput Phys 23:327–341.
Driving Forces of Peptide Folding
Shen-Shu Sung* Department of Pharmacology, College of Medicine, Penn State University, Hershey, Pennsylvania 17033 Received 20 December 2014; Accepted 13 May 2015 DOI: 10.1002/pro.2710 Published online 26 May 2015 proteinscience.org
Abstract: Contrary to the widespread view that hydrogen bonding and its entropy effect play a dominant role in protein folding, folding into helical and hairpin-like structures is observed in molecular dynamics (MD) simulations without hydrogen bonding in the peptide-solvent system. In the widely used point charge model, hydrogen bonding is calculated as part of the interaction between atomic partial charges. It is removed from these simulations by setting atomic charges of the peptide and water to zero. Because of the structural difference between the peptide and water, van der Waals (VDW) interactions favor peptide intramolecular interactions and are a major contributing factor to the structural compactness. These compact structures are amino acid sequence dependent and closely resemble standard secondary structures, as a consequence of VDW interactions and covalent bonding constraints. Hydrogen bonding is a short range interaction and it locks the approximate structure into the specific secondary structure when it is included in the simulation. In contrast to standard molecular simulations where the total energy is dominated by chargecharge interactions, these simulation results will give us a new view of the folding mechanism. Keywords: intramolecular interaction; compact structure; secondary structure; covalent bonding constraint; hydrogen bonding; folding mechanism Introduction Historically, protein intramolecular hydrogen bonding was a determining factor in the original prediction of protein secondary structures.1,2 Currently, the hydrophobic effect is considered to be the most important single factor in protein folding in our textbooks,3,4 and the hydrophobic effect is considered to arise from entropy changes in water hydrogen bonding around water-exposed protein nonpolar groups. However, the helical structure could be predicted based on the energy minimum of peptide intramolecular VDW interactions.5 The entropy and enthalpy effect of VDW interactions contribute importantly to protein folding in the VDW dynamic hydration shell proposal.6,7 Studies also showed that hydrophobic interactions are an important driving force for forming the transition state8 and for molecular recogniGrant sponsor: College of Medicine, Penn State University. *Correspondence to: Shen-Shu Sung, Department of Pharmacology, R130, Penn State College of Medicine, 500 University Drive, Hershey, PA 17033-0850. E-mail: [email protected] or [email protected]
C 2015 The Protein Society Published by Wiley-Blackwell. V
tion.9 It may be difficult to separate these contributing factors in laboratory tests, but it is feasible to change or remove hydrogen bonding effects in computations. In the widely applied atom-based point charge model, each atom is modeled by a point in space with a partial charge representing the molecular charge distribution and with a 6–12 VDW potential representing dispersion forces and electron shell repulsions. VDW interactions and covalent constraints of bond lengths, bond angles, and dihedral angles are necessary for keeping atomic distances realistic and their parameters are largely derived from experimental structural data. Atomic partial charges are mainly based on calculated values.10 Hydrogen bonding is calculated as part of Coulomb interactions between atomic charges. Setting the atomic charges to zero is a simple way to hypothetically remove hydrogen bonding and other charge– charge interactions from a computation, and it is tested in this study focusing on the qualitative nature of the peptide folding mechanism.
PROTEIN SCIENCE 2015 VOL 24:1383—1388
1383
Figure 1. From the left, the helix turn formation observed at 0.12 ns, a helix-turn-helix at 0.94 ns, an approximate helix at 6 ns, and a partial helix at 10 ns during the first simulation with AMBER force field parameters without atomic charges, and a helical structure at 17.75 ns during the simulation with CHARMM parameters without atomic charges. The carbon atoms are colored green, oxygen red, nitrogen blue, and polar hydrogen white. Water molecules and hydrogen atoms bound to carbon atoms are not shown. The molecular graphics was generated using Maestro display software from Schrodinger LLC.
Results Helix folding Early peptide folding simulations11,12 were carried out successfully on synthetic peptide sequences. Among them, a well characterized alanine-based ahelical peptide13 Ac(AAQAA)3Y-NH2 is studied here in helix folding simulations without atomic partial charges of the peptide and water molecules. The first simulation was carried out using the AMBER force field14 covalent bonding constraints and VDW parameters. The initial structure was a fully extended peptide in a solvent box of 4695 water molecules, built using the AMBER software package.15 This simulation was carried out at 273 K for 20 nanoseconds (ns). Helical turns formed at both ends of the peptide at 0.12 ns, and tertiary interactions between helical segments were observed between 0.6 and 1.1 ns, as shown in Figure 1. Helical structures with partially unfolded N- or C-termini were observed repeatedly during several periods of time, including 2.1-2.3 ns, 3.8-4.3 ns, and 4.9-5.6 ns. More complete helical structures were observed between 5.8 and 6.2 ns. The structure at 6.0 ns in Figure 1 is an example of an approximate helical structure. Partial unfolding and refolding were observed repeatedly. The structure at 10 ns in Figure 1 is a partially unfolded helical structure. Approximate helical structures were observed again at 11.4–11.6 ns, 16.3–16.7 ns, and 17.5–17.9 ns. Helix-turn-helixlike structures were observed in several periods of time, including 10.3–10.6 ns and 19.7–19.9 ns.
1384
PROTEINSCIENCE.ORG
Two views of a helical structure observed at 6.05 ns are shown in Figure 2, one perpendicular to the helix axis and the other along the axis from the C-terminus. Its backbone amide oxygen and hydro˚ distance are connected by gen atoms within 2.6A dotted lines. There are 10 (i, i14) oxygen-hydrogen pairs corresponding to the a-helix structure. There
Figure 2. Two views of the helical structure observed at 6.05 ns. One or two a-carbon atoms are labeled in each turn for the visualization of the number of residues per turn. The view perpendicular to the helix axis shows backbone oxygen and hydrogen atom pairs within 2.6 A˚, including ten (i, i14) pairs corresponding to a-turns and four (i, i13) pairs corresponding to b-turns, including 310 helix turns. For clarity, sidechains are not shown in the view along the helical axis from the Cterminus (lower right). Atom colors and other descriptions are the same as in Figure 1.
Driving Forces of Peptide Folding
are 4 (i, i13) pairs, 2 of them are b-turns near both ends of the helix, and the other two share atoms with (i, i14) pairs, indicating possible 310 helix turns. The oxygen-hydrogen atom pairs not connected by the dotted lines have their distances ˚ . One or two a-carbon atoms per larger than 2.6 A turn are labeled for visualizing the number of residues per turn, which is close to 3.6 residues per turn of the a-helix. The pitch of the helix is between ˚ , depending on the atom pairs used in 5.3 and 5.5 A the distance measurement. Without hydrogen bonds connecting the backbone hydrogen and oxygen atoms, amide groups tilt with oxygen atoms outwards, as shown in both views in Figure 2. This amide tilt makes the oxygen-hydrogen distance larger than that in a standard a-helix. During this simulation, atomic coordinates were recorded and the backbone /, w angles calculated every 10 ps. Several backbone dihedral angle criteria have been used in literatures16,17 for calculating the helical content or the helical ratio, such as / 5 2578 6 308 and w 5 2478 6 308, / 5 2658 6 358 and w 5 242.58 6 37.58, / 5 2658 6 358, and w 5 2378 6 308, etc. In the current study, when the /, w angles of two or more consecutive amino acid residues are within 308 of the standard a-helical angles (/ 5 2578 and w 5 2478), these residues are assumed to be in the helical conformation, which includes the majority of a-helix and 310 helix residues. The ratio of the number of helical residues to the total number of the residues (16 residues of the alanine based peptide) as a function of time in the simulation is shown in Figure 3. The average helical ratio during the 20 ns is 18.66%. In comparison, a standard MD simulation was carried out with atomic charges, which shows a higher helical ratio of 35.03%, indicating that hydrogen bonds stabilize helical structures. Amide tilting changes the backbone dihedral angles from those of the standard helix, which is another reason for the lower helical ratio without hydrogen bonding. To confirm these qualitative observations, multiple simulations with different initial structures were carried out. The initial structures were unfolded coil structures, generated from a high temperature (900 K) simulation at 10 ps time intervals. AMBER force field parameters were used without atomic charges of the peptide and water molecules. A group of 16 simulations was carried out for 20 ns each at 273 K with 1529 solvent water molecules. Two helical segments with tertiary interactions forming helix-turn-helix-like structures were observed in all these simulations, and the peptide formed single helices in 11 of the 16 simulations. The results of multiple simulations are qualitatively in agreement with those of the previous simulation with the extended initial structure. The helical ratio of each simulation averaged over the 20 ns ranges from
Sung
Figure 3. The helical ratio (the number of residues in helical conformations to the total number of residues of the peptide) as a function of time in simulations with AMBER force field parameters. The solid line shows the helical ratio in the simulation with the extended initial structure and without atomic charges, the dashed line shows the average helical ratio of 16 simulations with various unfolded initial structures and without atomic charges, and the dotted line shows the helical ratio in a standard MD simulation with the extended initial structure and with atomic charges.
10.15% to 20.21%. The helical ratio of the previous simulation (18.66%) is within this range. The average helical ratio of all 16 simulations as a function of time is shown in Figure 3. To see the force field dependence of these observations, CHARMM force field18 covalent bonding constraints and VDW parameters have been applied to a simulation of the same alanine-based peptide without atomic charges of the peptide and water molecules. This simulation was carried out at 273K for 20 ns using the NAMD program.19 The initial structure was a fully extended peptide in a solvent box of 4,427 water molecules, built using the VMD program.20 During the first 10 ns, the peptide was in random coil conformations with isolated helical turns. After the first 10 ns, various helical structures were observed, including partial helices and helix-turn-helix-like structures. Between 17.20 and 17.75 ns, the peptide formed a single helix. One of these structures observed at 17.75 ns is shown in Figure 1. Compared with the structure with AMBER parameters, the backbone amide tilt is more pronounced with CHARMM parameters. The periodicity of the turns is between 3.5 and 4 residues. The pitch ˚ , because of the larger of the helix is shorter at 5.2 A amide tilt. Helix folding occurred with different force field parameters without peptide intramolecular and water hydrogen bonding.
Hairpin folding Without hydrogen bonding, is the secondary structure amino acid sequence dependent? Besides the ahelix, the b-sheet is another common secondary structure. With limited computing resources, a small peptide is the first choice for a folding simulation. Blanco et al.21 have successfully designed a small b-
PROTEIN SCIENCE VOL 24:1383—1388
1385
Figure 4. A hairpin-like structure observed in the simulation using CHARMM force field parameters without atomic charges (left) and that using AMBER parameters without atomic charges (right), and a b–hairpin structure using CHARMM parameters with a fraction of atomic partial charges (center). Dotted lines connect the potential hydrogen bonding atom pairs with their distances shown in A˚. Atom colors and other descriptions are the same as in Figure 1.
hairpin peptide YQNPDGSQA, and Wu et al.22 have carried out a folding simulation of this b-hairpin peptide using the SGMD method23 with explicit solvent. Can this peptide fold into a hairpin structure without hydrogen bonds? A MD simulation was carried out with CHARMM force field parameters18 without atomic charges of the peptide and water molecules. The initial structure was a fully extended peptide in a solvent box of 1002 water molecules, built using the VMD program.20 The simulation was carried out at 300 K for 10 ns using the NAMD program.19 During the simulation, structures were recorded every 5 ps. Various transient hairpin-like structures were observed during the time period of 4,320 to 4,360 ps. A structure at 4340 ps is shown in Figure 4. Hydrogen bonding atom pairs indicated by NOE data21 are connected by dotted lines. Three out of the four distances between the hydrogen and the ˚ , about 0.5 A ˚ oxygen atoms are 2.93, 2.47, and 2.58 A larger than normal hydrogen bond distances. The backbone NH group of the amino acid A9 is pointing away from the Y1 carbonyl group. To test the effect of hydrogen bonding on the hairpin structure, a short simulation is carried out starting from the hairpin-like structure with a fraction (1/2) of atomic partial charges of the CHARMM18 force field. Within 50 ps, all four hydrogen bonds formed with the hydrogen-oxygen distan˚ , consistent with the NOE ces smaller than 2.6 A 21 data, as shown in Figure 4. These results show hairpin-like structures can form rapidly as a consequence of VDW interactions and covalent constraints, and hydrogen bonding locks them into the specific secondary structure. This observation may suggest a computational approach using mainly VDW interactions to direct folding into compact structures and including certain charge-charge interactions to stabilize the folded structure.
1386
PROTEINSCIENCE.ORG
AMBER force field14 parameters have also been tested in hairpin folding simulations. 20 unfolded structures were taken from a simulation at 900 K in 10 ps time intervals as initial structures of the following hairpin folding simulations. Each simulation was carried out for 20 ns at 300 K with 1,273 water molecules, using AMBER force field14 parameters without atomic charges of the peptide and water molecules. During these simulations b-turns were observed repeatedly, but they are not always flanked by antiparallel b-strands forming the hairpin structure. For example, in the 12th simulation hairpinlike structures were observed at 0.80–1.25 ns, but the two strands moved apart after 1.25 ns while the b-turn at the P4-D5 remained. A hairpin-like structure observed at 1.14 ns in this simulation is included in Figure 4. In this structure, three potential hydrogen bonding atom pairs are at distances ˚ , but the N-terminal NH3 between 2.5 and 3.5 A group does not point towards the oxygen atom in the C-terminus. Using the structures saved every 5 ps, the / and w angles were calculated. Standard type I bturn angles are /i11 5 2608, wi11 5 2308, / i12 5 2908, and w i12 5 08. The difference of these /, w angles between i11 and i12 positions is 308. The largest dihedral angle difference between the type I b-turn and the a-turn is 478 in wi12. To clearly distinguish between the a-turn and the b-turn and between the i11, i12 positions of the type I b-turn, the type I, II, I’, and II’ b-turns and the a-turn are assigned using a smaller variation range of 6158 from their standard /, w angles. The percentage of a type of turns at a (i11, i12) site in the total number of structures recorded during the simulation is the ratio of that type of turns at the given site. The type I b-turn ratios at i 1 1 5 2 to 7 (i 1 2 5 3 to 8) sites (the 20 simulation average) are 0.000, 0.117, 3.625,
Driving Forces of Peptide Folding
2.320, 0.870, and 1.420%, and helical ratios at the same sites are 0.052, 0.085, 0.695, 0.273, 0.140, and 0.442%. The first and the last residues do not have sufficient numbers of atoms for their dihedral angle calculations. The type II, I’, II’ b-turn ratios are all lower than 0.1%. The highest turn ratio is 3.625% of the type I b-turn at the P4-D5 site, consistent with the model structure based on NOE data.21 During the simulation, the b-hairpin conformation lasted shorter time than the b-turns, which is consistent with findings in previous studies that the hairpin conformation may not be the only structure of the peptide21 and that cooperative folding remains a rare event in the simulation.22 The helix and hairpin folding simulations show that without hydrogen bonding peptide structures are amino acid sequence dependent based on VDW interactions and covalent bonding constraints.
Discussion The electric polarity of water molecules and peptide groups has been the main focus of the protein folding theory. Without hydrogen bonding and other charge–charge interactions, what forces drive secondary structure folding in simulations? The covalently bound groups of the peptide chain have very different thermal motions from those of individual water molecules, including the translational motion, the rotational motion, and vibrational frequencies and amplitudes. Water molecules near the peptide chain are restricted by the covalently connected chain structure to move differently from those in bulk water, and consequently push the peptide chain toward intramolecular interactions forming compact structures. These intermolecular collisions in the thermal motion are repulsive VDW interactions. Furthermore, peptide covalent bonds are shorter than intermolecular distances. Within a given cutoff distance, the attractive VDW forces between two peptide segments include more atom pairs contributing to stronger attractions than those of each segment with water, which leads to peptide compact structure formation. Therefore, both repulsive and attractive VDW interactions contribute to peptide intramolecular interactions leading to compact structures. This rationalization agrees with the VDW dynamic hydration shell proposal.6,7 The compact structures are amino acid sequence dependent and closely resemble the a-helix and the b-hairpin, as a consequence of VDW interactions and covalent bond constraints, including torsional constraints. Both hydrogen bonding and VDW interactions have been considered important driving forces for compact packing.24 This MD study demonstrates that VDW interactions with covalent constraints are sufficient to fold these peptides into the proximity of their known structures. Hydrogen bonding is a short range interaction and it locks the
Sung
approximate structure into the specific secondary structure. In contrast to standard molecular simulations where the total interaction energy is dominated by atomic charge interactions, these results without atomic charges and hydrogen bonding provide the evidence for a new view of the protein folding mechanism.
Materials and Methods MD simulations were carried out without atomic charges and without any hydrogen bonding potential functions for all atoms of the peptide and water molecules. VDW parameters and covalent bonding constraints of AMBER14 and CHARMM18 force fields with the TIP3P water model were tested. Periodic boundary conditions were included in constant vol˚, ume MD simulations with the cutoff distance of 8 A the time step of 0.002 picoseconds (ps), and the bond length connecting to hydrogen atoms constrained using SHAKE algorithm.25
Acknowledgments The author is grateful to Dr. R. L. Baldwin and Dr. G. Matters for helpful discussions and valuable suggestions.
References 1. Pauling L, Corey RB, Branson HR (1951) The structure of proteins: two hydrogen-bonded helical configurations of the polypeptide chain. Proc Natl Acad Sci USA 37: 205–211. 2. Pauling L, Corey RB, Branson HR (1951) Configurations of polypeptide chains with favored orientations around single bonds: two new pleated sheets. Proc Natl Acad Sci USA 37:729–740. 3. Nelson DL, Cox MM (2008) Lehninger principles of biochemistry, 5th ed. New York: W. H. Freeman, pp 114– 115. 4. Berg JM. Tymoczko JL, Stryer L (2002) Biochemistry, 5th ed. New York: W. H. Freeman, pp 14. 5. De Santis P, Giglio E, Liquori AM, Ripamonti A (1965) Van der waals interaction and the stability of helical polypeptide chains. Nature 206:456–458. 6. Baldwin RL (2014) Dynamic hydration shell restores Kauzmann’s 1959 explanation of how the hydrophobic factor drives protein folding. Proc Natl Acad Sci USA 111:13052–13056. 7. Jorgensen WL, Gao J, Ravimohan C (1985) Monte Carlo simulations of alkanes in water: hydration numbers and the hydrophobic effect. J Phys Chem 89:3470– 3473. 8. Anil B, Sato S, Cho JH, Raleigh DP (2005) Fine structure analysis of a protein folding transition state: distinguishing between hydrophobic stabilization and specific packing. J Mol Biol 354:693–705. 9. Bachmann A, Wildemann D, Oraetorius F, Fischer G, Kiefhaber T (2011) Mapping backbone and side-chain interactions in the transition state of a coupled protein folding and binding reaction. Proc Natl Acad Sci USA 108:3952–3957. 10. Schulz GE, Schirmer RH (1990) Principles of protein structure. New York: Springer-Verlag, pp 30–33.
PROTEIN SCIENCE VOL 24:1383—1388
1387
11. Sung SS (1994) Helix folding simulations with various initial conformations. Biophys J 66:1796–1803. 12. Wang H, Sung SS (2000) Molecular dynamics simulations of three-strand b-sheet folding. J Am Chem Soc 122:1999–2009. 13. Scholtz JM, York EJ, Steward JM, Baldwin RL (1991) A natural water-soluble, a-helical peptide: the effect of ionic strength on the helix-coil equilibrium. J Am Chem Soc 113:5102–5104. 14. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Jr., Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA (1995) A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc 117:5179–5197. 15. Pearlman DA, Case DA, Caldwell JW, Ross WS, Cheatham TE, III, DeBolt S, Ferguson D, Seibel G, Kollman PA (1995) AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comp Phys Commun 91:1–41. 16. Daggett V, Levitt M (1992) Molecular dynamics simulations of helix denaturation. J Mol Biol 223:1121–1138. 17. Best RB, Hummer G (2009) Optimized molecular dynamics force fields applied to the helix-coil transition of polypeptides. J Phys Chem 113:9004–9015.
1388
PROTEINSCIENCE.ORG
18. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M (1983) CHARMM: a program for macromolecular energy, minimization, and dynamics calculations. J Comp Chem 4:187–217. 19. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kale L, Schulten K (2005) Scalable molecular dynamics with NAMD. J Comput Chem 26:1781–1802. 20. Humphrey W, Dalke A, Schulten K (1996) VMD—visual molecular dynamics. J Mol Graph 14:33–38. 21. Blanco FJ, Jimenez MA, Herranz J, Rico M, Santoro J, Nieto JL (1993) NMR evidence of a short linear peptide that folds into a b-hairpin in aqueous solution. J Am Chem Soc 115:5887–5888. 22. Wu X, Wang S, Brooks BR (2002) Direct observation of folding and unfolding of a b-hairpin in explicit water through computer simulation. J Am Chem Soc 124:5282–5283. 23. Wu X, Wang S (1998) Self-guided molecular dynamics simulation for efficient conformational search. J Phys Chem B 102:7238–7250. 24. Dill KA, Ozkan SB, Shell MS, Weikl TR (2008) The protein folding problem. Annu Rev Biophys 37:289–316. 25. Ryckaert JP, Ciccotti G, Berendsen HJC (1997) Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of nalkanes. J Comput Phys 23:327–341.
Driving Forces of Peptide Folding
