Temperature and pressure denaturation of chignolin: Folding and unfolding simulation by multibaric-multithermal molecular dynamics method

A multibaric‐multithermal molecular dynamics (MD) simulation of a 10‐residue protein, chignolin, was performed. All‐atom model with the Amber parm99SB force field was used for the protein and the TIP3P model was used for the explicit water molecules. This MD simulation covered wide ranges of temperature between 260 and 560 K and pressure between 0.1 and 600 MPa and sampled many conformations without getting trapped in local‐minimum free‐energy states. Folding events to the native β‐hairpin structure occurred five times and unfolding events were observed four times. As the temperature and/or pressure increases, fraction of folded chignolin decreases. The partial molar enthalpy change ΔH and partial molar volume change ΔV of unfolding were calculated as ΔH = 24.1 ± 4.9 kJ/mol and ΔV = ?5.6 ± 1.5 cm³/mol, respectively. These values agree well with recent experimental results. Illustrating typical local‐minimum free‐energy conformations, folding and unfolding pathways were revealed. When chignolin unfolds from the β‐hairpin structure, only the C terminus or both C and N termini open first. It may undergo an α‐helix or 3₁₀‐helix structure and finally unfolds to the extended structure. Difference of the mechanism between temperature denaturation and pressure denaturation is also discussed. Temperature denaturation is caused by making the protein transferred to a higher entropy state and making it move around more with larger space. The reason for pressure denaturation is that water molecules approach the hydrophobic residues, which are not well hydrated at the folded state, and some hydrophobic contacts are broken. Proteins 2012;. © 2012 Wiley Periodicals, Inc.     

Exploring the free energy landscape of a model β‐hairpin peptide and its isoform

Secondary structural transitions from α‐helix to β‐sheet conformations are observed in several misfolding diseases including Alzheimer's and Parkinson's. Determining factors contributing favorably to the formation of each of these secondary structures is therefore essential to better understand these disease states. β‐hairpin peptides form basic components of anti‐parallel β‐sheets and are suitable model systems for characterizing the fundamental forces stabilizing β‐sheets in fibrillar structures. In this study, we explore the free energy landscape of the model β‐hairpin peptide GB1 and its E2 isoform that preferentially adopts α‐helical conformations at ambient conditions. Umbrella sampling simulations using all‐atom models and explicit solvent are performed over a large range of end‐to‐end distances. Our results show the strong preference of GB1 and the E2 isoform for β‐hairpin and α‐helical conformations, respectively, consistent with previous studies. We show that the unfolded states of GB1 are largely populated by misfolded β‐hairpin structures which differ from each other in the position of the β‐turn. We discuss the energetic factors contributing favorably to the formation of α‐helix and β‐hairpin conformations in these peptides and highlight the energetic role of hydrogen bonds and non‐bonded interactions. Proteins 2014; 82:2394–2402. © 2014 Wiley Periodicals, Inc.     

Conformational transitions of a dipeptide in water: Effects of imposed pathways using umbrella sampling techniques

The free energy difference between two states of a molecular system separated by an energy barrier can generally be computed using the technique of umbrella sampling along a chosen reaction coordinate or pathway. The effect of a particular choice of pathway upon the obtained free energy difference is investigated by molecular dynamics simulation of a model system consisting of a glycine dipeptide in aqueous solution. Two different reaction coordinates connecting the so-called C₅ and C₇ conformations, one involving intramolecular hydrogen bonds and the other involving the peptide ?, ψ angles, are considered. The Gibbs free energy differences ΔG(C₅ – C₇) are small in both cases, 1.5 ± 1 kJ mol^?1 and 2.2 ± 1 kJ mol ^?1, respectively. The two different reaction coordinates yield free energy differences that are identical to within their statistical error. It is found that the exchange of solute–solute, solute–water, and water–water hydrogen bonds involves free energy changes of less than k_BT, which points at the existence of a multitutde of low free energy pathways connecting the C₅ and C₇ dipeptide conformations. © 1994 John Wiley & Sons, Inc.     

Josephin Domain Structural Conformations Explored by Metadynamics in Essential Coordinates

The Josephin Domain (JD), i.e. the N-terminal domain of Ataxin 3 (At3) protein, is an interesting example of competition between physiological function and aggregation risk. In fact, the fibrillogenesis of Ataxin 3, responsible for the spinocerebbellar ataxia 3, is strictly related to the JD thermodynamic stability. Whereas recent NMR studies have demonstrated that different JD conformations exist, the likelihood of JD achievable conformational states in solution is still an open issue. Marked differences in the available NMR models are located in the hairpin region, supporting the idea that JD has a flexible hairpin in dynamic equilibrium between open and closed states. In this work we have carried out an investigation on the JD conformational arrangement by means of both classical molecular dynamics (MD) and Metadynamics employing essential coordinates as collective variables. We provide a representation of the free energy landscape characterizing the transition pathway from a JD open-like structure to a closed-like conformation. Findings of our in silico study strongly point to the closed-like conformation as the most likely for a Josephin Domain in water.     

(Un)Folding Mechanisms of the FBP28 WW Domain in Explicit Solvent Revealed by Multiple Rare Event Simulation Methods

We report a numerical study of the (un)folding routes of the truncated FBP28 WW domain at ambient conditions using a combination of four advanced rare event molecular simulation techniques. We explore the free energy landscape of the native state, the unfolded state, and possible intermediates, with replica exchange molecular dynamics. Subsequent application of bias-exchange metadynamics yields three tentative unfolding pathways at room temperature. Using these paths to initiate a transition path sampling simulation reveals the existence of two major folding routes, differing in the formation order of the two main hairpins, and in hydrophobic side-chain interactions. Having established that the hairpin strand separation distances can act as reasonable reaction coordinates, we employ metadynamics to compute the unfolding barriers and find that the barrier with the lowest free energy corresponds with the most likely pathway found by transition path sampling. The unfolding barrier at 300 K is ∼17 k_BT ≈ 42 kJ/mol, in agreement with the experimental unfolding rate constant. This work shows that combining several powerful simulation techniques provides a more complete understanding of the kinetic mechanism of protein folding.     

Multiscale simulations of protein folding: application to formation of secondary structures

A multiscale simulation method of protein folding is proposed, using atomic representation of protein and solvent, combing genetic algorithms to determine the key protein structures from a global view, with molecular dynamic simulations to reveal the local folding pathways, thus providing an integrated landscape of protein folding. The method is found to be superior to previously investigated global search algorithms or dynamic simulations alone. For secondary structure formation of a selected peptide, RN24, the structures and dynamics produced by this method agree well with corresponding experimental results. Three most populated conformations are observed, including hairpin, β-sheet and α-helix. The energetic barriers separating these three structures are comparable to the kinetic energy of the atoms of the peptide, implying that the transition between these states can be easily triggered by kinetic perturbations, mainly through electrostatic interactions between charged atoms. Transitions between α-helix and β-sheet should jump over at least two energy barriers and may stay in the energetic trap of hairpin. It is proposed that the structure of proteins should be jointly governed by thermodynamic and dynamic factors; free energy is not the exclusive dominant for stability of proteins.     

A shorter peptide model from staphylococcal nuclease for the folding-unfolding equilibrium of a beta-hairpin shows that unfolded state has significant contribution from compact conformational states

It is important to understand the conformational features of the unfolded state in equilibrium with folded state under physiological conditions. In this paper, we consider a short peptide model LMYKGQPM from staphylococcal nuclease to model the conformational equilibrium between a hairpin conformation and its unfolded state using molecular dynamics simulation under NVT conditions at 300K using GROMOS96 force field. The free energy landscape has overall funnel-like shape with hairpin conformations sampling the minima. The "unfolded" state has a higher free energy of approximately 12kJ/mol with respect to native hairpin minimum and occupies a plateau region. We find that the unfolded state has significant contributions from compact conformations. Many of these conformations have hairpin-like topology. Further, these compact conformational forms are stabilized by hydrophobic interactions. Conversion between native and non-native hairpins occurs via unfolded states. Frequent conversions between folded and unfolded hairpins are observed with single exponential kinetics. We compare our results with the emerging picture of unfolded state from both experimental and theoretical studies.     

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.     

Nonlinear machine learning in simulations of soft and biological materials

Abstract

Interpretable parameterisations of free energy landscapes for soft and biological materials calculated from molecular simulation require the availability of ‘good’ collective variables (CVs) capable of discriminating the metastable states of the system and the barriers between them. If these CVs are coincident with the slow collective modes governing the long-time dynamical evolution, then they also furnish good coordinates in which to perform enhanced sampling to surmount high free energy barriers and efficiently explore and recover the landscape. Non-linear manifold learning techniques provide a means to systematically extract such CVs from molecular simulation trajectories by identifying and extracting low-dimensional manifolds lying latent within the high-dimensional coordinate space. We survey recent advances in data-driven CV discovery and enhanced sampling using non-linear manifold learning, describe the mathematical and theoretical underpinnings of these techniques, and present illustrative examples to molecular folding and colloidal self-assembly. We close with our outlook and perspective on future advances in this rapidly evolving field.     

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.     

Stability of polypeptide conformational states as determined by computer simulation of the free energy

A method is developed to extract the entropy of polypeptides and proteins from samples of conformations. It is based on techniques suggested previously by Meirovitch, and has the advantage that it can be applied not only to states in which the molecule undergoes harmonic or quasiharmonic conformational fluctuations, but also to the random coil, as well as to mixtures of these extreme states. In order to confine the search to a region of conformational space corresponding to a stable state, the transition probabilities are determined not by “looking to the future,” as in the previous method [H. Meirovitch and H. A. Scheraga (1986) J. Chem. Phys. 84 , 6369–6375], but by analyzing the previous steps in the generation of the chain. The method is applied to a model of decaglycine with rigid geometry, using the potential energy function ECEPP (Empirical Conformational Energy Program for Peptides). The model is simulated with the Metropolis Monte Carlo method to generate samples of conformations in the α-helical and hairpin regions, respectively, at T = 100 K. For the α-helix, the four dihedral angles of the N- and C-terminal residues are found to undergo full rotational variation. The results show that the α-helix is a more stable structure than the hairpin. Both its Helmholtz free energy F and energy E are lower than those of the hairpin by ΔF ～ 0.4 and ΔE ～ 0.3 kcal/mole/residue, respectively. It should be noted that the contribution of the entropy ΔS to ΔF is significant (TΔS ～ 0.1 kcal/mole/residue). Also, the entropy of the α-helix is found to be larger than that of the hairpin. This is a result of the extra entropy arising from the rotational freedom about the four terminal single bonds of the α-helix.     

Folding thermodynamics of β‐hairpins studied by replica‐exchange molecular dynamics simulations

We study the differences in folding stability of β‐hairpin peptides, including GB1 hairpin and a point mutant GB1 K10G, as well as tryptophan zippers (TrpZips): TrpZip1, TrpZip2, TrpZip3‐1, and TrpZip4. By performing replica‐exchange molecular dynamics simulations with Amber03* force field (a modified version of Amber ff03) in explicit solvent, we observe ab initio folding of all the peptides except TrpZip3‐1, which is experimentally known to be the least stable among the peptides studied here. By calculating the free energies of unfolding of the peptides at room temperature and folding midpoint temperatures for thermal unfolding of peptides, we find that TrpZip4 and GB1 K10G peptides are the most stable β‐hairpins followed by TrpZip1, GB1, and TrpZip2 in the given order. Hence, the proposed K10G mutation of GB1 peptide results in enhanced stability compared to wild‐type GB1. An important goal of our study is to test whether simulations with Amber 03* model can reproduce experimentally predicted folding stability differences between these peptides. While the stabilities of GB1 and TrpZip1 yield close agreement with experiment, TrpZip2 is found to be less stable than predicted by experiment. However, as heterogenous folding of TrpZip2 may yield divergent thermodynamic parameters by different spectroscopic methods, mismatching of results with previous experimental values are not conclusive of model shortcomings. For most of the cases, molecular simulations with Amber03* can successfully reproduce experimentally known differences between the mutated peptides, further highlighting the predictive capabilities of current state‐of‐the‐art all‐atom protein force fields. Proteins 2015; 83:1307–1315. © 2015 Wiley Periodicals, Inc.     

Energy landscape of a small peptide revealed by dihedral angle principal component analysis

A 100 ns molecular dynamics simulation of penta-alanine in explicit water is performed to study the reversible folding and unfolding of the peptide. Employing a standard principal component analysis (PCA) using Cartesian coordinates, the resulting free-energy landscape is found to have a single minimum, thus suggesting a simple, relatively smooth free-energy landscape. Introducing a novel PCA based on a transformation of the peptide dihedral angles, it is found, however, that there are numerous free energy minima of comparable energy (less than or approximately 1 kcal/mol), which correspond to well-defined structures with characteristic hydrogen-bonding patterns. That is, the true free-energy landscape is actually quite rugged and its smooth appearance in the Cartesian PCA represents an artifact of the mixing of internal and overall motion. Well-separated minima corresponding to specific conformational structures are also found in the unfolded part of the free energy landscape, revealing that the unfolded state of penta-alanine is structured rather than random. Performing a connectivity analysis, it is shown that neighboring states are connected by low barriers of similar height and that each state typically makes transitions to three or four neighbor states. Several principal pathways for helix nucleation are identified and discussed in some detail.     

Determination of low-energy structures of a small RNA hairpin using Monte Carlo–based techniques

The energy landscape of RNA is known to be extremely rugged, and hence finding low-energy structures starting from a random structure is a challenging task for any optimization algorithm. In the current work, we have investigated the ability of one Monte Carlo-based optimization algorithm, Temperature Basin Paving, to explore the energy landscape of a small RNA T-loop hairpin. In this method, the history of the simulation is used to increase the probability of states less visited in the simulation. It has been found that using both energy and end-to-end distance as the biasing parameters in the simulation, the partially folded structure of the hairpin starting from random structures could be obtained.     

Energetics of the hairpin to mismatched duplex transition of d(GCCGCAGC) on NaCl solution

We report Potential of Mean Force studies to describe the relative thermodynamic stabilities of d(GCCGCAGC) in a mismatched duplex and a hairpin monomer conformation in NaCl solution. The PMF calculations are combined with previous molecular mechanics and normal mode analysis in order to estimate the role of different components of the free energy in determining the relative stability of the duplex and hairpin structures. The high entropy associated with the loop region and the lack of minor groove phosphate-phosphate interactions in the hairpin compete against the gain in enthalpic contribution to the free energy due to base pairing in the mismatched duplex. The combined free energy calculations show that the hairpin is the most stable conformation at low salt and that a hairpin to duplex transition takes place at approximately 0.47 M NaCl. In addition, we studied the hairpin to partially stacked single helical conformation equilibrium at low salt. We found a small variation in transition temperature in salt concentration, delta Tm/delta log10(cs) approximately 2-3 degrees K/decade, in contrast to the duplex to hairpin or duplex to partially stacked single helix transition where the transition temperature exhibited marked dependence on salt concentration. This is in qualitative agreement with experimental data. Based on the Potential of Mean Force free energy calculation, the order of relative stability of the three-conformations studied varies with salt concentration. We observed the following orders of stability: stacked single helix greater than hairpin greater than duplex for cs less than 0.77 M NaCl; single helix greater than duplex greater than hairpin for 0.77 less than Cs less than 2.1 M; and duplex greater than hairpin greater than single strand for cs greater than 2.1 M. From the calculated PMF free energy curves in the NaCl concentration range, 0.012 less than cs less than 5.0 M, we can assign upper and lower bounds for the non-ionic differences in free energy between the duplex, hairpin, and stacked single helical states (at standard conditions: cs = 1.0 M, T = 25 degrees C, and 1 M oligomer concentration). We found that for delta G duplex single helix = G duplex - 2 x G single helix less than -7.38 Kcal/mol, the single helix is the least stable state. For the duplex-to-hairpin free energy difference in the range, -1.87 less than delta G duplex-hairpin less than 0.03 Kcal/mol, there will always be a salt-induced hairpin-to-duplex transition for 0.01 less than cs less than 1.6 M NaCl. If delta G duplex-hairpin less than -1.87, the duplex is always more stable than the hairpin; and for delta G duplex-hairpin greater than Kcal/mol, the hairpin state is always more stable than the duplex, for all salt concentrations.     

A molecular dynamics study of the effect of glycosidic linkage type in the hemicellulose backbone on the molecular chain flexibility

The macromolecular conformation of the constituent polysaccharides in lignocellulosic biomass influences their supramolecular interactions, and therefore their function in plants and their performance in technical products. The flexibility of glycosidic linkages from the backbone of hemicelluloses was studied by evaluating the conformational freedom of the φ and ψ dihedral angles using molecular dynamic simulations, additionally selected molecules were correlated with experimental data by nuclear magnetic resonance spectroscopy. Three types of β‐(1→4) glycosidic linkages involving the monosaccharides (Glcp, Xylp and Manp) present in the backbone of hemicelluloses were defined. Different di‐ and tetrasaccharides with combinations of such sugar monomers from hemicelluloses were simulated, and free energy maps of the φ – ψ space and hydrogen‐bonding patterns were obtained. The glycosidic linkage between Glc‐Glc or Glc‐Man (C‐type) was the stiffest with mainly one probable conformation; the linkage from Man‐Man or Man‐Glc (M‐type) was similar but with an increased probability for an alternative conformation making it more flexible, and the linkage between two Xyl‐units (X‐type) was the most flexible with two almost equally populated conformations. Glycosidic linkages of the same type showed essentially the same conformational space in both disaccharides and in the central region of tetrasaccharides. Different probabilities of glycosidic linkage conformations in the backbone of hemicelluloses can be directly estimated from the free energy maps, which to a large degree affect the overall macromolecular conformations of these polymers. The information gained contributes to an increased understanding of the function of hemicelluloses both in the cell wall and in technical products.     

Molecular simulation uncovers the conformational space of the λ Cro dimer in solution

The significant variation among solved structures of the λ Cro dimer suggests its flexibility. However, contacts in the crystal lattice could have stabilized a conformation which is unrepresentative of its dominant solution form. Here we report on the conformational space of the Cro dimer in solution using replica exchange molecular dynamics in explicit solvent. The simulated ensemble shows remarkable correlation with available x-ray structures. Network analysis and a free energy surface reveal the predominance of closed and semi-open dimers, with a modest barrier separating these two states. The fully open conformation lies higher in free energy, indicating that it requires stabilization by DNA or crystal contacts. Most NMR models are found to be unstable conformations in solution. Intersubunit salt bridging between Arg⁴ and Glu⁵³ during simulation stabilizes closed conformations. Because a semi-open state is among the low-energy conformations sampled in simulation, we propose that Cro-DNA binding may not entail a large conformational change relative to the dominant dimer forms in solution.     

Combinatorial pattern discovery approach for the folding trajectory analysis of a beta-hairpin

The study of protein folding mechanisms continues to be one of the most challenging problems in computational biology. Currently, the protein folding mechanism is often characterized by calculating the free energy landscape versus various reaction coordinates, such as the fraction of native contacts, the radius of gyration, RMSD from the native structure, and so on. In this paper, we present a combinatorial pattern discovery approach toward understanding the global state changes during the folding process. This is a first step toward an unsupervised (and perhaps eventually automated) approach toward identification of global states. The approach is based on computing biclusters (or patterned clusters)-each cluster is a combination of various reaction coordinates, and its signature pattern facilitates the computation of the Z-score for the cluster. For this discovery process, we present an algorithm of time complexity c in RO((N + nm) log n), where N is the size of the output patterns and (n x m) is the size of the input with n time frames and m reaction coordinates. To date, this is the best time complexity for this problem. We next apply this to a beta-hairpin folding trajectory and demonstrate that this approach extracts crucial information about protein folding intermediate states and mechanism. We make three observations about the approach: (1) The method recovers states previously obtained by visually analyzing free energy surfaces. (2) It also succeeds in extracting meaningful patterns and structures that had been overlooked in previous works, which provides a better understanding of the folding mechanism of the beta-hairpin. These new patterns also interconnect various states in existing free energy surfaces versus different reaction coordinates. (3) The approach does not require calculating the free energy values, yet it offers an analysis comparable to, and sometimes better than, the methods that use free energy landscapes, thus validating the choice of reaction coordinates. (An abstract version of this work was presented at the 2005 Asia Pacific Bioinformatics Conference [1].).