Using massively parallel simulation and Markovian models to study protein folding: examining the dynamics of the villin headpiece

We report on the use of large-scale distributed computing simulation and novel analysis techniques for examining the dynamics of a small protein. Matters addressed include folding rate, very long time scale kinetics, ensemble properties, and interaction with water. The target system for the study, the villin headpiece, has been of great interest to experimentalists and theorists both. Sampling totaled nearly 500 mus-the most extensive published to date for a system of villin's size in explicit solvent with all atom detail-and was in the form of tens of thousands of independent molecular dynamics trajectories, each several tens of nanoseconds in length. We report on kinetics sensitivity analyses that, using a set of short simulations, probed the role of water in villin's folding and sensitivity to the simulation's electrostatics treatment. By constructing Markovian state models (MSMs) from the collected data, we were able to propagate dynamics to times far beyond those directly simulated and to rapidly compute mean first passage times, long time kinetics (tens of microseconds), and evolution of ensemble property distributions over long times, otherwise currently impossible. We also tested our MSM by using it to predict the structure of villin de novo.     

Markov state models based on milestoning

Markov state models (MSMs) have become the tool of choice to analyze large amounts of molecular dynamics data by approximating them as a Markov jump process between suitably predefined states. Here we investigate "Core Set MSMs," a new type of MSMs that build on metastable core sets acting as milestones for tracing the rare event kinetics. We present a thorough analysis of Core Set MSMs based on the existing milestoning framework, Bayesian estimation methods and Transition Path Theory (TPT). We show that Core Set MSMs can be used to extract phenomenological rate constants between the metastable sets of the system and to approximate the evolution of certain key observables. The performance of Core Set MSMs in comparison to standard MSMs is analyzed and illustrated on a toy example and in the context of the torsion angle dynamics of alanine dipeptide.     

Probing molecular kinetics with Markov models: metastable states, transition pathways and spectroscopic observables

Markov (state) models (MSMs) have attracted a lot of interest recently as they (1) can probe long-term molecular kinetics based on short-time simulations, (2) offer a way to analyze great amounts of simulation data with relatively little subjectivity of the analyst, (3) provide insight into microscopic quantities such as the ensemble of transition pathways, and (4) allow simulation data to be reconciled with measurement data in a rigorous and explicit way. Here we sketch our current perspective of Markov models and explain in short their theoretical basis and assumptions. We describe transition path theory which allows the entire ensemble of protein folding pathways to be investigated and that combines naturally with Markov models. Experimental observations can be naturally linked to Markov models with the dynamical fingerprint theory, by which experimentally observable timescales can be equipped with an understanding of the structural rearrangement processes that take place at these timescales. The concepts of this paper are illustrated by a simple kinetic model of protein folding.     

Markov state model reveals folding and functional dynamics in ultra-long MD trajectories

Two strategies have been recently employed to push molecular simulation to long, biologically relevant time scales: projection-based analysis of results from specialized hardware producing a small number of ultralong trajectories and the statistical interpretation of massive parallel sampling performed with Markov state models (MSMs). Here, we assess the MSM as an analysis method by constructing a Markov model from ultralong trajectories, specifically two previously reported 100 μs trajectories of the FiP35 WW domain (Shaw, D. E. Science 2010, 330, 341-346). We find that the MSM approach yields novel insights. It discovers new statistically significant folding pathways, in which either beta-hairpin of the WW domain can form first. The rates of this process approach experimental values in a direct quantitative comparison (time scales of 5.0 μs and 100 ns), within a factor of ～2. Finally, the hub-like topology of the MSM and identification of a holo conformation predicts how WW domains may function through a conformational selection mechanism.     

Investigation on the folding mode of a polymer chain in a spiral crystal by single molecule force spectroscopy

Investigation on the folding mode of a single polymer chain in its crystal is significant to the understanding of the mechanism of the fundamental crystallization as well as the engineering of new polymer crystal-based materials. Herein, we use the combined techniques of atomic force microscopy （AFM） imaging and force spectroscopy to pull a single polyethylene oxide （PEO） chain out of its spiral crystal in amyl acetate. From these data, the folding mode of polymer chains in the spiral crystal has been reconstructed. We find that the stems tilt in the typical flat area, leading to the decrease in the apparent lamellar height. While in the area of screw dislocation, the lamellar height gradually increases in the range of several nanometers. These results indicate that the combined techniques present a novel tool to directly unravel the chain folding mode of spiral crystals at single-molecule level.     

Machine learning algorithms for predicting protein folding rates and stability of mutant proteins: comparison with statistical methods

Machine learning algorithms have wide range of applications in bioinformatics and computational biology such as prediction of protein secondary structures, solvent accessibility, binding site residues in protein complexes, protein folding rates, stability of mutant proteins, and discrimination of proteins based on their structure and function. In this work, we focus on two aspects of predictions: (i) protein folding rates and (ii) stability of proteins upon mutations. We briefly introduce the concepts of protein folding rates and stability along with available databases, features for prediction methods and measures for prediction performance. Subsequently, the development of structure based parameters and their relationship with protein folding rates will be outlined. The structure based parameters are helpful to understand the physical basis for protein folding and stability. Further, basic principles of major machine learning techniques will be mentioned and their applications for predicting protein folding rates and stability of mutant proteins will be illustrated. The machine learning techniques could achieve the highest accuracy of predicting protein folding rates and stability. In essence, statistical methods and machine learning algorithms are complimenting each other for understanding and predicting protein folding rates and the stability of protein mutants. The available online resources on protein folding rates and stability will be listed.     

Markov models of molecular kinetics: generation and validation

Markov state models of molecular kinetics (MSMs), in which the long-time statistical dynamics of a molecule is approximated by a Markov chain on a discrete partition of configuration space, have seen widespread use in recent years. This approach has many appealing characteristics compared to straightforward molecular dynamics simulation and analysis, including the potential to mitigate the sampling problem by extracting long-time kinetic information from short trajectories and the ability to straightforwardly calculate expectation values and statistical uncertainties of various stationary and dynamical molecular observables. In this paper, we summarize the current state of the art in generation and validation of MSMs and give some important new results. We describe an upper bound for the approximation error made by modeling molecular dynamics with a MSM and we show that this error can be made arbitrarily small with surprisingly little effort. In contrast to previous practice, it becomes clear that the best MSM is not obtained by the most metastable discretization, but the MSM can be much improved if non-metastable states are introduced near the transition states. Moreover, we show that it is not necessary to resolve all slow processes by the state space partitioning, but individual dynamical processes of interest can be resolved separately. We also present an efficient estimator for reversible transition matrices and a robust test to validate that a MSM reproduces the kinetics of the molecular dynamics data.     

Calculation of the distribution of eigenvalues and eigenvectors in Markovian state models for molecular dynamics

Markovian state models (MSMs) are a convenient and efficient means to compactly describe the kinetics of a molecular system as well as a formalism for using many short simulations to predict long time scale behavior. Building a MSM consists of grouping the conformations into states and estimating the transition probabilities between these states. In a previous paper, we described an efficient method for calculating the uncertainty due to finite sampling in the mean first passage time between two states. In this paper, we extend the uncertainty analysis to derive similar closed-form solutions for the distributions of the eigenvalues and eigenvectors of the transition matrix, quantities that have numerous applications when using the model. We demonstrate the accuracy of the distributions on a six-state model of the terminally blocked alanine peptide. We also show how to significantly reduce the total number of simulations necessary to build a model with a given precision using these uncertainty estimates for the blocked alanine system and for a 2454-state MSM for the dynamics of the villin headpiece.     

A semi-analytical description of protein folding that incorporates detailed geometrical information

Much has been done to study the interplay between geometric and energetic effects on the protein folding energy landscape. Numerical techniques such as molecular dynamics simulations are able to maintain a precise geometrical representation of the protein. Analytical approaches, however, often focus on the energetic aspects of folding, including geometrical information only in an average way. Here, we investigate a semi-analytical expression of folding that explicitly includes geometrical effects. We consider a Hamiltonian corresponding to a Gaussian filament with structure-based interactions. The model captures local features of protein folding often averaged over by mean-field theories, for example, loop contact formation and excluded volume. We explore the thermodynamics and folding mechanisms of beta-hairpin and alpha-helical structures as functions of temperature and Q, the fraction of native contacts formed. Excluded volume is shown to be an important component of a protein Hamiltonian, since it both dominates the cooperativity of the folding transition and alters folding mechanisms. Understanding geometrical effects in analytical formulae will help illuminate the consequences of the approximations required for the study of larger proteins.     

Assessment of protein folding potentials with an evolutionary method

Many different protein folding potentials have been developed in the last decades, based upon knowledge of experimentally determined protein structures. Decoy-based techniques are frequently used to assess these force fields, but other methods can explore different features in the performance of the interaction schemes, thus helping in their evaluation. Here, we propose an evolutionary strategy to efficiently assess folding potentials. We apply it to three potentials with different characteristics, taken from the bibliography. A search for minimum energy protein topologies, treated as arrangements of rigid protein fragments, is performed. The method, applied to a set of helix bundle proteins, shows the different behavior of the studied potentials, providing a reasonably fast tool to evaluate their advantages and limitations.     

Effects of long-range electrostatic forces on simulated protein folding kinetics

Molecular dynamics simulations are a useful tool for characterizing protein folding pathways. There are several methods of treating electrostatic forces in these simulations with varying degrees of physical fidelity and computational efficiency. In this article, we compare the reaction field (RF) algorithm, particle-mesh Ewald (PME), and tapered cutoffs with increasing cutoff radii to address the impact of the electrostatics method employed on the folding kinetics. We quantitatively compare different methods by a correlation of quantitative measures of protein folding kinetics. The results of these comparisons show that for protein folding kinetics, the RF algorithm can quantitatively reproduce the kinetics of the more costly PME algorithm. These results not only assist the selection of appropriate algorithms for future simulations, but also give insight on the role that long-range electrostatic forces have in protein folding.     

A Single‐Molecule Perspective on the Role of Solvent Hydrogen Bonds in Protein Folding and Chemical Reactions

We present an array of force spectroscopy experiments that aim to identify the role of solvent hydrogen bonds in protein folding and chemical reactions at the single‐molecule level. In our experiments we control the strength of hydrogen bonds in the solvent environment by substituting water (H₂O) with deuterium oxide (D₂O). Using a combination of force protocols, we demonstrate that protein unfolding, protein collapse, protein folding and a chemical reaction are affected in different ways by substituting H₂O with D₂O. We find that D₂O molecules form an integral part of the unfolding transition structure of the immunoglobulin module of human cardiac titin, I27. Strikingly, we find that D₂O is a worse solvent than H₂O for the protein I27, in direct contrast with the behaviour of simple hydrocarbons. We measure the effect of substituting H₂O with D₂O on the force dependent rate of reduction of a disulphide bond engineered within a single protein. Altogether, these experiments provide new information on the nature of the underlying interactions in protein folding and chemical reactions and demonstrate the power of single‐molecule techniques to identify the changes induced by a small change in hydrogen bond strength.     

Foldamer simulations: novel computational methods and applications to poly-phenylacetylene oligomers

We apply several methods to probe the ensemble kinetic and structural properties of a model system of poly-phenylacetylene (pPA) oligomer folding trajectories. The kinetic methods employed included a brute force accounting of conformations, a Markovian state matrix method, and a nonlinear least squares fit to a minimalist kinetic model used to extract the folding time. Each method gave similar measures for the folding time of the 12-mer chain, calculated to be on the order of 7 ns for the complete folding of the chain from an extended conformation. Utilizing both a linear and a nonlinear scaling relationship between the viscosity and the folding time to correct for a low simulation viscosity, we obtain an upper and a lower bound for the approximate folding time within the range 70 ns相似文献   

Effects of DNA secondary structure on oligonucleotide probe binding efficiency

Secondary structure motifs in nucleic acid probes generally impair intended hybridization reactions and so efforts to predict and avoid such structures are commonly employed in probe design schemes. Another key facet of probe design that has received much less attention, however, is that secondary structure at targeted probe binding site regions may also impair hybridization. Thus, evaluation of both probe and target site secondary structures together should improve hybridization prediction and design effectiveness. Several challenges confound this goal, including imperfect empirical rules and parameters underlying predictions and the fact that folding algorithms scale poorly with respect to sequence length. Here, we attempt to quantify the consequences of target site structure on predicted hybridization using sequences sampled from the human genome. We also provide a methodology for choosing a reasonable “window size” around target sites that is as small as possible without compromising folding algorithm prediction accuracy.     

Folding of human telomerase RNA pseudoknot using ion-jump and temperature-quench simulations

Globally RNA folding occurs in multiple stages involving chain compaction and subsequent rearrangement by a number of parallel routes to the folded state. However, the sequence-dependent details of the folding pathways and the link between collapse and folding are poorly understood. To obtain a comprehensive picture of the thermodynamics and folding kinetics we used molecular simulations of coarse-grained model of a pseudoknot found in the conserved core domain of the human telomerase (hTR) by varying both temperature (T) and ion concentration (C). The phase diagram in the [T,C] plane shows that the boundary separating the folded and unfolded state for the finite 47-nucleotide system is relatively sharp, implying that from a thermodynamic perspective hTR behaves as an apparent two-state system. However, the folding kinetics following single C-jump or T-quench is complicated, involving multiple channels to the native state. Although globally folding kinetics triggered by T-quench and C-jump are similar, the kinetics of chain compaction are vastly different, which reflects the role of initial conditions in directing folding and collapse. Remarkably, even after substantial reduction in the overall size of hTR, the ensemble of compact conformations are far from being nativelike, suggesting that the search for the folded state occurs among the ensemble of low-energy fluidlike globules. The rate of unfolding, which occurs in a single step, is faster upon C-decrease compared to a jump in temperature. To identify "hidden" states that are visited during the folding process we performed simulations by periodically interrupting the approach to the folded state by lowering C. These simulations show that hTR reaches the folded state through a small number of connected clusters that are repeatedly visited during the pulse sequence in which the folding or unfolding is interrupted. The results from interrupted folding simulations, which are in accord with non-equilibrium single-molecule folding of a large ribozyme, show that multiple probes are needed to reveal the invisible states that are sampled by RNA as it folds. Although we have illustrated the complexity of RNA folding using hTR as a case study, general arguments and qualitative comparisons to time-resolved scattering experiments on Azoarcus group I ribozyme and single-molecule non-equilibrium periodic ion-jump experiments establish the generality of our findings.     

Forced Folding of a Disordered Protein Accesses an Alternative Folding Landscape

Intrinsically disordered proteins (IDPs) are involved in diverse cellular functions. Many IDPs can interact with multiple binding partners, resulting in their folding into alternative ligand‐specific functional structures. For such multi‐structural IDPs, a key question is whether these multiple structures are fully encoded in the protein sequence, as is the case in many globular proteins. To answer this question, here we employed a combination of single‐molecule and ensemble techniques to compare ligand‐induced and osmolyte‐forced folding of α‐synuclein. Our results reveal context‐dependent modulation of the protein′s folding landscape, suggesting that the codes for the protein′s native folds are partially encoded in its primary sequence, and are completed only upon interaction with binding partners. Our findings suggest a critical role for cellular interactions in expanding the repertoire of folds and functions available to disordered proteins.     

Molecular latent space simulators

Small integration time steps limit molecular dynamics (MD) simulations to millisecond time scales. Markov state models (MSMs) and equation-free approaches learn low-dimensional kinetic models from MD simulation data by performing configurational or dynamical coarse-graining of the state space. The learned kinetic models enable the efficient generation of dynamical trajectories over vastly longer time scales than are accessible by MD, but the discretization of configurational space and/or absence of a means to reconstruct molecular configurations precludes the generation of continuous atomistic molecular trajectories. We propose latent space simulators (LSS) to learn kinetic models for continuous atomistic simulation trajectories by training three deep learning networks to (i) learn the slow collective variables of the molecular system, (ii) propagate the system dynamics within this slow latent space, and (iii) generatively reconstruct molecular configurations. We demonstrate the approach in an application to Trp-cage miniprotein to produce novel ultra-long synthetic folding trajectories that accurately reproduce atomistic molecular structure, thermodynamics, and kinetics at six orders of magnitude lower cost than MD. The dramatically lower cost of trajectory generation enables greatly improved sampling and greatly reduced statistical uncertainties in estimated thermodynamic averages and kinetic rates.

Latent space simulators learn kinetic models for atomistic simulations and generate novel trajectories at six orders of magnitude lower cost.     

Solvent viscosity dependence of the folding rate of a small protein: distributed computing study

By using distributed computing techniques and a supercluster of more than 20,000 processors we simulated folding of a 20-residue Trp Cage miniprotein in atomistic detail with implicit GB/SA solvent at a variety of solvent viscosities (gamma). This allowed us to analyze the dependence of folding rates on viscosity. In particular, we focused on the low-viscosity regime (values below the viscosity of water). In accordance with Kramers' theory, we observe approximately linear dependence of the folding rate on 1/gamma for values from 1-10(-1)x that of water viscosity. However, for the regime between 10(-4)-10(-1)x that of water viscosity we observe power-law dependence of the form k approximately gamma(-1/5). These results suggest that estimating folding rates from molecular simulations run at low viscosity under the assumption of linear dependence of rate on inverse viscosity may lead to erroneous results.     

Propensity for local folding induced by the urea fragment in short-chain oligomers

Examination of local folding and H-bonding patterns in model compounds can be extremely informative to gain insight into the propensity of longer-chain oligomers to adopt specific folding patterns (i.e. foldamers) based on remote interactions. Using a combination of experimental techniques (i.e. X-ray diffraction, FT-IR absorption and NMR spectroscopy) and theoretical calculations at the density functional theory (DFT) level, we have examined the local folding patterns induced by the urea fragment in short-chain aza analogues of beta- and gamma-amino acid derivatives. We found that the urea-turn, a robust C(8) conformation based on 1<--3 H-bond interaction, is largely populated in model ureidopeptides (I-IV) obtained by replacing the alpha-carbon of a beta-amino acid by a nitrogen. This H-bonding scheme is likely to compete with remote H-bond interactions, thus preventing the formation of secondary structures based on remote intrastrand interactions in longer oligomers. In related oligomers obtained by the addition of a methylene in the main chain (V-VIII), nearest-neighbour H-bonded interactions are unfavourable i.e. the corresponding C9 folding pattern is hardly populated. In this series, folding based on remote intrastrand interactions becomes possible for longer oligomers. We present spectroscopic evidence that tetraurea VIII is likely to be the smallest unit capable of reproducing the H-bonded motif found in 2.5-helical N,N'-linked oligoureas.