MSMBuilder: Statistical Models for Biomolecular Dynamics

Combined approaches from physics, statistics, and computer science for ab initio protein structure prediction: ex unitate vires (unity is strength)?

Connecting the dots among the amino acid sequence of a protein, its structure, and its function remains a central theme in molecular biology, as it would have many applications in the treatment of illnesses related to misfolding or protein instability. As a result of high-throughput sequencing methods, biologists currently live in a protein sequence-rich world. However, our knowledge of protein structure based on experimental data remains comparatively limited. As a consequence, protein structure prediction has established itself as a very active field of research to fill in this gap. This field, once thought to be reserved for theoretical biophysicists, is constantly reinventing itself, borrowing ideas informed by an ever-increasing assembly of scientific domains, from biology, chemistry, (statistical) physics, mathematics, computer science, statistics, bioinformatics, and more recently data sciences. We review the recent progress arising from this integration of knowledge, from the development of specific computer architecture to allow for longer timescales in physics-based simulations of protein folding to the recent advances in predicting contacts in proteins based on detection of coevolution using very large data sets of aligned protein sequences.

Variational encoding of complex dynamics

Often the analysis of time-dependent chemical and biophysical systems produces high-dimensional time-series data for which it can be difficult to interpret which individual features are most salient. While recent work from our group and others has demonstrated the utility of time-lagged covariate models to study such systems, linearity assumptions can limit the compression of inherently nonlinear dynamics into just a few characteristic components. Recent work in the field of deep learning has led to the development of the variational autoencoder (VAE), which is able to compress complex datasets into simpler manifolds. We present the use of a time-lagged VAE, or variational dynamics encoder (VDE), to reduce complex, nonlinear processes to a single embedding with high fidelity to the underlying dynamics. We demonstrate how the VDE is able to capture nontrivial dynamics in a variety of examples, including Brownian dynamics and atomistic protein folding. Additionally, we demonstrate a method for analyzing the VDE model, inspired by saliency mapping, to determine what features are selected by the VDE model to describe dynamics. The VDE presents an important step in applying techniques from deep learning to more accurately model and interpret complex biophysics.

MSM/RD: Coupling Markov state models of molecular kinetics with reaction-diffusion simulations

Molecular dynamics (MD) simulations can model the interactions between macromolecules with high spatiotemporal resolution but at a high computational cost. By combining high-throughput MD with Markov state models (MSMs), it is now possible to obtain long time-scale behavior of small to intermediate biomolecules and complexes. To model the interactions of many molecules at large length scales, particle-based reaction-diffusion (RD) simulations are more suitable but lack molecular detail. Thus, coupling MSMs and RD simulations (MSM/RD) would be highly desirable, as they could efficiently produce simulations at large time and length scales, while still conserving the characteristic features of the interactions observed at atomic detail. While such a coupling seems straightforward, fundamental questions are still open: Which definition of MSM states is suitable? Which protocol to merge and split RD particles in an association/dissociation reaction will conserve the correct bimolecular kinetics and thermodynamics? In this paper, we make the first step toward MSM/RD by laying out a general theory of coupling and proposing a first implementation for association/dissociation of a protein with a small ligand (A + B ⇌ C). Applications on a toy model and CO diffusion into the heme cavity of myoglobin are reported.

Examining a Thermodynamic Order Parameter of Protein Folding

Dimensionality reduction with a suitable choice of order parameters or reaction coordinates is commonly used for analyzing high-dimensional time-series data generated by atomistic biomolecular simulations. So far, geometric order parameters, such as the root mean square deviation, fraction of native amino acid contacts, and collective coordinates that best characterize rare or large conformational transitions, have been prevailing in protein folding studies. Here, we show that the solvent-averaged effective energy, which is a thermodynamic quantity but unambiguously defined for individual protein conformations, serves as a good order parameter of protein folding. This is illustrated through the application to the folding-unfolding simulation trajectory of villin headpiece subdomain. We rationalize the suitability of the effective energy as an order parameter by the funneledness of the underlying protein free energy landscape. We also demonstrate that an improved conformational space discretization is achieved by incorporating the effective energy. The most distinctive feature of this thermodynamic order parameter is that it works in pointing to near-native folded structures even when the knowledge of the native structure is lacking, and the use of the effective energy will also find applications in combination with methods of protein structure prediction.

Predicting the Conformational Variability of Abl Tyrosine Kinase using Molecular Dynamics Simulations and Markov State Models

Understanding protein conformational variability remains a challenge in drug discovery. The issue arises in protein kinases, whose multiple conformational states can affect the binding of small-molecule inhibitors. To overcome this challenge, we propose a comprehensive computational framework based on Markov state models (MSMs). Our framework integrates the information from explicit-solvent molecular dynamics simulations to accurately rank-order the accessible conformational variants of a target protein. We tested the methodology using Abl kinase with a reference and blind-test set. Only half of the Abl conformational variants discovered by our approach are present in the disclosed X-ray structures. The approach successfully identified a protein conformational state not previously observed in public structures but evident in a retrospective analysis of Lilly in-house structures: the X-ray structure of Abl with WHI-P154. Using a MSM-derived model, the free energy landscape and kinetic profile of Abl was analyzed in detail highlighting opportunities for targeting the unique metastable states.

Linking time-series of single-molecule experiments with molecular dynamics simulations by machine learning

Single-molecule experiments and molecular dynamics (MD) simulations are indispensable tools for investigating protein conformational dynamics. The former provide time-series data, such as donor-acceptor distances, whereas the latter give atomistic information, although this information is often biased by model parameters. Here, we devise a machine-learning method to combine the complementary information from the two approaches and construct a consistent model of conformational dynamics. It is applied to the folding dynamics of the formin-binding protein WW domain. MD simulations over 400 μs led to an initial Markov state model (MSM), which was then "refined" using single-molecule Förster resonance energy transfer (FRET) data through hidden Markov modeling. The refined or data-assimilated MSM reproduces the FRET data and features hairpin one in the transition-state ensemble, consistent with mutation experiments. The folding pathway in the data-assimilated MSM suggests interplay between hydrophobic contacts and turn formation. Our method provides a general framework for investigating conformational transitions in other proteins.

Recognition of protein allosteric states and residues: Machine learning approaches

Allostery is a process by which proteins transmit the effect of perturbation at one site to a distal functional site upon certain perturbation. As an intrinsically global effect of protein dynamics, it is difficult to associate protein allostery with individual residues, hindering effective selection of key residues for mutagenesis studies. The machine learning models including decision tree (DT) and artificial neural network (ANN) models were applied to develop classification model for a cell signaling allosteric protein with two states showing extremely similar tertiary structures in both crystallographic structures and molecular dynamics simulations. Both DT and ANN models were developed with 75% and 80% of predicting accuracy, respectively. Good agreement between machine learning models and previous experimental as well as computational studies of the same protein validates this approach as an alternative way to analyze protein dynamics simulations and allostery. In addition, the difference of distributions of key features in two allosteric states also underlies the population shift hypothesis of dynamics-driven allostery model. © 2018 Wiley Periodicals, Inc.

Ensemble Docking in Drug Discovery

Ensemble docking corresponds to the generation of an "ensemble" of drug target conformations in computational structure-based drug discovery, often obtained by using molecular dynamics simulation, that is used in docking candidate ligands. This approach is now well established in the field of early-stage drug discovery. This review gives a historical account of the development of ensemble docking and discusses some pertinent methodological advances in conformational sampling.

How structural elements evolving from bacterial to human SLC6 transporters enabled new functional properties

Background

Much of the structure-based mechanistic understandings of the function of SLC6A neurotransmitter transporters emerged from the study of their bacterial LeuT-fold homologs. It has become evident, however, that structural differences such as the long N- and C-termini of the eukaryotic neurotransmitter transporters are involved in an expanded set of functional properties to the eukaryotic transporters. These functional properties are not shared by the bacterial homologs, which lack the structural elements that appeared later in evolution. However, mechanistic insights into some of the measured functional properties of the eukaryotic transporters that have been suggested to involve these structural elements are sparse or merely descriptive.

Results

To learn how the structural elements added in evolution enable mechanisms of the eukaryotic transporters in ways not shared with their bacterial LeuT-like homologs, we focused on the human dopamine transporter (hDAT) as a prototype. We present the results of a study employing large-scale molecular dynamics simulations and comparative Markov state model analysis of experimentally determined properties of the wild-type and mutant hDAT constructs. These offer a quantitative outline of mechanisms in which a rich spectrum of interactions of the hDAT N-terminus and C-terminus contribute to the regulation of transporter function (e.g., by phosphorylation) and/or to entirely new phenotypes (e.g., reverse uptake (efflux)) that were added in evolution.

Conclusions

The findings are consistent with the proposal that the size of eukaryotic neurotransmitter transporter termini increased during evolution to enable more functions (e.g., efflux) not shared with the bacterial homologs. The mechanistic explanations for the experimental findings about the modulation of function in DAT, the serotonin transporter, and other eukaryotic transporters reveal separate roles for the distal and proximal segments of the much larger N-terminus in eukaryotic transporters compared to the bacterial ones. The involvement of the proximal and distal segments — such as the role of the proximal segment in sustaining transport in phosphatidylinositol 4,5-bisphosphate-depleted membranes and of the distal segment in modulating efflux — may represent an evolutionary adaptation required for the function of eukaryotic transporters expressed in various cell types of the same organism that differ in the lipid composition and protein complement of their membrane environment.

Electronic supplementary material

The online version of this article (10.1186/s12915-018-0495-6) contains supplementary material, which is available to authorized users.

Markov modeling of peptide folding in the presence of protein crowders

We use Markov state models (MSMs) to analyze the dynamics of a β-hairpin-forming peptide in Monte Carlo (MC) simulations with interacting protein crowders, for two different types of crowder proteins [bovine pancreatic trypsin inhibitor (BPTI) and GB1]. In these systems, at the temperature used, the peptide can be folded or unfolded and bound or unbound to crowder molecules. Four or five major free-energy minima can be identified. To estimate the dominant MC relaxation times of the peptide, we build MSMs using a range of different time resolutions or lag times. We show that stable relaxation-time estimates can be obtained from the MSM eigenfunctions through fits to autocorrelation data. The eigenfunctions remain sufficiently accurate to permit stable relaxation-time estimation down to small lag times, at which point simple estimates based on the corresponding eigenvalues have large systematic uncertainties. The presence of the crowders has a stabilizing effect on the peptide, especially with BPTI crowders, which can be attributed to a reduced unfolding rate ku, while the folding rate kf is left largely unchanged.

Molecular mechanisms of RNA polymerase II transcription elongation elucidated by kinetic network models

Transcription elongation cycle (TEC) of RNA polymerase II (Pol II) is a process of adding a nucleoside triphosphate to the growing messenger RNA chain. Due to the long timescale events in Pol II TEC, an advanced computational technique, such as Markov State Model (MSM), is needed to provide atomistic mechanism and reaction rates. The combination of MSM and experimental results can be used to build a kinetic network model (KNM) of the whole TEC. This review provides a brief protocol to build MSM and KNM of the whole TEC, along with the latest findings of MSM and other computational studies of Pol II TEC. Lastly, we offer a perspective on potentially using a sequence dependent KNM to predict genome-wide transcription error.

A Minimum Variance Clustering Approach Produces Robust and Interpretable Coarse-Grained Models

Markov state models (MSMs) are a powerful framework for the analysis of molecular dynamics data sets, such as protein folding simulations, because of their straightforward construction and statistical rigor. The coarse-graining of MSMs into an interpretable number of macrostates is a crucial step for connecting theoretical results with experimental observables. Here we present the minimum variance clustering approach (MVCA) for the coarse-graining of MSMs into macrostate models. The method utilizes agglomerative clustering with Ward's minimum variance objective function, and the similarity of the microstate dynamics is determined using the Jensen-Shannon divergence between the corresponding rows in the MSM transition probability matrix. We first show that MVCA produces intuitive results for a simple tripeptide system and is robust toward long-duration statistical artifacts. MVCA is then applied to two protein folding simulations of the same protein in different force fields to demonstrate that a different number of macrostates is appropriate for each model, revealing a misfolded state present in only one of the simulations. Finally, we show that the same method can be used to analyze a data set containing many MSMs from simulations in different force fields by aggregating them into groups and quantifying their dynamical similarity in the context of force field parameter choices. The minimum variance clustering approach with the Jensen-Shannon divergence provides a powerful tool to group dynamics by similarity, both among model states and among dynamical models themselves.

VAMPnets for deep learning of molecular kinetics

There is an increasing demand for computing the relevant structures, equilibria, and long-timescale kinetics of biomolecular processes, such as protein-drug binding, from high-throughput molecular dynamics simulations. Current methods employ transformation of simulated coordinates into structural features, dimension reduction, clustering the dimension-reduced data, and estimation of a Markov state model or related model of the interconversion rates between molecular structures. This handcrafted approach demands a substantial amount of modeling expertise, as poor decisions at any step will lead to large modeling errors. Here we employ the variational approach for Markov processes (VAMP) to develop a deep learning framework for molecular kinetics using neural networks, dubbed VAMPnets. A VAMPnet encodes the entire mapping from molecular coordinates to Markov states, thus combining the whole data processing pipeline in a single end-to-end framework. Our method performs equally or better than state-of-the-art Markov modeling methods and provides easily interpretable few-state kinetic models.

Extracting kinetic models from high-throughput molecular dynamics (MD) simulations is laborious and prone to human error. Here the authors introduce a deep learning framework that automates construction of Markov state models from MD simulation data.

Hybrid Methods for Modeling Protein Structures Using Molecular Dynamics Simulations and Small-Angle X-Ray Scattering Data

Small-angle X-ray scattering (SAXS) is an efficient experimental tool to measure the overall shape of macromolecular structures in solution. However, due to the low resolution of SAXS data, high-resolution data obtained from X-ray crystallography or NMR and computational methods such as molecular dynamics (MD) simulations are complementary to SAXS data for understanding protein functions based on their structures at atomic resolution. Because MD simulations provide a physicochemically proper structural ensemble for flexible proteins in solution and a precise description of solvent effects, the hybrid analysis of SAXS and MD simulations is a promising method to estimate reasonable solution structures and structural ensembles in solution. Here, we review typical and useful in silico methods for modeling three dimensional protein structures, calculating theoretical SAXS profiles, and analyzing ensemble structures consistent with experimental SAXS profiles. We also review two examples of the hybrid analysis, termed MD-SAXS method in which MD simulations are carried out without any knowledge of experimental SAXS data, and the experimental SAXS data are used only to assess the consistency of the solution model from MD simulations with those observed in experiments. One example is an investigation of the intrinsic dynamics of EcoO109I using the computational method to obtain a theoretical profile from the trajectory of an MD simulation. The other example is a structural investigation of the vitamin D receptor ligand-binding domain using snapshots generated by MD simulations and assessment of the snapshots by experimental SAXS data.

Mechanisms of Lipid Scrambling by the G Protein-Coupled Receptor Opsin

Several class-A G protein-coupled receptor (GPCR) proteins act as constitutive phospholipid scramblases catalyzing the transbilayer translocation of >10,000 phospholipids per second when reconstituted into synthetic vesicles. To address the molecular mechanism by which these proteins facilitate rapid lipid scrambling, we carried out large-scale ensemble atomistic molecular dynamics simulations of the opsin GPCR. We report that, in the process of scrambling, lipid head groups traverse a dynamically revealed hydrophilic pathway in the region between transmembrane helices 6 and 7 of the protein while their hydrophobic tails remain in the bilayer environment. We present quantitative kinetic models of the translocation process based on Markov State Model analysis. As key residues on the lipid translocation pathway are conserved within the class-A GPCR family, our results illuminate unique aspects of GPCR structure and dynamics while providing a rigorous basis for the design of variants of these proteins with defined scramblase activity.

Prediction of New Stabilizing Mutations Based on Mechanistic Insights from Markov State Models

Protein stabilization is fundamental to enzyme function and evolution, yet understanding the determinants of a protein's stability remains a challenge. This is largely due to a shortage of atomically detailed models for the ensemble of relevant protein conformations and their relative populations. For example, the M182T substitution in TEM β-lactamase, an enzyme that confers antibiotic resistance to bacteria, is stabilizing but the precise mechanism remains unclear. Here, we employ Markov state models (MSMs) to uncover how M182T shifts the distribution of different structures that TEM adopts. We find that M182T stabilizes a helix that is a key component of a domain interface. We then predict the effects of other mutations, including a novel stabilizing mutation, and experimentally test our predictions using a combination of stability measurements, crystallography, NMR, and in vivo measurements of bacterial fitness. We expect our insights and methodology to provide a valuable foundation for protein design.

Similar Active Sites and Mechanisms Do Not Lead to Cross-Promiscuity in Organophosphate Hydrolysis: Implications for Biotherapeutic Engineering

Organophosphate hydrolases are proficient catalysts of the breakdown of neurotoxic organophosphates and have great potential as both biotherapeutics for treating acute organophosphate toxicity and as bioremediation agents. However, proficient organophosphatases such as serum paraoxonase 1 (PON1) and the organophosphate-hydrolyzing lactonase SsoPox are unable to hydrolyze bulkyorganophosphates with challenging leaving groups such as diisopropyl fluorophosphate (DFP) or venomous agent X, creating a major challenge for enzyme design. Curiously, despite their mutually exclusive substrate specificities, PON1 and diisopropyl fluorophosphatase (DFPase) have essentially identical active sites and tertiary structures. In the present work, we use empirical valence bond simulations to probe the catalytic mechanism of DFPase as well as temperature, pH, and mutational effects, demonstrating that DFPase and PON1 also likely utilize identical catalytic mechanisms to hydrolyze their respective substrates. However, detailed examination of both static structures and dynamical simulations demonstrates subtle but significant differences in the electrostatic properties and solvent penetration of the two active sites and, most critically, the role of residues that make no direct contact with either substrate in acting as "specificity switches" between the two enzymes. Specifically, we demonstrate that key residues that are structurally and functionally critical for the paraoxonase activity of PON1 prevent it from being able to hydrolyze DFP with its fluoride leaving group. These insights expand our understanding of the drivers of the evolution of divergent substrate specificity in enzymes with identical active sites and guide the future design of organophosphate hydrolases that hydrolyze compounds with challenging leaving groups.

Millisecond dynamics of BTK reveal kinome-wide conformational plasticity within the apo kinase domain

Bruton tyrosine kinase (BTK) is a key enzyme in B-cell development whose improper regulation causes severe immunodeficiency diseases. Design of selective BTK therapeutics would benefit from improved, in-silico structural modeling of the kinase’s solution ensemble. However, this remains challenging due to the immense computational cost of sampling events on biological timescales. In this work, we combine multi-millisecond molecular dynamics (MD) simulations with Markov state models (MSMs) to report on the thermodynamics, kinetics, and accessible states of BTK’s kinase domain. Our conformational landscape links the active state to several inactive states, connected via a structurally diverse intermediate. Our calculations predict a kinome-wide conformational plasticity, and indicate the presence of several new potentially druggable BTK states. We further find that the population of these states and the kinetics of their inter-conversion are modulated by protonation of an aspartate residue, establishing the power of MD & MSMs in predicting effects of chemical perturbations.

Ab initio sampling of transition paths by conditioned Langevin dynamics

We propose a novel stochastic method to generate Brownian paths conditioned to start at an initial point and end at a given final point during a fixed time t_f under a given potential U(x). These paths are sampled with a probability given by the overdamped Langevin dynamics. We show that these paths can be exactly generated by a local stochastic partial differential equation. This equation cannot be solved in general but we present several approximations that are valid either in the low temperature regime or in the presence of barrier crossing. We show that this method warrants the generation of statistically independent transition paths. It is computationally very efficient. We illustrate the method first on two simple potentials, the two-dimensional Mueller potential and the Mexican hat potential, and then on the multi-dimensional problem of conformational transitions in proteins using the "Mixed Elastic Network Model" as a benchmark.

Modeling the mechanism of CLN025 beta-hairpin formation

Beta-hairpins are substructures found in proteins that can lend insight into more complex systems. Furthermore, the folding of beta-hairpins is a valuable test case for benchmarking experimental and theoretical methods. Here, we simulate the folding of CLN025, a miniprotein with a beta-hairpin structure, at its experimental melting temperature using a range of state-of-the-art protein force fields. We construct Markov state models in order to examine the thermodynamics, kinetics, mechanism, and rate-determining step of folding. Mechanistically, we find the folding process is rate-limited by the formation of the turn region hydrogen bonds, which occurs following the downhill hydrophobic collapse of the extended denatured protein. These results are presented in the context of established and contradictory theories of the beta-hairpin folding process. Furthermore, our analysis suggests that the AMBER-FB15 force field, at this temperature, best describes the characteristics of the full experimental CLN025 conformational ensemble, while the AMBER ff99SB-ILDN and CHARMM22* force fields display a tendency to overstabilize the native state.

A Kepler Workflow Tool for Reproducible AMBER GPU Molecular Dynamics

With the drive toward high throughput molecular dynamics (MD) simulations involving ever-greater numbers of simulation replicates run for longer, biologically relevant timescales (microseconds), the need for improved computational methods that facilitate fully automated MD workflows gains more importance. Here we report the development of an automated workflow tool to perform AMBER GPU MD simulations. Our workflow tool capitalizes on the capabilities of the Kepler platform to deliver a flexible, intuitive, and user-friendly environment and the AMBER GPU code for a robust and high-performance simulation engine. Additionally, the workflow tool reduces user input time by automating repetitive processes and facilitates access to GPU clusters, whose high-performance processing power makes simulations of large numerical scale possible. The presented workflow tool facilitates the management and deployment of large sets of MD simulations on heterogeneous computing resources. The workflow tool also performs systematic analysis on the simulation outputs and enhances simulation reproducibility, execution scalability, and MD method development including benchmarking and validation.

Direct Learning Hidden Excited State Interaction Patterns from ab initio Dynamics and Its Implication as Alternative Molecular Mechanism Models

The excited states of polyatomic systems are rather complex, and often exhibit meta-stable dynamical behaviors. Static analysis of reaction pathway often fails to sufficiently characterize excited state motions due to their highly non-equilibrium nature. Here, we proposed a time series guided clustering algorithm to generate most relevant meta-stable patterns directly from ab initio dynamic trajectories. Based on the knowledge of these meta-stable patterns, we suggested an interpolation scheme with only a concrete and finite set of known patterns to accurately predict the ground and excited state properties of the entire dynamics trajectories, namely, the prediction with ensemble models (PEM). As illustrated with the example of sinapic acids, The PEM method does not require any training data beyond the clustering algorithm, and the estimation error for both ground and excited state is very close, which indicates one could predict the ground and excited state molecular properties with similar accuracy. These results may provide us some insights to construct molecular mechanism models with compatible energy terms as traditional force fields.

RNA folding kinetics using Monte Carlo and Gillespie algorithms

RNA secondary structure folding kinetics is known to be important for the biological function of certain processes, such as the hok/sok system in E. coli. Although linear algebra provides an exact computational solution of secondary structure folding kinetics with respect to the Turner energy model for tiny ([Formula: see text]20 nt) RNA sequences, the folding kinetics for larger sequences can only be approximated by binning structures into macrostates in a coarse-grained model, or by repeatedly simulating secondary structure folding with either the Monte Carlo algorithm or the Gillespie algorithm. Here we investigate the relation between the Monte Carlo algorithm and the Gillespie algorithm. We prove that asymptotically, the expected time for a K-step trajectory of the Monte Carlo algorithm is equal to [Formula: see text] times that of the Gillespie algorithm, where [Formula: see text] denotes the Boltzmann expected network degree. If the network is regular (i.e. every node has the same degree), then the mean first passage time (MFPT) computed by the Monte Carlo algorithm is equal to MFPT computed by the Gillespie algorithm multiplied by [Formula: see text]; however, this is not true for non-regular networks. In particular, RNA secondary structure folding kinetics, as computed by the Monte Carlo algorithm, is not equal to the folding kinetics, as computed by the Gillespie algorithm, although the mean first passage times are roughly correlated. Simulation software for RNA secondary structure folding according to the Monte Carlo and Gillespie algorithms is publicly available, as is our software to compute the expected degree of the network of secondary structures of a given RNA sequence-see http://bioinformatics.bc.edu/clote/RNAexpNumNbors .

Computational and Experimental Evaluation of Designed β-Cap Hairpins Using Molecular Simulations and Kinetic Network Models

Molecular simulation has been used to model the detailed folding properties of peptides, yet prospective computational peptide design by such approaches remains challenging and nontrivial. To test the accuracy of simulation-based hairpin design, we characterized the folding properties of a series of so-called β-cap hairpin peptides designed to mimic a conserved hairpin of LapD, a bacterial intracellular signaling protein, both experimentally by NMR spectroscopy and computationally by implicit-solvent replica-exchange molecular dynamics using three different AMBER force fields (ff96, ff99sb-ildn, and ff99sb-ildn-NMR). A unique challenge presented by these designs is the presence of both a terminal Trp-Trp capping motif and a conserved GWxQ motif in the hairpin turn required for binding to LapG. Consistent with previous studies, we found AMBER ff96 to be the most accurate when used with the OBC GBSA implicit solvent model, despite its known bias toward β-sheet conformations when used in explicit-solvent simulations. To gain microscopic insight into the folding landscape of the hairpin designs, we additionally performed parallel simulations on the Folding@home distributed computing platform using AMBER ff99sb-ildn-NMR with TIP3P explicit solvent. Markov state models (MSMs) built from trajectory data reveal a number of non-native interactions between Trp and other amino acid side chains, creating potential problems in achieving well-folded hairpin structures in solution.

Ward Clustering Improves Cross-Validated Markov State Models of Protein Folding

Markov state models (MSMs) are a powerful framework for analyzing protein dynamics. MSMs require the decomposition of conformation space into states via clustering, which can be cross-validated when a prediction method is available for the clustering method. We present an algorithm for predicting cluster assignments of new data points with Ward's minimum variance method. We then show that clustering with Ward's method produces better or equivalent cross-validated MSMs for protein folding than other clustering algorithms.

Identification of simple reaction coordinates from complex dynamics

Reaction coordinates are widely used throughout chemical physics to model and understand complex chemical transformations. We introduce a definition of the natural reaction coordinate, suitable for condensed phase and biomolecular systems, as a maximally predictive one-dimensional projection. We then show that this criterion is uniquely satisfied by a dominant eigenfunction of an integral operator associated with the ensemble dynamics. We present a new sparse estimator for these eigenfunctions which can search through a large candidate pool of structural order parameters and build simple, interpretable approximations that employ only a small number of these order parameters. Example applications with a small molecule's rotational dynamics and simulations of protein conformational change and folding show that this approach can filter through statistical noise to identify simple reaction coordinates from complex dynamics.