Search tips
Search criteria 


Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
PLoS One. 2010; 5(4): e10056.
Published online 2010 April 7. doi:  10.1371/journal.pone.0010056
PMCID: PMC2850926

Evaluating Molecular Mechanical Potentials for Helical Peptides and Proteins

Jörg Langowski, Editor


Multiple variants of the AMBER all-atom force field were quantitatively evaluated with respect to their ability to accurately characterize helix-coil equilibria in explicit solvent simulations. Using a global distributed computing network, absolute conformational convergence was achieved for large ensembles of the capped A21 and Fs helical peptides. Further assessment of these AMBER variants was conducted via simulations of a flexible 164-residue five-helix-bundle protein, apolipophorin-III, on the 100 ns timescale. Of the contemporary potentials that had not been assessed previously, the AMBER-99SB force field showed significant helix-destabilizing tendencies, with beta bridge formation occurring in helical peptides, and unfolding of apolipophorin-III occurring on the tens of nanoseconds timescale. The AMBER-03 force field, while showing adequate helical propensities for both peptides and stabilizing apolipophorin-III, (i) predicts an unexpected decrease in helicity with ALA→ARG+ substitution, (ii) lacks experimentally observed 310 helical content, and (iii) deviates strongly from average apolipophorin-III NMR structural properties. As is observed for AMBER-99SB, AMBER-03 significantly overweighs the contribution of extended and polyproline backbone configurations to the conformational equilibrium. In contrast, the AMBER-99[var phi] force field, which was previously shown to best reproduce experimental measurements of the helix-coil transition in model helical peptides, adequately stabilizes apolipophorin-III and yields both an average gyration radius and polar solvent exposed surface area that are in excellent agreement with the NMR ensemble.


Simulating protein dynamics remains a daunting task in computational chemistry and biophysics. While many advances have been made in the last decade, processor speeds limit the timescales on which biomolecules can be simulated using explicit representations of aqueous solvent to the sub-microsecond timescale–far below the threshold of most interesting biomolecular events. More importantly, the ability of contemporary molecular models to accurately characterize the energetics of biopolymers depends greatly on the computational methods used to assess the models in question and the systems studied during such assessment. It has been shown that very subtle modifications to commonly used molecular mechanical potentials, such as changes to the scaling factors applied to non-bonded interactions [1], can significantly alter the behavior of those potentials with respect to stabilizing, or destabilizing, protein substructure. Most often, these modifications have aimed to improve upon the torsional potentials around the [var phi] and ψ dihedrals in protein backbones [2], which were fit to the relative quantum mechanical energies of alternate rotamers of small GLY and ALA peptides in Cornell's seminal AMBER-94 force field [3]: (a) AMBER-96 recalibrates these parameters to accurately predict energy differences between constrained and extended α-helical conformations of ALA peptides [4]; (b) AMBER-99 refits the parameterization using calculations on representative ALA tetrapeptides [5]; and (c) the Garcia-Sanbonmatsu variant of AMBER-94 (commonly referred to as AMBER-GS) zeroes both torsional potentials [6].

In an effort to properly assess the conformational preferences and equilibria of model systems simulated under the various force fields in use by the computational community, we thus set out to systematically study contemporary potentials in their application to one of the most ubiquitous and fundamental of protein substructures: the α-helix [1], [7], [8], [9]. These studies included the simulation of large ensembles of model helical peptides starting from fully helical and fully unfolded states to convergence of conformational equilibrium on the hundreds of nanoseconds timescale. Our results allowed absolute characterization of the equilibrium helical content and dynamics of many published molecular models including the AMBER-94, AMBER-96, AMBER-99, and AMBER-GS force fields. At that time, a new variant dubbed AMBER-99[var phi], which replaces the [var phi] potential in the helix-destabilizing AMBER-99 with that of the helix-friendly AMBER-94 force field, was shown to yield the best agreement with numerous experimental characterizations of the helical Fs peptide including the helix folding time, mean 310 helicity, mean residue dwell time in the coil state, mean radius of gyration, and Lifson-Roig (LR) parameters. AMBER-99[var phi] also gave the best agreement with quantum mechanical sampling of the alanine dimer and a survey of alanine conformations in the Protein Data Bank (PDB) [7].

In recent years, Duan and coworkers have published their third generation AMBER-03 force field [10], which fits the dihedral potentials to new quantum mechanical calculations using a low-dielectric continuum model, and Simmerling and coworkers have reported on their AMBER-99 variant, commonly referred to as AMBER-99SB, which includes reparameterized [var phi]/ψ backbone torsions that were fit to ab initio calculations of ALA and GLY tetrapeptides [2]. The backbone torsional potentials for each of these force fields, as well as the previously characterized AMBER-94 and AMBER-99[var phi] force fields, are shown in Figure 1, which highlights both the relative magnitude of the respective [var phi]/ψ potentials and the locations of energetic minima and maxima. Regions of the [var phi]/ψ space corresponding to alpha helix, beta strand, and polyproline backbone configurations are also indicated based on the definitions of Garcia for polyalanine peptides [11].

Figure 1
The [var phi]/ψ potentials for the four AMBER force fields compared in this work.

In tandem with non-bonded interactions, these torsional potentials comprise one of the most dominant factors in determining the configurational preferences of each molecular mechanical force field. This article therefore reports on our continued and broadened assessment of contemporary potentials in simulating helical peptides and proteins, which includes evaluation of the conformational equilibrium of large ensembles of small alanine-based peptides, and examines the ability of these force fields to stabilize the native structure of a large and flexible five-helix-bundle lipid-transport protein, apolipophorin-III (apoLp-III) [12], on the 100 ns timescale.


Helical peptide simulations

The capped A21 (Ace-A21-NMe) and Fs (Ace-A5[AAAR+A]3A-NMe) peptides were simulated using the AMBER-99[var phi] [7], AMBER-99SB [2], and AMBER-03 [10] all-atom potentials, which were ported into the GROMACS molecular dynamics (MD) suite [13] as part of our ffAMBER distribution of force field ports (available at, and quantitatively validated against the AMBER 8 software. The GROMACS suite, known to be one of the fastest MD packages available for biomolecular simulation, was modified for the emoH@gnidloF [14] infrastructure ( Following our previous methodology [1], [7], [15], a canonical helix ([var phi] = −57°, ψ = −47°) and a random coil configuration with no helical content were generated and centered in 40 Å cubic boxes. For Fs peptide simulations, electroneutrality was achieved by placing three Cl ions randomly around the solute with minimum ion-ion and ion-solute separations of 5 Å. As the degree of helicity of the Fs peptide has been shown to be nearly constant between pH 1 and pH 8 for NaCl concentrations up to 4 M [16], [17], ionic strength should not influence secondary structure in our simulations. The helix and coil conformations were then solvated with 2,075 and 2,065 TIP3P water molecules [18], respectively. After energy minimization using a steepest descent algorithm and solvent annealing for 500 ps of MD with the peptide held fixed, each of these conformations served as the starting point for 1,000 independent MD trajectories in each AMBER potential listed above, which were simulated on ~10,000 CPUs within the emoH@gnidloF supercluster.

All simulations reported herein were conducted under constant NPT conditions [19] at 1 atm and 305 K, the approximate Fs midpoint temperature detected by circular dichroism [20] and ultraviolet resonance Raman [21]. As Ewald approaches are known to overstabilize helices in periodic cells [22], long range electrostatic interactions were treated using the reaction field method [6], [7], [23] with a dielectric constant of 80, and 9 Å cutoffs imposed on all non-bonded interactions, which has proven highly successful in our previous studies of helical peptides [1], [7], [8], [15], [24]. Nonbonded pair lists were updated every 10 steps, and covalent bonds involving hydrogen atoms were constrained with the LINCS algorithm [25]. An integration step size of 2 fs was used, with coordinates stored every 100 ps. As shown in Table 1, the cumulative sampling achieved for the Fs peptide totals over 350 µs, with equilibrium sampling of over 150 µs–orders of magnitude longer than the ~16 ns folding time of Fs. Similar sampling of the A21 peptide was also collected.

Table 1
Simulated ensemble statistics for the Fs and A21 peptides.

Apolipoprotein simulations

The 164-residue apolipophorin-III protein (PDB 1LS4, model 1) [26] was simulated according to the methods outlined above with minimal changes. Electroneutrality was achieved by placing four Cl ions and twelve Na+ ions randomly around the solute with minimum ion-ion and ion-solute separations of 5 Å, and the temperature was set at 300 K. The simulation box size was set at 66 Å×67 Å×88 Å, and the system was solvated with 12,552 TIP3P water molecules. Simulations were then run in parallel on a local computing cluster using the AMBER-94, AMBER-99[var phi], AMBER-99SB, and AMBER-03 force fields with 12 Å cutoffs imposed on non-bonded interactions to accurately account for tertiary contacts within the helical bundle. Eight 100 ns simulations per force field were conducted, yielding a total apoLp-III sampling time of 3.2 µs. Each 100 ns simulation of this system, composed of 40,118 atoms in all, required approximately 6 months of wall-clock time on a single 2.5 GHz Xeon processing core.


For the helical peptides and proteins studied herein, structural content per residue was assessed using the Dictionary of Secondary Structure in Proteins (DSSP) [27], which has gained general acceptance among the biophysical community. As in our previous analyses employing DSSP [7], total helical content (H) includes α, 310, and π helical types. Beta structure (B) has been defined as consisting of both β-sheet and β-bridge conformations, and turn regions (T) are distinguished from random coil (C) configurations. In addition, we have assigned polyproline type II (P) structure within our data based on backbone torsions according to the region of the Ramachandran map defined by [var phi] = −75°±50° and ψ = 150°±50°. While seemingly broad, this window was determined via quantitative sampling of the native collagen triple-helix sequences (POG)10 and (PPG)10 in various AMBER force fields.

Other structural metrics used for comparison of these force fields include the all-atom root-mean-squared-deviation (RMSD) and gyration radius (Rg), which were calculated using GROMACS analysis tools, as well as the solvent accessible surface area (SASA). We note that analysis tools within molecular simulation packages that calculate SASA, such as the g_sas module in GROMACS, often base definitions of polar and non-polar surface areas on atomic partial charges taken from the molecular potential being used. To avoid this force field dependence due to significant differences in point charges between AMBER variants, and obtain SASA quantities that were based solely on sequence and structure, reported SASA values were calculated using the VEGA package [28], with a probe sphere radius of 1.4 Å and a point density of 24 points/Å2. Reported SASA values, including polar and non-polar components, are thus directly comparable between force fields and with the 21 energy-minimized NMR models of apoLp-III [26].

Results and Discussion

Fs and A21 helical peptides

Ensemble convergence of mean helical structure on the 40 ns timescale when using the AMBER-99[var phi], AMBER-99SB, and AMBER-03 force fields, as demonstrated in Figure 2, is in agreement with our previous observations [1], [7]. Additional structural metrics, all of which showed fully convergent behavior on this timescale, are reported in Table 2. Helical metrics include the mean total helicity (Nhelix), mean 310 helicity (N310), mean number of helical segments (Nseg), and mean longest contiguous helical segment (Ncont). As our previous efforts have demonstrated that the Lifson-Roig two-state helix-coil model [29] does not adequately capture the complex character of helix-coil equilibria, and LR counting is not an intuitive, nor necessarily accurate, method of assigning helical status to residues in a given conformation, the values shown in Figure 2 and Table 2 come from the generally-accepted DSSP.

Figure 2
Convergence of mean helical content for the (a) Fs and (b) A21 ensembles.
Table 2
Ensemble averaged equilibrium structural properties for the Fs and A21 peptides.

Most noticeable in both Figure 2 and Table 2 is the lack of any significant helical content for both the Fs and A21 peptides under the AMBER-99SB potential at equilibrium. Additionally, the ensemble average RMSD from the ideal helix is significantly larger than predicted by the AMBER-99[var phi] and AMBER-03 potentials. We note that for these small helical peptides the loss of helicity takes 10–20 ns, and simulations of larger proteins, particularly those that include tertiary stabilization of helical regions, should thus reach significantly longer timescales than this to adequately evaluate force field behavior. Moreover, our AMBER-99SB ensembles showed a slight tendency (~4% of residues in our equilibrium data set) to form beta structure, composed primarily of bridges of one to two β residues, which was not seen in the AMBER-03 or AMBER-99[var phi] data. Our simulations of polyalanine-based helical peptides thus suggest that this force field not only destabilizes helical structure, as does the original AMBER-99 force field [7], but also exhibits β−structure forming tendencies.

Unlike AMBER-99SB, the AMBER-03 force field of Duan et al. clearly exhibits significant helical content in a manner that is qualitatively similar to our AMBER-99[var phi] potential, and a quantitative comparison of these models must therefore focus on the fine details of the established conformational equilibria. Notably, while the AMBER-99[var phi] ensembles converge to nearly identical equilibria for the Fs and A21 peptides, this is not the case for the AMBER-03 potential, which predicts a significant increase in overall helicity of ~30% when moving from the ARG+ substituted Fs peptide to the polyalanine 21-mer. In tandem with this increase in overall helicity, the AMBER-03 potential predicts an increase in average contiguous helix length of ~20%.

It is well known that polyalanine peptides are insoluble in aqueous solution, which led the Baldwin laboratory to develop many model helical peptides by substituting a limited number of ALA residues with polar amino acids [30], such as the ARG+ substitutions made in Fs. It has been postulated that these sidechains serve not only to make these peptides soluble, but also to stabilize the helical structure within substituted polyalanine peptides [6], [31]. Herein lies an important distinction between these force fields: while AMBER-99[var phi] predicts essentially identical helical propensities between the insoluble polyalanine peptide and its soluble ARG+ substituted analog, the AMBER-03 potential predicts increased helicity for the insoluble peptide in comparison to its soluble ARG+ substituted analog, suggesting that the polar ARG+ sidechains actually serve to destabilize helical structure in the AMBER-03 potential. While it is difficult to compare this aspect of these potentials when applied to an insoluble peptide, which is a purely non-physical model for which no experimental data can be obtained, the resulting decrease in helicity that accompanies ALA→ARG+ substitutions when using AMBER-03 seems to diverge significantly from current knowledge of substituted polyalanine peptides.

We can better understand the differences between these force fields by considering the equilibrium sampling of Fs peptide backbone torsions and the resulting Ramachandran maps shown in Figure 3. For visual comparison to Figure 1, and discussion of our apoLp-III simulations in the section below, a map for the helix-stabilizing AMBER-94 force field is included (Figure 3a). In our previous work [7], we stated that of the many AMBER force fields studied at that time, “the best agreement with the Protein Data Bank and quantum mechanical sampling is achieved by the AMBER-99[var phi] variant, which captures distributions that are underweighted by other force fields without overweighting other regions of the phase space.” We thus compare both AMBER-99SB and AMBER-03 to our original data set (Figure 3b).

Figure 3
Free energy landscapes projected onto the Ramachandran map.

While exhibiting an energetic minimum in the appropriate helical region of the [var phi]/ψ map, the AMBER-99SB variant (Figure 3c) shows a strong bias toward extended β and polyproline conformations ([var phi]<0° and ψ>90°) that is qualitatively similar to that observed for the AMBER-03 potential. As demonstrated in Figure 1, both of these potentials inherently include local minima in this region, with AMBER-99SB including several local minima that span all high ψ values. In addition, AMBER-99SB heavily overweighs the prevalence of left-handed helices in the ([var phi], ψ) = (57°, 47°) region of the map. Together, these characteristics of AMBER-99SB act to overwhelmingly minimize the helical content of these model helix-forming peptides.

It is certainly of import and interest that our observations regarding the AMBER-99SB force field clearly contradict recent reports that apply this force field to polyalanine peptides and larger proteins, such as the simulational study of ubiquitin by Showalter and Brüschweiler [32], and suggests that AMBER-99SB outperforms its predecessor, AMBER-99, in accurately predicting N-H NMR S2 order parameters. Their limited sampling, however, on the 20 ns timescale for ubiquitin, which is composed of predominantly beta structure with only a small helical region, is not indicative of the ability of this force field to adequately stabilize helical structure on longer timescales. In addition, recent reports by Best and Hummer have shown that the AMBER-99SB force field performs poorly in reproducing experimental helical content for small peptides [33], [34]. Contradicting this report, a recent study by Simmerling and co-workers, developers of the AMBER-99SB force field, reported that this potential performs among the best of currently available models in reproducing experimentally observed J-coupling constants [35]. While this concise review of recent studies on force field behavior only scratches the surface of the heavy recent activity within the literature, it is clear that contradictions abound, and we can only report here our observations of the behavior of this system with respect to the well characterized helical Fs peptide.

In contrast to the observed behavior of AMBER-99SB, the [var phi]/ψ map representing the AMBER-03 force field shown in Figure 3d is qualitatively similar to that derived using AMBER-99[var phi]. However, some differences are readily apparent. Most importantly, the energetic minimum corresponding to left-handed helix conformations is slightly less favored in the AMBER-03 force field, while the polyproline region of the map is significantly more favorable than observed for the AMBER-99[var phi] potential. In addition, AMBER-03 overweighs the region of the map that corresponds to extended β conformations and, in general, samples a much broader portion of the Ramachandran map than any of the other force fields studied. The lack of proper helical stability for Fs, as well as the decreased occurrence of 310 helix and mean number of helical segments per peptide (consistently 22% less than observed in AMBER-99[var phi]), are thus attributed to these aspects of the AMBER-03 potential.

To be fair, we must acknowledge recent studies that have made arguments against the use of the AMBER-99[var phi] force field in similar fashion to criticisms of AMBER-99SB [33], [36]. The recent increase in force field assessment has included a number of approaches, a variety of molecular systems ranging from short polyalanine peptides (generally ALA3 to ALA5) to larger protein systems, and various inconsistencies in methodology, such as the equilibration method followed and the solvation model used. We stress that our analysis of the AMBER-99[var phi] force field, and the comparisons made to other force fields, follows a previously outlined goal [7] of comparing a model helical system, the Fs peptide, to known experimental observables on a quantitative basis. We have previously shown that this force field does well at reproducing experimental thermodynamic and kinetic quantities for Fs [7], but made no claims about the application of this potential to the highly varying chemical moieties upon which the many recent reports described above have been based.

Indeed, following our initial publication of the AMBER-99[var phi] force field, we also reported on the inability of this potential and other AMBER potentials to adequately characterize polyproline type II structure in the blocked ALA7 peptide [24], suggesting in that report that force field behavior should be dependent on peptide length, as highly diverging structural character has been observed in varying lengths of polyalanine peptides. While AMBER-99[var phi] was shown to best predict the experimental radius of gyration and alanine J-coupling constant for this small polyalanine peptide compared to the AMBER-94, AMBER-GS, AMBER-96, AMBER-99, and GROMOS 53A6 force fields, a lack of experimental evidence for significant helix formation in that peptide suggests that this success was ambiguous. We thus look forward to a more systematic and consistent approach being adopted by the biosimulation community in assessing specific force fields with respect to specific peptide and protein sequences and lengths, and move forward in an effort to assess these force fields with respect to a much larger helical system: the flexible apolipophorin-III helical bundle protein.

Apolipophorin-III helical bundle protein

The NMR structure of the apoLp-III lipid-transport protein used to initiate our simulations is shown in Figure 4, with turn regions colored in bright green and helices 1 thru 5 colored in shades ranging from blue to green, respectively. This protein participates in lipid transport via an apparent “unhinging” or “bundle opening” mechanism [12], which exposes hydrophobic residues within the core of the native fold for lipid binding. It is known to have a relatively low structural stability of ~2 kcal/mol when folded and a midpoint temperature of T1/2~325 K [37]. Assuming similar folding kinetics to its relative, apolipoprotein A-I, this protein should fold on the seconds timescale [38]. We therefore expect the protein to be structurally stable, but also flexible, at 300 K on the 100 ns timescale. While the data from these 100 ns runs cannot be used to predict thermodynamic stability on longer timescales, they are adequate for assessing the ability of these force fields to model such systems on timescales currently accessible in silico.

Figure 4
Ribbon views of the NMR model of apolipophorin-III.

Figure 5 shows the time-averaged structural properties of the protein when using the force fields studied herein. For comparative purposes, the seminal Cornell force field, AMBER-94 [3], was also employed in this study. The AMBER-94 force field, while being quantitatively the closest to our AMBER-99[var phi] potential in accurately reproducing experimental values for the Fs peptide in our previous study [7], was also shown to significantly overstabilize helical structure and thus serves as a good benchmark by which the helicity of more recently reported potentials can be considered. We note that the AMBER-94 force field maintains the most helical and rigid of structures seen in our 100 ns trials, as would be expected based on our previous results and similar observations reported in the literature.

Figure 5
Mean structural properties of simulated apolipophorin-III ensembles.

In contrast, the AMBER-03 and AMBER-99[var phi] potentials both maintain slightly larger RMSD values and slightly lower helicity on average, as was observed for the helical peptides studied herein. The AMBER-03 ensemble does show a slight continued increase in both RMSD and Rg over this window, and it is thus not possible to draw conclusions about the behavior of this force field on longer time scales. The AMBER-99[var phi] ensemble shows a more constant behavior across the 100 ns trial window with respect to these structural metrics, with one exception: late in this simulation window a single apoLp-III simulation has the protein undergoing an unhinging motion of helix 1 in which the protein bundle opens, resulting in larger RMSD and Rg values, with no significant change in helicity. Motion of this sort has been suggested as a possible means by which lipid moieties are bound for transport [12], and it is thus assumed that such motion is not an indication of loss of structure, but rather a native dynamics of this protein. Further study of this behavior, particularly on much longer time scales, is currently underway.

As is apparent in Figure 5, the lack of helical content observed for small helical peptides under the AMBER-99SB force field was again observed in our simulations of apoLp-III. Unfolding of the protein occurs on the tens of nanoseconds timescale, with the entire ensemble of eight simulations showing distinct unfolding behavior throughout the 100 ns simulation window and the mean number of helical residues decreasing ~40% alongside a mean RMSD reaching ~20 Å. This unfolding is also clearly demonstrated in Table 3, which details the mean solvent-accessible surface areas observed in our apoLp-III simulations and amongst the 21 NMR models. The values in Table 3 were generated only after discarding the first 10.0 ns of simulation time in each force field to allow for structural relaxation. Unlike the other AMBER force fields studied, the AMBER-99SB potential shows very large divergence toward greater polar and, much more notably, non-polar surface areas. These large values correspond to the unfolding described in Figure 5.

Table 3
Mean SASA* for apolipophorin-III simulations.

Corresponding structural probabilities per residue are shown in Figure 6 for each of these force fields. The helix schematic at the top of that figure depicts the NMR model used to initiate our simulations, including turn and coil regions. Figure 6 demonstrates a significant extension of helical regions in apoLp-III when employing AMBER-94, primarily into regions that are expected to act as turns between helices. In contrast, AMBER-99[var phi] yields the most flexible molecular shape and size with standard deviations in RMSD and Rg of 2.8 and 1.5 Å, respectively, while also preserving a proper molecular size. In comparison, the AMBER-03 potential yields the highest RMSD and Rg of the three helix-friendly force fields. AMBER-99SB shows a much lower helical content over this timescale than these helix-friendly force fields, along with substantial preference for turn and coil conformations. Interestingly, the AMBER-03 and AMBER-99SB force fields show similar polyproline type II conformational content, which is lacking from the AMBER-94 and AMBER-99[var phi] potentials.

Figure 6
Structural sampling of apolipophorin-III per residue.

We note that AMBER-99[var phi] and AMBER-03 both yield significantly lower mean helical content than AMBER-94, and slightly larger mean RMSD values. While the NMR models provide a picture of the average structure of apoLp-III, it is expected that helical regions will fluctuate on the nanosecond timescale, as observed for smaller helical peptides. Moreover, as this protein is expected to undergo some large scale fluctuations to function as a lipid transport protein [12], the observed flexibility of this protein is not unexpected as it is a primary aspect of what makes the apolipoprotein family interesting from both biochemical and computational vantages. Indeed, accurate simulations of these proteins should not maintain the overstabilized helical propensities exhibited in AMBER-94 simulations, nor a tendency to extend helical regions. Both AMBER-99[var phi] and AMBER-03 exhibit significantly less helical propensity in regions that are initially in turn or coil conformations and show global flexibility that should be expected for such proteins.

As noted above, Table 3 lists the mean total SASA, polar SASA, and non-polar SASA predicted by each force field. While it should be noted that the three helix-friendly force fields are all in agreement with the NMR ensemble within the errors reported in Table 3, the mean values observed for each force field do show significant differences on the timescale of our study. The NMR models yield a consistently lower total SASA than observed in any of our simulated ensembles. In agreement with the largest RMSD and Rg, AMBER-03 predicts a total SASA that is ~10% larger than NMR models, with most of this discrepancy resulting from the exposure of non-polar groups to solvent during simulation. With a total SASA that is closest to the NMR average, AMBER-94 sampling demonstrates an unexpected trend: a decrease in polar SASA below the NMR average and an increase in non-polar SASA well above the NMR average. We postulate that this strikingly counterintuitive behavior results from the extension of helices into turn and coil regions as described above, which would require the rearrangement of sidechains away from their NMR positions, and we are currently investigating this behavior in order to offer a more definitive description of this phenomenon. Of the three force fields, AMBER-99[var phi] yields a reasonably low total SASA and a polar SASA that is nearly identical to the NMR average. As is the case for both AMBER-94 and AMBER-03, the AMBER-99[var phi] potential also exposes significant non-polar SASA to solvent, suggesting a force field independent, systematic trend in all-atom simulations of this protein, which is also currently under investigation.

Finally, the AMBER-99[var phi] and AMBER-03 sampling shown in Figure 6 demonstrates distinct structural trends that should be noted. First and foremost, while these two models do not suffer from the overstabilization of helical structure observed in AMBER-94 simulations, each exhibits clear trends with regard to the turn and coil regions connecting helices 1 thru 5. In our AMBER-03 sampling this includes a significantly larger probability of being in polyproline type II and random coil conformations between helices, with greatly diminished likelihood of turn regions. In stark contrast, AMBER-99[var phi] favors more compact turn conformations over coil and polyproline configurations, thereby explaining the accurate Rg and more accurate mean SASA values observed in AMBER-99[var phi] sampling.

Of these three models which is “best” for simulating large helical proteins? As is the case in other contexts, proteins that are inherently flexible cannot be well characterized by a single quantitative metric, such as RMSD. We thus address this question by considering other quantitative metrics of structural integrity. While three of these four potentials prove to stabilize apoLp-III simulations, some significant disparities are observed. AMBER-94 overstabilizes helical regions, extending them into the turn and coil linker regions between helices, and strays from the expected trend in moving toward higher polar SASA and lower non-polar SASA. AMBER-03 yields a very reasonable total helical content, yet also overestimates the molecular size, significantly overestimates both components of the SASA, and favors coil and polyproline conformations between helices where turns should be prevalent. Of the three, AMBER-99[var phi] maintains proper molecular size, overall helical content, and polar SASA, while favoring turns in regions between helices rather than coil or extended configurations.


We have employed several molecular mechanical potentials and evaluated these models to assess their relative accuracy in simulations of helical peptides and proteins. The AMBER-99SB and AMBER-03 potentials were compared to the AMBER-99[var phi] helix-coil force field as applied to polyalanine-based helical peptides and several observations were made. Most notably, the AMBER-99SB potential is decidedly helix-destabilizing, as demonstrated by the rapid unfolding of both small helical peptides and the apoLp-III helix bundle on the tens of nanoseconds timescale. We also found that the AMBER-03 potential, while yielding reasonable and qualitatively similar results for these helical peptides when compared to our AMBER-99[var phi] data, also suffers from multiple shortcomings. These include a tendency to show a significant decrease in helical content with ALA→ARG+ substitutions, a lack of expected 310 helical content, and an overweighing of the extended β and polyproline portions of the Ramachandran map.

These force fields were also employed in simulations of the 164-residue five-helix-bundle apolipophorin-III protein, as was the seminal AMBER-94 force field of Cornell et al. This latter potential strongly overstabilized the helical content in apoLp-III, showing a tendency to expand helical regions into the turn and coil linkers separating helices 1 thru 5, while maintaining a very rigid gyration radius of 16.9±0.2 Å, very close to the NMR average of 17.2 Å. Surprisingly, the AMBER-94 data showed distinct shifting to lower polar SASA and higher non-polar SASA, a trend that warrants further investigation. AMBER-03 yielded the largest mean RMSD, Rg, and SASA values of the three helix-friendly potentials, demonstrating a definitive preference for extended and polyproline conformations, as was observed in our equilibrium sampling of small helical peptides.

Of the four force fields discussed above, AMBER-99[var phi] yields the best quantitative results in simulating helical peptides. This includes the occurrence of a low but appreciable 310 helical content, as well as consistent results for the insoluble polyalanine 21-mer and its ARG+ substituted analog. In addition, this potential yields nearly perfect agreement with the apoLp-III NMR average gyration radius and polar SASA, while also exhibiting flexibility and fluctuation in helical regions expected of a helix bundle protein with low stability.

In assessing our apoLp-III data, a trend toward increasing non-polar SASA was observed with all AMBER variants tested, the most significant of which (AMBER-03) showed a nearly 14% increase in mean non-polar solvent exposed surface area within 100 ns. There is no doubt that simple point-charge molecular mechanical models, such as the AMBER force fields employed in this study, suffer from imperfections that vary from model to model. We find it striking that these three AMBER variants, which differ significantly both in terms of their respective point-charges and their torsional parameterizations, would consistently favor an increasing non-polar surface area while showing such dramatically different trends in polar SASA. We are thus investigating this phenomenon and look forward to providing a follow-up report on the cause for such trends.


We thank Carlos Simmerling for his communication regarding our findings, members of the Sorin lab for their thoughtful discussion of this manuscript, and the worldwide emoH@gnidloF volunteers who contributed invaluable processor time to this effort (


Competing Interests: The authors have declared that no competing interests exist.

Funding: E.J.T. was supported by a Provost's Undergraduate Student Summer Research Stipend and A.J.D. was supported by a Women & Philanthropy Undergraduate Research and Creative Activity Scholarship. The computation was supported by a Research Corporation Cottrell College Science Award ( and support from the CSULB Academic Technologies Division. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


1. Sorin EJ, Pande VS. Empirical Force Field Assessment: the Interplay Between Backbone Torsions and Non-covalent Term Scaling. Journal of Computational Chemistry. 2005;26:682–690. [PubMed]
2. Hornak V, Abel R, Okur A, Strockbine B, Roitberg A, et al. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Structure, Function, and Bioinformatics. 2006;65:712–725. [PubMed]
3. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, et al. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal Of The American Chemical Society. 1995;117:5179–5197.
4. Kollman PA. Advances and Continuing Challenges in Achieving Realistic and Predictive Simulations of the Properties of Organic and Biological Molecules. Accounts of Chemical Research. 1996;29:461–469.
5. Wang J, Cieplak P, Kollman PA. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? Journal of Computational Chemistry. 2000;21:1049–1074.
6. Garcia AE, Sanbonmatsu KY. α-Helical stabilization by side chain shielding of backbone hydrogen bonds. Proceedings Of The National Academy Of Sciences Of The United States Of America. 2001;99:2782–2787. [PubMed]
7. Sorin EJ, Pande VS. Exploring the Helix-coil Transition via All-atom Equilibrium Ensemble Simulations. Biophysical Journal. 2005;88:2472–2493. [PubMed]
8. Sorin EJ, Rhee YM, Shirts MR, Pande VS. The solvation interface is a determining factor in peptide conformational preferences. Journal Of Molecular Biology. 2006;356:248–256. [PubMed]
9. Monticelli L, Sorin EJ, Tieleman DP, Pande VS, Colombo G. Molecular Simulation of Multistate Peptide Dynamics: A Comparison Between Microsecond Timescale Sampling and Multiple Shorter Trajectories. Journal of Computational Chemistry. 2008;29:1740–1752. [PubMed]
10. Duan Y, Wu C, Chowdhury S, Lee MC, Xiong G, et al. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. Journal of Computational Chemistry. 2003;24:1999–2012. [PubMed]
11. Garcia AE. Characterization of non-alpha helical conformations in Ala peptides. Polymer. 2004;45:669–676.
12. Leon LJ, Idangodage H, Wan C-PL, Weers PMM. Apolipophorin III: Lipopolysaccharide binding requires helix bundle opening. Biochemical and Biophysical Research Communications. 2006;348:1328–1333. [PMC free article] [PubMed]
13. Lindahl E, Hess B, van der Spoel D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. Journal Of Molecular Modeling. 2001;7:306–317.
14. Zagrovic B, Sorin EJ, Pande V. β-hairpin folding simulations in atomistic detail using an implicit solvent model. Journal Of Molecular Biology. 2001;313:151–169. [PubMed]
15. Sorin EJ, Pande VS. Nanotube confinement denatures protein helices. Journal Of The American Chemical Society. 2006;128:6316–6317. [PubMed]
16. Williams S, Causgrove TP, Gilmanshin R, Fang KS, Callender RH, et al. Fast Events in Protein Folding: Helix Melting and Formation in a Small Peptide. Biochemistry. 1996;35:691–697. [PubMed]
17. Lockhart DJ, Kim PS. Electrostatic screening of charge and dipole interactions with the helix backbone. Science. 1993;260:198–202. [PubMed]
18. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. Journal Of Chemical Physics. 1983;79:926–935.
19. Berendsen HJC, Postma JPM, Van Gunsteren WF, Dinola A, Haak J. Molecular-dynamics with coupling to an external bath. Journal Of Chemical Physics. 1984;81:3684–3690.
20. Thompson PA, Eaton WA, Hofrichter J. Laser Temperature Jump Study of the Helix-Coil Kinetics of an Alanine Peptide Interpreted with a ‘Kinetic Zipper’ Model. Biochemistry. 1997;36:9200–9210. [PubMed]
21. Ianoul A, Mikhonin A, Lednev IK, Asher SA. UV Resonance Raman Study of the Spatial Dependence of α-Helix Unfolding. Journal of Physical Chemistry A. 2002;106:3621–3624.
22. Weber W, Hunenberger PH, McCammon JA. Molecular dynamics simulations of a polyalanine octapeptide under Ewald boundary conditions: Influence of artificial periodicity on peptide conformation. Journal of Physical Chemistry B. 2000;104:3668–3675.
23. Tironi IG, Sperb R, Smith PE, van Gunsteren WF. A generalized reaction field method for molecular dynamics simulations. Journal of Chemical Physics. 1995;102:5451–5459.
24. Zagrovic B, Lipfert J, Sorin EJ, Millett IS, Van Gunsteren WF, et al. Unusual compactness of a polyproline type II structure. Proceedings Of The National Academy Of Sciences Of The United States Of America. 2005;102:11698–11703. [PubMed]
25. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM. LINCS: A Linear Constraint Solver for molecular simulations. Journal of Computational Chemistry. 1997;18:1463–1472.
26. Fan D, Zheng Y, Yang D, Wang J. NMR Solution Structure and Dynamics of an Exchangeable Apolipoprotein, Locusta migratoria Apolipophorin III. Journal Of Biological Chemistry. 2003;278:21212–21220. [PubMed]
27. Kabsch W, Sander C. Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features. Biopolymers. 1983;22:2577–2637. [PubMed]
28. Pedretti A, Villa L, Vistoli G. VEGA - An open platform to develop chemo-bio-informatics applications, using plug-in architecture and script programming. Journal of Computer-Aided Molecular Design. 2004;18:167–173. [PubMed]
29. Lifson S, Roig A. Theory of Helix-Coil Transition in Polypeptides. Journal Of Chemical Physics. 1961;34:1963–1974.
30. Marqusee S, Robbins VH, Baldwin RL. Unusually stable Helix Formation in Short Alanine Based Peptides. Proceedings Of The National Academy Of Sciences Of The United States Of America. 1989;86:5286–5290. [PubMed]
31. Shepherd CM, van der Spoel D, Vogel HJ. Molecular dynamics simulations of peptides from the central domain of smooth muscle caldesmon. Journal of Biomolecular Structure & Dynamics. 2004;21:555–565. [PubMed]
32. Showalter SA, Bruschweiler R. Validation of Molecular Dynamics Simulations of Biomolecules Using NMR Spin Relaxation as Benchmarks: Application to the AMBER99SB Force Field. Journal of Chemical Theory and Computation. 2007;3:961–975.
33. Best RB, Buchete N-V, Hummer G. Are Current Molecular Dynamics Force Fields too Helical? Biophysical Journal. 2008;95:L07–L09. [PubMed]
34. Best RB, Hummer G. Optimized Molecular Dynamics Force Fields Applied to the Helix-Coil Transition of Polypeptides. The Journal of Physical Chemistry B. 2009;113:9004–9015. [PMC free article] [PubMed]
35. Wickstrom L, Okur A, Simmerling C. Evaluating the Performance of the ff99SB Force Field Based on NMR Scalar Coupling Data. Biophysical Journal. 2009;97:853–856. [PubMed]
36. Viktor H, Robert A, Asim O, Bentley S, Adrian R, et al. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Structure, Function, and Bioinformatics. 2006;65:712–725. [PubMed]
37. Weers PMM, Abdullahi WE, Cabrera JM, Hsu T-C. Role of Buried Polar Residues in Helix Bundle Stability and Lipid Binding of Apolipophorin III: Destabilization by Threonine 31. Biochemistry. 2005;44:8810–8816. [PubMed]
38. Gross E, Peng D-Q, Hazen SL, Smith JD. A Novel Folding Intermediate State for Apolipoprotein A-I: Role of the Amino and Carboxy Termini. Biophysical Journal. 2006;90:1362–1370. [PubMed]

Articles from PLoS ONE are provided here courtesy of Public Library of Science