An accurate potential energy model is crucial for biomolecular simulations. Despite many recent improvements of classical protein force fields, there are remaining key issues: much weaker temperature dependence of folding/unfolding equilibrium and overly collapsed unfolded or disordered states. For the latter problem, a new water model (TIP4P-D) has been proposed to correct the significantly underestimated water dispersion interactions. Here, using TIP4P-D, we reveal problems in current force fields through failures in folding model systems (a polyalanine peptide, Trp-cage, and the GB1 hairpin). By using residue-specific parameters to achieve better match between amino acid sequences and native structures and adding a small H-bond correction to partially compensate the missing many-body effects in α-helix formation, the new RSFF2+ force field with the TIP4P-D water model can excellently reproduce experimental melting curves of both α-helical and β-hairpin systems. The RSFF2+/TIP4P-D method also gives less collapsed unfolded structures and describes well folded proteins simultaneously.

Site-specific phosphorylation of an intrinsically disordered protein, eIF4E-binding protein isoform 2 (4E-BP2), can suppress its native function by folding it into a four-stranded β-sheet, but the mechanism of this phosphorylation-induced folding is unclear. In this work, we use all-atom molecular dynamics simulations to investigate both the folded and unfolded states of 4E-BP2 under different phosphorylation states of T37 and T46. The results show that the phosphorylated forms of both T37 and T46 play important roles in stabilizing the folded structure, especially for the β-turns and the sequestered binding motif. The phosphorylated residues not only guide the folding of the protein through several intermediate states but also affect the conformational distribution of the unfolded ensemble. Significantly, the phosphorylated residues can function as nucleation sites for the folding of the protein by forming certain local structures that are stabilized by hydrogen bonding involving the phosphate group. The region around phosphorylated T46 appears to fold before that around phosphorylated T37. These findings provide new insight into the intricate effects of protein phosphorylation.

RfaH protein functions as both transcription anti-terminator and translation enhancer in bacteria. Recent studies have shown that the C-terminal domain (CTD) is an α-helical hairpin (two-helix bundle) in full-length RfaH, despite the intrinsically favored β-barrel structure. Here, we carried out μs-timescale molecular dynamics (MD) simulations for the wild-type (WT) RfaH, its E48S mutant and an established model without the intrinsically disordered region (IDR1) linking the CTD and the N-terminal domain (NTD). Our simulations showed that the WT can be well stabilized by our RSFF1 force field, while the E48S mutant and the model without IDR1 undergo considerable structural change, which is in good agreement with experimental observations. The IDR1 plays important roles in stabilizing the hydrophobic environment near the crucial E48–R138 salt-bridge as well as in tethering α4 helix in CTD to α3 helix in NTD. In the absence of the IDR1, destabilization of key interdomain contacts and unfolding of the CTD α5 helix were observed in the simulation. In addition, the intrinsically disordered tail of the CTD (IDR2) is also of great significance to stabilize the bound conformation of CTD. These findings provide important implications for consideration of simulations in revealing the functions of residues invisible in a crystal structure.Figure optionsDownload full-size imageDownload as PowerPoint slide

Cyclic peptides (CPs) are promising candidates for drugs, chemical biology tools, and self-assembling nanomaterials. However, the development of reliable and accurate computational methods for their structure prediction has been challenging. Here, 20 all-trans CPs of 5–12 residues selected from Cambridge Structure Database have been simulated using replica-exchange molecular dynamics with four different force fields. Our recently developed residue-specific force fields RSFF1 and RSFF2 can correctly identify the crystal-like conformations of more than half CPs as the most populated conformation. The RSFF2 performs the best, which consistently predicts the crystal structures of 17 out of 20 CPs with rmsd < 1.1 Å. We also compared the backbone (ϕ, ψ) sampling of residues in CPs with those in short linear peptides and in globular proteins. In general, unlike linear peptides, CPs have local conformational free energies and entropies quite similar to globular proteins.

α-Helical hairpin (two-helix bundle) is a structure motif composed of two interacting helices connected by a turn or a short loop. It is an important model for protein folding studies, filling the gap between isolated α-helix and larger all-α domains. Here, we present, for the first time, successful folding simulations of an α-helical hairpin. Our RSFF1 and RSFF2 force fields give very similar predicted structures of this αtα peptide, which is in good agreement with its NMR structure. Our simulations also give site-specific stability of α-helix formation in good agreement with amide hydrogen exchange experiments. Combining the folding free energy landscapes and analyses of structures sampled in five different ranges of the fraction of native contacts (Q), a folding mechanism of αtα is proposed. The most stable sites of Q9-E15 in helix-1 and E24-A30 in helix-2 close to the loop region act as the folding initiation sites. The formation of interhelix side-chain contacts also initiates near the loop region, but some residues in the central parts of the two helices also form contacts quite early. The two termini fold at a final stage, and the loop region remains flexible during the whole folding process. This mechanism is similar to the “zipping out” pathway of β-hairpin folding.

To test whether our recently developed residue-specific force field RSFF2 can reproduce the mutational effect on the thermal stability of Trp-cage mini-protein and decipher its detailed folding mechanism, we carried out long-time replica-exchange molecular dynamics (REMD) simulations on five Trp-cage variants, including TC5b and TC10b. Initiated from their unfolded structures, the simulations not only well-reproduce their experimental structures but also their melting temperatures and folding enthalpies reasonably well. For each Trp-cage variant, the overall folding free energy landscape is apparently two-state, but some intermediate states can be observed when projected on more detailed coordinates. We also found different variants have the same major folding pathway, including the well formed PII-helix in the unfolded state, the formation of W6-P12/P18/P19 contacts and the α-helix before the transition state, the following formation of most native contacts, and the final native loop formation. The folding mechanism derived here is consistent with many previous simulations and experiments.

An important application of all-atom explicit-solvent molecular dynamics (MD) simulations is the refinement of protein structures from low-resolution experiments or template-based modeling. A critical requirement is that the native structure is stable with the force field. We have applied a recently developed residue-specific force field, RSFF1, to a set of 30 refinement targets from recent CASP experiments. Starting from their experimental structures, 1.0 μs unrestrained simulations at 298 K retain most of the native structures quite well except for a few flexible terminals and long internal loops. Starting from each homology model, a 150 ns MD simulation at 380 K generates the best RMSD improvement of 0.85 Å on average. The structural improvements roughly correlate with the RMSD of the initial homology models, indicating possible consistent structure refinement. Finally, targets TR614 and TR624 have been subjected to long-time replica-exchange MD simulations. Significant structural improvements are generated, with RMSD of 1.91 and 1.36 Å with respect to their crystal structures. Thus, it is possible to achieve realistic refinement of protein structure models to near-experimental accuracy, using accurate force field with sufficient conformational sampling.

Recently, we developed a residue-specific force field (RSFF1) based on conformational free-energy distributions of the 20 amino acid residues from a protein coil library. Most parameters in RSFF1 were adopted from the OPLS-AA/L force field, but some van der Waals and torsional parameters that effectively affect local conformational preferences were introduced specifically for individual residues to fit the coil library distributions. Here a similar strategy has been applied to modify the Amber ff99SB force field, and a new force field named RSFF2 is developed. It can successfully fold α-helical structures such as polyalanine peptides, Trp-cage miniprotein, and villin headpiece subdomain and β-sheet structures such as Trpzip-2, GB1 β-hairpins, and the WW domain, simultaneously. The properties of various popular force fields in balancing between α-helix and β-sheet are analyzed based on their descriptions of local conformational features of various residues, and the analysis reveals the importance of accurate local free-energy distributions. Unlike the RSFF1, which overestimates the stability of both α-helix and β-sheet, RSFF2 gives melting curves of α-helical peptides and Trp-cage in good agreement with experimental data. Fitting to the two-state model, RSFF2 gives folding enthalpies and entropies in reasonably good agreement with available experimental results.

Ab initio protein folding via physical-based all-atom simulation is still quite challenging. Using a recently developed residue-specific force field (RSFF1) in explicit solvent, we are able to fold a diverse set of 14 model proteins. The obtained structural features of unfolded state are in good agreement with previous observations. The replica-exchange molecular dynamics simulation is found to be efficient, resulting in multiple folding events for each protein. Transition path time is found to be significantly reduced under elevated temperature.

Traditional protein force fields use one set of parameters for most of the 20 amino acids (AAs), allowing transferability of the parameters. However, a significant shortcoming is the difficulty to fit the Ramachandran plots of all AA residues simultaneously, affecting the accuracy of the force field. In this Feature Article, we report a new strategy for protein force field parametrization. Backbone and side-chain conformational distributions of all 20 AA residues obtained from protein coil library were used as the target data. The dihedral angle (torsion) potentials and some local nonbonded (1-4/1-5/1-6) interactions in OPLS-AA/L force field were modified such that the target data can be excellently reproduced by molecular dynamics simulations of dipeptides (blocked AAs) in explicit water, resulting in a new force field with AA-specific parameters, RSFF1. An efficient free energy decomposition approach was developed to separate the corrections on ϕ and ψ from the two-dimensional Ramachandran plots. RSFF1 is shown to reproduce the experimental NMR 3J-coupling constants of AA dipeptides better than other force fields. It has a good balance between α-helical and β-sheet secondary structures. It can successfully fold a set of α-helix proteins (Trp-cage and Homeodomain) and β-hairpins (Trpzip-2, GB1 hairpin), which cannot be consistently stabilized by other state-of-the-art force fields. Interestingly, the RSFF1 force field systematically overestimates the melting temperature (and the stability of native state) of these peptides/proteins. It has a potential application in the simulation of protein folding and protein structure refinement.