Folding Kinetics and Unfolded State Dynamics of the GB1 Hairpin from Molecular Simulation

Uncertainties in Markov State Models of Small Proteins

Markov state models are widely used to describe and analyze protein dynamics based on molecular dynamics simulations, specifically to extract functionally relevant characteristic time scales and motions. Particularly for larger biomolecules such as proteins, however, insufficient sampling is a notorious concern and often the source of large uncertainties that are difficult to quantify. Furthermore, there are several other sources of uncertainty, such as choice of the number of Markov states and lag time, choice and parameters of dimension reduction preprocessing step, and uncertainty due to the limited number of observed transitions; the latter is often estimated via a Bayesian approach. Here, we quantified and ranked all of these uncertainties for four small globular test proteins. We found that the largest uncertainty is due to insufficient sampling and initially increases with the total trajectory length T up to a critical tipping point, after which it decreases as 1 / T , thus providing guidelines for how much sampling is required for given accuracy. We also found that single long trajectories yielded better sampling accuracy than many shorter trajectories starting from the same structure. In comparison, the remaining sources of the above uncertainties are generally smaller by a factor of about 5, rendering them less of a concern but certainly not negligible. Importantly, the Bayes uncertainty, commonly used as the only uncertainty estimate, captures only a relatively small part of the true uncertainty, which is thus often drastically underestimated.

GB1 hairpin kinetics: capturing the folding pathway with molecular dynamics, replica exchange and optimal dimensionality reduction

We have performed molecular dynamics (MD) and replica-exchange (REMD) simulations of folding of the GB1 hairpin peptide in aqueous solution. REMD results were consistent with a cooperative zipper folding model. 120 μ s MD trajectories at 320 K yielded relaxation times of 1.8 μ s and 100 ns, with the slower assigned to global folding. The MD folding/unfolding transitions also followed the cooperative zipper model, specifying nucleation at the central turn followed by consecutive hydrogen bond formation. Formation of hydrogen bonds and hydrophobic contacts were highly correlated. Coarse-grained kinetic models constructed with the Optimal Dimensionality Reduction (ODR) approach found a folding time of 3.3 μ s and unfolding time of 4.0 μ s . Additionally, relaxation times in the 130-170 ns range could be assigned to formation of the transition state and off-path intermediates. The unfolded state was the most highly populated and, significantly, most heterogenous, assembling the largest number of microstates, primarily composed of extended and turn structures. The folded state was also heterogenous, but a to a lesser degree, involving the fully folded and partially folded in-register hairpins at early stages of the zipper pathway. The transition state corresponded to the nucleated hairpin, with central turn and first beta-sheet hydrogen bond, while the off-path intermediates were off-register partial hairpins. Our simulation results were in excellent agreement with experimental data on folded fraction, relaxation time and folding mechanism. The new findings from this work suggest a highly cooperative zipper folding mechanism, nascent hairpin transition state and ∼100 ns relaxation related to intermediate formation.Communicated by Ramaswamy H. Sarma.

Kirkwood-Buff-Derived Force Field for Peptides and Proteins: Applications of KBFF20

Here, we perform structural, thermodynamic, and kinetics tests of the Kirkwood-Buff-derived force field, KBFF20, for peptides and proteins developed in the previous article. The physical/structural tests measure the ability of KBFF20 to capture the experimental J-couplings for small peptides, to keep globular monomeric and oligomeric proteins folded, and to produce the experimentally relevant expanded conformational ensembles of intrinsically disordered proteins. The thermodynamic-based tests probe KBFF20's ability to quantify the preferential interactions of sodium chloride around native β-lactoglobulin and urea around native lysozyme, to reproduce the melting curves for small helix- and sheet-based peptides, and to fold the small proteins Trp-cage and Villin. The kinetics-based tests quantify how well KBFF20 can match the experimental contact formation rates of small, repeat-sequence peptides of variable lengths and the rotational diffusion coefficients of globular proteins. The results suggest that KBFF20 is naturally able to reproduce properties of both folded and disordered proteins, which we attribute to the use of the Kirkwood-Buff theory as the foundation of the force field's development. However, we show that KBFF20 tends to lose some well-defined secondary structural elements and increases the percentage of coil regions, indicating that the perfect balance of all interactions remains elusive. Nevertheless, we argue that KBFF20 is an improvement over recently modified force fields that require ad hoc interventions to prevent the collapse of intrinsically disordered proteins.

The nature of persistent interactions in two model β-grasp proteins reveals the advantage of symmetry in stability

Two proteins within the β-grasp superfamily, the B1-domain of protein G and the small archaeal modifier protein 1, were investigated to elucidate the key determinants of structural stability at the level of individual interactions. These symmetrical proteins both contain two β-hairpins which form a sheet flanked by a central α-helix. They were subjected to high temperature molecular dynamics simulations and the detailed behavior of each long-range interaction was characterized. The results revealed that in GB1 the most stable region was the C-terminal hairpin and in SAMP1 it was the opposite, the N-terminal hairpin. Experimental results for GB1 support this finding. In conclusion, it appears that the difference in the location and number of hydrophobic interactions dictate the differential stability which is accommodated due to structural symmetry of the β-grasp fold. Thus, the hairpins are interchangeable and in nature this lends itself to adaptability and flexibility.

Effect of heterochiral inversions on the structure of a β-hairpin peptide

We study computationally a family of β-hairpin peptides with systematically introduced chiral inversions, in explicit water, and we investigate the extent to which the backbone structure is able to fold in the presence of heterochiral perturbations. In contrast to the recently investigated case of a helical peptide, we do not find a monotonic change in secondary structure content as a function of the number of L- to D-inversions. The effects of L- to D-inversions are instead found to be highly position-specific. Additionally, in contrast to the helical peptide, some inversions increase the stability of the folded peptide: in such cases, we compute an increase in β-sheet content in the aqueous solution equilibrium ensemble. However, the tertiary structures of the stable (folded) configurations for peptides for which inversions cause an increase in β-sheet content show differences from one another, as well as from the native fold of the nonchirally perturbed β-hairpin. Our results suggest that although some chiral perturbations can increase folding stability, chirally perturbed proteins may still underperform functionally, given the relationship between structure and function.

Dynamics of Structural Elements of GB1 β-Hairpin Revealed by Tryptophan-Cysteine Contact Formation Experiments

Quenching of the triplet state of tryptophan by close contact with cysteine provides a tool for measuring the rate of intramolecular contact formation, one of the most elementary events in the folding process, in peptides and proteins using only natural probes. Here we present a study performed on a stabilized mutant of the second β-hairpin of the GB1 domain, where we combine steady-state fluorescence, laser-induced temperature-jump, and contact formation measurements to unveil the role of elementary structural components on hairpin dynamics and overall stability. In particular, our methodology provides access to the conformational dynamics of both the folded and unfolded state of the hairpin under native conditions, revealing the presence of extremely slow dynamics on the microsecond time scale in the unfolded state and coexistence of structures with partial pairing of the tails in the folded state. Comparing model peptides that mimic the turn sequence, we found that both ion pairing and hydrogen bonding due to the threonine side chain contribute to the propensity of turn formation but not to the much slower dynamics of the hydrophobic core formation. Interestingly, the dynamics of the turn region in isolation are significantly faster than the dynamics measured for the unfolded state of the complete hairpin, suggesting that non-native hydrophobic contacts slow down the reconfiguration dynamics of the unfolded state. Overall, the information extracted from these experiments provides kinetic limits on interconversions among conformational populations, hence enabling a simplified multistate free-energy landscape for the GB1 hairpin to be drawn.

Assessment and optimization of collective variables for protein conformational landscape: GB1 β-hairpin as a case study

Collective variables (CVs), when chosen judiciously, can play an important role in recognizing rate-limiting processes and rare events in any biomolecular systems. However, high dimensionality and inherent complexities associated with such biochemical systems render the identification of an optimal CV a challenging task, which in turn precludes the elucidation of an underlying conformational landscape in sufficient details. In this context, a relevant model system is presented by a 16-residue β-hairpin of GB1 protein. Despite being the target of numerous theoretical and computational studies for understanding the protein folding, the set of CVs optimally characterizing the conformational landscape of the β-hairpin of GB1 protein has remained elusive, resulting in a lack of consensus on its folding mechanism. Here we address this by proposing a pair of optimal CVs which can resolve the underlying free energy landscape of the GB1 hairpin quite efficiently. Expressed as a linear combination of a number of traditional CVs, the optimal CV for this system is derived by employing the recently introduced time-structured independent component analysis approach on a large number of independent unbiased simulations. By projecting the replica-exchange simulated trajectories along these pair of optimized CVs, the resulting free energy landscape of this system is able to resolve four distinct well-separated metastable states encompassing the extensive ensembles of folded, unfolded, and molten globule states. Importantly, the optimized CVs were found to be capable of automatically recovering a novel partial helical state of this protein, without needing to explicitly invoke helicity as a constituent CV. Furthermore, a quantitative sensitivity analysis of each constituent in the optimized CV provided key insights on the relative contributions of the constituent CVs in the overall free energy landscapes. Finally, the kinetic pathways connecting these metastable states, constructed using a Markov state model, provide an optimum description of the underlying folding mechanism of the peptide. Taken together, this work offers a quantitatively robust approach toward comprehensive mapping of the underlying folding landscape of a quintessential model system along its optimized CV.

Microscopic interpretation of folding ϕ-values using the transition path ensemble

All-atom molecular dynamics simulations now allow us to create movies of proteins folding and unfolding. However, it is difficult to assess the accuracy of the folding mechanisms observed because experiments cannot yet directly resolve events occurring along the transition paths between unfolded and folded states. Protein folding ϕ-values provide residue-resolved information about folding mechanisms by comparing the effects of mutations on folding rates and stability, but determining ϕ-values by separately simulating mutant proteins would be computationally demanding and prone to large statistical errors. Here we use transition path theory to develop a method for computing ϕ-values directly from the transition path ensemble, without the need for additional simulations. This path-based approach uses the full transition path information available from equilibrium folding and unfolding trajectories, or from transition path sampling, and does not require identification of folding transition states. Applying our approach to a set of simulations of 10 small proteins by Shaw and coworkers [Lindorff-Larsen K, Piana S, Dror RO, Shaw DE (2011) Science 334(6055):517-520; Piana S, Lindorff-Larsen K, Shaw DE (2011) Biophys J100(9):L47-L49; and Piana S, Lindorff-Larsen K, Shaw DE (2013) Proc Natl Acad Sci USA 110(15):5915-5920], we find good agreement with experiments in most cases where data are available. We can further resolve the contributions to fractional ϕ-values coming from partial contact formation versus transition path heterogeneity. Although in some cases, there is substantial heterogeneity of folding mechanism, in others, such as Ubiquitin, the mechanism is strongly conserved.

Coarse kMC-based replica exchange algorithms for the accelerated simulation of protein folding in explicit solvent

In this paper, we present a coarse replica exchange molecular dynamics (REMD) approach, based on kinetic Monte Carlo (kMC).

Transient β-hairpin formation in α-synuclein monomer revealed by coarse-grained molecular dynamics simulation

Parkinson's disease, originating from the intrinsically disordered peptide α-synuclein, is a common neurodegenerative disorder that affects more than 5% of the population above age 85. It remains unclear how α-synuclein monomers undergo conformational changes leading to aggregation and formation of fibrils characteristic for the disease. In the present study, we perform molecular dynamics simulations (over 180 μs in aggregated time) using a hybrid-resolution model, Proteins with Atomic details in Coarse-grained Environment (PACE), to characterize in atomic detail structural ensembles of wild type and mutant monomeric α-synuclein in aqueous solution. The simulations reproduce structural properties of α-synuclein characterized in experiments, such as secondary structure content, long-range contacts, chemical shifts, and (3)J(HNHCα )-coupling constants. Most notably, the simulations reveal that a short fragment encompassing region 38-53, adjacent to the non-amyloid-β component region, exhibits a high probability of forming a β-hairpin; this fragment, when isolated from the remainder of α-synuclein, fluctuates frequently into its β-hairpin conformation. Two disease-prone mutations, namely, A30P and A53T, significantly accelerate the formation of a β-hairpin in the stated fragment. We conclude that the formation of a β-hairpin in region 38-53 is a key event during α-synuclein aggregation. We predict further that the G47V mutation impedes the formation of a turn in the β-hairpin and slows down β-hairpin formation, thereby retarding α-synuclein aggregation.

Identification of Mutational Hot Spots for Substrate Diffusion: Application to Myoglobin

The pathways by which small molecules (substrates or inhibitors) access active sites are a key aspect of the function of enzymes and other proteins. A key problem in designing or altering such proteins is to identify sites for mutation that will have the desired effect on the substrate transport properties. While specific access channels have been invoked in the past, molecular simulations suggest that multiple routes are possible, complicating the analysis. This complexity, however, can be captured by a Markov State Model (MSM) of the ligand diffusion process. We have developed a sensitivity analysis of the resulting rate matrix, which identifies the locations where mutations should have the largest effect on the diffusive on rate. We apply this method to myoglobin, which is the best characterized example both from experiment and simulation. We validate the approach by translating the sensitivity parameter obtained from this method into the CO binding rates in myoglobin upon mutation, resulting in a semi-quantitative correlation with experiments. The model is further validated against an explicit simulation for one of the experimental mutants.

Dependence of internal friction on folding mechanism

An outstanding challenge in protein folding is understanding the origin of "internal friction" in folding dynamics, experimentally identified from the dependence of folding rates on solvent viscosity. A possible origin suggested by simulation is the crossing of local torsion barriers. However, it was unclear why internal friction varied from protein to protein or for different folding barriers of the same protein. Using all-atom simulations with variable solvent viscosity, in conjunction with transition-path sampling to obtain reaction rates and analysis via Markov state models, we are able to determine the internal friction in the folding of several peptides and miniproteins. In agreement with experiment, we find that the folding events with greatest internal friction are those that mainly involve helix formation, while hairpin formation exhibits little or no evidence of friction. Via a careful analysis of folding transition paths, we show that internal friction arises when torsion angle changes are an important part of the folding mechanism near the folding free energy barrier. These results suggest an explanation for the variation of internal friction effects from protein to protein and across the energy landscape of the same protein.

Balanced Protein-Water Interactions Improve Properties of Disordered Proteins and Non-Specific Protein Association

Some frequently encountered deficiencies in all-atom molecular simulations, such as nonspecific protein-protein interactions being too strong, and unfolded or disordered states being too collapsed, suggest that proteins are insufficiently well solvated in simulations using current state-of-the-art force fields. To address these issues, we make the simplest possible change, by modifying the short-range protein-water pair interactions, and leaving all the water-water and protein-protein parameters unchanged. We find that a modest strengthening of protein-water interactions is sufficient to recover the correct dimensions of intrinsically disordered or unfolded proteins, as determined by direct comparison with small-angle X-ray scattering (SAXS) and Förster resonance energy transfer (FRET) data. The modification also results in more realistic protein-protein affinities, and average solvation free energies of model compounds which are more consistent with experiment. Most importantly, we show that this scaling is small enough not to affect adversely the stability of the folded state, with only a modest effect on the stability of model peptides forming α-helix and β-sheet structures. The proposed adjustment opens the way to more accurate atomistic simulations of proteins, particularly for intrinsically disordered proteins, protein-protein association, and crowded cellular environments.

Detailed microscopic unfolding pathways of an α-helix and a β-hairpin: direct observation and molecular dynamics

We present a combined experimental and computational study of unfolding pathways of a model 21-residue α-helical heteropeptide (W1H5-21) and a 16-residue β-hairpin (GB41-56). Experimentally, we measured fluorescence energy transfer efficiency as a function of temperature, employing natural tryptophans as donors and dansylated lysines as acceptors. Secondary structural analysis was performed with circular dichroism and Fourier transform infrared spectroscopy. Our studies present markedly different unfolding pathways of the two elementary secondary structural elements. During thermal denaturation, the helical peptide exhibits an initial decrease in length, followed by an increase, while the hairpin undergoes a systematic increase in length. In the complementary computational part of the project, we performed microsecond length replica-exchange molecular dynamics simulations of the peptides in explicit solvent, yielding a detailed microscopic picture of the unfolding processes. For the α-helical peptide, we found a large heterogeneous population of intermediates that are primarily frayed single helices or helix-turn-helix motifs. Unfolding starts at the termini and proceeds through a stable helical region in the interior of the peptide but shifted off-center toward the C-terminus. The simulations explain the experimentally observed non-monotonic variation of helix length with temperature as due primarily to the presence of frayed-end single-helix intermediate structures. For the β-hairpin peptide, our simulations indicate that folding is initiated at the turn, followed by formation of the hairpin in zipper-like fashion, with Cα···Cα contacts propagating from the turn to termini and hairpin hydrogen bonds forming in parallel with these contacts. In the early stages of hairpin formation, the hydrophobic side-chain contacts are only partly populated. Intermediate structures with low numbers of β-hairpin hydrogen bonds have very low populations. This is in accord with the "broken zipper" model of Scheraga. The monotonic increase in length with temperature may be explained by the zipper-like breaking of the hairpin hydrogen bonds and backbone contacts.

Markov state models of biomolecular conformational dynamics

It has recently become practical to construct Markov state models (MSMs) that reproduce the long-time statistical conformational dynamics of biomolecules using data from molecular dynamics simulations. MSMs can predict both stationary and kinetic quantities on long timescales (e.g. milliseconds) using a set of atomistic molecular dynamics simulations that are individually much shorter, thus addressing the well-known sampling problem in molecular dynamics simulation. In addition to providing predictive quantitative models, MSMs greatly facilitate both the extraction of insight into biomolecular mechanism (such as folding and functional dynamics) and quantitative comparison with single-molecule and ensemble kinetics experiments. A variety of methodological advances and software packages now bring the construction of these models closer to routine practice. Here, we review recent progress in this field, considering theoretical and methodological advances, new software tools, and recent applications of these approaches in several domains of biochemistry and biophysics, commenting on remaining challenges.

Energy landscape and dynamics of proteins: an exact analysis of a simplified lattice model

We present the results of exact numerical studies of the energy landscape and the dynamics of a 12-monomer chain with contact interactions encoding the ground state on a square lattice. In spite of its simplicity, the model is shown to exhibit behavior at odds with the standard picture of proteins.

How long does it take to equilibrate the unfolded state of a protein?

How long does it take to equilibrate the unfolded state of a protein? The answer to this question has important implications for our understanding of why many small proteins fold with two state kinetics. When the equilibration within the unfolded state U is much faster than the folding, the folding kinetics will be two state even if there are many folding pathways with different barriers. Yet the mean first passage times (MFPTs) between different regions of the unfolded state can be much longer than the folding time. This seems to imply that the equilibration within U is much slower than the folding. In this communication we resolve this paradox. We present a formula for estimating the time to equilibrate the unfolded state of a protein. We also present a formula for the MFPT to any state within U, which is proportional to the average lifetime of that state divided by the state population. This relation is valid when the equilibration within U is very fast as compared with folding as it often is for small proteins. To illustrate the concepts, we apply the formulas to estimate the time to equilibrate the unfolded state of Trp-cage and MFPTs within the unfolded state based on a Markov State Model using an ultra-long 208 microsecond trajectory of the miniprotein to parameterize the model. The time to equilibrate the unfolded state of Trp-cage is ∼100 ns while the typical MFPTs within U are tens of microseconds or longer.

Importance of the ion-pair interactions in the OPEP coarse-grained force field: parametrization and validation

We have derived new effective interactions that improve the description of ion-pairs in the OPEP coarse-grained force field without introducing explicit electrostatic terms. The iterative Boltzmann inversion method was used to extract these potentials from all atom simulations by targeting the radial distribution function of the distance between the center of mass of the side-chains. The new potentials have been tested on several systems that differ in structural properties, thermodynamic stabilities and number of ion-pairs. Our modeling, by refining the packing of the charged amino-acids, impacts the stability of secondary structure motifs and the population of intermediate states during temperature folding/unfolding; it also improves the aggregation propensity of peptides. The new version of the OPEP force field has the potentiality to describe more realistically a large spectrum of situations where salt-bridges are key interactions.