Article pubs.acs.org/jcim

Force-Field Induced Bias in the Structure of Aβ21−30: A Comparison of OPLS, AMBER, CHARMM, and GROMOS Force Fields Micholas Dean Smith,†,§ J. Srinivasa Rao,†,‡ Elizabeth Segelken,† and Luis Cruz*,† †

Department of Physics, Drexel University, 3141 Chestnut Street, Philadelphia, Pennsylvania 19104, United States Department of Physics, New Jersey Institute of Technology, University Heights, Newark, New Jersey 07102-1982, United States

‡

S Supporting Information *

ABSTRACT: In this work we examine the dynamics of an intrinsically disordered protein fragment of the amyloid β, the Aβ21−30, under seven commonly used molecular dynamics force fields (OPLS-AA, CHARMM27-CMAP, AMBER99, AMBER99SB, AMBER99SB-ILDN, AMBER03, and GROMOS53A6), and three water models (TIP3P, TIP4P, and SPC/E). We find that the tested force fields and water models have little effect on the measures of radii of gyration and solvent accessible surface area (SASA); however, secondary structure measures and intrapeptide hydrogen-bonding are significantly modified, with AMBER (99, 99SB, 99SB-ILDN, and 03) and CHARMM22/27 force-fields readily increasing helical content and the variety of intrapeptide hydrogen bonds. On the basis of a comparison between the population of helical and β structures found in experiments, our data suggest that force fields that suppress the formation of helical structure might be a better choice to model the Aβ21−30 peptide.

■

INTRODUCTION All-atom molecular dynamics (MD) is a powerful computational technique for studying the dynamics of proteins in the submicrosecond regime; a regime important, among other things, for early stage peptide/protein folding events. This technique, however, is only as accurate as the force fields and accompanying water models (energy/solvent models) used to describe the bonded, dihedral, torsional, nonbonded, and other energy terms of particle interaction. To obtain the parameters necessary to match experimental results, force fields are typically developed by fitting parameters describing these molecular interactions to the results of highly detailed quantum-mechanics (QM) calculations of small peptide fragments and/or common empirical measures (e.g., heat of vaporization, heat capacities, liquid densities) of experimentally well-characterized short peptides (or both organic liquids and peptides, as in the case of OPLS).1−12 Of the large variety of force fields, four of the most commonly used in simulations of biomolecules are the GROMOS,1 OPLS,3,4 AMBER,2,7−9 and CHARMM5,6,10,11 families which were developed for the simulation of biological systems (and/or organic liquids) and typically differ in the systems and protocols used at the parametrization phase of development. Over time, the AMBER, CHARMM, GROMOS, and OPLS force fields were reparametrized (indeed reparameterization is ongoing, in the last year a new AMBER parametrization has been released13), with AMBER99, CHARMM22/27-CMAP, and OPLS-AA being improvements to their respective force fields. Following AMBER99’s initial © 2015 American Chemical Society

reparametrization, three additional forms were made bearing the AMBER moniker: AMBER99SB,8 AMBER99SB-ILDN,9 and AMBER03;7 where the first two were corrections to the AMBER99 force-field and AMBER03 was developed as its own separate entity. External validation of these force fields come from comparisons of simulation results to well-described experimental systems and other computational techniques. Example comparison systems include: dipeptides, 14 polyalanine chains,15−17 the Trp-cage,18 and Chignolin,18 though others are also used.18−20 The results of these external evaluations, overall, have demonstrated the strength of the AMBER, OPLS, GROMOS, and CHARMM families; however, minor biases have been found. For instance, helical content in some peptides has been overestimated in some AMBER18,19 and CHARMM force fields10,18,19 while both GROMOS and (to a lesser extent) OPLS may overestimate the populations of β secondary structures.15,18 Interestingly, the appearance of these biases may be peptide dependent,18 which suggests that the choice of systems (short peptides and organic liquids) used during the initial parametrization may introduce artifacts in simulations of unusual systems. Here, we address the question of whether biases in these force fields are a dominant factor in the simulated behavior of an intrinsically disordered peptide (IDP), the peptide fragment Aβ21−30. Briefly, Aβ21−30 is a protease resistant core of the fullReceived: May 24, 2015 Published: December 2, 2015 2587

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling

using 1000 steepest decent steps followed by a short 20 ps position restraint simulation used to relax the system. Once these preliminary steps were complete, production runs using a 2 fs time step and a frame recording rate of 1 ns were performed. In both production and position restraint trajectories, bonds were fixed with the use of SETTLE49 and LINCS.50 The total amount of production run simulation time (the sum of all trajectories performed) was 50 μs. To determine the convergence of our simulations, measurements of the mean fraction/percentage of ϕ/Ψ space mapped over time was also computed, as described in Scott et al.51 and is provided in the Supporting Information (SI), (Figures S2 provides the measure of convergence, while Figure S3 provides the Ramachandran plots constructed from the 2.5 μs of simulated time generated for each force field/water pair). Analysis of the trajectories included measurements of solvent accessible surface area (SASA), radius of gyration (Rg), intrapeptide and water−peptide hydrogen bonds (HBs, with a distance cutoff of 3 Å and a donor−acceptor angle cutoff of 20°), and secondary structure. Measurements of intrapeptide HBs and secondary structure were performed with the VMD software package52 and the STRIDE plug-in.53,54 All other measurements were made with internal GROMACS utilities. Although secondary structure measurements were performed with the STRIDE, a second categorization system was used to further reduce the secondary structure measurements at each frame for the entire chain into five classes (β, Bridge, Helix, Turn, and Coil).34 In this secondary categorization system, β structures are any STRIDE sequence with at least two residues in the “E” (extended β) state separated by at least one residue (e.g., CEETTEECCC). Bridges are taken as any sequence not already identified as a β-structure and contains two or more “B” (β-bridge) state residues. Helix labeled sequences are sequences not already defined as a β or Bridge and has three adjacent “H” (α-helix), “G” (3−10 helix), or “I” (π-helix) classified residues (e.g., TCCHHHCCTC). Turns are taken as sequences not already defined as a β, Bridge, or Helix and have residues identified as “T” in positions 1 and 10, 2−5, 6−9, or 4−6 (eg.: CCCTTTCCCC or TTTCCCCCTTT). Coils are taken to be any sequence not identified as β, Bridge, Helix, or Turn. Sample representations of the secondary structure classifications are provided in the SI as Figure S1. In addition to these categorized secondary structure measures, the traditional measurements of per-residue secondary structure were also performed and are provided in the SI. To determine the nature of differences found between the OPLS-AA and AMBER99SB-ILDN force fields (see Results), we swapped a particular set of parameters corresponding to the nonbonded terms (Lennard-Jones) in the OPLS-AA by those found in the AMBER99SB-ILDN (i.e we have generated two modified force fields, one being constructed from OPLS-AA’s bonded terms with AMBER99SB-ILDN’s nonbonded terms and the other having AMBER99SB-ILDN’s bonded terms and OPLS-AA’s nonbonded terms). As OPLS-AA and AMBER99SB-ILDN use different rules for combining LennardJones terms, it should be noted that the adapted Lennard-Jones parameters are approximations. The specific parameters for this modified OPLS-AA force field along with the new values, are listed in Table S1. Five trajectories for the modified OPLS-AA force-field were simulated each with a total time of 500 ns, while ten trajectories of the modified AMBER99SB-ILDN force-field were generated and simulated for 500 ns. In both

length Alzheimer’s Amyloid-β protein and has been the subject of both experimental studies21−25 and previous molecular dynamics studies21,23,25−34 that have shown it to exist in an ensemble of coil and turn structures. Furthermore, these studies have consistently shown signals, from solution 2D-NMR (including ROESY, and TOCSY), limited protolysis, and amide proton exchange observations (to quantify hydrogenbonding), the existence of salt-bridges between the decapeptide’s second and/or third (Glu22/Asp23) residues with its eighth (Lys28) residues (noted by cross-peaks in the NMR spectra between Glu22/Asp23-Lys28 along with a lack of with a lack of hydrogen-bonding of Lys28 found from proton amine exchange measurements, strongly suggesting salt-bridge formation), important backbone−backbone hydrogen bonds between Asp23 and residues within the turn, particularly Ser26 (determined by structural reconstruction generated from 2D-NMR data), and a metastable β-hairpin motif, along with being predominately in a turn or coil state.21−25 It should also be noted that even when this decapeptide region is not in a free-peptide form, but rather a component of the full-length Amyloid-β protein, the Aβ21−30 region is known, from circular dicrhoism (CD) spectra, infrared (IR) spectra, and 2D-NMR spectra (through the turn rich region is shifted slightly to be between residue 20−26) to have a secondary structure that lacks helices (though the nearby N-terminal residues from 10− 20 may exhibit helices) and is made up of random coil, turn, or β structures.35−39

■

METHODS To examine the effects of various force fields on the dynamics of the Aβ21−30, decapeptide (sequence AEDGVSNKGA) the following force fields were used: OPLS-AA,3,4 AMBER99,2 AMBER99SB,8 AMBER99SB-ILDN,9 AMBER03,7 CHARMM22/27-CMAP 5,6,10,11 (referred to simply as CHARMM27 in all figures), and GROMOS53A6.1 To test the preceding list of force fields for biases, the GROMACS 4.5.6 software package40−43 was used to generate a total of 100 independent, 500 ns long, trajectories of the uncapped decapeptide’s dynamics from an initial stretched coil, i.e. a conformation with no initial hydrogen bonds and no secondary-structure features. Although the different force fields were originally parametrized for use with specific water models (e.g CHARMM with TIP3P), it is also interesting to examine possible solvent model-dependent biases. For this reason, for each force-field 15 trajectories were generated with (excluding GROMOS53A6), of which 5 trajectories used the TIP3P water model, 5 used the TIP4P water-model,44 and the remaining 5 used the SPC/E45 water model. As the implementation of GROMOS53A6 in GROMACS 4.5.6 lacked the parameters for TIP3P and TIP4P, only five simulations under this force-field, with the SPC/E model, were used. Further, simulations using in-house modified OPLS-AA and AMBER99SB-ILDN force fields (described below) were also performed. In all production trajectories the pressure was taken as 1 atm, at a temperature of 300 K, and inside a periodic box of ∼120 nm3 (chosen such that at least a 0.8 nm thick layer of water was placed between the peptide and the periodic walls), with the temperature and pressure kept fixed by using the Parrinello− Rahman46,47 and V-rescale48 baro/thermostats, respectively. Each system was prepared with the decapeptide centered in the box and solvated, followed by the replacement of one water molecule with a Na+ cation to render the system a zero net charge. After this, the energy of each system was minimized 2588

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling

Figure 1. Average (a) radius of gyration (Rg) and (b) solvent accessible surface area (SASA). Error bars are standard error of the mean.

Figure 2. Percentage of each trajectory exhibiting a given secondary structure. Error bars are standard error of the mean.

found between the AMBER99 and OPLS-AA and the remaining force fields. It is interesting to note, however, that even the largest differences are only on the order of ∼1 nm2 SASA and ∼0.1 nm for the Rg measures. Secondary Structure. Figure 2 presents the population percentages of secondary structures for each force-field and water model combination using our STRIDE-based classification scheme (see Methods or our previous work34). Secondary structure propensity per amino acid measured directly from STRIDE are provided in the SI, Figure S4. The data presented in these figures (Figures 2 and S4) show that the more recent parametrizations of AMBER (AMBER99SB-ILDN), CHARMM27-CMAP, and GROMOS53A6 find β structures (regardless of water model) in roughly ∼1% of total the simulation time, while the OPLS-AA and AMBER99SB force fields find β structures between 5% and 23% of the time, depending on the water model. An interesting aside is to note that the maximum number of β structures under the OPLS-AA and AMBER99SB force fields were found with the TIP4P water

cases, after the simulations completed, the secondary structure was determined using the scheme described above.

■

RESULTS Measurements of the average SASA and Rg values of the decapeptide under each force field/water model combination are presented along with secondary structure population fractions, average numbers of water−peptide HBs per frame, average number of intrapeptide HBs per frame, and average intrapeptide HB contact maps. These measures are used to compare and contrast force-field effects on secondary structure distributions, HB behavior, and additional conformational changes to the decapeptide. Solvent Accessible Surface Area (SASA) and Radius of Gyration (Rg). Figure 1 shows the mean Rg and SASA values for each of the force field/water combinations. It is evident from this data that in both measures there are only minor differences between each of the force field/water models with the largest difference between SASA and Rg measures being 2589

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling model and that unlike AMBER99SB, OPLS-AA still has a significant number of β structures with the use of the older TIP3P. For bridges (Figure 2; upper-right) AMBER99, AMBER99SB, AMBER99SB-ILDN, AMBER03, and CHARMM27-CMAP exhibit population percentages in the neighborhood of 5% of the simulation time, with AMBER99 exhibiting the least (less than 1% of simulation time). GROMOS53A6 and OPLS-AA, both exhibit between ∼3% and 9%. OPLS-AA, interestingly, displayed a minor water model dependence for this measure, with decreased bridges under TIP3P and SPC/E. Further comparison between the water models and force fields reveals that unlike OPLS, both the AMBER family and CHARMM27-CMAP force fields have negligible water model dependence for this structure. For the case of helices, OPLS-AA and GROMOS53A6 force fields have negligible numbers, while the remaining force fields do exhibit helical structures. AMBER03 and CHARMM27CMAP exhibit these structures the most, followed by AMBER99SB-ILDN, AMBER99, and AMBER99SB. Unlike the relative insensitivity of bridge formation to the choice of water model, for the AMBER and CHARMM force fields, substantial differences with water model choice are evident in helix formation. Indeed, the largest variation over helical appearance across the water models, was found with the CHARMM27-CMAP force-field. Closer examination of the water dependent variations, however, show that aside from AMBER99SB and AMBER03, error bars between force-field and water model combinations typically overlap, indicating that the water dependence, although apparent in this measure (helix formation), is overall weak. Unlike the β, Bridge, and Helix secondary structure categories, Turn structures were exhibited by all force-field and water model combinations. Besides OPLS-AA/TIP4P (and to a lesser extent AMBER99SB), dependence on water model choice was minor for this structural metric. Examination of this difference and comparison to the other secondary structure categories revealed that a driving factor for this difference was that the production of Helix, β, and/or Bridge structures diminished the trajectory percentages of the Turn structures. This is also apparent from the data shown in Figure S3. Hydrogen-Bonding. Figure 3 shows the average number of all intrapeptide HBs per frame. It is clear from this figure AMBER99 and OPLS-AA simulations had substantially higher numbers of intrapeptide HBs. Further, comparison between the water models for each force-field shows a reduction in intrapeptide HB (compared to TIP3P and SPC/E) with TIP4P water. By examining the number of HBs between the water and the decapeptide (Figure 4), it is clear that the reduction of intrapeptide HBs due to the TIP4P water is associated with an increased number of water−peptide HBs, suggesting that the TIP4P water molecules may be preventing formation of intrapeptide contacts by forming HBs with the decapeptide. This screening-like behavior is also found, to a lesser (but still significant) extent, with the SPC/E water-model in combination with the AMBER99, AMBER99SB, OPLS-AA, and CHARMM27-CMAP force fields. To further elucidate differences in the structure of the decapeptide, two different measures of HB contacts were implemented: (i) average HBs per frame of the simulations (frequency of HB formation) and (ii) the fractions of the total intrapeptide HBs contact maps. These are displayed in Figures 5 and 6, respectively. Focusing first on Figure 5, the HB

Figure 3. Average number of intrapeptide hydrogen bonds per frame. Error bars are standard error of the mean.

Figure 4. Average number of water−peptide hydrogen bonds per frame. Error bars are standard error of the mean.

frequency contact map, it is clear that, as with the results from Figure 3, OPLS-AA and AMBER99 have substantially more intrapeptide HBs. Interestingly, the largest differences in the frequency contact maps of OPLS-AA and AMBER99 are (i) the absence of A21-E22 contacts in OPLS-AA simulations, (ii) the absence of E22-D23 contacts in AMBER99, (iii) increased frequency of contacts between E22-(G25/S26/N27) with AMBER99, and for simulations using SPC/E and TIP3P, and (iv) AMBER99 having substantially more D23-G29 contacts than OPLS-AA. Turning to the contact maps of the fraction of total HBs per amino-acid pairs shown in Figure 6, the results show that in all of the force-field combinations tested, the D23-S26 HB was ubiquitous, with OPLS-AA having a predominant fraction of its total contacts associated with this HB and CHARMM27CMAP and GROMOS53A6 having the least amount of this contact. Further, occurrence of the A21-E22 HB was common for all but the OPLS-AA and CHARMM27-CMAP force fields. 2590

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling

Figure 5. Contact maps of the average HBs per frame (formation frequency) independent of other HBs.

Figure 6. Contact maps of the fraction of the total number of HBs per amino-acid pairs. The sum of all contacts per simulation add to 1.

Of particular interest from this figure, however, is that in most of the AMBER force fields, CHARMM27-CMAP, and GROMOS53A6 force fields there was a wide range of intrapeptide HBs including many not previously established by experimental methods, while OPLS-AA primarily sampled established HBs (as noted by both Lazo et al.23 and Fawzi et al.25). Further, the E22 residue is found to be “flexible” in its HB “partner” in all force fields, except OPLS-AA. On Helicity in AMBER99 Derivative Force Fields. To test whether the increased helicity found in the AMBER force fields was due to the value of their nonbonded parameters, the secondary structure of both a modified OPLS-AA and a modified AMBER99SB-ILDN force fields were determined, as described in the Methods. The results are listed in Tables S2 and S3. These data indicate a dramatic increase in the occurrence of helical structures in the modified OPLS-AA

force fields using the TIP4P water model. Where before there was only ∼1% occurrence, now there is a 47% occurrence in helical content. On the opposite end, by replacing the AMBER99SB-ILDN nonbonded terms with those of OPLSAA, the secondary structure is only slightly modified, with a minor decrease in the helical content (from ∼17% to 13%). It is interesting that although the helical content is not entirely lost in the modified AMBER99SB-ILDN force field, it is the only structure that is significantly changed by swapping nonbonded terms with OPLS; this would suggest that although the bonded parameters in AMBER99SB-ILDN have a preference for helical structures, the nonbonded contribute as well, though to a lesser extent than the bonded terms. Recalling (see Introduction) that revisions to AMBER to correct for an overestimation of helicity were made to bonded terms (indeed a recent evaluation of the newest AMBER (AMBER14SB) mentions that bonded terms 2591

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling

the previous experimentally suggested metastable β-hairpin motifs. Of the water models, the overall largest effect is on the frequency of intrapeptide and water−peptide HBs. This is in contrast to small modifications in secondary structure measures (in particular the weak effect on finding β/bridge structures and moderate influence on finding helical structures), though the direction of the influence is force-field dependent. Comparing SASA and Rg measures under different water models further reveals that the water models do not pose a discernible influence. From these observations it is clear, that outside of significant changes in intrapeptide HBs in particular force fields (OPLS-AA and AMBER99) and moderate modification to helical secondary structure, the results of the simulations are largely invariant to the choice of water. When considering our results, it is interesting to compare to the work of Matthes et al.18 and others56−59,61,62 on other peptides in order to provide a more clear picture of the dependency of secondary structure sampling on force-field choice. Recently, Matthes et al.18 demonstrated that in the de novo protein Chignolin, the folding dynamics are readily accelerated by AMBER force fields, while the rate of the folding dynamics under the OPLS-AA force-field are closer to expectations.18 What is particularly interesting about this observation is that similar to Aβ21−30, Chignolin exhibit turn and β-hairpin structures (with the latter defined as its folded conformation as Chignolin is not intrinsically disordered), yet AMBER force fields readily sample these states and have limited, if any, sampling of helical structures. By comparing the hydropathy (using the Kyte−Doolittle59 and Hopp−Woods60 scales, where values near 1 are hydrophilic for Hopp−Woods, and more negative are more hydrophobic in the Kyte-doolittle) of the primary structures of Aβ21−30 (AEDGVSNKGA; Hopp− Woods 0.7; Kyte−Doolittle −1.2) and Chignolin (GYDPETGTWG; Hopp−Woods −0.05; Kyte−Doolittle −1.6) along with previous studies,20−33,59,60 it is clear that while Aβ21−30’s sampling of turn (and β-hairpin motifs) is driven by intrapeptide salt-bridging and hydrogen-bonding, Chignolin’s collapse to turn (and eventually β-hairpin) structures is driven by hydrophobic collapse and cooperative hydrogen bonding. This suggests that an enhanced hydrophobic effect may be at play in AMBER force fields. To address this possibility, we performed a preliminary set of simulations (see the SI) comparing the residue−residue interaction energies of two cystine dimers in vacu and in water (TIP3P) with a fixed position restraint of 0.35 nm between the centers of mass of the two dimers. Table S4 reports our preliminary results and demonstrates that short-range interaction energies of the AMBER family of force fields are significantly enhanced by the inclusion of water, while GROMOS53A, OPLS-AA, and CHARMM are only slightly changed. This preliminary data strengthens our view that an enhanced hydrophobic effect may be at play; however, a more detailed study (that is beyond the scope of this work) is necessary to fully explore this possible hydrophobic enhancement. It should also be noted that this possible enhanced hydrophobic effect may not exhibit itself as an across-the-board increase in helicity. Indeed, this proposed effect may result in increased helicity for Aβ21−30 while at the same time be responsible for the accelerated folding of Chignolin, a peptide with minimum traces of helices. In addition to providing insights into a possible general force field bias in MD simulations, we are also able to address the long-standing question within the Aβ21−30 community of which

are used to correct for possibly deficient nonbonded interactions55), it is interesting that the intrinsic helical component remains. We would suggest that in future parametrizations of the AMBER force fields, it may be beneficial to pause the sole improvement of the bonded terms and instead re-evaluate the nonbonded terms in tandem with the bonded terms.

■

DISCUSSION This study compared molecular dynamics simulations of the Aβ21−30 decapeptide under seven common force fields (AMBER99, AMBER99SB, AMBER99SB-ILDN, AMBER03, CHARMM22/27-CMAP, GROMOS53A6, and OPLS-AA) and three different water models (TIP3P, TIP4P, and SPC/E). The comparison between the force fields and water models was achieved through the use of SASA, Rg, secondary structure distributions, intrapeptide and water−peptide HBs. Focusing first on the influence of the force field, it is interesting to note that measures associated with the “size” and solvent exposure of the decapeptide (SASA and Rg) revealed little differences between the force fields, suggesting that the range of configurations sampled by each force field is, roughly, the same. This is further supported by our tests of simulation convergence (SI Figure S2) and the associated Ramachandran plots (SI Figure S3), where, excluding the case of GROMOS53A6, the amount of unique ϕ/Ψ space sampled over 500 ns is roughly equal, and that the location of minima (along with the general sampling itself) of the ϕ/Ψ space is remarkably similar, with the exception of GROMOS53A6, between the force fields and water models. We should note that although the same minima regions are reached, the force fields do vary in the amount of sampling at each minima (i.e., the depth of the minima), which corresponds to differences in both in secondary structures (see below) and the sampling of the nonminima ϕ/Ψ space. Regarding the outlier, GROMOS53A6, it is instructive to note that it is a united-atom force field, where nonpolar hydrogens are “lumped” into heavy atoms, and as a result, this likely reduces some steric interference, allowing a larger fraction of the ϕ/Ψ space to be sampled compared to the other force fields, which may (and in this case does) increase the width of each minimum sampled at the cost of the depth of each minimum (though as we do not see any new/additional minima not already found in the other force fields, it is likely that the force-field only introduces many more short-lived conformations). The key differences between the force fields reveal themselves in the measures of secondary structure and intrapeptide HBs where rates of sampling of different portions of the conformational space are modified. In particular, some AMBER-Family and CHARMM27-CMAP force fields, as previously reported,10,18,19 enhance the occurrence of helical structures. Also, the overall number of intrapeptide HB are substantially greater in the OPLS-AA and AMBER99 than in the other force fields. An interesting secondary note is that the variety of intrapeptide HB types occurring in all but OPLS-AA and AMBER99 was large, with >45% of all possible HB combinations sampled above a threshold of 0.5% (i.e., the percentage of a given HB map in Figure 6 which had a measured HB combination greater than 0.5%). This wide range of HBs in combination with the secondary structure measures indicates that the AMBER-family, (excluding AMBER99), CHARMM27-CMAP, and GROMOS53A6 (though to a lesser extent) are hindered in finding the local minima associated with 2592

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling force field and water model should be used to produce accurate simulations of Aβ21−30. Previous experimental work21−27 has given qualitative insight regarding secondary structure and intrapeptide HBs (see the Introduction); further work by Fawzi et al.25 has previously attempted to address this same question by combining NMR derived J-coupling constants with the same measures from MD and argued that AMBER99SB is a reasonable force-field choice. We alternatively use the metrics of intrapeptide HBs and secondary structures, which are qualitatively comparable to the structural details reported in structural characterization experiments and find that the CHARMM27-CMAP and AMBER-family force fields are not the best choices as they have substantial sampling of helical structures, which as noted in the Introduction is not observed experimentally for this region.21−25,35−39 Furthermore, by noting that OPLS-AA and GROMOS53A6 lack helical structures and do exhibit the characteristic E22/D23-K28 and D23-S26 salt-bridge/HBs, it can be suggested that either is a reasonably better force-field choice; though, OPLS may be slightly preferred as its sampling of ramachandra space is more attune to steric restrictions (as seen by the differences in Figure S2). Of the choice of water model, our results suggest that TIP3P is a reasonable choice as it has substantial sampling of the expected Turn structure, lacks Helical structures, and does a significant, though still limited, sampling of metastable β structures. To conclude, this work has compared MD simulations with OPLS-AA, AMBER99, AMBER99SB, AMBER99SB-ILDN, AMBER03, CHARMM27-CMAP, and GROMOS53A6 force fields with TIP3P, TIP4P, and SPC/E water models using the Aβ21−30 decapeptide as a probe and found that force-field biases are prevalent in the distribution of intrapeptide HBs and secondary structures and that the extent of the differences in secondary structure measures and frequency of intrapeptide HBs is also weakly dependent on water model. By comparing our results to similar work by Matthes et al.,18 we suggest that these biases may be associated with the parametrization of hydrophobic amino acids in the Amber family of force fields. Despite there being no explicit hydrophobicity parameter in these force fields, it would appear that the full-range of behaviors for these amino acids is not entirely taken into account resulting in a possible enhancement of the hydrophobic effect; though further investigation is necessary to quantify the sequence dependence of this effect and to rule out other possible interactions. Further, comparison of the results of our simulations with qualitative measures from experiments on the Aβ21−30 decapeptide also allow us, for this particular peptide, to recommend the use of OPLS-AA or GROMOS53A6.

■

Present Address §

M.D.S.: Center for Molecular Biophysics, University of Tennessee/Oak Ridge National Laboratory, Oak Ridge, TN, 37830-6309 USA.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Computational time was partially provided by XSEDE computational resources grant TG-MCB110142.

■

(1) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656−1676. (2) Wang, J. M.; Cieplak, P.; Kollman, P. A. How Well Does A Restrained Electrostatic Potential (RESP) Model Perform In Calculating Conformational Energies Of Organic And Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (3) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (4) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105, 6474−6487. (5) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential For Molecular Modeling And Dynamics Studies Of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (6) Mackerell, A. D.; Feig, M.; Brooks, C. L. Extending The Treatment Of Backbone Energetics In Protein Force Fields: Limitations Of Gas-Phase Quantum Mechanics In Reproducing Protein Conformational Distributions In Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25, 1400−1415. (7) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J. M.; Kollman, P. A Point-Charge Force Field For Molecular Mechanics Simulations Of Proteins Based On Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999− 2012. (8) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison Of Multiple Amber Force Fields And Development Of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (9) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials For The Amber ff99SB Protein Force Field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (10) Bjelkmar, P.; Larsson, P.; Cuendet, M. A.; Hess, B.; Lindahl, E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models. J. Chem. Theory Comput. 2010, 6, 459−466. (11) Mackerell, A. D.; Wiorkiewiczkuczera, J.; Karplus, M. An AllAtom Empirical Energy Function for the Simulation of Nucleic Acids. J. Am. Chem. Soc. 1995, 117, 11946−11975. (12) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. K. A New Force Field For Molecular Mechanical Simulation Of Nucleic Acids And Proteins. J. Am. Chem. Soc. 1984, 106, 765−784.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.5b00308. Details regarding modified force-field simulations, simulation convergence measures, secondary structure measures at a residue-level resolution, and cystine dimer water/vacu energy comparisons (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*Tel.: 215-895-2739. E-mail: [email protected]. 2593

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling

Electrostatic And Hydrophobic Interactions. J. Mol. Biol. 2008, 379, 815−829. (34) Smith, M. D.; Cruz, L. Effect of Ionic Aqueous Environments on the Structure and Dynamics of the Aβ21−30 Fragment: A MolecularDynamics Study. J. Phys. Chem. B 2013, 117, 6614−6624. (35) Soto, C.; Castano, E. M.; Kumar, R. A.; Beavis, R. C.; Frangione, B. Fibrillogenesis of Synthetic Amyloid-β Peptides Is Dependent On Their Initial Secondary Structure. Neurosci. Lett. 1995, 200, 105−108. (36) Hilbich, C.; Kisters-Woike, B.; Reed, J.; Masters, C. L.; Beyreuther, K. Aggregation and Secondary Structure Of Synthetic Amyloid βA4 Peptides Of Alzheimer’s Disease. J. Mol. Biol. 1991, 218, 149−163. (37) Soto, C.; Castano, E. M.; Frangione, B.; Inestrosa, N. C. The αHelical to β-Strand Transition in the Amino-terminal Fragment of the Amyloid β-Peptide Modulates Amyloid Formation. J. Biol. Chem. 1995, 270, 3063−3067. (38) Lim, K. H.; Collver, H. H.; Le, Y. T. H; Nagchowdhuri, P.; Kenney, J. M. Characterizations of Distinct Amyloidogenic Conformations of the Aβ(1−40) and (1−42) Peptides. Biochem. Biophys. Res. Commun. 2007, 353, 443−449. (39) Hou, L. M.; Shao, H. Y.; Zhang, Y. B.; Li, H.; Menon, N. K.; Neuhaus, E. B.; Brewer, J. M.; Byeon, I. J. L; Ray, D. G.; Vitek, M. P.; Iwashita, T.; Makula, R. A.; Przybyla, A. B.; Zagorski, M. G. Solution NMR Studies of the Aβ1−40 and Aβ(1−42) Peptides Establish That The Met35 Oxidation State Affects The Mechanim Of Amyloid Formation. J. Am. Chem. Soc. 2004, 126, 1992−2005. (40) Berendsen, H. J. C.; Vanderspoel, D.; Vandrunen, R. GromacsA Message-Passing Parallel Molecular-Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (41) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (42) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A High-Throughput And Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (43) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C GROMACS: Fast, Flexible, And Free. J. Comput. Chem. 2005, 26, 1701−1718. (44) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (45) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular Forces, The Jerusalem Symposia on Quantum Chemistry and Biochemistry 1981, 14, 331−342. (46) Parrinello, M.; Rahman, A. Polymorphic Transitions In Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (47) Nose, S.; Klein, M. L. Constant Pressure Molecular Dynamics For Molecular Dynamics. Mol. Phys. 1983, 50, 1055−1076. (48) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (49) Miyamoto, S.; Kollman, P. A. SETTLE: An Analytical Version of the SHAKE and RATTLEAlgorithms For Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (50) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. A Linear Constraint Solver For Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (51) Scott, K. A.; Alonso, D. O. V.; Sato, S.; Fersht, A. R.; Daggett, V. Conformational Entropy of Alanine versus Glycine in Protein Denatured States. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 2661−2666. (52) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33. (53) Heinig, M.; Frishman, D. STRIDE: A Web Server For Secondary Structure Assignment From Known Atomic Coordinates Of Proteins. Nucleic Acids Res. 2004, 32, W500−W502.

(13) Cerutti, D. S.; Swope, W. C.; Rice, J. E.; Case, D. A. ff14ipq: A Self Consistent Force Field for Condensed-Phase Simulations of Proteins. J. Chem. Theory Comput. 2014, 10, 4515−4534. (14) Vymetal, J.; Vondrasek, J. Critical Assessment of Current Force Fields. Short Peptide Test Case. J. Chem. Theory Comput. 2013, 9, 441−451. (15) Mu, Y. G.; Kosov, D. S.; Stock, G. Conformational Dynamics of Trialanine in Water. 2. Comparison of AMBER, CHARMM, GROMOS, and OPLS Force Fields to NMR and Infrared Experiments. J. Phys. Chem. B 2003, 107, 5064−5073. (16) Raucci, R.; Colonna, G.; Castello, G.; Costantini, S. Peptide Folding Problem: A Molecular Dynamics Study on Polyalanines Using Different Force Fields. Int. J. Pept. Res. Ther. 2013, 19, 117−123. (17) Best, R. B.; Buchete, N. V.; Hummer, G. Are Current Molecular Dynamics Force Fields Too Helical? Biophys. J. 2008, 95, L7−L09. (18) Matthes, D.; de Groot, B. L. Secondary Structure Propensities in Peptide Folding Simulations: A Systematic Comparison of Molecular Mechanics Interaction Schemes. Biophys. J. 2009, 97, 599−608. (19) Sakae, Y.; Okamoto, Y. Folding Simulations of Three Proteins Having All Alpha-Helix, All Beta-Stand And Alpha/Beta Structures. Mol. Simul. 2010, 36, 302−310. (20) Freddolino, P. L.; Park, S.; Roux, B.; Schulten, K. Force Field Bias In Protein Folding Simulations. Biophys. J. 2009, 96, 3772−3780. (21) Teplow, D. B.; Lazo, N. D.; Bitan, G.; Bernstein, S.; Wyttenbach, T.; Bowers, M. T.; Baumketner, A.; Shea, J. E.; Urbanc, B.; Cruz, L.; Borreguero, J.; Stanley, H. E. Elucidating Amyloid βProtein Folding and Assembly: A Multidisciplinary Approach. Acc. Chem. Res. 2006, 39, 635−645. (22) Roychaudhuri, R.; Yang, M. F.; Condron, M. M.; Teplow, D. B. Structural Dynamics of the Amyloid β-protein Monomer Folding Nucleus. Biochemistry 2012, 51, 3957−3959. (23) Murray, M. M.; Krone, M. G.; Bernstein, S. L.; Baumketner, A.; Condron, M. M.; Lazo, N. D.; Teplow, D. B.; Wyttenbach, T.; Shea, J. E.; Bowers, M. T. Amyloid-β protein: Experiment and Theory on the 21−30 Fragment. J. Phys. Chem. B 2009, 113, 6041−6046. (24) Lazo, N. D.; Grant, M. A.; Condron, M. C.; Rigby, A. C.; Teplow, D. B, On the Nucleation of Amyloid β-protein Monomer Folding. Protein Sci. 2005, 14, 1581−1596. (25) Fawzi, N. L.; Phillips, A. H.; Ruscio, J. Z.; Doucleff, M.; Wemmer, D. E.; Head-Gordon, T. Structure and Dynamics of the Alzheimer’s Aβ21−30 Peptide from the Interplay of NMR Experiments and Simulation. J. Am. Chem. Soc. 2008, 130, 6145−6158. (26) Rao, J. S.; Cruz, L. Effects of Confinement on the Structure and Dynamics of an Intrinsically Disordered Peptide: A MolecularDynamics Study. J. Phys. Chem. B 2013, 117, 3707−3719. (27) Krone, M. G.; Baumketner, A.; Bernstein, S. L.; Wyttenbach, T.; Lazo, N. D.; Teplow, D. B.; Bowers, M. T.; Shea, J. E. Effects of Familial Alzheimer’s Disease Mutations on the Folding Nucleation of the Amyloid-ß protein. J. Mol. Biol. 2008, 381, 221−228. (28) Cruz, L.; Rao, J. S.; Teplow, D. B.; Urbanc, B. Dynamics of Metastable β-Hairpin Structures in the Folding Nucleus of Amyloid βProtein. J. Phys. Chem. B 2012, 116, 6311−6325. (29) Cruz, L.; Urbanc, B.; Borreguero, J. M.; Lazo, N. D.; Teplow, D. B.; Stanley, H. E. Solvent and Mutation Effects on the Nucleation of Amyloid Beta-Protein Folding. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 18258−18263. (30) Chen, W.; Mousseau, N.; Derreumaux, P. The Conformations of the Amyloid-beta (21−30) Fragment Can Be Described By Three Families In Solution. J. Chem. Phys. 2006, 125, 084911. (31) Borreguero, J. M.; Urbanc, B.; Lazo, N. D.; Buldyrev, S. V.; Teplow, D. B.; Stanley, H. E. Folding Events in the 21−30 Region of Amyloid-Beta-Protein (Aβ) Studied In Silico. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6015−6020. (32) Baumketner, A.; Bernstein, S. L.; Wyttenbach, T.; Lazo, N. D.; Teplow, D. B.; Bowers, M. T.; Shea, J. E. Structure of the 21−30 Fragment of Amyloid β-protein. Protein Sci. 2006, 15, 1239−1247. (33) Tarus, B.; Straub, J. E.; Thirumalai, D. Structures and the Free Energy Landscapes Of The Wild Type And Mutants of the Aβ21−30peptide Are Determined By An Interplay Between Intra-Peptide 2594

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595

Article

Journal of Chemical Information and Modeling (54) Frishman, D.; Argos, P. Knowledge-Based Protein Secondary Structure Assignment. Proteins: Struct., Funct., Genet. 1995, 23, 566− 579. (55) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696−3713. (56) Yoda, T.; Sugita, Y.; Okamoto, Y. Secondary-Structure Preferences Of Force Fields For Proteins Evaluated By GeneralizedEnsemble Simulations. Chem. Phys. 2004, 307, 269−283. (57) Duan, L. L.; Mei, Y.; Li, Y. L.; Zhang, Q. G.; Zhang, D. W.; Zhang, J. Z. Simulation Of The Thermodynamics Of Folding And Unfolding Of The Trp-Cage Mini-Protein TC5b Using Different Combinations Of Force Fields And Solvation Models. Sci. China: Chem. 2010, 53, 196−201. (58) Nguyen, P. H.; Li, M. S.; Derreumaux, P. Effects of All-Atom Force Fields On Amyloid Oligomerization: Replica Exchange Molecular Dynamics Simulations of the Aβ(16−22) Dimer And Trimer. Phys. Chem. Chem. Phys. 2011, 13, 9778−9788. (59) Kyte, J.; Doolittle, R. F. A Simple Method For Displaying The Hydropathic Character Of A Protein. J. Mol. Biol. 1982, 157, 105−132. (60) Hopp, T. P.; Woods, K. R. Prediction of Protein Antigenic Determinants From Amino Acid Sequences. Proc. Natl. Acad. Sci. U. S. A. 1981, 78, 3824−3828. (61) Florova, P.; Sklenovsky, P.; Banas, P.; Otyepka, M. Explicit Water Models Affect the Specific Solvation and Dynamics of Unfolded Peptides While the Conformational Behavior and Flexibility of Folded Peptides Remain Intact. J. Chem. Theory Comput. 2010, 6, 3569−3579. (62) Honda, S.; Akiba, T.; Kato, Y. S.; Sawada, Y.; Sekijima, M.; Ishimura, M.; Ooishi, A.; Watanabe, H.; Odahara, T.; Harata, K. Crystal Structure of a Ten-Amino Acid Protein. J. Am. Chem. Soc. 2008, 130, 15327−15331.

2595

DOI: 10.1021/acs.jcim.5b00308 J. Chem. Inf. Model. 2015, 55, 2587−2595