Microsecond simulations of the folding/unfolding thermodynamics of the Trp‐cage miniprotein

We study the unbiased folding/unfolding thermodynamics of the Trp‐cage miniprotein using detailed molecular dynamics simulations of an all‐atom model of the protein in explicit solvent using the Amberff99SB force field. Replica‐exchange molecular dynamics simulations are used to sample the protein ensembles over a broad range of temperatures covering the folded and unfolded states at two densities. The obtained ensembles are shown to reach equilibrium in the 1 μs/replica timescale. The total simulation time used in the calculations exceeds 100 μs. Ensemble averages of the fraction folded, pressure, and energy differences between the folded and unfolded states as a function of temperature are used to model the free energy of the folding transition, ΔG(P, T), over the whole region of temperatures and pressures sampled in the simulations. The ΔG(P, T) diagram describes an ellipse over the range of temperatures and pressures sampled, predicting that the system can undergo pressure‐induced unfolding and cold denaturation at low temperatures and high pressures, and unfolding at low pressures and high temperatures. The calculated free energy function exhibits remarkably good agreement with the experimental folding transition temperature (T_f = 321 K), free energy, and specific heat changes. However, changes in enthalpy and entropy are significantly different than the experimental values. We speculate that these differences may be due to the simplicity of the semiempirical force field used in the simulations and that more elaborate force fields may be required to describe appropriately the thermodynamics of proteins. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

Effective energy function for proteins in lipid membranes

A simple extension of the EEF1 energy function to heterogeneous membrane-aqueous media is proposed. The extension consists of (a) development of solvation parameters for a nonpolar phase using experimental data for the transfer of amino acid side-chains from water to cyclohexane, (b) introduction of a heterogeneous membrane-aqueous system by making the reference solvation free energy of each atom dependent on the vertical coordinate, (c) a modification of the distance-dependent dielectric model to account for reduced screening of electrostatic interactions in the membrane, and (d) an adjustment of the EEF1 aqueous model in light of recent calculations of the potential of mean force between amino acid side-chains in water. The electrostatic model is adjusted to match experimental observations for polyalanine, polyleucine, and the glycophorin A dimer. The resulting energy function (IMM1) reproduces the preference of Trp and Tyr for the membrane interface, gives reasonable energies of insertion into or adsorption onto a membrane, and allows stable 1-ns MD simulations of the glycophorin A dimer. We find that the lowest-energy orientation of melittin in bilayers varies, depending on the thickness of the hydrocarbon layer.     

Linear response theory: An alternative to PB and GB methods for the analysis of molecular dynamics trajectories?

We explore the use of classical Linear Response Theory (LRT) as an alternative strategy to the use of Molecular Mechanics/Poisson-Boltzmann strategies to compute the solvation free energy of macromolecules from molecular dynamics simulations using an explicit representation of solvent. The method reproduces well the free energy of solvation of standard amino acid side chains, small peptides, and proteins. The use of a fully discrete representation of solvent avoids the possible problems of continuum models to represent the solvation of systems containing tightly bound water molecules.     

An Effective Solvation Term Based on Atomic Occupancies for Use in Protein Simulations

Abstract

Several approaches to the treatment of solvent effects based on continuum models are reviewed and a new method based on occupied atomic volumes (occupancies) is proposed and tested. The new method describes protein-water interactions in terms of atomic solvation parameters, which represent the solvation free energy per unit of volume. These parameters were determined for six different atoms types, using experimental free energies of solvation. The method was implemented in the GROMOS and PRESTO molecular simulation program suites. Simulations with the solvation term require 20-50% more CPU time than the corresponding vacuum simulations and are approximately 20 times faster than explicit water simulations. The method and parameters were tested by carrying out 200 ps simulations of BPTI in water, in vacuo, and with the solvation term. The performance of the solvation term was assessed by comparing the structures and energies from the solvation simulations with the equivalent quantities derived from several BPTI crystal structures and from the explicit water and vacuum simulations. The model structures were evaluated in terms of exposed total surface, buried and exposed polar surfaces, secondary structure preservation, number of hydrogen bonds, energy contributions, and positional deviations from BPTI crystal structures. Vacuum simulations produced unrealistic structures with respect to all criteria applied. The structures resulting from the simulations with explicit water were closer to the 5PTI crystal structure, although part of the secondary structure dissolved. The simulations with the effective solvation term produce structures that are normal according to all evaluations and in most respects are remarkably similar to the 5PTI crystal structure despite considerable positional fluctuations during the simulations. The segments where the model and crystal structures differ are known to be flexible and the observed difference may be physically realistic. The effective solvation term based on occupancies is not only very efficient in terms of computer time but also results in meaningful structural properties for BPTI. It may therefore be generally useful in molecular dynamics of macromolecules.     

Methods for molecular dynamics simulations of protein folding/unfolding in solution

All atom molecular dynamics simulations have become a standard method for mapping equilibrium protein dynamics and non-equilibrium events like folding and unfolding. Here, we present detailed methods for performing such simulations. Generic protocols for minimization, solvation, simulation, and analysis derived from previous studies are also presented. As a measure of validation, our water model is compared with experiment. An example of current applications of these methods, simulations of the ultrafast folding protein Engrailed Homeodomain are presented including the experimental evidence used to verify their results. Ultrafast folders are an invaluable tool for studying protein behavior as folding and unfolding events measured by experiment occur on timescales accessible with the high-resolution molecular dynamics methods we describe. Finally, to demonstrate the prospect of these methods for folding proteins, a temperature quench simulation of a thermal unfolding intermediate of the Engrailed Homeodomain is described.     

Discrimination of the native from misfolded protein models with an energy function including implicit solvation.

An essential requirement for theoretical protein structure prediction is an energy function that can discriminate the native from non-native protein conformations. To date most of the energy functions used for this purpose have been extracted from a statistical analysis of the protein structure database, without explicit reference to the physical interactions responsible for protein stability. The use of the statistical functions has been supported by the widespread belief that they are superior for such discrimination to physics-based energy functions. An effective energy function which combined the CHARMM vacuum potential with a Gaussian model for the solvation free energy is tested for its ability to discriminate the native structure of a protein from misfolded conformations; the results are compared with those obtained with the vacuum CHARMM potential. The test is performed on several sets of misfolded structures prepared by others, including sets of about 650 good decoys for six proteins, as well as on misfolded structures of chymotrypsin inhibitor 2. The vacuum CHARMM potential is successful in most cases when energy minimized conformations are considered, but fails when applied to structures relaxed by molecular dynamics. With the effective energy function the native state is always more stable than grossly misfolded conformations both in energy minimized and molecular dynamics-relaxed structures. The present results suggest that molecular mechanics (physics-based) energy functions, complemented by a simple model for the solvation free energy, should be tested for use in the inverse folding problem, and supports their use in studies of the effective energy surface of proteins in solution. Moreover, the study suggests that the belief in the superiority of statistical functions for these purposes may be ill founded.     

Validating Solution Ensembles from Molecular Dynamics Simulation by Wide-Angle X-ray Scattering Data

Wide-angle x-ray scattering (WAXS) experiments of biomolecules in solution have become increasingly popular because of technical advances in light sources and detectors. However, the structural interpretation of WAXS profiles is problematic, partly because accurate calculations of WAXS profiles from structural models have remained challenging. In this work, we present the calculation of WAXS profiles from explicit-solvent molecular dynamics (MD) simulations of five different proteins. Using only a single fitting parameter that accounts for experimental uncertainties because of the buffer subtraction and dark currents, we find excellent agreement to experimental profiles both at small and wide angles. Because explicit solvation eliminates free parameters associated with the solvation layer or the excluded solvent, which would require fitting to experimental data, we minimize the risk of overfitting. We further find that the influence from water models and protein force fields on calculated profiles are insignificant up to q ≈ 15 nm^?1. Using a series of simulations that allow increasing flexibility of the proteins, we show that incorporating thermal fluctuations into the calculations significantly improves agreement with experimental data, demonstrating the importance of protein dynamics in the interpretation of WAXS profiles. In addition, free MD simulations up to one microsecond suggest that the calculated profiles are highly sensitive with respect to minor conformational rearrangements of proteins, such as an increased flexibility of a loop or an increase of the radius of gyration by  <  1%. The present study suggests that quantitative comparison between MD simulations and experimental WAXS profiles emerges as an accurate tool to validate solution ensembles of biomolecules.     

Molecular dynamics simulations of the unfolding of beta(2)-microglobulin and its variants

In this study, we examined the unfolding processes of native beta(2)-microglobulin and two related variants, one with an N-terminal hexapeptide deletion DeltaN6 and another with Lys57-Asp58 cleavage, by high-temperature molecular dynamics simulations. Three simulation models were used, molecular dynamics (MD) simulations with explicit water solvation, MD simulations with the CHARMM EEF1 force field and Langevin dynamics with the CHARMM EEF1 force field. Our simulations reproduce many of the experimentally observed structural changes. The most striking agreement is in the beta-strands to alpha-helix transition. In our simulations, strands beta(3), beta(4) and beta(5) consistently change to alpha-helix, whereas beta(8) changes to an alpha-helix only briefly. Through comparisons of the conformational behavior of the native, the DeltaN6 and the Lys-cut beta(2)-m, using the three simulation methods, we identified the consensus conformational changes that differentiate between the native beta(2)-m and its two variants. We found that the main effect of the removal of the N-terminal hexapeptide is to increase the separation between strands beta(2) and beta(6) and to facilitate the beta to alpha transition. On the other hand, the lysine cleavage only increases the flexibility of strand beta(5) and does not affect the interactions between strands beta(2) and beta(6). These conformational changes may relate to polymerization tendencies of these variants.     

Solid-State NMR-Restrained Ensemble Dynamics of a Membrane Protein in Explicit Membranes

Solid-state NMR has been used to determine the structures of membrane proteins in native-like lipid bilayer environments. Most structure calculations based on solid-state NMR observables are performed using simulated annealing with restrained molecular dynamics and an energy function, where all nonbonded interactions are represented by a single, purely repulsive term with no contributions from van der Waals attractive, electrostatic, or solvation energy. To our knowledge, this is the first application of an ensemble dynamics technique performed in explicit membranes that uses experimental solid-state NMR observables to obtain the refined structure of a membrane protein together with information about its dynamics and its interactions with lipids. Using the membrane-bound form of the fd coat protein as a model membrane protein and its experimental solid-state NMR data, we performed restrained ensemble dynamics simulations with different ensemble sizes in explicit membranes. For comparison, a molecular dynamics simulation of fd coat protein was also performed without any restraints. The average orientation of each protein helix is similar to a structure determined by traditional single-conformer approaches. However, their variations are limited in the resulting ensemble of structures with one or two replicas, as they are under the strong influence of solid-state NMR restraints. Although highly consistent with all solid-state NMR observables, the ensembles of more than two replicas show larger orientational variations similar to those observed in the molecular dynamics simulation without restraints. In particular, in these explicit membrane simulations, Lys⁴⁰, residing at the C-terminal side of the transmembrane helix, is observed to cause local membrane curvature. Therefore, compared to traditional single-conformer approaches in implicit environments, solid-state NMR restrained ensemble simulations in explicit membranes readily characterize not only protein dynamics but also protein-lipid interactions in detail.     

Free-energy calculations highlight differences in accuracy between X-ray and NMR structures and add value to protein structure prediction

BACKGROUND: While X-ray crystallography structures of proteins are considerably more reliable than those from NMR spectroscopy, it has been difficult to assess the inherent accuracy of NMR structures, particularly the side chains. RESULTS: For 15 small single-domain proteins, we used a molecular mechanics-/dynamics-based free-energy approach to investigate native, decoy, and fully extended alpha conformations. Decoys were all less energetically favorable than native conformations in nine of the ten X-ray structures and in none of the five NMR structures, but short 150 ps molecular dynamics simulations on the experimental structures caused them to have the lowest predicted free energy in all 15 proteins. In addition, a strong correlation exists (r(2) = 0.86) between the predicted free energy of unfolding, from native to fully extended conformations, and the number of residues. CONCLUSIONS: This work suggests that the approximate treatment of solvent used in solving NMR structures can lead NMR model conformations to be less reliable than crystal structures. This conclusion was reached because of the considerably higher calculated free energies and the extent of structural deviation during aqueous dynamics simulations of NMR models compared to those determined by X-ray crystallography. Also, the strong correlation found between protein length and predicted free energy of unfolding in this work suggests, for the first time, that a free-energy function can allow for identification of the native state based on calculations on an extended state and in the absence of an experimental structure.     

Empirical solvation models in the context of conformational energy searches: application to bovine pancreatic trypsin inhibitor.

Continuum solvation models that estimate free energies of solvation as a function of solvent accessible surface area are computationally simple enough to be useful for predicting protein conformation. The behavior of three such solvation models has been examined by applying them to the minimization of the conformational energy of bovine pancreatic trypsin inhibitor. The models differ only with regard to how the constants of proportionality between free energy and surface area were derived. Each model was derived by fitting to experimentally measured equilibrium solution properties. For two models, the solution property was free energy of hydration. For the third, the property was NMR coupling constants. The purpose of this study is to determine the effect of applying these solvation models to the nonequilibrium conformations of a protein arising in the course of global searches for conformational energy minima. Two approaches were used: (1) local energy minimization of an ensemble of conformations similar to the equilibrium conformation and (2) global search trajectories using Monte Carlo plus minimization starting from a single conformation similar to the equilibrium conformation. For the two models derived from free energy measurements, it was found that both the global searches and local minimizations yielded conformations more similar to the X-ray crystallographic structures than did searches or local minimizations carried out in the absence of a solvation component of the conformational energy. The model derived from NMR coupling constants behaved similarly to the other models in the context of a global search trajectory. For one of the models derived from measured free energies of hydration, it was found that minimization of an ensemble of near-equilibrium conformations yielded a new ensemble in which the conformation most similar to the X-ray determined structure PTI4 had the lowest total free energy. Despite the simplicity of the continuum solvation models, the final conformation generated in the trajectories for each of the models exhibited some of the characteristics that have been reported for conformations obtained from molecular dynamics simulations in the presence of a bath of explicit water molecules. They have smaller root mean square (rms) deviations from the experimentally determined conformation, fewer incorrect hydrogen bonds, and slightly larger radii of gyration than do conformations derived from search trajectories carried out in the absence of solvent.     

A critical analysis of continuum electrostatics: the screened Coulomb potential--implicit solvent model and the study of the alanine dipeptide and discrimination of misfolded structures of proteins

An analysis of the screened Coulomb potential--implicit solvent model (SCP--ISM) is presented showing that general equations for both the electrostatic and solvation free energy can be derived in a continuum approach, using statistical averaging of the polarization field created by the solvent around the molecule. The derivation clearly shows how the concept of boundary, usually found in macroscopic approaches, is eliminated when the continuum model is obtained from a microscopic treatment using appropriate averaging techniques. The model is used to study the alanine dipeptide in aqueous solution, as well as the discrimination of native protein structures from misfolded conformations. For the alanine dipeptide the free energy surface in the phi--psi space is calculated and compared with recently reported results of a detailed molecular dynamics simulation using an explicit representation of the solvent, and with other available data. The study showed that the results obtained using the SCP--ISM are comparable to those of the explicit water calculation and compares favorably to the FDPB approach. Both transition states and energy minima show a high correlation (r > 0.98) with the results obtained in the explicit water analysis. The study of the misfolded structures of proteins comprised the analysis of three standard decoy sets, namely, the EMBL, Park and Levitt, and Baker's CASP3 sets. In all cases the SCP--ISM discriminated well the native structures of the proteins, and the best-predicted structures were always near-native (cRMSD approximately 2 A).     

Optimized atomic radii for protein continuum electrostatics solvation forces

Recently, we presented a Green's function approach for the calculation of analytic continuum electrostatic solvation forces based on numerical solutions of the finite-difference Poisson-Botzmann (FDPB) equation [Im et al., Comp. Phys. Comm. 111 (1998) 59]. In this treatment the analytic forces were explicitly defined as the first derivative of the FDPB continuum electrostatic free energy with respect to the coordinates of the solute atoms. A smooth intermediate region for the solute-solvent dielectric boundary needed to be introduced to avoid abrupt discontinuous variations in the solvation free energy and forces as a function of the atomic positions. In the present paper we extend the set of optimized radii, which was previously parametrized from molecular dynamics free energy simulations of the 20 standard amino acids with explicit solvent molecules [Nina et al., J. Phys. Chem. 101 (1997) 5239], to yield accurate solvation free energy by taking the influence of the smoothed dielectric region into account.     

A Practical Implicit Membrane Potential for NMR Structure Calculations of Membrane Proteins

The highly anisotropic environment of the lipid bilayer membrane imposes significant constraints on the structures and functions of membrane proteins. However, NMR structure calculations typically use a simple repulsive potential that neglects the effects of solvation and electrostatics, because explicit atomic representation of the solvent and lipid molecules is computationally expensive and impractical for routine NMR-restrained calculations that start from completely extended polypeptide templates. Here, we describe the extension of a previously described implicit solvation potential, eefxPot, to include a membrane model for NMR-restrained calculations of membrane protein structures in XPLOR-NIH. The key components of eefxPot are an energy term for solvation free energy that works together with other nonbonded energy functions, a dedicated force field for conformational and nonbonded protein interaction parameters, and a membrane function that modulates the solvation free energy and dielectric screening as a function of the atomic distance from the membrane center, relative to the membrane thickness. Initial results obtained for membrane proteins with structures determined experimentally in lipid bilayer membranes show that eefxPot affords significant improvements in structural quality, accuracy, and precision. Calculations with eefxPot are straightforward to implement and can be used to both fold and refine structures, as well as to run unrestrained molecular-dynamics simulations. The potential is entirely compatible with the full range of experimental restraints measured by various techniques. Overall, it provides a useful and practical way to calculate membrane protein structures in a physically realistic environment.     

Correlation between knowledge-based and detailed atomic potentials: application to the unfolding of the GCN4 leucine zipper.

The relationship between the unfolding pseudo free energies of reduced and detailed atomic models of the GCN4 leucine zipper is examined. Starting from the native crystal structure, a large number of conformations ranging from folded to unfolded were generated by all-atom molecular dynamics unfolding simulations in an aqueous environment at elevated temperatures. For the detailed atomic model, the pseudo free energies are obtained by combining the CHARMM all-atom potential with a solvation component from the generalized Born, surface accessibility, GB/SA, model. Reduced model energies were evaluated using a knowledge-based potential. Both energies are highly correlated. In addition, both show a good correlation with the root mean square deviation, RMSD, of the backbone from native. These results suggest that knowledge-based potentials are capable of describing at least some of the properties of the folded as well as the unfolded states of proteins, even though they are derived from a database of native protein structures. Since only conformations generated from an unfolding simulation are used, we cannot assess whether these potentials can discriminate the native conformation from the manifold of alternative, low-energy misfolded states. Nevertheless, these results also have significant implications for the development of a methodology for multiscale modeling of proteins that combines reduced and detailed atomic models.     

A Gaussian statistical mechanical model for the equilibrium thermodynamics of barnase folding

Given an all non-hydrogen-atom potential function that implicitly includes solvation effects, it is possible to adjust its parameters to favor the correct native structure for several proteins over decoys produced by ungapped threading. It is also possible to further train it to reproduce the experimental free energy of unfolding in aqueous solution at 298 K for wild-type barnase and 66 mutants. For this, the native state is represented by the crystal structure at a single energy level with a calculated low degeneracy; the denatured state is represented by the extended conformation and a high calculated degeneracy. The same two-state model can be extended to account for the stability of all 67 sequences toward urea denaturation at 298 K by building in a solvation term that depends on urea concentration. With the addition of one more parameter set to give the correct heat capacity of unfolded barnase in solution, it is possible to approximate the experimental thermodynamics of barnase thermal denaturation: melting temperature, width of thermal transition, deltaG, deltaH, deltaS, and deltaCp. This requires a novel sort of statistical mechanical model where the two states each have a Gaussian density of microscopic state distribution as a function of energy.     

Solvation Free Energy of the Peptide Group: Its Model Dependence and Implications for the Additive-Transfer Free-Energy Model of Protein Stability

The group-additive decomposition of the unfolding free energy of a protein in an osmolyte solution relative to that in water poses a fundamental paradox: whereas the decomposition describes the experimental results rather well, theory suggests that a group-additive decomposition of free energies is, in general, not valid. In a step toward resolving this paradox, here we study the peptide-group transfer free energy. We calculate the vacuum-to-solvent (solvation) free energies of (Gly)_n and cyclic diglycine (cGG) and analyze the data according to experimental protocol. The solvation free energies of (Gly)_n are linear in n, suggesting group additivity. However, the slope interpreted as the free energy of a peptide unit differs from that for cGG scaled by a factor of half, emphasizing the context dependence of solvation. However, the water-to-osmolyte transfer free energies of the peptide unit are relatively independent of the peptide model, as observed experimentally. To understand these observations, a way to assess the contribution to the solvation free energy of solvent-mediated correlation between distinct groups is developed. We show that linearity of solvation free energy with n is a consequence of uniformity of the correlation contributions, with apparent group-additive behavior in the water-to-osmolyte transfer arising due to their cancellation. Implications for inferring molecular mechanisms of solvent effects on protein stability on the basis of the group-additive transfer model are suggested.     

Solvation Free Energy of the Peptide Group: Its Model Dependence and Implications for the Additive-Transfer Free-Energy Model of Protein Stability

The group-additive decomposition of the unfolding free energy of a protein in an osmolyte solution relative to that in water poses a fundamental paradox: whereas the decomposition describes the experimental results rather well, theory suggests that a group-additive decomposition of free energies is, in general, not valid. In a step toward resolving this paradox, here we study the peptide-group transfer free energy. We calculate the vacuum-to-solvent (solvation) free energies of (Gly)_n and cyclic diglycine (cGG) and analyze the data according to experimental protocol. The solvation free energies of (Gly)_n are linear in n, suggesting group additivity. However, the slope interpreted as the free energy of a peptide unit differs from that for cGG scaled by a factor of half, emphasizing the context dependence of solvation. However, the water-to-osmolyte transfer free energies of the peptide unit are relatively independent of the peptide model, as observed experimentally. To understand these observations, a way to assess the contribution to the solvation free energy of solvent-mediated correlation between distinct groups is developed. We show that linearity of solvation free energy with n is a consequence of uniformity of the correlation contributions, with apparent group-additive behavior in the water-to-osmolyte transfer arising due to their cancellation. Implications for inferring molecular mechanisms of solvent effects on protein stability on the basis of the group-additive transfer model are suggested.     

Thermal unfolding simulations of a multimeric protein--transition state and unfolding pathways

The folding of an oligomeric protein poses an extra challenge to the folding problem because the protein not only has to fold correctly; it has to avoid nonproductive aggregation. We have carried out over 100 molecular dynamics simulations using an implicit solvation model at different temperatures to study the unfolding of one of the smallest known tetramers, p53 tetramerization domain (p53tet). We found that unfolding started with disruption of the native tetrameric hydrophobic core. The transition state for the tetramer to dimer transition was characterized as a diverse ensemble of different structures using Phi value analysis in quantitative agreement with experimental data. Despite the diversity, the ensemble was still native-like with common features such as partially exposed tetramer hydrophobic core and shifts in the dimer-dimer arrangements. After passing the transition state, the secondary and tertiary structures continued to unfold until the primary dimers broke free. The free dimer had little secondary structure left and the final free monomers were random-coil like. Both the transition states and the unfolding pathways from these trajectories were very diverse, in agreement with the new view of protein folding. The multiple simulations showed that the folding of p53tet is a mixture of the framework and nucleation-condensation mechanisms and the folding is coupled to the complex formation. We have also calculated the entropy and effective energy for the different states along the unfolding pathway and found that the tetramerization is stabilized by hydrophobic interactions.