Replica exchange simulation of reversible folding/unfolding of the Trp-cage miniprotein in explicit solvent: on the structure and possible role of internal water

We simulate the folding/unfolding equilibrium of the 20-residue miniprotein Trp-cage. We use replica exchange molecular dynamics simulations of the AMBER94 atomic detail model of the protein explicitly solvated by water, starting from a completely unfolded configuration. We employ a total of 40 replicas, covering the temperature range between 280 and 538 K. Individual simulation lengths of 100 ns sum up to a total simulation time of about 4 micros. Without any bias, we observe the folding of the protein into the native state with an unfolding-transition temperature of about 440 K. The native state is characterized by a distribution of root mean square distances (RMSD) from the NMR data that peaks at 1.8A, and is as low as 0.4A. We show that equilibration times of about 40 ns are required to yield convergence. A folded configuration in the entire extended ensemble is found to have a lifetime of about 31 ns. In a clamp-like motion, the Trp-cage opens up during thermal denaturation. In line with fluorescence quenching experiments, the Trp-residue sidechain gets hydrated when the protein opens up, roughly doubling the number of water molecules in the first solvation shell. We find the helical propensity of the helical domain of Trp-cage rather well preserved even at very high temperatures. In the folded state, we can identify states with one and two buried internal water molecules interconnecting parts of the Trp-cage molecule by hydrogen bonds. The loss of hydrogen bonds of these buried water molecules in the folded state with increasing temperature is likely to destabilize the folded state at elevated temperatures.     

Folding Trp-cage to NMR resolution native structure using a coarse-grained protein model

We develop a coarse-grained protein model with a simplified amino acid interaction potential. Using this model, we perform discrete molecular dynamics folding simulations of a small 20-residue protein--Trp-cage--from a fully extended conformation. We demonstrate the ability of the Trp-cage model to consistently reach conformations within 2-angstroms backbone root-mean-square distance from the corresponding NMR structures. The minimum root-mean-square distance of Trp-cage conformations in simulations can be <1 angstroms. Our findings suggest that, at least in the case of Trp-cage, a detailed all-atom protein model with a molecular mechanics force field is not necessary to reach the native state of a protein. Our results also suggest that the success of folding Trp-cage in our simulations and in the reported all-atom molecular mechanics simulation studies may be mainly due to the special stabilizing features specific to this miniprotein.     

Possible locally driven folding pathways of TC5b,a 20-residue protein

A novel computational procedure for modeling possible locally driven folding pathways by stepwise elongations of the peptide chain was successfully applied to TC5b, a 20-residue miniprotein. Systematic exploration of the possible locally driven pathways showed that the Trp-cage structure of TC5b could be obtained by stepwise elongation starting from the noncentral local nucleation centers preexisting in the unfolded state of TC5b. The probable locally driven folding pathway starts with folding of alpha-helical fragment 4-9, followed by formation of the proper three-dimensional structure of fragment 4-12, and then 4-18. Accordingly, the Trp-cage-forming interactions emerge successively, first Trp(6)-Pro(12), then Trp(6)-Pro(18), and then Trp(6)-Tyr(3). The Trp-cage-like structures of TC5b found in this study by independent energy calculations are in excellent agreement with the NMR experimental data. The same procedure rationalizes the incomplete Trp-cage formation observed for two analogs of TC5b. Generally, the success of this novel approach is encouraging and provides some justification for the use of computational simulations of locally driven protein folding.     

Folding pathways for initiator and effector procaspases from computer simulations

The folding pathways of procaspases 3, 7, and 8 have been studied using a Go-like Hamiltonian and molecular dynamics simulations coupled with a parallel tempering scheme. The folding pathways and the overall structures of procaspases 3 and 7 are similar, and are characterized by monomeric as well as dimeric folding intermediates in agreement with the available structural and thermochemical data. The folding pathway of procaspase 8, on the other hand, is characterized by a larger population of monomers and partially folded dimer intermediates, and only a relatively small population of folded dimer species. The most stable structure predicted for procaspase 8 is a dimer, in which the position of the linker is remarkably different from the one observed in procaspases 3 and 7, leading to the fact that all the contacts that stabilize the active site are essentially formed. This novel and unexpected structure provides a rationale for the observed activity of the procaspase 8 dimer, and thus could be highly relevant for the initiation of FAS-mediated apoptosis.     

Kinetic studies of folding of the B-domain of staphylococcal protein A with molecular dynamics and a united-residue (UNRES) model of polypeptide chains

Langevin dynamics is used with our physics-based united-residue (UNRES) force field to study the folding pathways of the B-domain of staphylococcal protein A (1BDD (alpha; 46 residues)). With 400 trajectories of protein A started from the extended state (to gather meaningful statistics), and simulated for more than 35 ns each, 380 of them folded to the native structure. The simulations were carried out at the optimal folding temperature of protein A with this force field. To the best of our knowledge, this is the first simulation study of protein-folding kinetics with a physics-based force field in which reliable statistics can be gathered. In all the simulations, the C-terminal alpha-helix forms first. The ensemble of the native basin has an average RMSD value of 4 A from the native structure. There is a stable intermediate along the folding pathway, in which the N-terminal alpha-helix is unfolded; this intermediate appears on the way to the native structure in less than one-fourth of the folding pathways, while the remaining ones proceed directly to the native state. Non-native structures persist until the end of the simulations, but the native-like structures dominate. To express the kinetics of protein A folding quantitatively, two observables were used: (i) the average alpha-helix content (averaged over all trajectories within a given time window); and (ii) the fraction of conformations (averaged over all trajectories within a given time window) with Calpha RMSD values from the native structure less than 5 A (fraction of completely folded structures). The alpha-helix content grows quickly with time, and its variation fits well to a single-exponential term, suggesting fast two-state kinetics. On the other hand, the fraction of folded structures changes more slowly with time and fits to a sum of two exponentials, in agreement with the appearance of the intermediate, found when analyzing the folding pathways. This observation demonstrates that different qualitative and quantitative conclusions about folding kinetics can be drawn depending on which observable is monitored.     

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.     

The high-resolution NMR structure of the early folding intermediate of the Thermus thermophilus ribonuclease H

Elucidation of the high-resolution structures of folding intermediates is a necessary but difficult step toward the ultimate understanding of the mechanism of protein folding. Here, using hydrogen-exchange-directed protein engineering, we populated the folding intermediate of the Thermus thermophilus ribonuclease H, which forms before the rate-limiting transition state, by removing the unfolded regions of the intermediate, including an α-helix and two β-strands (51 folded residues). Using multidimensional NMR, we solved the structure of this intermediate mimic to an atomic resolution (backbone rmsd, 0.51 Å). It has a native-like backbone topology and shows some local deviations from the native structure, revealing that the structure of the folded region of an early folding intermediate can be as well defined as the native structure. The topological parameters calculated from the structures of the intermediate mimic and the native state predict that the intermediate should fold on a millisecond time scale or less and form much faster than the native state. Other factors that may lead to the slow folding of the native state and the accumulation of the intermediate before the rate-limiting transition state are also discussed.     

Simulation and experiment at high temperatures: ultrafast folding of a thermophilic protein by nucleation-condensation

We used Phi-value analysis to characterise the transition state for folding of a thermophilic protein at the relatively high temperature of 325 K. PhiF values for the folding of the three-helix bundle, peripheral subunit binding domain from Bacillus stearothermophilus (E3BD) were determined by temperature-jump experiments in the absence of chemical denaturants. E3BD folded in microseconds through a highly diffuse transition state. Excellent agreement was observed between experiment and the results from eight (independent) molecular dynamics simulations of unfolding at 373 K. We used a combination of heteronuclear NMR experiments and molecular dynamics simulations to characterise the denatured ensemble, and found that it contained very little persistent, residual structure. However, those regions that adopt helical structure in the native state were found by simulation to be poised for helix formation in the denatured state. These regions also had significant structure in the transition state for folding. The overall folding pathway appears to be nucleation-condensation.     

Three force fields' views of the 3(10) helix

Slowly but steadily bibliographic evidence is accumulating that the apparent convergence of the various biomolecular force fields as evidenced from simulations of proteins in the folded state does not hold true for folding simulations. Here we add one more example to the growing list of peptides and proteins for which different force fields show irreconcilable differences in their folding predictions, even at such a fundamental level as that of a peptide's secondary structure. We show that for an undecamer peptide that is known from two independent NMR structure determinations to have a mainly 3₁₀-helical structure in solution, three mainstream biomolecular force fields give completely disparate predictions: The CHARMM force field (with the CMAP correction) predicts an outstandingly stable α-helical structure, in disagreement not only with the experimental structures, but also with experimental evidence obtained from circular dichroism. OPLS-AA shows an almost totally disordered peptide with the most frequently observed folded conformation corresponding to a β-hairpin-like structure, again in disagreement with all available experimental evidence. Only the AMBER99SB force field appears to qualitatively agree with not only the general structural characteristics of the peptide (on the account of both NMR- and CD-based experiments), but to also correctly predict some of the experimentally observed interactions at the level of side chains. Possible interpretations of these findings are discussed.     

A Compact Native 24-Residue Supersecondary Structure Derived from the Villin Headpiece Subdomain

Many small proteins fold highly cooperatively in an all-or-none fashion and thus their native states are well protected from thermal fluctuations by an extensive network of interactions across the folded structure. Because protein structures are stabilized by local and nonlocal interactions among distal residues, dissecting individual substructures from the context of folded proteins results in large destabilization and loss of unique three-dimensional structure. Thus, mini-protein substructures can only rarely be derived from natural templates. Here, we describe a compact native 24-residues-long supersecondary structure derived from the hyperstable villin headpiece subdomain consisting of helices 2 and 3 (HP24). Using a combination of experimental techniques, including NMR and small-angle x-ray scattering, as well as all-atom replica exchange molecular-dynamics simulations, we show that a variant with stabilizing substitutions (HP24stab) forms a densely packed and compact conformation. In HP24stab, interactions between helices 2 and 3 are similar to those observed in native folded HP35, and the two helices cooperatively stabilize each other by completing the hydrophobic core lining the central part of HP35. Interestingly, even though the HP24wt fragment shows a more expanded and less structured conformation, NMR and simulations demonstrate a preference for a native-like topology. Thus, the two stabilizing residues in HP24stab shift the energy balance toward the native state, leading to a minimal folding motif.     

Study of the Villin headpiece folding dynamics by combining coarse-grained Monte Carlo evolution and all-atom molecular dynamics

The folding mechanism of the Villin headpiece (HP36) is studied by means of a novel approach which entails an initial coarse-grained Monte Carlo (MC) scheme followed by all-atom molecular dynamics (MD) simulations in explicit solvent. The MC evolution occurs in a simplified free-energy landscape and allows an efficient selection of marginally-compact structures which are taken as viable initial conformations for the MD. The coarse-grained MC structural representation is connected to the one with atomic resolution through a "fine-graining" reconstruction algorithm. This two-stage strategy is used to select and follow the dynamics of seven different unrelated conformations of HP36. In a notable case the MD trajectory rapidly evolves towards the folded state, yielding a typical root-mean-square deviation (RMSD) of the core region of only 2.4 A from the closest NMR model (the typical RMSD over the whole structure being 4.0 A). The analysis of the various MC-MD trajectories provides valuable insight into the details of the folding and mis-folding mechanisms and particularly about the delicate influence of local and nonlocal interactions in steering the folding process.     

The nature of folded states of globular proteins.

We suggest, using dynamical simulations of a simple heteropolymer modelling the alpha-carbon sequence in a protein, that generically the folded states of globular proteins correspond to statistically well-defined metastable states. This hypothesis, called the metastability hypothesis, states that there are several free energy minima separated by barriers of various heights such that the folded conformations of a polypeptide chain in each of the minima have similar structural characteristics but have different energies from one another. The calculated structural characteristics, such as bond angle and dihedral angle distribution functions, are assumed to arise from only those configurations belonging to a given minimum. The validity of this hypothesis is illustrated by simulations of a continuum model of a heteropolymer whose low temperature state is a well-defined beta-barrel structure. The simulations were done using a molecular dynamics algorithm (referred to as the "noisy" molecular dynamics method) containing both friction and noise terms. It is shown that for this model there are several distinct metastable minima in which the structural features are similar. Several new methods of analyzing fluctuations in structures belonging to two distinct minima are introduced. The most notable one is a dynamic measure of compactness that can in principle provide the time required for maximal compactness to be achieved. The analysis shows that for a given metastable state in which the protein has a well-defined folded structure the transition to a state of higher compactness occurs very slowly, lending credence to the notion that the system encounters a late barrier in the process of folding to the most compact structure. The examination of the fluctuations in the structures near the unfolding----folding transition temperature indicates that the transition state for the unfolding to folding process occurs closer to the folded state.     

Folding and unfolding of gammaTIM monomers and dimers

Kinetic simulations of the folding and unfolding of triosephosphate isomerase (TIM) from yeast were conducted using a single monomer gammaTIM polypeptide chain that folds as a monomer and two gammaTIM chains that fold to the native dimer structure. The basic protein model used was a minimalist Gō model using the native structure to determine attractive energies in the protein chain. For each simulation type--monomer unfolding, monomer refolding, dimer unfolding, and dimer refolding--thirty simulations were conducted, successfully capturing each reaction in full. Analysis of the simulations demonstrates four main conclusions. First, all four simulation types have a similar "folding order", i.e., they have similar structures in intermediate stages of folding between the unfolded and folded state. Second, despite this similarity, different intermediate stages are more or less populated in the four different simulations, with 1), no intermediates populated in monomer unfolding; 2), two intermediates populated with beta(2)-beta(4) and beta(1)-beta(5) regions folded in monomer refolding; 3), two intermediates populated with beta(2)-beta(3) and beta(2)-beta(4) regions folded in dimer unfolding; and 4), two intermediates populated with beta(1)-beta(5) and beta(1)-beta(5) + beta(6) + beta(7) + beta(8) regions folded in dimer refolding. Third, simulations demonstrate that dimer binding and unbinding can occur early in the folding process before complete monomer-chain folding. Fourth, excellent agreement is found between the simulations and MPAX (misincorporation proton alkyl exchange) experiments. In total, this agreement demonstrates that the computational Gō model is accurate for gammaTIM and that the energy landscape of gammaTIM appears funneled to the native state.     

Multiplexed-Replica Exchange Molecular Dynamics Method for Protein Folding Simulation

Simulating protein folding thermodynamics starting purely from a protein sequence is a grand challenge of computational biology. Here, we present an algorithm to calculate a canonical distribution from molecular dynamics simulation of protein folding. This algorithm is based on the replica exchange method where the kinetic trapping problem is overcome by exchanging noninteracting replicas simulated at different temperatures. Our algorithm uses multiplexed-replicas with a number of independent molecular dynamics runs at each temperature. Exchanges of configurations between these multiplexed-replicas are also tried, rendering the algorithm applicable to large-scale distributed computing (i.e., highly heterogeneous parallel computers with processors having different computational power). We demonstrate the enhanced sampling of this algorithm by simulating the folding thermodynamics of a 23 amino acid miniprotein. We show that better convergence is achieved compared to constant temperature molecular dynamics simulation, with an efficient scaling to large number of computer processors. Indeed, this enhanced sampling results in (to our knowledge) the first example of a replica exchange algorithm that samples a folded structure starting from a completely unfolded state.     

Folding of a DNA hairpin loop structure in explicit solvent using replica-exchange molecular dynamics simulations

Hairpin loop structures are common motifs in folded nucleic acids. The 5'-GCGCAGC sequence in DNA forms a characteristic and stable trinucleotide hairpin loop flanked by a two basepair stem helix. To better understand the structure formation of this hairpin loop motif in atomic detail, we employed replica-exchange molecular dynamics (RexMD) simulations starting from a single-stranded DNA conformation. In two independent 36 ns RexMD simulations, conformations in very close agreement with the experimental hairpin structure were sampled as dominant conformations (lowest free energy state) during the final phase of the RexMDs ( approximately 35% at the lowest temperature replica). Simultaneous compaction and accumulation of folded structures were observed. Comparison of the GCA trinucleotides from early stages of the simulations with the folded topology indicated a variety of central loop conformations, but arrangements close to experiment that are sampled before the fully folded structure also appeared. Most of these intermediates included a stacking of the C(2) and G(3) bases, which was further stabilized by hydrogen bonding to the A(5) base and a strongly bound water molecule bridging the C(2) and A(5) in the DNA minor groove. The simulations suggest a folding mechanism where these intermediates can rapidly proceed toward the fully folded hairpin and emphasize the importance of loop and stem nucleotide interactions for hairpin folding. In one simulation, a loop motif with G(3) in syn conformation (dihedral flip at N-glycosidic bond) accumulated, resulting in a misfolded hairpin. Such conformations may correspond to long-lived trapped states that have been postulated to account for the folding kinetics of nucleic acid hairpins that are slower than expected for a semiflexible polymer of the same size.     

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.

     

Theory of protein folding

Protein folding should be complex. Proteins organize themselves into specific three-dimensional structures, through a myriad of conformational changes. The classical view of protein folding describes this process as a nearly sequential series of discrete intermediates. In contrast, the energy landscape theory of folding considers folding as the progressive organization of an ensemble of partially folded structures through which the protein passes on its way to the natively folded structure. As a result of evolution, proteins have a rugged funnel-like landscape biased toward the native structure. Connecting theory and simulations of minimalist models with experiments has completely revolutionized our understanding of the underlying mechanisms that control protein folding.     

All-atom Monte Carlo simulation of GCAA RNA folding

We report a detailed all-atom simulation of the folding of the GCAA RNA tetraloop. The GCAA tetraloop motif is a very common and thermodynamically stable secondary structure in natural RNAs. We use our simulation methods to study the folding behavior of a 12-base GCAA tetraloop structure with a four-base helix adjacent to the tetraloop proper. We implement an all-atom Monte Carlo (MC) simulation of RNA structural dynamics using a Go potential. Molecular dynamics (MD) simulation of RNA and protein has realistic energetics and sterics, but is extremely expensive in terms of computational time. By coarsely treating non-covalent energetics, but retaining all-atom sterics and entropic effects, all-atom MC techniques are a useful method for the study of protein and now RNA. We observe a sharp folding transition for this structure, and in simulations at room temperature the state histogram shows three distinct minima: an unfolded state (U), a more narrow intermediated state (I), and a narrow folded state (F). The intermediate consists primarily of structures with the GCAA loop and some helix hydrogen bonds formed. Repeated kinetic folding simulations reveal that the number of helix base-pairs forms a simple 1D reaction coordinate for the I-->N transition.