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.     

Decoding Structural Properties of a Partially Unfolded Protein Substrate: En Route to Chaperone Binding

Many proteins comprising of complex topologies require molecular chaperones to achieve their unique three-dimensional folded structure. The E.coli chaperone, GroEL binds with a large number of unfolded and partially folded proteins, to facilitate proper folding and prevent misfolding and aggregation. Although the major structural components of GroEL are well defined, scaffolds of the non-native substrates that determine chaperone-mediated folding have been difficult to recognize. Here we performed all-atomistic and replica-exchange molecular dynamics simulations to dissect non-native ensemble of an obligate GroEL folder, DapA. Thermodynamics analyses of unfolding simulations revealed populated intermediates with distinct structural characteristics. We found that surface exposed hydrophobic patches are significantly increased, primarily contributed from native and non-native β-sheet elements. We validate the structural properties of these conformers using experimental data, including circular dichroism (CD), 1-anilinonaphthalene-8-sulfonic acid (ANS) binding measurements and previously reported hydrogen-deutrium exchange coupled to mass spectrometry (HDX-MS). Further, we constructed network graphs to elucidate long-range intra-protein connectivity of native and intermediate topologies, demonstrating regions that serve as central “hubs”. Overall, our results implicate that genomic variations (or mutations) in the distinct regions of protein structures might disrupt these topological signatures disabling chaperone-mediated folding, leading to formation of aggregates.     

Protein refolding versus aggregation: computer simulations on an intermediate-resolution protein model.

Computer simulations are performed on a system of eight model peptide chains to study how the competition between protein refolding and aggregation affects the optimal conditions for refolding of four-helix bundles. The discontinuous molecular dynamics algorithm is utilized along with an intermediate-resolution protein model that we developed for this work. Physically, the model is much more detailed than any model used to date for simulations of protein aggregation. Each model residue consists of a detailed, three-bead backbone and a simplified, single-bead side-chain. Excluded volume, hydrogen bond, and hydrophobic interactions are modeled with discontinuous (i.e. hard-sphere and square-well) potentials. Simulations efficiently sample conformational space, and complete folding trajectories from random initial configurations to two four-helix bundles are possible within two days on a single processor workstation. Folding of the bundles follows two main pathways, one through a trimeric intermediate and the other through an intermediate with two dimers. The proportion of trajectories that follow each route is significantly different for the eight-peptide system in this work than in a previously studied four-peptide system, which yields one four-helix bundle, suggesting, as our previous simulations have, that protein folding properties are strongly influenced by the presence of other proteins. Folding of the bundles is optimal within a fixed temperature range, with the high-temperature boundary a function of the complexity of the protein (or oligomer) to be folded and the low-temperature boundary a function of the complexity of the protein's environment. Above the optimal temperature range for folding, the model chains tend to unfold; below the optimal range, the model chains tend to aggregate. As has been seen previously, aggregates have substantial levels of native secondary structure, suggesting that aggregates are composed largely of partially folded intermediates, not denatured chains.     

Secondary Structure Propensities in Peptide Folding Simulations: A Systematic Comparison of Molecular Mechanics Interaction Schemes

We present a systematic study directed toward the secondary structure propensity and sampling behavior in peptide folding simulations with eight different molecular dynamics force-field variants in explicit solvent. We report on the combinational result of force field, water model, and electrostatic interaction schemes and compare to available experimental characterization of five studied model peptides in terms of reproduced structure and dynamics. The total simulation time exceeded 18 μs and included simulations that started from both folded and extended conformations. Despite remaining sampling issues, a number of distinct trends in the folding behavior of the peptides emerged. Pronounced differences in the propensity of finding prominent secondary structure motifs in the different applied force fields suggest that problems point in particular to the balance of the relative stabilities of helical and extended conformations.     

Structural analysis of kinetic folding intermediates for a TIM barrel protein, indole-3-glycerol phosphate synthase, by hydrogen exchange mass spectrometry and Gō model simulation

The structures of partially folded states appearing during the folding of a (βα)₈ TIM barrel protein, the indole-3-glycerol phosphate synthase from Sulfolobus solfataricus (sIGPS), was assessed by hydrogen exchange mass spectrometry (HX-MS) and Gō model simulations. HX-MS analysis of the peptic peptides derived from the pulse-labeled product of the sub-millisecond folding reaction from the urea-denatured state revealed strong protection in the (βα)₄ region, modest protection in the neighboring (βα)_1-3 and (βα)₅β₆ segments and no significant protection in the remaining N and C-terminal segments. These results demonstrate that this species is not a collapsed form of the unfolded state under native-favoring conditions nor is it the native state formed via fast-track folding. However, the striking contrast of these results with the strong protection observed in the (βα)_2-5β₆ region after 5 s of folding demonstrates that these species represent kinetically distinct folding intermediates that are not identical as previously thought. A re-examination of the kinetic folding mechanism by chevron analysis of fluorescence data confirmed distinct roles for these two species: the burst-phase intermediate is predicted to be a misfolded, off-pathway intermediate, while the subsequent 5 s intermediate corresponds to an on-pathway equilibrium intermediate. Comparison with the predictions using a C^α Gō model simulation of the kinetic folding reaction for sIGPS shows good agreement with the core of the structure offering protection against exchange in the on-pathway intermediate(s). Because the native-centric Gō model simulations do not explicitly include sequence-specific information, the simulation results support the hypothesis that the topology of TIM barrel proteins is a primary determinant of the folding free energy surface for the productive folding reaction. The early misfolding reaction must involve aspects of non-native structure not detected by the Gō model simulation.     

The Folding Trajectory of RNase H Is Dominated by Its Topology and Not Local Stability: A Protein Engineering Study of Variants that Fold via Two-State and Three-State Mechanisms

Proteins can sample a variety of partially folded conformations during the transition between the unfolded and native states. Some proteins never significantly populate these high-energy states and fold by an apparently two-state process. However, many proteins populate detectable, partially folded forms during the folding process. The role of such intermediates is a matter of considerable debate. A single amino acid change can convert Escherichia coli ribonuclease H from a three-state folder that populates a kinetic intermediate to one that folds in an apparent two-state fashion. We have compared the folding trajectories of the three-state RNase H and the two-state RNase H, proteins with the same native-state topology but altered regional stability, using a protein engineering approach. Our data suggest that both versions of RNase H fold through a similar trajectory with similar high-energy conformations. Mutations in the core and the periphery of the protein affect similar aspects of folding for both variants, suggesting a common trajectory with folding of the core region followed by the folding of the periphery. Our results suggest that formation of specific partially folded conformations may be a general feature of protein folding that can promote, rather than hinder, efficient folding.     

Separability between overall and internal motion: A protein folding problem

The separability between overall and internal motions is evaluated over multiple folding trajectories of the villin headpiece subdomain. The analysis, which relies on the Prompers‐Brüschweiler separability index, offers a potentially useful perspective on protein folding. The protein is considered folded in this study, not when it reaches some static target, but rather when it tumbles as a dynamically constrained object. The analysis also demonstrates how the separability index, when applied to protein folding simulations, can facilitate the analysis of NMR relaxation data. Proteins 2012;. © 2012 Wiley Periodicals, Inc.     

Folding Simulations of a De Novo Designed Protein with a βαβ Fold

βαβ structural motifs are commonly used building blocks in protein structures containing parallel β-sheets. However, to our knowledge, no stand-alone βαβ structure has been observed in nature to date. Recently, for the first time that we know of, a small protein with an independent βαβ structure (DS119) was successfully designed in our laboratory. To understand the folding mechanism of DS119, in the study described here, we carried out all-atom molecular dynamics and coarse-grained simulations to investigate its folding pathways and energy landscape. From all-atom simulations, we successfully observed the folding event and got a stable folded structure with a minimal root mean-square deviation of 2.6 Å with respect to the NMR structure. The folding process can be described as a fast collapse phase followed by rapid formation of the central helix, and then slow formation of a parallel β-sheet. By using a native-centric Gō-like model, the cooperativity of the system was characterized in terms of the calorimetric criterion, sigmoidal transitions, conformation distribution shifts, and free-energy profiles. DS119 was found to be an incipient downhill folder that folds more cooperatively than a downhill folder, but less cooperatively than a two-state folder. This may reflect the balance between the two structural elements of DS119: the rapidly formed α-helix and the slowly formed parallel β-sheet. Folding times estimated from both the all-atom simulations and the coarse-grained model were at microsecond level, making DS119 another fast folder. Compared to fast folders reported previously, DS119 is, to the best of our knowledge, the first that exhibits a parallel β-sheet.     

Hairpins participating in folding of human telomeric sequence quadruplexes studied by standard and T-REMD simulations

DNA G-hairpins are potential key structures participating in folding of human telomeric guanine quadruplexes (GQ). We examined their properties by standard MD simulations starting from the folded state and long T-REMD starting from the unfolded state, accumulating ∼130 μs of atomistic simulations. Antiparallel G-hairpins should spontaneously form in all stages of the folding to support lateral and diagonal loops, with sub-μs scale rearrangements between them. We found no clear predisposition for direct folding into specific GQ topologies with specific syn/anti patterns. Our key prediction stemming from the T-REMD is that an ideal unfolded ensemble of the full GQ sequence populates all 4096 syn/anti combinations of its four G-stretches. The simulations can propose idealized folding pathways but we explain that such few-state pathways may be misleading. In the context of the available experimental data, the simulations strongly suggest that the GQ folding could be best understood by the kinetic partitioning mechanism with a set of deep competing minima on the folding landscape, with only a small fraction of molecules directly folding to the native fold. The landscape should further include non-specific collapse processes where the molecules move via diffusion and consecutive random rare transitions, which could, e.g. structure the propeller loops.     

Insights into nucleic acid conformational dynamics from massively parallel stochastic simulations

The helical hairpin is one of the most ubiquitous and elementary secondary structural motifs in nucleic acids, capable of serving functional roles and participating in long-range tertiary contacts. Yet the self-assembly of these structures has not been well-characterized at the atomic level. With this in mind, the dynamics of nucleic acid hairpin formation and disruption have been studied using a novel computational tool: large-scale, parallel, atomistic molecular dynamics simulation employing an inhomogeneous distributed computer consisting of more than 40,000 processors. Using multiple methodologies, over 500 micro s of atomistic simulation time has been collected for a large ensemble of hairpins (sequence 5'-GGGC[GCAA]GCCU-3'), allowing characterization of rare events not previously observable in simulation. From uncoupled ensemble dynamics simulations in unperturbed folding conditions, we report on 1), competing pathways between the folded and unfolded regions of the conformational space; 2), observed nonnative stacking and basepairing traps; and 3), a helix unwinding-rewinding mode that is differentiated from the unfolding and folding dynamics. A heterogeneous transition state ensemble is characterized structurally through calculations of conformer-specific folding probabilities and a multiplexed replica exchange stochastic dynamics algorithm is used to derive an approximate folding landscape. A comparison between the observed folding mechanism and that of a peptide beta-hairpin analog suggests that although native topology defines the character of the folding landscape, the statistical weighting of potential folding pathways is determined by the chemical nature of the polymer.     

Beta-hairpin folding simulations in atomistic detail using an implicit solvent model

We have used distributed computing techniques and a supercluster of thousands of computer processors to study folding of the C-terminal beta-hairpin from protein G in atomistic detail using the GB/SA implicit solvent model at 300 K. We have simulated a total of nearly 38 micros of folding time and obtained eight complete and independent folding trajectories. Starting from an extended state, we observe relaxation to an unfolded state characterized by non-specific, temporary hydrogen bonding. This is followed by the appearance of interactions between hydrophobic residues that stabilize a bent intermediate. Final formation of the complete hydrophobic core occurs cooperatively at the same time that the final hydrogen bonding pattern appears. The folded hairpin structures we observe all contain a closely packed hydrophobic core and proper beta-sheet backbone dihedral angles, but they differ in backbone hydrogen bonding pattern. We show that this is consistent with the existing experimental data on the hairpin alone in solution. Our analysis also reveals short-lived semi-helical intermediates which define a thermodynamic trap. Our results are consistent with a three-state mechanism with a single rate-limiting step in which a varying final hydrogen bond pattern is apparent, and semi-helical off-pathway intermediates may appear early in the folding process. We include details of the ensemble dynamics methodology and a discussion of our achievements using this new computational device for studying dynamics at the atomic level.     

Force Field Bias in Protein Folding Simulations

Long timescale (>1 μs) molecular dynamics simulations of protein folding offer a powerful tool for understanding the atomic-scale interactions that determine a protein's folding pathway and stabilize its native state. Unfortunately, when the simulated protein fails to fold, it is often unclear whether the failure is due to a deficiency in the underlying force fields or simply a lack of sufficient simulation time. We examine one such case, the human Pin1 WW domain, using the recently developed deactivated morphing method to calculate free energy differences between misfolded and folded states. We find that the force field we used favors the misfolded states, explaining the failure of the folding simulations. Possible further applications of deactivated morphing and implications for force field development are discussed.     

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.     

Assembly of a tetrameric alpha-helical bundle: computer simulations on an intermediate-resolution protein model.

Discontinuous molecular dynamics (DMD) simulation on an intermediate-resolution protein model is used to study the folding of an isolated, small model peptide to an amphipathic alpha-helix and the assembly of four of these model peptides into a four-helix bundle. A total of 129 simulations were performed on the isolated peptide, and 50 simulations were performed on the four-peptide system. Simulations efficiently sample conformational space allowing complete folding trajectories from random initial configurations to be observed within 15 min for the one-peptide system and within 15 h for the four-peptide system on a 500-MHz workstation. The native structures of both the alpha-helix and the four-helix bundle are consistent with experimental characterization studies and with results from previous simulations on these model peptides. In both the one- and four-peptide systems, the native state is achieved during simulations within an optimal temperature range, a phenomenon also observed experimentally. The ease with which our simulations yield reasonable estimates of folded structures demonstrates the power of the intermediate-resolution model developed for this work and the DMD algorithm and suggests that simulations of very long times and of multiprotein systems may be possible with this model.     

Heterogeneity even at the speed limit of folding: large-scale molecular dynamics study of a fast-folding variant of the villin headpiece

We have performed molecular dynamics simulations on a set of nine unfolded conformations of the fastest-folding protein yet discovered, a variant of the villin headpiece subdomain (HP-35 NleNle). The simulations were generated using a new distributed computing method, yielding hundreds of trajectories each on a time scale comparable to the experimental folding time, despite the large (10,000 atom) size of the simulation system. This strategy eliminates the need to assume a two-state kinetic model or to build a Markov state model. The relaxation to the folded state at 300 K from the unfolded configurations (generated by simulation at 373 K) was monitored by a method intended to reflect the experimental observable (quenching of tryptophan by histidine). We also monitored the relaxation to the native state by directly comparing structural snapshots with the native state. The rate of relaxation to the native state and the number of resolvable kinetic time scales both depend upon starting structure. Moreover, starting structures with folding rates most similar to experiment show some native-like structure in the N-terminal helix (helix 1) and the phenylalanine residues constituting the hydrophobic core, suggesting that these elements may exist in the experimentally relevant unfolded state. Our large-scale simulation data reveal kinetic complexity not resolved in the experimental data. Based on these findings, we propose additional experiments to further probe the kinetics of villin folding.     

Probing Position-Dependent Diffusion in Folding Reactions Using Single-Molecule Force Spectroscopy

Folding of proteins and nucleic acids involves a diffusive search over a multidimensional conformational energy landscape for the minimal-energy structure. When examining the projection of conformational motions onto a one-dimensional reaction coordinate, as done in most experiments, the diffusion coefficient D is generally position dependent. However, it has proven challenging to measure such position-dependence experimentally. We investigated the position-dependence of D in the folding of DNA hairpins as a simple model system in two ways: first, by analyzing the round-trip time to return to a given extension in constant-force extension trajectories measured by force spectroscopy, and second, by analyzing the fall time required to reach a given extension in force jump measurements. These methods yielded conflicting results: the fall time implied a fairly constant D, but the round-trip time implied variations of over an order of magnitude. Comparison of experiments with computational simulations revealed that both methods were strongly affected by experimental artifacts inherent to force spectroscopy measurements, which obscured the intrinsic position-dependence of D. Lastly, we applied Kramers’s theory to the kinetics of hairpins with energy barriers located at different positions along the hairpin stem, as a crude probe of D at different stem positions, and we found that D did not vary much as the barrier was moved along the reaction coordinate. This work underlines the difficulties faced when trying to deduce position-dependent diffusion coefficients from experimental folding trajectories.