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.     

Simulation of folding of a small alpha-helical protein in atomistic detail using worldwide-distributed computing

By employing thousands of PCs and new worldwide-distributed computing techniques, we have simulated in atomistic detail the folding of a fast-folding 36-residue alpha-helical protein from the villin headpiece. The total simulated time exceeds 300 micros, orders of magnitude more than previous simulations of a molecule of this size. Starting from an extended state, we obtained an ensemble of folded structures, which is on average 1.7A and 1.9A away from the native state in C(alpha) distance-based root-mean-square deviation (dRMS) and C(beta) dRMS sense, respectively. The folding mechanism of villin is most consistent with the hydrophobic collapse view of folding: the molecule collapses non-specifically very quickly ( approximately 20ns), which greatly reduces the size of the conformational space that needs to be explored in search of the native state. The conformational search in the collapsed state appears to be rate-limited by the formation of the aromatic core: in a significant fraction of our simulations, the C-terminal phenylalanine residue packs improperly with the rest of the hydrophobic core. We suggest that the breaking of this interaction may be the rate-determining step in the course of folding. On the basis of our simulations we estimate the folding rate of villin to be approximately 5micros. By analyzing the average features of the folded ensemble obtained by simulation, we see that the mean folded structure is more similar to the native fold than any individual folded structure. This finding highlights the need for simulating ensembles of molecules and averaging the results in an experiment-like fashion if meaningful comparison between simulation and experiment is to be attempted. Moreover, our results demonstrate that (1) the computational methodology exists to simulate the multi-microsecond regime using distributed computing and (2) that potential sets used to describe interatomic interactions may be sufficiently accurate to reach the folded state, at least for small proteins. We conclude with a comparison between our results and current protein-folding theory.     

Principles of carbopeptoid folding: a molecular dynamics simulation study.

The conformational spaces of five oligomers of tetrahydrofuran-based carbopeptoids in chloroform and dimethyl sulfoxide were investigated through nine molecular dynamics simulations. Prompted by nuclear magnetic resonance experiments that indicated various stable folds for some but not all of these carbopeptoids, their folding behaviour was investigated as a function of stereochemistry, chain length and solvent. The conformational distributions of these molecules were analysed in terms of occurrence of hydrogen bonds, backbone torsional-angle distributions, conformational clustering and solute configurational entropy. While a cis-linkage across the tetrahydrofuran ring favours right-handed helical structures, a trans-linkage results in a larger conformational variability. Intra-solute hydrogen bonding is reduced with increasing chain length and with increasing solvent polarity. Solute configurational entropies confirm the picture obtained: they are smaller for cis- than for trans-linked peptides, for chloroform than for dimethyl sulfoxide as solvent and for shorter peptide chains. The simulations provide an atomic picture of molecular conformational variability that is consistent with the available experimental data.     

The folding mechanism of collagen-like model peptides explored through detailed molecular simulations

Collagen has a unique folding mechanism that begins with the formation of a triple-helical structure near its C terminus followed by propagation of this structure to the N terminus. To elucidate factors that affect the folding of collagen, we explored the folding pathway of collagen-like model peptides using detailed molecular simulations with explicit solvent. Using biased molecular dynamics we examined the latter stages of folding of a peptide model of native collagen, (Pro-Hyp-Gly)10, and a peptide that models a Gly --> Ser mutation found in several forms of osteogenesis imperfecta, (Pro-Hyp-Gly)3-Pro-Hyp-Ser-(Pro-Hyp-Gly)6. Starting from an unfolded state that contains a C-terminal nucleated trimer, (Pro-Hyp-Gly)10 folds to a structure where two of the three chains associate through water-mediated hydrogen bonds and the third is relatively separated from this dimer. Calculated free-energy profiles for folding from this intermediate to the final triple-helical structure suggest that further folding occurs at a rate of approximately one Pro-Hyp-Gly triplet per msec. In contrast, after 6 nsec of biased dynamics, the region N-terminal to the Ser residue in (Pro-Hyp-Gly)3-Pro-Hyp-Ser-(Pro-Hyp-Gly)6 folds to a structure where the three chains form close contacts near the N terminus, away from the mutation site. Further folding to an ideal triple-helical structure at the site of the mutation is unfavorable as the free energy of a triple-helical conformation at this position is more than 20 kcal/mol higher than that of a structure with unassociated chains. These data provide insights into the folding pathway of native collagen and the events underlying the formation of misfolded structures.     

Reproducible polypeptide folding and structure prediction using molecular dynamics simulations

The folding of a polypeptide from an extended state to a well-defined conformation is studied using microsecond classical molecular dynamics (MD) simulations and replica exchange molecular dynamics (REMD) simulations in explicit solvent and in vacuo. It is shown that the solvated peptide folds many times in the REMD simulations but only a few times in the conventional simulations. From the folding events in the classical simulations we estimate an approximate folding time of 1-2 micros. The REMD simulations allow enough sampling to deduce a detailed Gibbs free energy landscape in three dimensions. The global minimum of the energy landscape corresponds to the native state of the peptide as determined previously by nuclear magnetic resonance (NMR) experiments. Starting from an extended state it takes about 50 ns before the native structure appears in the REMD simulations, about an order of magnitude faster than conventional MD. The calculated melting curve is in good qualitative agreement with experiment. In vacuo, the peptide collapses rapidly to a conformation that is substantially different from the native state in solvent.     

Understanding the folding and stability of a designed WW domain protein with replica exchange molecular dynamics simulations

WW domain proteins are usually regarded as simple models for understanding the folding mechanism of β-sheet. CC45 is an artificial protein that is capable of folding into the same structure as WW domain. In this article, the replica exchange molecular dynamics simulations are performed to investigate the folding mechanism of CC45. The analysis of thermal stability shows that β-hairpin 1 is more stable than β-hairpin 2 during the unfolding process. Free energy analysis shows that the unfolding of this protein substantially proceeds through solvating the smaller β-hairpin 2, followed by the unfolding of β-hairpin 1. We further propose the unfolding process of CC45 and the folding mechanism of two β-hairpins. These results are similar to the previous folding studies of formin binding protein 28 (FBP28). Compared with FBP28, it is found that CC45 has more aromatic residues in N-terminal loop, and these residues contact with C-terminal loop to form the outer hydrophobic core, which increases the stability of CC45. Knowledge about the stability and folding behaviour of CC45 may help in understanding the folding mechanisms of the β-sheet and in designing new WW domains.     

Integrating experimental data with molecular simulations to investigate RNA structural dynamics

Conformational dynamics is crucial for ribonucleic acid (RNA) function. Techniques such as nuclear magnetic resonance, cryo-electron microscopy, small- and wide-angle X-ray scattering, chemical probing, single-molecule Förster resonance energy transfer, or even thermal or mechanical denaturation experiments probe RNA dynamics at different time and space resolutions. Their combination with accurate atomistic molecular dynamics (MD) simulations paves the way for quantitative and detailed studies of RNA dynamics. First, experiments provide a quantitative validation tool for MD simulations. Second, available data can be used to refine simulated structural ensembles to match experiments. Finally, comparison with experiments allows for improving MD force fields that are transferable to new systems for which data is not available. Here we review the recent literature and provide our perspective on this field.     

Characterizing a partially folded intermediate of the villin headpiece domain under non-denaturing conditions: contribution of His41 to the pH-dependent stability of the N-terminal subdomain

The contribution of interactions involving the imidazole ring of His41 to the pH-dependent stability of the villin headpiece (HP67) N-terminal subdomain has been investigated by nuclear magnetic resonance (NMR) spin relaxation. NMR-derived backbone N-H order parameters (S2) for wild-type (WT) HP67 and H41Y HP67 indicate that reduced conformational flexibility of the N-terminal subdomain in WT HP67 is due to intramolecular interactions with the His41 imidazole ring. These interactions, together with desolvation effects, contribute to significantly depress the pKa of the buried imidazole ring in the native state. 15N R1rho relaxation dispersion data indicate that WT HP67 populates a partially folded intermediate state that is 10.9 kJ mol(-1) higher in free energy than the native state under non-denaturing conditions at neutral pH. The partially folded intermediate is characterized as having an unfolded N-terminal subdomain while the C-terminal subdomain retains a native-like fold. Although the majority of the residues in the N-terminal subdomain sample a random-coil distribution of conformations, deviations of backbone amide 1H and 15N chemical shifts from canonical random-coil values for residues within 5A of the His41 imidazole ring indicate that a significant degree of residual structure is maintained in the partially folded ensemble. The pH-dependence of exchange broadening is consistent with a linear three-state exchange model whereby unfolding of the N-terminal subdomain is coupled to titration of His41 in the partially folded intermediate with a pKa,I=5.69+/-0.07. Although maintenance of residual interactions with the imidazole ring in the unfolded N-terminal subdomain appears to reduce pKa,I compared to model histidine compounds, protonation of His41 disrupts these interactions and reduces the difference in free energy between the native state and partially folded intermediate under acidic conditions. In addition, chemical shift changes for residues Lys70-Phe76 in the C-terminal subdomain suggest that the HP67 actin binding site is disrupted upon unfolding of the N-terminal subdomain, providing a potential mechanism for regulating the villin-dependent bundling of actin filaments.     

Kinetics of cytochrome C folding: atomically detailed simulations

The vast range of time scales (from nanoseconds to seconds) during protein folding is a challenge for experiments and computations. To make concrete predictions on folding mechanisms, atomically detailed simulations of protein folding, using potentials derived from chemical physics principles, are desired. However, due to their computational complexity, straightforward molecular dynamics simulations of protein folding are impossible today. An alternative algorithm is used that makes it possible to compute approximate atomically detailed long time trajectories (the Stochastic Difference Equation in Length). This algorithm is used to compute 26 atomically detailed folding trajectories of cytochrome c (a millisecond process). The early collapse of the protein chain (with marginal formation of secondary structure), and the earlier formation of the N and C helices (compare to the 60's helix) are consistent with the experiment. The existence of an energy barrier upon entry to the molten globule is examined as well. In addition to (favorable) comparison to experiments, we show that non-native contacts drive the formation of the molten globule. In contrast to popular folding models, the non-native contacts do not form off-pathway kinetic traps in cytochrome c.     

Ab initio RNA folding by discrete molecular dynamics: from structure prediction to folding mechanisms

RNA molecules with novel functions have revived interest in the accurate prediction of RNA three-dimensional (3D) structure and folding dynamics. However, existing methods are inefficient in automated 3D structure prediction. Here, we report a robust computational approach for rapid folding of RNA molecules. We develop a simplified RNA model for discrete molecular dynamics (DMD) simulations, incorporating base-pairing and base-stacking interactions. We demonstrate correct folding of 150 structurally diverse RNA sequences. The majority of DMD-predicted 3D structures have <4 A deviations from experimental structures. The secondary structures corresponding to the predicted 3D structures consist of 94% native base-pair interactions. Folding thermodynamics and kinetics of tRNA(Phe), pseudoknots, and mRNA fragments in DMD simulations are in agreement with previous experimental findings. Folding of RNA molecules features transient, non-native conformations, suggesting non-hierarchical RNA folding. Our method allows rapid conformational sampling of RNA folding, with computational time increasing linearly with RNA length. We envision this approach as a promising tool for RNA structural and functional analyses.     

Folding simulations of Trp‐cage mini protein in explicit solvent using biasing potential replica‐exchange molecular dynamics simulations

Replica exchange molecular dynamics (RexMD) simulations are frequently used for studying structure formation and dynamics of peptides and proteins. A significant drawback of standard temperature RexMD is, however, the rapid increase of the replica number with increasing system size to cover a desired temperature range. A recently developed Hamiltonian RexMD method has been used to study folding of the Trp‐cage protein. It employs a biasing potential that lowers the backbone dihedral barriers and promotes peptide backbone transitions along the replica coordinate. In two independent applications of the biasing potential RexMD method including explicit solvent and starting from a completely unfolded structure the formation of near‐native conformations was observed after 30–40 ns simulation time. The conformation representing the most populated cluster at the final simulation stage had a backbone root mean square deviation of ～1.3 Å from the experimental structure. This was achieved with a very modest number of five replicas making it well suited for peptide and protein folding and refinement studies including explicit solvent. In contrast, during five independent continuous 70 ns molecular dynamics simulations formation of collapsed states but no near native structure formation was observed. The simulations predict a largely collapsed state with a significant helical propensity for the helical domain of the Trp‐cage protein already in the unfolded state. Hydrogen bonded bridging water molecules were identified that could play an active role by stabilizing the arrangement of the helical domain with respect to the rest of the chain already in intermediate states of the protein. Proteins 2009. © 2008 Wiley‐Liss, Inc.     

Dimerization of the p53 oligomerization domain: identification of a folding nucleus by molecular dynamics simulations

Dimerization of the p53 oligomerization domain involves coupled folding and binding of monomers. To examine the dimerization, we have performed molecular dynamics (MD) simulations of dimer folding from the rate-limiting transition state ensemble (TSE). Among 799 putative transition state structures that were selected from a large ensemble of high-temperature unfolding trajectories, 129 were identified as members of the TSE via calculation of a 50% transmission coefficient from at least 20 room-temperature simulations. This study is the first to examine the refolding of a protein dimer using MD simulations in explicit water, revealing a folding nucleus for dimerization. Our atomistic simulations are consistent with experiment and offer insight that was previously unobtainable.     

Understanding the folding and stability of a zinc finger-based full sequence design protein with replica exchange molecular dynamics simulations

Full sequence design protein FSD-1 is a designed protein based on the motif of zinc finger protein. In this work, its folding mechanism and thermal stability are investigated using the replica exchange molecular dynamics model with the water molecules being treated explicitly. The results show that the folding of the FSD-1 is initiated by the hydrophobic collapse, which is accompanied with the formation of the C-terminal alpha-helix. Then the folding proceeds with the formation of the beta-hairpin and the further package of the hydrophobic core. Compared with the beta-hairpin, the alpha-helix has much higher stability. It is also found that the N-capping motif adopted by the FSD-1 contributes to the stability of the alpha-helix dramatically. The hydrophobic contacts made by the side chain of Tyr3 in the native state are essential for the stabilization of the beta-hairpin. It is also found that the folding of the N-terminal beta-hairpin and the C-terminal alpha-helix exhibits weak cooperativity, which is consistent with the experimental data. Meanwhile, the folding pathway is compared between the FSD-1 and the target zinc finger peptide, and the possible role of the zinc ion on the folding pathway of zinc finger is proposed. Proteins 2007. (c) 2007 Wiley-Liss, Inc.     

Reconstruction of the src-SH3 protein domain transition state ensemble using multiscale molecular dynamics simulations

We use an integrated computational approach to reconstruct accurately the transition state ensemble (TSE) for folding of the src-SH3 protein domain. We first identify putative TSE conformations from free energy surfaces generated by importance sampling molecular dynamics for a fully atomic, solvated model of the src-SH3 protein domain. These putative TSE conformations are then subjected to a folding analysis using a coarse-grained representation of the protein and rapid discrete molecular dynamics simulations. Those conformations that fold to the native conformation with a probability (P(fold)) of approximately 0.5, constitute the true transition state. Approximately 20% of the putative TSE structures were found to have a P(fold) near 0.5, indicating that, although correct TSE conformations are populated at the free energy barrier, there is a critical need to refine this ensemble. Our simulations indicate that the true TSE conformations are compact, with a well-defined central beta sheet, in good agreement with previous experimental and theoretical studies. A structured central beta sheet was found to be present in a number of pre-TSE conformations, however, indicating that this element, although required in the transition state, does not define it uniquely. An additional tight cluster of contacts between highly conserved residues belonging to the diverging turn and second beta-sheet of the protein emerged as being critical elements of the folding nucleus. A number of commonly used order parameters to identify the transition state for folding were investigated, with the number of native Cbeta contacts displaying the most satisfactory correlation with P(fold) values.     

Motional timescale predictions by molecular dynamics simulations: Case study using proline and hydroxyproline sidechain dynamics

We propose a new approach for force field optimizations which aims at reproducing dynamics characteristics using biomolecular MD simulations, in addition to improved prediction of motionally averaged structural properties available from experiment. As the source of experimental data for dynamics fittings, we use ¹³C NMR spin‐lattice relaxation times T₁ of backbone and sidechain carbons, which allow to determine correlation times of both overall molecular and intramolecular motions. For structural fittings, we use motionally averaged experimental values of NMR J couplings. The proline residue and its derivative 4‐hydroxyproline with relatively simple cyclic structure and sidechain dynamics were chosen for the assessment of the new approach in this work. Initially, grid search and simplexed MD simulations identified large number of parameter sets which fit equally well experimental J couplings. Using the Arrhenius‐type relationship between the force constant and the correlation time, the available MD data for a series of parameter sets were analyzed to predict the value of the force constant that best reproduces experimental timescale of the sidechain dynamics. Verification of the new force‐field (termed as AMBER99SB‐ILDNP) against NMR J couplings and correlation times showed consistent and significant improvements compared to the original force field in reproducing both structural and dynamics properties. The results suggest that matching experimental timescales of motions together with motionally averaged characteristics is the valid approach for force field parameter optimization. Such a comprehensive approach is not restricted to cyclic residues and can be extended to other amino acid residues, as well as to the backbone. Proteins 2014; 82:195–215. © 2013 Wiley Periodicals, Inc.     

β‐Hairpin folding and stability: molecular dynamics simulations of designed peptides in aqueous solution

The structural properties of a 10‐residue and a 15‐residue peptide in aqueous solution were investigated by molecular dynamics simulation. The two designed peptides, SYINSDGTWT and SESYINSDGTWTVTE, had been studied previously by NMR at 278 K and the resulting model structures were classified as 3:5 β‐hairpins with a type I + G1 β‐bulge turn. In simulations at 278 K, starting from the NMR model structure, the 3:5 β‐hairpin conformers proved to be stable over the time period evaluated (30 ns). Starting from an extended conformation, simulations of the decapeptide at 278 K, 323 K and 353 K were also performed to study folding. Over the relatively short time scales explored (30 ns at 278 K and 323 K, 56 ns at 353 K), folding to the 3:5 β‐hairpin could only be observed at 353 K. At this temperature, the collapse to β‐hairpin‐like conformations is very fast. The conformational space accessible to the peptide is entirely dominated by loop structures with different degrees of β‐hairpin character. The transitions between different types of ordered loops and β‐hairpins occur through two unstructured loop conformations stabilized by a single side‐chain interaction between Tyr2 and Trp9, which facilitates the changes of the hydrogen‐bond register. In agreement with previous experimental results, β‐hairpin formation is initially driven by the bending propensity of the turn segment. Nevertheless, the fine organization of the turn region appears to be a late event in the folding process. Copyright © 2004 European Peptide Society and John Wiley & Sons, Ltd.