Direct folding simulation of alpha-helices and beta-hairpins based on a single all-atom force field with an implicit solvation model

Recently, we have shown that a modified energy model based on the param99 force field with the generalized Born (GB) solvation model produces reliable free energy landscapes of mini-proteins with a betabetaalpha motif (BBA5, 1FSD, and 1PSV), with the native structures of the mini-proteins located in their lowest free energy minimum states. One of the main features in the modified energy model is a significant improvement for more balanced treatments of alpha and beta strands in proteins. In this study, using the replica exchange molecular dynamics (REMD) simulation method with this new force field, we have carried out extensive ab initio folding studies of several well-known peptides with alpha or beta strands (C-peptide, EK-peptide, le0q, and gbl). Starting from fully extended conformations as the initial conditions, all of the native-like structures of the target peptides were successfully identified by REMD, with reasonable representations of free energy surfaces. The present simulation results with the modified energy model are consistent with experiments, demonstrating an extended applicability of the energy model to folding studies of a variety of alpha-helices, beta-strands, and alpha/beta proteins.     

Free energy landscape of protein folding in water: explicit vs. implicit solvent

The Generalized Born (GB) continuum solvent model is arguably the most widely used implicit solvent model in protein folding and protein structure prediction simulations; however, it still remains an open question on how well the model behaves in these large-scale simulations. The current study uses the beta-hairpin from C-terminus of protein G as an example to explore the folding free energy landscape with various GB models, and the results are compared to the explicit solvent simulations and experiments. All free energy landscapes are obtained from extensive conformation space sampling with a highly parallel replica exchange method. Because solvation model parameters are strongly coupled with force fields, five different force field/solvation model combinations are examined and compared in this study, namely the explicit solvent model: OPLSAA/SPC model, and the implicit solvent models: OPLSAA/SGB (Surface GB), AMBER94/GBSA (GB with Solvent Accessible Surface Area), AMBER96/GBSA, and AMBER99/GBSA. Surprisingly, we find that the free energy landscapes from implicit solvent models are quite different from that of the explicit solvent model. Except for AMBER96/GBSA, all other implicit solvent models find the lowest free energy state not the native state. All implicit solvent models show erroneous salt-bridge effects between charged residues, particularly in OPLSAA/SGB model, where the overly strong salt-bridge effect results in an overweighting of a non-native structure with one hydrophobic residue F52 expelled from the hydrophobic core in order to make better salt bridges. On the other hand, both AMBER94/GBSA and AMBER99/GBSA models turn the beta-hairpin in to an alpha-helix, and the alpha-helical content is much higher than the previously reported alpha-helices in an explicit solvent simulation with AMBER94 (AMBER94/TIP3P). Only AMBER96/GBSA shows a reasonable free energy landscape with the lowest free energy structure the native one despite an erroneous salt-bridge between D47 and K50. Detailed results on free energy contour maps, lowest free energy structures, distribution of native contacts, alpha-helical content during the folding process, NOE comparison with NMR, and temperature dependences are reported and discussed for all five models.     

Exploring protein native states and large-scale conformational changes with a modified generalized born model

Implicit solvation models provide, for many applications, a reasonably accurate and computationally effective way to describe the electrostatics of aqueous solvation. Here, a popular analytical Generalized Born (GB) solvation model is modified to improve its accuracy in calculating the solvent polarization part of free energy changes in large-scale conformational transitions, such as protein folding. In contrast to an earlier GB model (implemented in the AMBER-6 program), the improved version does not overstabilize the native structures relative to the finite-difference Poisson-Boltzmann continuum treatment. In addition to improving the energy balance between folded and unfolded conformers, the algorithm (available in the AMBER-7 and NAB molecular modeling packages) is shown to perform well in more than 50 ns of native-state molecular dynamics (MD) simulations of thioredoxin, protein-A, and ubiquitin, as well as in a simulation of Barnase/Barstar complex formation. For thioredoxin, various combinations of input parameters have been explored, such as the underlying gas-phase force fields and the atomic radii. The best performance is achieved with a previously proposed modification to the torsional potential in the Amber ff99 force field, which yields stable native trajectories for all of the tested proteins, with backbone root-mean-square deviations from the native structures being approximately 1.5 A after 6 ns of simulation time. The structure of Barnase/Barstar complex is regenerated, starting from an unbound state, to within 1.9 A relative to the crystal structure of the complex.     

Ab initio folding of proteins with all-atom discrete molecular dynamics

Discrete molecular dynamics (DMD) is a rapid sampling method used in protein folding and aggregation studies. Until now, DMD was used to perform simulations of simplified protein models in conjunction with structure-based force fields. Here, we develop an all-atom protein model and a transferable force field featuring packing, solvation, and environment-dependent hydrogen bond interactions. We performed folding simulations of six small proteins (20-60 residues) with distinct native structures by the replica exchange method. In all cases, native or near-native states were reached in simulations. For three small proteins, multiple folding transitions are observed, and the computationally characterized thermodynamics are in qualitative agreement with experiments. The predictive power of all-atom DMD highlights the importance of environment-dependent hydrogen bond interactions in modeling protein folding. The developed approach can be used for accurate and rapid sampling of conformational spaces of proteins and protein-protein complexes and applied to protein engineering and design of protein-protein interactions.     

Recent advances in implicit solvent-based methods for biomolecular simulations

Implicit solvent-based methods play an increasingly important role in molecular modeling of biomolecular structure and dynamics. Recent methodological developments have mainly focused on the extension of the generalized Born (GB) formalism for variable dielectric environments and accurate treatment of nonpolar solvation. Extensive efforts in parameterization of GB models and implicit solvent force fields have enabled ab initio simulation of protein folding to native or near-native structures. Another exciting area that has benefited from the advances in implicit solvent models is the development of constant pH molecular dynamics methods, which have recently been applied to the calculations of protein pK(a) values and the studies of pH-dependent peptide and protein folding.     

Blind Test of Physics-Based Prediction of Protein Structures

We report here a multiprotein blind test of a computer method to predict native protein structures based solely on an all-atom physics-based force field. We use the AMBER 96 potential function with an implicit (GB/SA) model of solvation, combined with replica-exchange molecular-dynamics simulations. Coarse conformational sampling is performed using the zipping and assembly method (ZAM), an approach that is designed to mimic the putative physical routes of protein folding. ZAM was applied to the folding of six proteins, from 76 to 112 monomers in length, in CASP7, a community-wide blind test of protein structure prediction. Because these predictions have about the same level of accuracy as typical bioinformatics methods, and do not utilize information from databases of known native structures, this work opens up the possibility of predicting the structures of membrane proteins, synthetic peptides, or other foldable polymers, for which there is little prior knowledge of native structures. This approach may also be useful for predicting physical protein folding routes, non-native conformations, and other physical properties from amino acid sequences.     

Force field influences in beta-hairpin folding simulations

All-atom force fields are now routinely used for more detailed understanding of protein folding mechanisms. However, it has been pointed out that use of all-atom force fields does not guarantee more accurate representations of proteins; in fact, sometimes it even leads to biased structural distributions. Indeed, several issues remain to be solved in force field developments, such as accurate treatment of implicit solvation for efficient conformational sampling and proper treatment of backbone interactions for secondary structure propensities. In this study, we first investigate the quality of several recently improved backbone interaction schemes in AMBER for folding simulations of a beta-hairpin peptide, and further study their influences on the peptide's folding mechanism. Due to the significant number of simulations needed for a thorough analysis of tested force fields, the implicit Poisson-Boltzmann solvent was used in all simulations. The chosen implicit solvent was found to be reasonable for studies of secondary structures based on a set of simulations of both alpha-helical and beta-hairpin peptides with the TIP3P explicit solvent as benchmark. Replica exchange molecular dynamics was also utilized for further efficient conformational sampling. Among the tested AMBER force fields, ff03 and a revised ff99 force field were found to produce structural and thermodynamic data in comparably good agreement with the experiment. However, detailed folding pathways, such as the order of backbone hydrogen bond zipping and the existence of intermediate states, are different between the two force fields, leading to force field-dependent folding mechanisms.     

Free energies of protein decoys provide insight into determinants of protein stability

We have calculated the stability of decoy structures of several proteins (from the CASP3 models and the Park and Levitt decoy set) relative to the native structures. The calculations were performed with the force field-consistent ES/IS method, in which an implicit solvent (IS) model is used to calculate the average solvation free energy for snapshots from explicit simulations (ESs). The conformational free energy is obtained by adding the internal energy of the solute from the ESs and an entropic term estimated from the covariance positional fluctuation matrix. The set of atomic Born radii and the cavity-surface free energy coefficient used in the implicit model has been optimized to be consistent with the all-atom force field used in the ESs (cedar/gromos with simple point charge (SPC) water model). The decoys are found to have a consistently higher free energy than that of the native structure; the gap between the native structure and the best decoy varies between 10 and 15 kcal/mole, on the order of the free energy difference that typically separates the native state of a protein from the unfolded state. The correlation between the free energy and the extent to which the decoy structures differ from the native (as root mean square deviation) is very weak; hence, the free energy is not an accurate measure for ranking the structurally most native-like structures from among a set of models. Analysis of the energy components shows that stability is attained as a result of three major driving forces: (1) minimum size of the protein-water surface interface; (2) minimum total electrostatic energy, which includes solvent polarization; and (3) minimum protein packing energy. The detailed fit required to optimize the last term may underlie difficulties encountered in recovering the native fold from an approximate decoy or model structure.     

Exploring the folding energy landscapes of heme proteins using a hybrid AWSEM-heme model

Heme is an active center in many proteins. Here we explore computationally the role of heme in protein folding and protein structure. We model heme proteins using a hybrid model employing the AWSEM Hamiltonian, a coarse-grained forcefield for the protein chain along with AMBER, an all-atom forcefield for the heme. We carefully designed transferable force fields that model the interactions between the protein and the heme. The types of protein–ligand interactions in the hybrid model include thioester covalent bonds, coordinated covalent bonds, hydrogen bonds, and electrostatics. We explore the influence of different types of hemes (heme b and heme c) on folding and structure prediction. Including both types of heme improves the quality of protein structure predictions. The free energy landscape shows that both types of heme can act as nucleation sites for protein folding and stabilize the protein folded state. In binding the heme, coordinated covalent bonds and thioester covalent bonds for heme c drive the heme toward the native pocket. The electrostatics also facilitates the search for the binding site.

     

Role of hydrophilic and hydrophobic contacts in folding of the second beta-hairpin fragment of protein G: molecular dynamics simulation studies of an all-atom model

Predicting the folding mechanism of the second beta-hairpin fragment of the Ig-binding domain B of streptococcal protein G is unexpectedly challenging for simplified reduced models because the models developed so far indicated a different folding mechanism from what was suggested from high-temperature unfolding and equilibrium free-energy surface analysis based on established all-atom empirical force fields in explicit or implicit solvent. This happened despite the use of empirical residue-based interactions, multibody hydrophobic interactions, and inclusions of hydrogen bonding effects in the simplified models. This article employs a recently developed all-atom (except nonpolar hydrogens) model interacting with simple square-well potentials to fold the peptide fragment by molecular dynamics simulation methods. In this study, 193 out of 200 trajectories are folded at two reduced temperatures (3.5 and 3.7) close to the transition temperature T* approximately 4.0. Each simulation takes <7 h of CPU time on a Pentium 800-MHz PC. Folding of the new all-atom model is found to be initiated by collapse before the formation of main-chain hydrogen bonds. This verifies the mechanism proposed from previous all-atom unfolding and equilibrium simulations. The new model further predicts that the collapse is initiated by two nucleation contacts (a hydrophilic contact between D46 and T49 and a hydrophobic contact between Y45 and F52), in agreement with recent NMR measurements. The results suggest that atomic packing and native contact interactions play a dominant role in folding mechanism.     

Multiscale Coarse-Graining of the Protein Energy Landscape

A variety of coarse-grained (CG) models exists for simulation of proteins. An outstanding problem is the construction of a CG model with physically accurate conformational energetics rivaling all-atom force fields. In the present work, atomistic simulations of peptide folding and aggregation equilibria are force-matched using multiscale coarse-graining to develop and test a CG interaction potential of general utility for the simulation of proteins of arbitrary sequence. The reduced representation relies on multiple interaction sites to maintain the anisotropic packing and polarity of individual sidechains. CG energy landscapes computed from replica exchange simulations of the folding of Trpzip, Trp-cage and adenylate kinase resemble those of other reduced representations; non-native structures are observed with energies similar to those of the native state. The artifactual stabilization of misfolded states implies that non-native interactions play a deciding role in deviations from ideal funnel-like cooperative folding. The role of surface tension, backbone hydrogen bonding and the smooth pairwise CG landscape is discussed. Ab initio folding aside, the improved treatment of sidechain rotamers results in stability of the native state in constant temperature simulations of Trpzip, Trp-cage, and the open to closed conformational transition of adenylate kinase, illustrating the potential value of the CG force field for simulating protein complexes and transitions between well-defined structural states.     

Free-energy landscape of the GB1 hairpin in all-atom explicit solvent simulations with different force fields: Similarities and differences

Although it is now possible to fold peptides and miniproteins in molecular dynamics simulations, it is well appreciated that force fields are not all transferable to different proteins. Here, we investigate the influence of the protein force field and the solvent model on the folding energy landscape of a prototypical two‐state folder, the GB1 hairpin. We use extensive replica‐exchange molecular dynamics simulations to characterize the free‐energy surface as a function of temperature. Most of these force fields appear similar at a global level, giving a fraction folded at 300 K between 0.2 and 0.8 in all cases, which is a difference in stability of 2.8 kT, and are generally consistent with experimental data at this temperature. The most significant differences appear in the unfolded state, where there are different residual secondary structures which are populated, and the overall dimensions of the unfolded states, which in most of the force fields are too collapsed relative to experimental Förster Resonance Energy Transfer (FRET) data. Proteins 2011. © 2010 Wiley‐Liss, Inc.     

A combined steepest descent and genetic algorithm (SD/GA) approach for the optimization of solvation parameters

An accurate solvation model is essential for computer modeling of protein folding and other biomolecular self-assembly processes. Compared to explicit solvent models, implicit solvent models, such as the Poisson-Boltzmann (PB) with solvent accessible surface area model (PB/SA), offer a much faster speed—the most compelling reason for the popularity of these implicit solvent models. Since these implicit solvent models typically use empirical parameters, such as atomic radii and the surface tensions, an optimal fit of these parameters is crucial for the final accuracy of properties such as solvation free energy and folding free energy. In this paper, we proposed a combined approach, namely SD/GA, which takes the advantage of both local optimization with the steepest descent (SD), and global optimization with the genetic algorithm (GA), for parameters optimization in multi-dimensional space. The SD/GA method is then applied to the optimization of solvation parameters in the non-polar cavity term of the PB/SA model. The results show that the newly optimized parameters from SD/GA not only increase the accuracy in the solvation free energies for ～200 organic molecules, but also significantly improve the free energy landscape of a β-hairpin folding. The current SD/GA method can be readily applied to other multi-dimensional parameter space optimization as well.     

Modeling loop reorganization free energies of acetylcholinesterase: a comparison of explicit and implicit solvent models

The treatment of hydration effects in protein dynamics simulations varies in model complexity and spans the range from the computationally intensive microscopic evaluation to simple dielectric screening of charge-charge interactions. This paper compares different solvent models applied to the problem of estimating the free-energy difference between two loop conformations in acetylcholinesterase. Molecular dynamics (MD) simulations were used to sample potential energy surfaces of the two basins with solvent treated by means of explicit and implicit methods. Implicit solvent methods studied include the generalized Born (GB) model, atomic solvation potential (ASP), and the distance-dependent dieletric constant. By using the linear response approximation (LRA), the explicit solvent calculations determined a free-energy difference that is in excellent agreement with the experimental estimate, while rescoring the protein conformations with GB or the Poisson equation showed inconsistent and inferior results. While the approach of rescoring conformations from explicit water simulations with implicit solvent models is popular among many applications, it perturbs the energy landscape by changing the solvent contribution to microstates without conformational relaxation, thus leading to non-optimal solvation free energies. Calculations applying MD with a GB solvent model produced results of comparable accuracy as observed with LRA, yet the electrostatic free-energy terms were significantly different due to optimization on a potential energy surface favored by an implicit solvent reaction field. The simpler methods of ASP and the distance-dependent scaling of the dielectric constant both produced considerable distortions in the protein internal free-energy terms and are consequently unreliable.     

Stabilization of coiled-coil peptide domains by introduction of trifluoroleucine

Substitution of leucine residues by 5,5,5-trifluoroleucine at the d-positions of the leucine zipper peptide GCN4-p1d increases the thermal stability of the coiled-coil structure. The midpoint thermal unfolding temperature of the fluorinated peptide is elevated by 13 degrees C at 30 microM peptide concentration. The modified peptide is more resistant to chaotropic denaturants, and the free energy of folding of the fluorinated peptide is 0.5-1.2 kcal/mol larger than that of the hydrogenated form. A similarly fluorinated form of the DNA-binding peptide GCN4-bZip binds to target DNA sequences with affinity and specificity identical to those of the hydrogenated form, while demonstrating enhanced thermal stability. Molecular dynamics simulation on the fluorinated GCN4-p1d peptide using the Surface Generalized Born implicit solvation model revealed that the coiled-coil binding energy is 55% more favorable upon fluorination. These results suggest that fluorination of hydrophobic substructures in peptides and proteins may provide new means of increasing protein stability, enhancing protein assembly, and strengthening receptor-ligand interactions.     

An effective solvent theory connecting the underlying mechanisms of osmolytes and denaturants for protein stability

An all-atom Gō model of Trp-cage protein is simulated using discontinuous molecular dynamics in an explicit minimal solvent, using a single, contact-based interaction energy between protein and solvent particles. An effective denaturant or osmolyte solution can be constructed by making the interaction energy attractive or repulsive. A statistical mechanical equivalence is demonstrated between this effective solvent model and models in which proteins are immersed in solutions consisting of water and osmolytes or denaturants. Analysis of these studies yields the following conclusions: 1), Osmolytes impart extra stability to the protein by reducing the entropy of the unfolded state. 2), Unfolded states in the presence of osmolyte are more collapsed than in water. 3), The folding transition in osmolyte solutions tends to be less cooperative than in water, as determined by the ratio of van 't Hoff to calorimetric enthalpy changes. The decrease in cooperativity arises from an increase in native structure in the unfolded state, and thus a lower thermodynamic barrier at the transition midpoint. 4), Weak denaturants were observed to destabilize small proteins not by lowering the unfolded enthalpy, but primarily by swelling the unfolded state and raising its entropy. However, adding a strong denaturant destabilizes proteins enthalpically. 5), The folding transition in denaturant-containing solutions is more cooperative than in water. 6), Transfer to a concentrated osmolyte solution with purely hard-sphere steric repulsion significantly stabilizes the protein, due to excluded volume interactions not present in the canonical Tanford transfer model. 7), Although a solution with hard-sphere interactions adds a solvation barrier to native contacts, the folding is nevertheless less cooperative for reasons 1–3 above, because a hard-sphere solvent acts as a protecting osmolyte.     

Computer simulations of protein folding by targeted molecular dynamics

We have performed 128 folding and 45 unfolding molecular dynamics runs of chymotrypsin inhibitor 2 (CI2) with an implicit solvation model for a total simulation time of 0.4 microseconds. Folding requires that the three-dimensional structure of the native state is known. It was simulated at 300 K by supplementing the force field with a harmonic restraint which acts on the root-mean-square deviation and allows to decrease the distance to the target conformation. High temperature and/or the harmonic restraint were used to induce unfolding. Of the 62 folding simulations started from random conformations, 31 reached the native structure, while the success rate was 83% for the 66 trajectories which began from conformations unfolded by high-temperature dynamics. A funnel-like energy landscape is observed for unfolding at 475 K, while the unfolding runs at 300 K and 375 K as well as most of the folding trajectories have an almost flat energy landscape for conformations with less than about 50% of native contacts formed. The sequence of events, i.e., secondary and tertiary structure formation, is similar in all folding and unfolding simulations, despite the diversity of the pathways. Previous unfolding simulations of CI2 performed with different force fields showed a similar sequence of events. These results suggest that the topology of the native state plays an important role in the folding process.     

Structural details,pathways, and energetics of unfolding apomyoglobin

Protein folding is often difficult to characterize experimentally because of the transience of intermediate states, and the complexity of the protein-solvent system. Atomistic simulations, which could provide more detailed information, have had to employ highly simplified models or high temperatures, to cope with the long time scales of unfolding; direct simulation of folding is even more problematic. We report a fully atomistic simulation of the acid-induced unfolding of apomyoglobin in which the protonation of acidic side-chains to simulate low pH is sufficient to induce unfolding at room temperature with no added biasing forces or other unusual conditions; and the trajectory is validated by comparison to experimental characterization of intermediate states. Novel insights provided by their analysis include: characterization of a dry swollen globule state forming a barrier to initial unfolding or final folding; observation of cooperativity in secondary and tertiary structure formation and its explanation in terms of dielectric environments; and structural details of the intermediate and the completely unfolded states. These insights involve time scales and levels of structural detail that are presently beyond the range of experiment, but come within reach through the simulation methods described here. An implicit solvation model is used to analyze the energetics of protein folding at various pH and ionic strength values, and a reasonable estimate of folding free energy is obtained. Electrostatic interactions are found to disfavor folding.     

An implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteins

Exploiting recent developments in generalized Born (GB) electrostatics theory, we have reformulated the calculation of the self-electrostatic solvation energy to account for the influence of biological membranes. Consistent with continuum Poisson-Boltzmann (PB) electrostatics, the membrane is approximated as an solvent-inaccessible infinite planar low-dielectric slab. The present membrane GB model closely reproduces the PB electrostatic solvation energy profile across the membrane. The nonpolar contribution to the solvation energy is taken to be proportional to the solvent-exposed surface area (SA) with a phenomenological surface tension coefficient. The proposed membrane GB/SA model requires minor modifications of the pre-existing GB model and appears to be quite efficient. By combining this implicit model for the solvent/bilayer environment with advanced computational sampling methods, like replica-exchange molecular dynamics, we are able to fold and assemble helical membrane peptides. We examine the reliability of this model and approach by applications to three membrane peptides: melittin from bee venom, the transmembrane domain of the M2 protein from Influenza A (M2-TMP), and the transmembrane domain of glycophorin A (GpA). In the context of these proteins, we explore the role of biological membranes (represented as a low-dielectric medium) in affecting the conformational changes in melittin, the tilt of transmembrane peptides with respect to the membrane normal (M2-TMP), helix-to-helix interactions in membranes (GpA), and the prediction of the configuration of transmembrane helical bundles (GpA). The present method is found to perform well in each of these cases and is anticipated to be useful in the study of folding and assembly of membrane proteins as well as in structure refinement and modeling of membrane proteins where a limited number of experimental observables are available.