Christof Schütte

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.

Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations

Characterizing the equilibrium ensemble of folding pathways, including their relative probability, is one of the major challenges in protein folding theory today. Although this information is in principle accessible via all-atom molecular dynamics simulations, it is difficult to compute in practice because protein folding is a rare event and the affordable simulation length is typically not sufficient to observe an appreciable number of folding events, unless very simplified protein models are used. Here we present an approach that allows for the reconstruction of the full ensemble of folding pathways from simulations that are much shorter than the folding time. This approach can be applied to all-atom protein simulations in explicit solvent. It does not use a predefined reaction coordinate but is based on partitioning the state space into small conformational states and constructing a Markov model between them. A theory is presented that allows for the extraction of the full ensemble of transition pathways from the unfolded to the folded configurations. The approach is applied to the folding of a PinWW domain in explicit solvent where the folding time is two orders of magnitude larger than the length of individual simulations. The results are in good agreement with kinetic experimental data and give detailed insights about the nature of the folding process which is shown to be surprisingly complex and parallel. The analysis reveals the existence of misfolded trap states outside the network of efficient folding intermediates that significantly reduce the folding speed.

Hierarchical analysis of conformational dynamics in biomolecules: Transition networks of metastable states

Molecular dynamics simulation generates large quantities of data that must be interpreted using physically meaningful analysis. A common approach is to describe the system dynamics in terms of transitions between coarse partitions of conformational space. In contrast to previous work that partitions the space according to geometric proximity, the authors examine here clustering based on kinetics, merging configurational microstates together so as to identify long-lived, i.e., dynamically metastable, states. As test systems microsecond molecular dynamics simulations of the polyalanines Ala(8) and Ala(12) are analyzed. Both systems clearly exhibit metastability, with some kinetically distinct metastable states being geometrically very similar. Using the backbone torsion rotamer pattern to define the microstates, a definition is obtained of metastable states whose lifetimes considerably exceed the memory associated with interstate dynamics, thus allowing the kinetics to be described by a Markov model. This model is shown to be valid by comparison of its predictions with the kinetics obtained directly from the molecular dynamics simulations. In contrast, clustering based on the hydrogen-bonding pattern fails to identify long-lived metastable states or a reliable Markov model. Finally, an approach is proposed to generate a hierarchical model of networks, each having a different number of metastable states. The model hierarchy yields a qualitative understanding of the multiple time and length scales in the dynamics of biomolecules.

Identification of Almost Invariant Aggregates in Reversible Nearly Uncoupled Markov Chains

The topic of the present paper has been motivated by a recent computational approach to identify metastable chemical conformations and patterns of conformational changes within molecular systems. After proper discretization, such conformations show up as almost invariant aggregates in reversible, nearly uncoupled Markov chains (NUMCs). Most of the former work on this subject treated the direct problem: given the aggregates, analyze the loose coupling in connection with the computation of the stationary distribution (aggregation/disaggregation techniques). In contrast to that, the present paper focuses on the inverse problem: given the system as a whole, identify the almost invariant aggregates together with the (small) probabilities of transitions between them. A robust algorithm is worked out on the basis of some detailed perturbation analysis and illustrated at a simple molecular system.

Transition Path Theory for Markov Jump Processes

The framework of transition path theory (TPT) is developed in the context of continuous-time Markov chains on discrete state-spaces. Under assumption of ergodicity, TPT singles out any two subsets in the state-space and analyzes the statistical properties of the associated reactive trajectories, i.e., those trajectories by which the random walker transits from one subset to another. TPT gives properties such as the probability distribution of the reactive trajectories, their probability current and flux, and their rate of occurrence and the dominant reaction pathways. In this paper the framework of TPT for Markov chains is developed in detail, and the relation of the theory to electric resistor network theory and data analysis tools such as Laplacian eigenmaps and diffusion maps is discussed as well. Various algorithms for the numerical calculation of the various objects in TPT are also introduced. Finally, the theory and the algorithms are illustrated in several examples.

On the Approximation Quality of Markov State Models

We consider a continuous-time Markov process on a large continuous or discrete state space. The process is assumed to have strong enough ergodicity properties and to exhibit a number of metastable sets. Markov state models (MSMs) are designed to represent the effective dynamics of such a process by a Markov chain that jumps between the metastable sets with the transition rates of the original process. MSMs have been used for a number of applications, including molecular dynamics, for more than a decade. Their approximation quality, however, has not yet been fully understood. In particular, it would be desirable to have a sharp error bound for the difference in propagation of probability densities between the MSM and the original process on long timescales. Here, we provide such a bound for a rather general class of Markov processes ranging from diffusions in energy landscapes to Markov jump processes on large discrete spaces. Furthermore, we discuss how this result provides formal support or shows the limitations of algorithmic strategies that have been found to be useful for the construction of MSMs. Our findings are illustrated by numerical experiments.

Illustration of transition path theory on a collection of simple examples

Transition path theory (TPT) has been recently introduced as a theoretical framework to describe the reaction pathways of rare events between long lived states in complex systems. TPT gives detailed statistical information about the reactive trajectories involved in these rare events, which are beyond the realm of transition state theory or transition path sampling. In this paper the TPT approach is outlined, its distinction from other approaches is discussed, and, most importantly, the main insights and objects provided by TPT are illustrated in detail via a series of low dimensional test problems.

EMMA: A Software Package for Markov Model Building and Analysis

The study of folding and conformational changes of macromolecules by molecular dynamics simulations often requires the generation of large amounts of simulation data that are difficult to analyze. Markov (state) models (MSMs) address this challenge by providing a systematic way to decompose the state space of the molecular system into substates and to estimate a transition matrix containing the transition probabilities between these substates. This transition matrix can be analyzed to reveal the metastable, i.e., long-living, states of the system, its slowest relaxation time scales, and transition pathways and rates, e.g., from unfolded to folded, or from dissociated to bound states. Markov models can also be used to calculate spectroscopic data and thus serve as a way to reconcile experimental and simulation data. To reduce the technical burden of constructing, validating, and analyzing such MSMs, we provide the software framework EMMA that is freely available at https://simtk.org/home/emma .

Quantum-classical Liouville approach to molecular dynamics: Surface hopping Gaussian phase-space packets

In mixed quantum-classical molecular dynamics few but important degrees of freedom of a molecular system are modeled quantum-mechanically while the remaining degrees of freedom are treated within the classical approximation. Such models can be systematically derived as a first order approximation to the partial Wigner transform of the quantum Liouville-von Neumann equation. The resulting adiabatic quantum-classical Liouville equation (QCLE) can be decomposed into three individual propagators by means of a Trotter splitting: Phase oscillations of the coherences resulting from the time evolution of the quantum-mechanical subsystem. Exchange of densities and coherences reflecting non-adiabatic effects in quantum-classical dynamics. Classical Liouvillian transport of densities and coherences along adiabatic potential energy surfaces or arithmetic means thereof. A novel stochastic implementation of the QCLE is proposed in the present work. In order to substantially improve the traditional algorithm based on surface hopping trajectories [J. C. Tully, J. Chem. Phys. 93 (2), 1061 (1990)], we model the evolution of densities and coherences by a set of surface hopping Gaussian phase-space packets (GPPs) with variable width and with adjustable real or complex amplitudes, respectively. The dense sampling of phase-space offers two main advantages over other numerical schemes to solve the QCLE. First, it allows to perform a quantum-classical simulation employing a constant number of particles, i. e. the generation of new trajectories at each surface hop is avoided. Second, the effect of non-local operators in the exchange of densities and coherences can be treated without having to invoke the momentum jump approximation. For the example of a single avoided crossing we demonstrate that convergence towards fully quantum-mechanical dynamics is much faster for surface hopping GPPs than for trajectory-based methods. For dual avoided crossings the Gaussian-based dynamics correctly reproduces the quantum-mechanical result even when trajectory-based methods not accounting for the transport of coherences fail qualitatively.

Serum Peptidome Profiling Revealed Platelet Factor 4 as a Potential Discriminating Peptide Associated With Pancreatic Cancer

Purpose: Mass spectrometry–based serum peptidome profiling is a promising tool to identify novel disease-associated biomarkers, but is limited by preanalytic factors and the intricacies of complex data processing. Therefore, we investigated whether standardized sample protocols and new bioinformatic tools combined with external data validation improve the validity of peptidome profiling for the discovery of pancreatic cancer–associated serum markers. Experimental Design: For the discovery study, two sets of sera from patients with pancreatic cancer (n = 40) and healthy controls (n = 40) were obtained from two different clinical centers. For external data validation, we collected an independent set of samples from patients (n = 20) and healthy controls (n = 20). Magnetic beads with different surface functionalities were used for peptidome fractionation followed by matrix-assisted laser desorption/ionization time-of-flight (MALDI-TOF) mass spectrometry (MS). Data evaluation was carried out by comparing two different bioinformatic strategies. Following proteome database search, the matching candidate peptide was verified by MALDI-TOF MS after specific antibody-based immunoaffinity chromatography and independently confirmed by an ELISA assay. Results: Two significant peaks (m/z 3884; 5959) achieved a sensitivity of 86.3% and a specificity of 97.6% for the discrimination of patients and healthy controls in the external validation set. Adding peak m/z 3884 to conventional clinical tumor markers (CA 19-9 and CEA) improved sensitivity and specificity, as shown by receiver operator characteristics curve analysis (AUROCcombined = 1.00). Mass spectrometry–based m/z 3884 peak identification and following immunologic quantitation revealed platelet factor 4 as the corresponding peptide. Conclusions: MALDI-TOF MS-based serum peptidome profiling allowed the discovery and validation of platelet factor 4 as a new discriminating marker in pancreatic cancer.