Force field bias in protein folding simulations

Protective V127 prion variant prevents prion disease by interrupting the formation of dimer and fibril from molecular dynamics simulations

Recent studies uncovered a novel protective prion protein variant: V127 variant, which was reported intrinsically resistant to prion conversion and propagation. However, the structural basis of its protective effect is still unknown. To uncover the origin of the protective role of V127 variant, molecular dynamics simulations were performed to explore the influence of G127V mutation on two key processes of prion propagation: dimerization and fibril formation. The simulation results indicate V127 variant is unfavorable to form dimer by reducing the main-chain H-bond interactions. The simulations of formed fibrils consisting of β1 strand prove V127 variant will make the formed fibril become unstable and disorder. The weaker interaction energies between layers and reduced H-bonds number for V127 variant reveal this mutation is unfavorable to the formation of stable fibril. Consequently, we find V127 variant is not only unfavorable to the formation of dimer but also unfavorable to the formation of stable core and fibril, which can explain the mechanism on the protective role of V127 variant from the molecular level. Our findings can deepen the understanding of prion disease and may guide the design of peptide mimetics or small molecule to mimic the protective effect of V127 variant.

Accelerating the Conformational Sampling of Intrinsically Disordered Proteins

Intrinsically disordered proteins (IDPs) are a class of proteins lacking a well-defined secondary structure. Instead, they are able to attain multiple conformations, bind to multiple targets, and respond to changes in their surroundings. Functionally, IDPs have been associated with molecular recognition, cell regulation, and signal transduction. The dynamic conformational ensemble of IDPs is highly environmental and binding partner dependent, rendering the characterization of IDPs extremely challenging. Here, we compare the sampling efficiencies of conventional molecular dynamics (MD), well-tempered metadynamics (WT-META), and bias-exchange metadynamics (BE-META). The total simulation time was over 10 μs, and a 20-mer peptide derived from the Neh2 domain of the Nuclear factor erythroid 2-related factor 2 (Nrf2) protein was simulated. BE-META, with a neutral replica and seven biased replicas employing a set of seven relevant collective variables (CVs), provided the most reliable and efficient sampling. Finally, we propose a free-energy reconstruction method based on the probability distribution of the secondary structure contents. This postprocessing analysis confirms the presence of not only the β-hairpin conformation of the free Neh2 peptide but also its rare bound-state-like conformation, both of that have been experimentally observed. In addition, our simulations also predict other possible conformations to be verified with future experiments.

Force-Field Induced Bias in the Structure of Aβ21-30: A Comparison of OPLS, AMBER, CHARMM, and GROMOS Force Fields

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.

Multiple Simulated Annealing-Molecular Dynamics (MSA-MD) for Conformational Space Search of Peptide and Miniprotein

Protein and peptide structure predictions are of paramount importance for understanding their functions, as well as the interactions with other molecules. However, the use of molecular simulation techniques to directly predict the peptide structure from the primary amino acid sequence is always hindered by the rough topology of the conformational space and the limited simulation time scale. We developed here a new strategy, named Multiple Simulated Annealing-Molecular Dynamics (MSA-MD) to identify the native states of a peptide and miniprotein. A cluster of near native structures could be obtained by using the MSA-MD method, which turned out to be significantly more efficient in reaching the native structure compared to continuous MD and conventional SA-MD simulation.

Selective refinement and selection of near-native models in protein structure prediction

In recent years in silico protein structure prediction reached a level where fully automated servers can generate large pools of near-native structures. However, the identification and further refinement of the best structures from the pool of models remain problematic. To address these issues, we have developed (i) a target-specific selective refinement (SR) protocol; and (ii) molecular dynamics (MD) simulation based ranking (SMDR) method. In SR the all-atom refinement of structures is accomplished via the Rosetta Relax protocol, subject to specific constraints determined by the size and complexity of the target. The best-refined models are selected with SMDR by testing their relative stability against gradual heating through all-atom MD simulations. Through extensive testing we have found that Mufold-MD, our fully automated protein structure prediction server updated with the SR and SMDR modules consistently outperformed its previous versions.

An Analysis of Biomolecular Force Fields for Simulations of Polyglutamine in Solution

Polyglutamine (polyQ) peptides are a useful model system for biophysical studies of protein folding and aggregation, both for their intriguing aggregation properties and their own relevance to human disease. The genetic expansion of a polyQ tract triggers the formation of amyloid aggregates associated with nine neurodegenerative diseases. Several clearly identifiable and separable factors, notably the length of the polyQ tract, influence the mechanism of aggregation, its associated kinetics, and the ensemble of structures formed. Atomistic simulations are well positioned to answer open questions regarding the thermodynamics and kinetics of polyQ folding and aggregation. The additional, explicit representation of water permits deeper investigation of the role of solvent dynamics, and it permits a direct comparison of simulation results with infrared spectroscopy experiments. The generation of meaningful simulation results hinges on satisfying two essential criteria: achieving sufficient conformational sampling to draw statistically valid conclusions, and accurately reproducing the intermolecular forces that govern system structure and dynamics. In this work, we examine the ability of 12 biomolecular force fields to reproduce the properties of a simple, 30-residue polyQ peptide (Q30) in explicit water. In addition to secondary and tertiary structure, we consider generic structural properties of polymers that provide additional dimensions for analysis of the highly degenerate disordered states of the molecule. We find that the 12 force fields produce a wide range of predictions. We identify AMBER ff99SB, AMBER ff99SB*, and OPLS-AA/L to be most suitable for studies of polyQ folding and aggregation.

Assessment of the utility of contact-based restraints in accelerating the prediction of protein structure using molecular dynamics simulations

Molecular dynamics (MD) simulation is a well-established tool for the computational study of protein structure and dynamics, but its application to the important problem of protein structure prediction remains challenging, in part because extremely long timescales can be required to reach the native structure. Here, we examine the extent to which the use of low-resolution information in the form of residue-residue contacts, which can often be inferred from bioinformatics or experimental studies, can accelerate the determination of protein structure in simulation. We incorporated sets of 62, 31, or 15 contact-based restraints in MD simulations of ubiquitin, a benchmark system known to fold to the native state on the millisecond timescale in unrestrained simulations. One-third of the restrained simulations folded to the native state within a few tens of microseconds-a speedup of over an order of magnitude compared with unrestrained simulations and a demonstration of the potential for limited amounts of structural information to accelerate structure determination. Almost all of the remaining ubiquitin simulations reached near-native conformations within a few tens of microseconds, but remained trapped there, apparently due to the restraints. We discuss potential methodological improvements that would facilitate escape from these near-native traps and allow more simulations to quickly reach the native state. Finally, using a target from the Critical Assessment of protein Structure Prediction (CASP) experiment, we show that distance restraints can improve simulation accuracy: In our simulations, restraints stabilized the native state of the protein, enabling a reasonable structural model to be inferred.

Secondary Structure of Rat and Human Amylin across Force Fields

The aggregation of human amylin has been strongly implicated in the progression of Type II diabetes. This 37-residue peptide forms a variety of secondary structures, including random coils, α-helices, and β-hairpins. The balance between these structures depends on the chemical environment, making amylin an ideal candidate to examine inherent biases in force fields. Rat amylin differs from human amylin by only 6 residues; however, it does not form fibrils. Therefore it provides a useful complement to human amylin in studies of the key events along the aggregation pathway. In this work, the free energy of rat and human amylin was determined as a function of α-helix and β-hairpin content for the Gromos96 53a6, OPLS-AA/L, CHARMM22/CMAP, CHARMM22*, Amberff99sb*-ILDN, and Amberff03w force fields using advanced sampling techniques, specifically bias exchange metadynamics. This work represents a first systematic attempt to evaluate the conformations and the corresponding free energy of a large, clinically relevant disordered peptide in solution across force fields. The NMR chemical shifts of rIAPP were calculated for each of the force fields using their respective free energy maps, allowing us to quantitatively assess their predictions. We show that the predicted distribution of secondary structures is sensitive to the choice of force-field: Gromos53a6 is biased towards β-hairpins, while CHARMM22/CMAP predicts structures that are overly α-helical. OPLS-AA/L favors disordered structures. Amberff99sb*-ILDN, AmberFF03w and CHARMM22* provide the balance between secondary structures that is most consistent with available experimental data. In contrast to previous reports, our findings suggest that the equilibrium conformations of human and rat amylin are remarkably similar, but that subtle differences arise in transient alpha-helical and beta-strand containing structures that the human peptide can more readily adopt. We hypothesize that these transient states enable dynamic pathways that facilitate the formation of aggregates and, eventually, amyloid fibrils.

Force field-dependent solution properties of glycine oligomers

Molecular simulations can be used to study disordered polypeptide systems and to generate hypotheses on the underlying structural and thermodynamic mechanisms that govern their function. As the number of disordered protein systems investigated with simulations increase, it is important to understand how particular force fields affect the structural properties of disordered polypeptides in solution. To this end, we performed a comparative structural analysis of Gly(3) and Gly(10) in aqueous solution from all atom, microsecond molecular dynamics (MD) simulations using the CHARMM 27 (C27), CHARMM 36 (C36), and Amber ff12SB force fields. For each force field, Gly(3) and Gly(10) were simulated for at least 300 ns and 1 μs, respectively. Simulating oligoglycines of two different lengths allows us to evaluate how force field effects depend on polypeptide length. Using a variety of structural metrics (e.g., end-to-end distance, radius of gyration, dihedral angle distributions), we characterize the distribution of oligoglycine conformers for each force field and show that each sample conformation space differently, yielding considerably different structural tendencies of the same oligoglycine model in solution. Notably, we find that C36 samples more extended oligoglycine structures than both C27 and ff12SB.

Conformational Dynamics of Two Natively Unfolded Fragment Peptides: Comparison of the AMBER and CHARMM Force Fields

Physics-based force fields are the backbone of molecular dynamics simulations. In recent years, significant progress has been made in the assessment and improvement of commonly used force fields for describing conformational dynamics of folded proteins. However, the accuracy for the unfolded states remains unclear. The latter is however important for detailed studies of protein folding pathways, conformational transitions involving unfolded states, and dynamics of intrinsically disordered proteins. In this work, we compare the three commonly used force fields, AMBER ff99SB-ILDN, CHARMM22/CMAP, and CHARMM36, for modeling the natively unfolded fragment peptides, NTL9(1-22) and NTL9(6-17), using explicit-solvent replica-exchange molecular dynamics simulations. All three simulations show that NTL9(6-17) is completely unstructured, while NTL9(1-22) transiently samples various β-hairpin states, reminiscent of the first β-hairpin in the structure of the intact NTL9 protein. The radius of gyration of the two peptides is force field independent but likely underestimated due to the current deficiency of additive force fields. Compared to the CHARMM force fields, ff99SB-ILDN gives slightly higher β-sheet propensity and more native-like residual structures for NTL9(1-22), which may be attributed to its known β preference. Surprisingly, only two sequence-local pairs of charged residues make appreciable ionic contacts in the simulations of NTL9(1-22), which are sampled slightly more by the CHARMM force fields. Taken together, these data suggest that the current CHARMM and AMBER force fields are globally in agreement in modeling the unfolded states corresponding to β-sheet in the folded structure, while differing in details such as the native-likeness of the residual structures and interactions.

Residue-Specific Force Field (RSFF2) Improves the Modeling of Conformational Behavior of Peptides and Proteins

A recent report of (3)J(HNHα) scalar coupling constants for hundreds of two-residue peptides has provided an important opportunity to test simulation force fields for proteins. Here, we compare the abilities of three derivatives of the Amber ff99SB force field to reproduce these data. We report molecular dynamics (MD) simulations of 256 two-residue peptides and show that the recently developed residue-specific force field (RSFF2) produces a dramatic improvement in the agreement with experimental (3)J(HNHα) coupling constants. We further show that RSFF2 also appears to produce a modest improvement in reproducing the (3)J(HNHα) coupling constants of five model proteins. Perhaps surprisingly, an analysis of neighboring residue effects (NREs) on the (3)J(HNHα) coupling constants of the two-residue peptides indicates little difference between the force fields' abilities to reproduce experimental NREs. We speculate that this might indicate limitations in the force fields' descriptions of nonbonded interactions between adjacent side chains or with terminal capping groups.

CHARMM additive and polarizable force fields for biophysics and computer-aided drug design

BACKGROUND

Molecular Mechanics (MM) is the method of choice for computational studies of biomolecular systems owing to its modest computational cost, which makes it possible to routinely perform molecular dynamics (MD) simulations on chemical systems of biophysical and biomedical relevance.

SCOPE OF REVIEW

As one of the main factors limiting the accuracy of MD results is the empirical force field used, the present paper offers a review of recent developments in the CHARMM additive force field, one of the most popular biomolecular force fields. Additionally, we present a detailed discussion of the CHARMM Drude polarizable force field, anticipating a growth in the importance and utilization of polarizable force fields in the near future. Throughout the discussion emphasis is placed on the force fields' parametrization philosophy and methodology.

MAJOR CONCLUSIONS

Recent improvements in the CHARMM additive force field are mostly related to newly found weaknesses in the previous generation of additive force fields. Beyond the additive approximation is the newly available CHARMM Drude polarizable force field, which allows for MD simulations of up to 1μs on proteins, DNA, lipids and carbohydrates.

GENERAL SIGNIFICANCE

Addressing the limitations ensures the reliability of the new CHARMM36 additive force field for the types of calculations that are presently coming into routine computational reach while the availability of the Drude polarizable force fields offers an inherently more accurate model of the underlying physical forces driving macromolecular structures and dynamics. This article is part of a Special Issue entitled "Recent developments of molecular dynamics".

Termini capping of metal-poly-His peptide complexes induces the formation of α-helix

Capping of both N- and C-terminal induce α-helix formation in Cu²⁺-His₆ peptide.

Current status of protein force fields for molecular dynamics simulations

The current status of classical force fields for proteins is reviewed. These include additive force fields as well as the latest developments in the Drude and AMOEBA polarizable force fields. Parametrization strategies developed specifically for the Drude force field are described and compared with the additive CHARMM36 force field. Results from molecular simulations of proteins and small peptides are summarized to illustrate the performance of the Drude and AMOEBA force fields.

Orientations of residues along the β-arch of self-assembled amylin fibril-like structures lead to polymorphism

Amylin is an endocrine hormone peptide that consists of 37 residues and is the main component of extracellular amyloid deposits found in the pancreas of most type 2 diabetes patients. Amylin peptides are self-assembled to form oligomers and fibrils. So far, four different molecular structures of the self-assembled amylin fibrils have been observed experimentally: two ssNMR models and two crystal models. This study reveals, for the first time, that there are four self-assembled amylin forms that differ in the orientations of the side chains along the β-arch and are all derived from the two ssNMR models. The two ssNMR models are composed of these four different self-assembled forms of amylin, and the two crystal models are composed of two different self-assembled forms of amylin. This study illustrates at the atomic level the differences among the four experimental models and proposes eight new models of self-assembled amylin that are also composed of the four different self-assembled forms of amylin. Our results show polymorphism of the self-assembled fibril-like amylin, with a slight preference of some of the newly constructed models over the experimental models. Finally, we propose that two different self-assembled fibril-like forms of amylin can interact to form a new fibril-like amylin. We investigated this argument and found that some fibril-like amylin prefers to interact to form stable fibril-like structures, whereas others disfavor it. Our work provides new insights that may suggest strategies for future pharmacological studies that aim to find ways to ameliorate the interactions between polymorphic oligomers and fibrils of amylin.

Some practical approaches to treating electrostatic polarization of proteins

Conspectus Electrostatic interaction plays a significant role in determining many properties of biomolecules, which exist and function in aqueous solution, a highly polar environment. For example, proteins are composed of amino acids with charged, polar, and nonpolar side chains and their specific electrostatic properties are fundamental to the structure and function of proteins. An important issue that arises in computational study of biomolecular interaction and dynamics based on classical force field is lack of polarization. Polarization is a phenomenon in which the charge distribution of an isolated molecule will be distorted when interacting with another molecule or presented in an external electric field. The distortion of charge distribution is intended to lower the overall energy of the molecular system, which is counter balanced by the increased internal energy of individual molecules due to the distorted charge distributions. The amount of the charge redistribution, which characterizes the polarizability of a molecule, is determined by the level of the charge distortion. Polarization is inherently quantum mechanical, and therefore classical force fields with fixed atomic charges are incapable of capturing this important effect. As a result, simulation studies based on popular force fields, AMBER, CHARMM, etc., lack the polarization effect, which is a widely known deficiency in most computational studies of biomolecules today. Many efforts have been devoted to remedy this deficiency, such as adding additional movable charge on the atom, allowing atomic charges to fluctuate, or including induced multipoles. Although various successes have been achieved and progress at various levels has been reported over the past decades, the issue of lacking polarization in force field based simulations is far from over. For example, some of these methods do not always give converged results, and other methods require huge computational cost. This Account reviews recent work on developing polarized and polarizable force fields based on fragment quantum mechanical calculations for proteins. The methods described here are based on quantum mechanical calculations of proteins in solution, but with a different level of rigor and different computational efficiency for the molecular dynamics applications. In the general approach, a fragment quantum mechanical calculation for protein with implicit solvation is carried out to derive a polarized protein-specific charge (PPC) for any given protein structure. The PPC correctly reflects the polarization state of the protein in a given conformation, and it can also be dynamically changed as the protein changes conformation in dynamics simulations. Another approach that is computationally more efficient is the effective polarizable bond method in which only polar bonds or groups can be polarized and their polarizabilities are predetermined from quantum mechanical calculations of these groups in external electric fields. Both methods can be employed for applications in various situations by taking advantage of their unique features.

Stability of amyloid oligomers

Molecular simulations are now commonly used to complement experimental techniques in investigating amyloids and their role in human diseases. In this chapter, we will summarize techniques and approaches often used in amyloid simulations and will present recent success stories. Our examples will be focused on lessons learned from molecular dynamics simulations in aqueous environments that start from preformed aggregates. These studies explore the limitations that arise from the choice of force field, the role of mutations in the growth of amyloid aggregates, segmental polymorphism, and the importance of cross-seeding. Furthermore, they give evidence for potential toxicity mechanisms. We finally discuss the role of molecular simulations in the search for aggregation inhibitors.

Statistical mechanics of the denatured state of a protein using replica-averaged metadynamics

The characterization of denatured states of proteins is challenging because the lack of permanent structure in these states makes it difficult to apply to them standard methods of structural biology. In this work we use all-atom replica-averaged metadynamics (RAM) simulations with NMR chemical shift restraints to determine an ensemble of structures representing an acid-denatured state of the 86-residue protein ACBP. This approach has enabled us to reach convergence in the free energy landscape calculations, obtaining an ensemble of structures in relatively accurate agreement with independent experimental data used for validation. By observing at atomistic resolution the transient formation of native and non-native structures in this acid-denatured state of ACBP, we rationalize the effects of single-point mutations on the folding rate, stability, and transition-state structures of this protein, thus characterizing the role of the unfolded state in determining the folding process.

Calculation of the "absolute" free energy of a β-hairpin in an all-atom force field

We propose a new approach to calculate the conformational free energy of a macromolecule in a compact stable state in implicit solvent. The free energy is evaluated with respect to an artificial reference system without internal interactions that is confined to a small well-defined multidimensional volume of a regular shape occupying approximately the same part of the conformational space as the macromolecule of interest. We present a practical implementation of our method, successfully apply it to a β-hairpin in all-atom representation, verify the results with direct parallel tempering simulations, and discuss the possibilities of further improvements.

Perspective: Coarse-grained models for biomolecular systems

By focusing on essential features, while averaging over less important details, coarse-grained (CG) models provide significant computational and conceptual advantages with respect to more detailed models. Consequently, despite dramatic advances in computational methodologies and resources, CG models enjoy surging popularity and are becoming increasingly equal partners to atomically detailed models. This perspective surveys the rapidly developing landscape of CG models for biomolecular systems. In particular, this review seeks to provide a balanced, coherent, and unified presentation of several distinct approaches for developing CG models, including top-down, network-based, native-centric, knowledge-based, and bottom-up modeling strategies. The review summarizes their basic philosophies, theoretical foundations, typical applications, and recent developments. Additionally, the review identifies fundamental inter-relationships among the diverse approaches and discusses outstanding challenges in the field. When carefully applied and assessed, current CG models provide highly efficient means for investigating the biological consequences of basic physicochemical principles. Moreover, rigorous bottom-up approaches hold great promise for further improving the accuracy and scope of CG models for biomolecular systems.

Observation of complete pressure-jump protein refolding in molecular dynamics simulation and experiment

Density is an easily adjusted variable in molecular dynamics (MD) simulations. Thus, pressure-jump (P-jump)-induced protein refolding, if it could be made fast enough, would be ideally suited for comparison with MD. Although pressure denaturation perturbs secondary structure less than temperature denaturation, protein refolding after a fast P-jump is not necessarily faster than that after a temperature jump. Recent P-jump refolding experiments on the helix bundle λ-repressor have shown evidence of a <3 μs burst phase, but also of a ~1.5 ms "slow" phase of refolding, attributed to non-native helical structure frustrating microsecond refolding. Here we show that a λ-repressor mutant is nonetheless capable of refolding in a single explicit solvent MD trajectory in about 19 μs, indicating that the burst phase observed in experiments on the same mutant could produce native protein. The simulation reveals that after about 18.5 μs of conformational sampling, the productive structural rearrangement to the native state does not occur in a single swift step but is spread out over a brief series of helix and loop rearrangements that take about 0.9 μs. Our results support the molecular time scale inferred for λ-repressor from near-downhill folding experiments, where transition-state population can be seen experimentally, and also agrees with the transition-state transit time observed in slower folding proteins by single-molecule spectroscopy.

Fast protein folding kinetics

Fast-folding proteins have been a major focus of computational and experimental study because they are accessible to both techniques: they are small and fast enough to be reasonably simulated with current computational power, but have dynamics slow enough to be observed with specially developed experimental techniques. This coupled study of fast-folding proteins has provided insight into the mechanisms, which allow some proteins to find their native conformation well <1 ms and has uncovered examples of theoretically predicted phenomena such as downhill folding. The study of fast folders also informs our understanding of even 'slow' folding processes: fast folders are small; relatively simple protein domains and the principles that govern their folding also govern the folding of more complex systems. This review summarizes the major theoretical and experimental techniques used to study fast-folding proteins and provides an overview of the major findings of fast-folding research. Finally, we examine the themes that have emerged from studying fast folders and briefly summarize their application to protein folding in general, as well as some work that is left to do.

Structural diversity and initial oligomerization of PrP106-126 studied by replica-exchange and conventional molecular dynamics simulations

Prion diseases are marked by cerebral accumulation of the abnormal isoform of the prion protein. A fragment of prion protein composed of residues 106–126 (PrP106–126) exhibits similar properties to full length prion and plays a key role in the conformational conversion from cellular prion to its pathogenic pattern. Soluble oligomers of PrP106–126 have been proposed to be responsible for neurotoxicity. However, the monomeric conformational space and initial oligomerization of PrP106–126 are still obscure, which are very important for understanding the conformational conversion of PrP106–126. In this study, replica exchange molecular dynamics simulations were performed to investigate monomeric and dimeric states of PrP106–126 in implicit solvent. The structural diversity of PrP106–126 was observed and this peptide did not acquire stable structure. The dimeric PrP106–126 also displayed structural diversity and hydrophobic interaction drove the dimerization. To further study initial oligomerization of PrP106–126, 1 µs conventional molecular dynamics simulations of trimer and tetramer formation were carried out in implicit solvent. We have observed the spontaneous formation of several basic oligomers and stable oligomers with high β-sheet contents were sampled in the simulations of trimer and tetramer formation. The β-hairpin formed in hydrophobic tail of PrP106–126 with residues 118–120 in turn may stabilize these oligomers and seed the formation oligomers. This study can provide insight into the detailed information about the structure of PrP106–126 and the dynamics of aggregation of monomeric PrP106–126 into oligomers in atomic level.

Estimation of Interaction Potentials through the Configurational Temperature Formalism

Molecular interaction potentials are difficult to measure experimentally and hard to compute from first principles, especially for large systems such as proteins. It is therefore desirable to estimate the potential energy that underlies a thermodynamic ensemble from simulated or experimentally determined configurations. This inverse problem of statistical mechanics is challenging because the various potential energy terms can exhibit subtle indirect and correlated effects on the resulting ensemble. A direct approach would try to adapt the force field parameters such that the given configurations are highly probable in the resulting ensemble. But this would require a full simulation of the system whenever a parameter changes. We introduce an extension of the configurational temperature formalism that allows us to circumvent these difficulties and efficiently estimate interaction potentials from molecular configurations. We illustrate the approach for various systems including fluids and a coarse-grained protein model.

Predicting the Thermodynamics and Kinetics of Helix Formation in a Cyclic Peptide Model

The peptide Ac-(cyclo-2,6)-R[KAAAD]-NH2 (cyc-RKAAAD) is a short cyclic peptide known to adopt a remarkably stable single turn α-helix in water. Due to its simplicity and the availability of thermodynamic and kinetic experimental data, cyc-RKAAAD poses as an ideal model for evaluating the aptness of current molecular dynamics (MD) simulation methodologies to accurately sample conformations that reproduce experimentally observed properties. In this work, we extensively sample the conformational space of cyc-RKAAAD using microsecond-timescale MD simulations. We characterize the peptide conformational preferences in terms of secondary structure propensities and, using Cartesian-coordinate principal component analysis (cPCA), construct its free energy landscape, thus obtaining a detailed weighted discrimination between the helical and nonhelical subensembles. The cPCA state discrimination, together with a Markov model built from it, allowed us to estimate the free energy of unfolding (-0.57 kJ/mol) and the relaxation time (∼0.435 μs) at 298.15 K, which are in excellent agreement with the experimentally reported values (-0.22 kJ/mol and 0.42 μs, Serrano, A. L.; Tucker, M. J.; Gai, F. J. Phys. Chem. B, 2011, 115, 7472-7478.). Additionally, we present simulations conducted using two enhanced sampling methods: replica-exchange molecular dynamics (REMD) and bias-exchange metadynamics (BE-MetaD). We compare the free energy landscape obtained by these two methods with the results from MD simulations and discuss the sampling and computational gains achieved. Overall, the results obtained attest to the suitability of modern simulation methods to explore the conformational behavior of peptide systems with a high level of realism.

The intrinsic conformational features of amino acids from a protein coil library and their applications in force field development

The local conformational (φ, ψ, χ) preferences of amino acid residues remain an active research area, which are important for the development of protein force fields. In this perspective article, we first summarize spectroscopic studies of alanine-based short peptides in aqueous solution. While most studies indicate a preference for the P(II) conformation in the unfolded state over α and β conformations, significant variations are also observed. A statistical analysis from various coil libraries of high-resolution protein structures is then summarized, which gives a more coherent view of the local conformational features. The φ, ψ, χ distributions of the 20 amino acids have been obtained from a protein coil library, considering both backbone and side-chain conformational preferences. The intrinsic side-chain χ(1) rotamer preference and χ(1)-dependent Ramachandran plot can be generally understood by combining the interaction of the side-chain Cγ/Oγ atom with two neighboring backbone peptide groups. Current all-atom force fields such as AMBER ff99sb-ILDN, ff03 and OPLS-AA/L do not reproduce these distributions well. A method has been developed by combining the φ, ψ plot of alanine with the influence of side-chain χ(1) rotamers to derive the local conformational features of various amino acids. It has been further applied to improve the OPLS-AA force field. The modified force field (OPLS-AA/C) reproduces experimental (3)J coupling constants for various short peptides quite well. It also better reproduces the temperature-dependence of the helix-coil transition for alanine-based peptides. The new force field can fold a series of peptides and proteins with various secondary structures to their experimental structures. MD simulations of several globular proteins using the improved force field give significantly less deviation (RMSD) to experimental structures. The results indicate that the local conformational features from coil libraries are valuable for the development of balanced protein force fields.

Insights into the Binding of Intrinsically Disordered Proteins from Molecular Dynamics Simulation

Intrinsically disordered proteins (IDPs) are a class of protein that, in the native state, possess no well-defined secondary or tertiary structure, existing instead as dynamic ensembles of conformations. They are biologically important, with approximately 20% of all eukaryotic proteins disordered, and found at the heart of many biochemical networks. To fulfil their biological roles, many IDPs need to bind to proteins and/or nucleic acids. And while unstructured in solution, IDPs typically fold into a well-defined three-dimensional structure upon interaction with a binding partner. The flexibility and structural diversity inherent to IDPs makes this coupled folding and binding difficult to study at atomic resolution by experiment alone, and computer simulation currently offers perhaps the best opportunity to understand this process. But simulation of coupled folding and binding is itself extremely challenging; these molecules are large and highly flexible, and their binding partners, such as DNA or cyclins, are also often large. Therefore, their study requires either or both simplified representations and advanced enhanced sampling schemes. It is not always clear that existing simulation techniques, optimized for studying folded proteins, are well-suited to IDPs. In this article, we examine the progress that has been made in the study of coupled folding and binding using molecular dynamics simulation. We summarise what has been learnt, and examine the state of the art in terms of both methodologies and models. We also consider the lessons to be learnt from advances in other areas of simulation and highlight the issues that remain of be addressed.

Relative free enthalpies for point mutations in two proteins with highly similar sequences but different folds

Enveloping distribution sampling was used to calculate free-enthalpy changes associated with single amino acid mutations for a pair of proteins, GA95 and GB95, that show 95% sequence identity yet fold into topologically different structures. Of the L → A, I → F, and L → Y mutations at positions 20, 30, and 45, respectively, of the 56-residue sequence, the first and the last contribute the most to the free-enthalpy difference between the native and non-native sequence-structure combinations, in agreement with the experimental findings for this protein pair. The individual free-enthalpy changes are almost sequence-independent in the four-strand/one-helix structure, the stable form of GB95, while in the three-helix bundle structure, the stable form of GA95, an interplay between residues 20 and 45 is observed.

A coupled two-dimensional main chain torsional potential for protein dynamics: generation and implementation

Main chain torsions of alanine dipeptide are parameterized into coupled 2-dimensional Fourier expansions based on quantum mechanical (QM) calculations at M06 2X/aug-cc-pvtz//HF/6-31G** level. Solvation effect is considered by employing polarizable continuum model. Utilization of the M06 2X functional leads to precise potential energy surface that is comparable to or even better than MP2 level, but with much less computational demand. Parameterization of the 2D expansions is against the full main chain torsion space instead of just a few low energy conformations. This procedure is similar to that for the development of AMBER03 force field, except unique weighting factor was assigned to all the grid points. To avoid inconsistency between quantum mechanical calculations and molecular modeling, the model peptide is further optimized at molecular mechanics level with main chain dihedral angles fixed before the calculation of the conformational energy on molecular mechanical level at each grid point, during which generalized Born model is employed. Difference in solvation models at quantum mechanics and molecular mechanics levels makes this parameterization procedure less straightforward. All force field parameters other than main chain torsions are taken from existing AMBER force field. With this new main chain torsion terms, we have studied the main chain dihedral distributions of ALA dipeptide and pentapeptide in aqueous solution. The results demonstrate that 2D main chain torsion is effective in delineating the energy variation associated with rotations along main chain dihedrals. This work is an implication for the necessity of more accurate description of main chain torsions in the future development of ab initio force field and it also raises a challenge to the development of quantum mechanical methods, especially the quantum mechanical solvation models.

Molecular dynamics simulations of a new branched antimicrobial peptide: a comparison of force fields

Branched antimicrobial peptides are promising as a new class of antibiotics displaying high activity and low toxicity and appear to work through a unique mechanism of action. We explore the structural dynamics of a covalently branched 18 amino acid peptide (referred to as B2088) in aqueous and membrane mimicking environments through molecular dynamics (MD) simulations. Towards this, we carry out conventional MD simulations and supplement these with replica exchange simulations. The simulations are carried out using four different force fields that are commonly employed for simulating biomolecular systems. These force fields are GROMOS53a6, CHARMM27 with cMAP, CHARMM27 without cMAP and AMBER99sb. The force fields are benchmarked against experimental data available from circular dichroism and nuclear magnetic resonance spectroscopies, and show that CHARMM27 without cMAP correction is the most successful in reproducing the structural dynamics of B2088 both in water and in the presence of micelles. Although the four force fields predict different structures of B2088, they all show that B2088 stabilizes against the head group of the lipid through hydrogen bonding of its Lys and Arg side chains. This leads us to hypothesize that B2088 is unlikely to penetrate into the hydrophobic region of the membrane owing to the high free energy costs of transfer from water, and possibly acts by carpeting and thus disrupting the membrane.

Capturing native/native like structures with a physico-chemical metric (pcSM) in protein folding

Specification of the three dimensional structure of a protein from its amino acid sequence, also called a "Grand Challenge" problem, has eluded a solution for over six decades. A modestly successful strategy has evolved over the last couple of decades based on development of scoring functions (e.g. mimicking free energy) that can capture native or native-like structures from an ensemble of decoys generated as plausible candidates for the native structure. A scoring function must be fast enough in discriminating the native from unfolded/misfolded structures, and requires validation on a large data set(s) to generate sufficient confidence in the score. Here we develop a scoring function called pcSM that detects true native structure in the top 5 with 93% accuracy from an ensemble of candidate structures. If we eliminate the native from ensemble of decoys then pcSM is able to capture near native structure (RMSD<=5Ǻ) in top 10 with 86% accuracy. The parameters considered in pcSM are a C-alpha Euclidean metric, secondary structural propensity, surface areas and an intramolecular energy function. pcSM has been tested on 415 systems consisting 142,698 decoys (public and CASP-largest reported hitherto in literature). The average rank for the native is 2.38, a significant improvement over that existing in literature. In-silico protein structure prediction requires robust scoring technique(s). Therefore, pcSM is easily amenable to integration into a successful protein structure prediction strategy. The tool is freely available at http://www.scfbio-iitd.res.in/software/pcsm.jsp.

On the foldability of tryptophan-containing tetra- and pentapeptides: an exhaustive molecular dynamics study

Short peptides serve as minimal model systems to decipher the determinants of foldability due to their simplicity arising from their smaller size, their ability to echo protein-like structural characteristics, and their direct implication in force field validation. Here, we describe an effort to identify small peptides that can still form stable structures in aqueous solutions. We followed the in silico folding of a selected set of 8640 tryptophan-containing tetra- and pentapeptides through 15 210 molecular dynamics simulations amounting to a total of 272.46 μs using explicit representation of the solute and full treatment of the electrostatics. The evaluation and sorting of peptides is achieved through scoring functions, which include terms based on interatomic vector distances, atomic fluctuations, and rmsd matrices between successive frames of a trajectory. Highly scored peptides are studied further via successive simulation rounds of increasing simulation length and using different empirical force fields. Our method suggested only a handful of peptides with strong foldability prognosis. The discrepancies between the predictions of the various force fields for such short sequences are also extensively discussed. We conclude that the vast majority of such short peptides do not adopt significantly stable structures in water solutions, at least based on our computational predictions. The present work can be utilized in the rational design and engineering of bioactive peptides with desired molecular properties.

Misplaced helix slows down ultrafast pressure-jump protein folding

Using a newly developed microsecond pressure-jump apparatus, we monitor the refolding kinetics of the helix-stabilized five-helix bundle protein λ*YA, the Y22W/Q33Y/G46,48A mutant of λ-repressor fragment 6-85, from 3 μs to 5 ms after a 1,200-bar P-drop. In addition to a microsecond phase, we observe a slower 1.4-ms phase during refolding to the native state. Unlike temperature denaturation, pressure denaturation produces a highly reversible helix-coil-rich state. This difference highlights the importance of the denatured initial condition in folding experiments and leads us to assign a compact nonnative helical trap as the reason for slower P-jump-induced refolding. To complement the experiments, we performed over 50 μs of all-atom molecular dynamics P-drop refolding simulations with four different force fields. Two of the force fields yield compact nonnative states with misplaced α-helix content within a few microseconds of the P-drop. Our overall conclusion from experiment and simulation is that the pressure-denatured state of λ*YA contains mainly residual helix and little β-sheet; following a fast P-drop, at least some λ*YA forms misplaced helical structure within microseconds. We hypothesize that nonnative helix at helix-turn interfaces traps the protein in compact nonnative conformations. These traps delay the folding of at least some of the population for 1.4 ms en route to the native state. Based on molecular dynamics, we predict specific mutations at the helix-turn interfaces that should speed up refolding from the pressure-denatured state, if this hypothesis is correct.

Transient helical structure during PI3K and Fyn SH3 domain folding

A growing list of proteins, including the β-sheet-rich SH3 domain, is known to transiently populate a compact α-helical intermediate before settling into the native structure. Examples have been discovered in cryogenic solvent as well as by pressure jumps. Earlier studies of λ repressor mutants showed that transient states with excess helix are robust in an all-α protein. Here we extend a previous study of src SH3 domain to two new SH3 sequences, phosphatidylinositol 3-kinase (PI3K) and a Fyn mutant, to see how robust such helix-rich transients are to sequence variations in this β-sheet fold. We quantify helical structure by circular dichroism (CD), protein compactness by small-angle X-ray scattering (SAXS), and transient helical populations by cryo-stopped-flow CD. Our results show that transient compact helix-rich intermediates are easily accessible on the folding landscape of different SH3 domains. In molecular dynamics simulations, force field errors are often blamed for transient non-native structure. We suggest that experimental examples of very fast α-rich transient misfolding could become a more subtle test for further force field improvements than observation of the native state alone.

Microsecond folding experiments and simulations: a match is made

For the past two decades, protein folding experiments have been speeding up from the second or millisecond time scale to the microsecond time scale, and full-atom simulations have been extended from the nanosecond to the microsecond and even millisecond time scale. Where the two meet, it is now possible to compare results directly, allowing force fields to be validated and refined, and allowing experimental data to be interpreted in atomistic detail. In this perspective we compare recent experiments and simulations on the microsecond time scale, pointing out the progress that has been made in determining native structures from physics-based simulations, refining experiments and simulations to provide more quantitative underlying mechanisms, and tackling the problems of multiple reaction coordinates, downhill folding, and complex underlying structure of unfolded or misfolded states.

Derivation of fixed partial charges for amino acids accommodating a specific water model and implicit polarization

We have developed the IPolQ method for fitting nonpolarizable point charges to implicitly represent the energy of polarization for systems in pure water. The method involves iterative cycles of molecular dynamics simulations to estimate the water charge density around the solute of interest, followed by quantum mechanical calculations at the MP2/cc-pV(T+d)Z level to determine updated solute charges. Lennard-Jones parameters are updated starting from the Amber FF99SB nonbonded parameter set to accommodate the new charge model, guided by the comparisons to experimental hydration free energies (HFEs) of neutral amino acid side chain analogs and assumptions about the computed HFEs for charged side chains. These Lennard-Jones parameter adjustments for side-chain analogs are assumed to be transferable to amino acids generally, and new charges for all standard amino acids are then derived in the presence of water modeled by TIP4P-Ew. Overall, the new charges depict substantially more polarized amino acids, particularly in the backbone moieties, than previous Amber charge sets. Efforts to complete a new force field with appropriate torsion parameters for this charge model are underway. The IPolQ method is general and applicable to arbitrary solutes.

Folding Kinetics and Unfolded State Dynamics of the GB1 Hairpin from Molecular Simulation

The C-terminal β-hairpin of protein G is a 16-residue peptide that folds in a two-state fashion akin to many larger proteins. However, with an experimental folding time of ∼6 μs, it remains a challenging system for all-atom, explicitly solvated, molecular dynamics simulations. Here, we use a large simulation data set (0.7 ms total) of the hairpin at 300 and 350 K to interpret its folding via a master equation approach. We find a separation of over an order of magnitude between the longest and second longest relaxation times, with the slowest relaxation corresponding to folding. However, in spite of this apparent two-state dynamics, the folding rate determined based on a first-passage time analysis depends on the initial conditions chosen, with a nonexponential distribution of first passage times being obtained in some cases. Using the master equation model, we are now able to account quantitatively for the observed distribution of first passage times. The deviation from the expected exponential distribution for a two-state system arises from slow dynamics in the unfolded state, associated with formation and melting of helical structures. Our results help to reconcile recent findings of slow dynamics in unfolded proteins with observed two-state folding kinetics. At the same time, they indicate that care is required in estimating folding kinetics from many short folding simulations. Last, we are able to use the master equation model to obtain details of the folding mechanism and folding transition state, which appear consistent with the "zipper" mechanism inferred from the experiment.

Identification of minimally interacting modules in an intrinsically disordered protein

The conformational characterization of intrinsically disordered proteins (IDPs) is complicated by their conformational heterogeneity and flexibility. If an IDP could somehow be divided into smaller fragments and reconstructed later, theoretical and spectroscopic studies could probe its conformational variability in detail. Here, we used replica molecular-dynamics simulations and network theory to explore whether such a divide-and-conquer strategy is feasible for α-synuclein, a prototypical IDP. We characterized the conformational variability of α-synuclein by conducting >100 unbiased all-atom molecular-dynamics simulations, for a total of >10 μs of trajectories. In these simulations, α-synuclein formed a heterogeneous ensemble of collapsed coil states in an aqueous environment. These states were stabilized by heterogeneous contacts between sequentially distant regions. We find that α-synuclein contains residual secondary structures in the collapsed states, and the heterogeneity in the collapsed state makes it feasible to split α-synuclein into sequentially contiguous minimally interacting fragments. This study reveals previously unknown characteristics of α-synuclein and provides a new (to our knowledge) approach for studying other IDPs.

To milliseconds and beyond: challenges in the simulation of protein folding

Quantitatively accurate all-atom molecular dynamics (MD) simulations of protein folding have long been considered a holy grail of computational biology. Due to the large system sizes and long timescales involved, such a pursuit was for many years computationally intractable. Further, sufficiently accurate forcefields needed to be developed in order to realistically model folding. This decade, however, saw the first reports of folding simulations describing kinetics on the order of milliseconds, placing many proteins firmly within reach of these methods. Progress in sampling and forcefield accuracy, however, presents a new challenge: how to turn huge MD datasets into scientific understanding. Here, we review recent progress in MD simulation techniques and show how the vast datasets generated by such techniques present new challenges for analysis. We critically discuss the state of the art, including reaction coordinate and Markov state model (MSM) methods, and provide a perspective for the future.