A. Performance of Force Fields and Solvent Models
Our assessment of each of the force fields is based on the extent to which they maintain correct basic hairpin and helical tendencies, under extensive conformational sampling. Because experimental measurements on these short peptides are difficult to convert directly into structural ensembles, we do not believe a more detailed analysis–comparing the fractional stabilities, for example–is particularly informative here. Instead, we focus primarily on the degree of balance between helical and β motifs, by examining the dominant structures produced in each simulation and comparing to a single, presumed native structure. Though certainly the peptides experience fluctuations away from these presumed structures in reality, these conformations still serve as a useful milepost for indicating proximity to the correct topology.
Our main results are shown in – and in . shows the dominant conformations of the EK peptide after REMD simulation. These conformations come from the most populous cluster after the clustering algorithm is performed on the trajectory; of the three run repeats for each condition, the one whose top cluster contained the most secondary structure is shown. The simulations with ff96, ff99, and ff99SB with igb1 or igb5 perform well, all retaining a near-perfect helix. The strong helical tendency of ff99 has been seen before.45
On the other hand, whereas ff96 was previously found to give poor helices in explicit solvent,17,50
here we find that ff96 works well with these implicit solvent models. The difference may be due to the fact that the implicit models increase the helical propensity of ff96, relative to explicit water. Another difference may be greater conformational sampling in the 10 ns simulations using implicit solvent, which lacks viscosity as compared to previous explicit water simulations.
Figure 1 Top cluster conformations for the EK peptide from replica exchange simulations. The conformation with the most secondary structure among the dominant clusters of the three repeat runs for each force field/solvent combination is shown. Colors are residue-specific: (more ...)
Figure 3 (left) Double-error plot of average backbone rmsd for the EK peptide and G hairpin for the last 1 ns of the 10 ns REMD simulations. The points represent the average over three independent REMD runs; the error bars show the maximum and minimum averaged (more ...)
Summary of Simulation Data Calculated over the Final 1 ns of Each REMD Simulationa
Other force field/solvent combinations are problematic. The ff03 force field unfolds the helix completely, forming a coil-dominated structure with igb5 and a hairpin with igb7. These results are consistent among all trials, indicating at least minor difficulties in achieving α-helical structures in ff03 with implicit solvent. A second weakness appears to be the igb7 solvent model. The results for igb7 consistently show significant deviations from helix: starting from the native structure, the helix essentially unfolds when igb7 is used with all force fields except ff99SB, where it partially unfolds. The variety of structures produced with different force fields suggests that igb7 does not adequately stabilize helical structures.
shows our tests on the protein G hairpin. Here, successful predictions appear to depend more on the force field than on the solvation model. It is clear that ff99 erroneously tends toward helical conformations, as previous researchers have seen. ff99SB gives slightly better structures than ff99, although the hairpin never fully develops in the former, and some helical tendency is still evident for igb5. However, ff96 performs well for all three solvent models, consistent with findings of Okamoto et al.17
The ff03 parameters also fare well; the particular combination of ff03 and igb5 generates an extremely stable hairpin with near-perfect alignment. Compared to its poor performance on the helical EK peptide, this might indicate that ff03’s dihedral torsions are slightly biased toward β
-structure when used with implicit solvent.
Top cluster conformations for the protein G hairpin from replica exchange simulations. See for further details.
In these visual assessments of structure, force field ff96 appears the most consistent and balanced. We find that igb5 performs well with it, and that igb1 is almost as good, but igb7 is problematic. presents a more quantitative comparison using a type of plot that we have found to be informative. This double-error plot shows the rmsd error in predicting the helix vs the rmsd error in predicting the β-hairpin. The rmsd is time-averaged over the last 1 ns of the simulation trajectories. This kind of diagnostic does not give detailed information about the ensemble of structures sampled in each simulation, but we believe it is at least a good indication of the basic extent to which helical and hairpin structures are sampled. Therefore, the plot in gives two simple types of information. First, points closer to the main diagonal, the y = x line, indicate that a force field is well-balanced, which is to say that errors in predicting β-hairpins are roughly equal to errors in predicting α-helical structures. Second, points closer to the origin indicate smaller overall errors of the model relative to the presumed native structures compared to points more distant from the origin. Hence the best force fields and solvent models give points on this graph that are closest to the main diagonal and closest to the origin.
By these criteria, the best overall force field/solvation pairs are ff96 + igb1, ff96 + igb5, and ff03 + igb5. We note that it is difficult to discriminate further, using small differences in average rmsd values (~1 Å), due to both the imperfect nature of the rmsd measurement to a presumed native structure and the run-to-run variations. However, we find that globally the most consistent force field, with multiple solvents, appears to be ff96. Similarly, we find igb5 to be the most consistent solvent when paired with multiple force fields. Therefore, we consider the choice of ff96 + igb5 to be a reasonable overall choice, with rmsd’s centered around 2.9 Å for the helix and 3.9 Å for the hairpin.
The error bars in show the minimum and maximum rmsd values of the three run repeats. Much of this variance is likely due to structural fluctuations around the stable conformation in each case, on the order of ~ 1–2 Å. On the other hand, two runs in particular exhibit fairly pronounced fluctuations: those involving ff96/igb1 and a modified version of ff96/igb5 that is described later in the paper. For those two cases, there is significant run-to-run variation among the three trials, and examining these further, we also found that each individual run could be characterized by significant structural fluctuations as the simulations progressed with time. This kind of behavior is typical of a rough folding landscape.
As a reference comparison, shows the correlation between averaged rmsd (as reported above) and the time-averaged fraction of native backbone hydrogen bonds formed in each peptide, compiled over all simulation test cases in our study. Both the EK peptide and protein G hairpin were reported to remain only partially stable in solution. Though a detailed comparison of structural ensembles with experiment is difficult (since it is not clear whether the reported fractional stabilities correspond to partial unfolding or temporal fluctuations), we can ask what kinds of fluctuations emerge from partially stable structures in the simulations. For the instructive case in which 50% of the native hydrogen bonding is present, suggests rmsd averages of 2.0 Å for the EK helix and 3.0 Å for the G hairpin. Thus, it is not unreasonable to expect that partially stable peptides experience rmsd fluctuations around a native structure that are on the order 2–3 Å, and that the fluctuations might be greater for the hairpin than the helix.
Due to the finite length of these simulations, it is indeed possible that equilibrium convergence is not fully attained. However, we believe that the dramatic differences among the systems studied in the 10 ns REMD runs provides a useful force field diagnostic. To illustrate some of these differences, shows some of the convergence properties of the simulations, showing a well- and poor-performing force field/solvent pair for each of the two peptides. A time progression of the rmsd from presumed native structure shows that the poor-performing pairs quickly depart the native basin, in under 1 ns of REMD simulation time. For these, the populations of the dominant clusters increase with time, indicating that the peptides find increasingly stable alternate conformations, away from the native. In contrast, the best-performing pairs largely remain a short distance from the native basin, although with time it appears that the rmsd fluctuations for these slightly increase. Moreover, in the well-performing simulations, the dominant cluster populations decrease during the run relative to the initial 2 ns. This could indicate that the systems locate slightly more diverse conformational ensembles than the native structure alone. On the time scale of the run, it appears these conformational populations start stabilizing around the 8 ns time point.
Figure 4 Comparison of simulation progress between the best- and worst-performing force fields. The top panel corresponds to the EK peptide and the bottom to the protein G hairpin. Continuous thin lines depict the rmsd to the presumed native structure, whereas (more ...)
B. Salt-Bridge Effects on Peptide Stability
Several research groups have previously found that implicit solvent models give non-native salt-bridge-trap conformations.51
We looked for salt-bridge traps in our simulations. The G hairpin contains a single lysine residue that can pair with either of two aspartate or two glutamate residues. One of these four potential salt bridges, Lys7 and Asp10, does not disrupt the native hairpin conformation. compares the potential of mean force (PMF) for this native salt bridge with different solvent models for ff96. These PMF curves were extracted from the last 2 ns of each simulation, and averaged over three runs. In the strictest sense, these curves must be taken as approximate given the finite simulation durations and the limitations of the 2 ns data window, and are best taken as a measure of probabilities during the sampling procedure. However, error bars measured from differences among the independent runs do give an indication as to the variation of these computations, and they show that there are some statistically distinct features in these measurements for different force fields. Note that each force field/solvent pair produces a minimum around 2.8 Å corresponding to the formation of the salt bridge and a shallower minimum around 4.5 Å due to the second charged aspartate oxygen.
Figure 5 Potential of mean force between native ion-pairing side-chain atoms in the protein G hairpin. The distance is calculated between the charged N in Lys7 and one charged O in Asp10; the two PMFs for each of the aspartate oxygens are averaged. The PMFs were (more ...)
An analysis of the PMF curves shows the relative populations of ion pairs in each of the three solvents tested. The main minima in the PMFs in igb1 and igb5 are roughly equal in depth, indicating that the salt bridge is formed roughly the same fraction of time in both solvents. The heights of the barriers between igb1 and igb5 are different, however, which suggests that salt bridges in igb1 less readily dissociate once formed. In the igb7 case, the higher absolute position of the minimum means that salt bridges are less populated overall. Qualitatively, our simulations show a ranking of igb1, igb5, and igb7 from strongest to weakest salt-bridge stabilization in the protein G hairpin, in terms of populations of different salt-bridge states.
It is interesting to consider these results in light of the recent work of Lwin and Luo24
that found that salt-bridge populations are also significantly affected by the choice of force field. Indeed, our results show that the other force fields, beyond ff96, alter the strength of salt-bridge formation; in many of these cases, however, the effect is indirect, through a major altering of the secondary structure preference of the peptide. Still, it is challenging to rigorously separate force field and implicit solvent effects on ion pairing in these simulations. Lwin and Luo also suggest that the nonpolar surface area term contributes significantly to ion pairing.25
Our simulations clearly show an effect of the electrostatic part alone, since each of the three GB models is paired with the same surface area term. A useful future study might examine the effect of turning off the surface area component of the implicit models.
We find in the ff96 + igb5 REMD trajectories that non-native salt bridges are convergence traps, persisting for 100–1000 ps in run time at the lowest temperature replica, but not stable over longer times. In the protein G hairpin, the most frequently seen such decoy configuration entailed a one-residue shift in the strand alignment to accommodate the non-native Lys7–Asp11 pairing. This shifted alignment was more frequent in the ff96 + igb1 case, consistent with the relative salt-bridge stabilization of igb1 and igb5. For the EK peptide, we found non-native salt-bridge pairing occasionally led to a hairpin conformation from which it was difficult to reach the native helical state.
We performed a simple test to establish whether the salt bridges contributed significantly to the stability of the native structures in the force field. We excluded the possibility of salt-bridge formation during the simulations by removing the associated configurations from the sampling space. In our molecular dynamics runs, we introduced a repulsive potential felt only upon close approach between any two ion-pairing atoms, steep enough to prevent salt bridges but not so steep that it introduced numerical instabilities. The form of the potential is shown in , chosen because it is an option in the AMBER 9 package. We chose parameters to create a steep repulsion that would mimic the effects of a hard overlap term: the potential is zero at large distances, rises harmonically from 6 to 4 Å with a force constant of 0.5 kcal/mol/Å2, and is continuously linear thereafter. For comparison, salt bridges typically form when the pairwise distance is less than 4 Å. We stress that this modification is not intended to be a reparametrization of the force field, but simply an easy way to exclude those portions of configuration space that correspond to salt-bridge formation. We have also verified that increasing the steepness of the repulsion does not have a statistically significant effect (data not shown).
Figure 6 Repulsive energy function used between potential ion-pairing heavy atoms to prevent salt-bridge formation. Decreasing from large distances, the potential is zero until r0, harmonic with a force constant k until r1, and continuously linear thereafter. (more ...)
We used the salt-bridge-free sampling on the best-performing force field, ff96 + igb5. shows the effect of removing the possibility of salt bridges on the PMF calculations for the protein G hairpin. The modification has two effects: as expected, it pushes the potential ion-pairing heavy atoms apart, but interestingly, it also increases the probability that the native contact is formed relative to the non-native salt bridge, as shown by the lower PMF for the former.
Figure 7 Comparison of salt-bridge interactions in the protein G hairpin for ff96 + igb5 with and without the salt-bridge-free sampling (using a repulsive potential to prevent salt bridges). The native interaction is Lys7–Asp10, while the non-native is (more ...)
In all of the peptides we examined, we found that removing the possibility of salt bridges from the sampling did not significantly destabilize the native secondary structures. shows a comparison of the structures for the EK helix, protein G hairpin, Baldwin’s C peptide, and the trpzip2 mini protein. All four peptides contain potential salt bridges. Even without the repulsion, the ff96 + igb5 force field/solvent model combination yields reasonable results across all four systems. With the salt-bridge sampling, similar cluster structures were obtained, and in some cases, we found that convergence was accelerated. One minor exception is the case of the protein G hairpin, which exhibited a large run-to-run variation in the kinds of structures sampled when the salt bridge was not allowed to form (): this could be an indication that the native basin was broadened by the removal of salt bridges, since it resulted in both structures closer to and farther from the normal ff96/igb5 case. Overall, however, the fact that the conformations of the most populated cluster for each peptide did not change significantly during salt-bridge-free sampling suggests that the very close, strong association between charged ions is not critical for stability in these peptides.
Figure 8 Top cluster conformations for four peptides from replica exchange simulations using the best-performing force field/solvent pair, 96 + igb5. The conformation with the most secondary structure among the dominant clusters of the three repeat runs for each (more ...)
We do caution that this salt-bridge-free sampling is a diagnostic, and not a general fix for ion-pairing overstabilization. In particular, it is possible that proteins with buried salt bridges, not investigated here, could be significantly destabilized by removal of salt-bridge interactions. Here, due to the small size of the systems studied, essentially every residue is exposed to solvent, and hence, all possible ion-pairing interactions are surface interactions. In contrast to this work, Simmerling and co-workers have attempted to address the ion-pairing overstabilization problem by rebalancing the force field, reducing the intrinsic (low-dieletric) radii of hydrogen atoms bound to charged nitrogen atoms.52
We also conducted a brief test of salt-bridge-free sampling on three larger proteins with surface salt-bridge pairs, to ensure it did not lead to unphysical artifacts in larger systems than peptides. We performed REMD simulations of three larger proteins, the 56-residue α/β protein G (PDB code 1GB1), 56-residue all-β src SH3 (1SRL), and the 46-residue all-α protein A (1BDD). Two simulations were performed for each, one with and one without the salt-bridge repulsion, and each run was initialized from the experimental native structure. We found that the rmsd-to-native of each protein after 10 ns of REMD remained the same in the salt-bridge-free case. The backbone rmsd’s for the final dominant cluster configurations in the normal/repulsion-applied cases were 1.2 Å/1.2 Å for protein G, 2.7 Å/2.6 Å for src SH3, and 2.2 Å/2.0 Å for protein A. Although such short REMD simulations for proteins of this size certainly do not reach equilibrium, these results do show that the added repulsion does not destabilize the native structure on this time scale.