Christoph Wehmeyer

PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models

Markov (state) models (MSMs) and related models of molecular kinetics have recently received a surge of interest as they can systematically reconcile simulation data from either a few long or many short simulations and allow us to analyze the essential metastable structures, thermodynamics, and kinetics of the molecular system under investigation. However, the estimation, validation, and analysis of such models is far from trivial and involves sophisticated and often numerically sensitive methods. In this work we present the open-source Python package PyEMMA ( http://pyemma.org ) that provides accurate and efficient algorithms for kinetic model construction. PyEMMA can read all common molecular dynamics data formats, helps in the selection of input features, provides easy access to dimension reduction algorithms such as principal component analysis (PCA) and time-lagged independent component analysis (TICA) and clustering algorithms such as k-means, and contains estimators for MSMs, hidden Markov models, and several other models. Systematic model validation and error calculation methods are provided. PyEMMA offers a wealth of analysis functions such that the user can conveniently compute molecular observables of interest. We have derived a systematic and accurate way to coarse-grain MSMs to few states and to illustrate the structures of the metastable states of the system. Plotting functions to produce a manuscript-ready presentation of the results are available. In this work, we demonstrate the features of the software and show new methodological concepts and results produced by PyEMMA.

Multiensemble Markov models of molecular thermodynamics and kinetics

Significance Molecular dynamics simulations can provide mechanistic understanding of biomolecular processes. However, direct simulation of slow transitions such as protein conformational transitions or protein–ligand dissociation are unfeasible with commonly available computational resources. Two typical strategies are (i) conducting large ensembles of short simulations and estimating the long-term kinetics with a Markov state model, and (ii) speeding up rare events by bias potentials or higher temperatures and estimating the unbiased thermodynamics with reweighting estimators. In this work, we introduce the transition-based reweighting analysis method (TRAM), a statistically optimal approach that combines the best of both worlds and estimates a multiensemble Markov model (MEMM) with full thermodynamic and kinetic information at all simulated ensembles. We introduce the general transition-based reweighting analysis method (TRAM), a statistically optimal approach to integrate both unbiased and biased molecular dynamics simulations, such as umbrella sampling or replica exchange. TRAM estimates a multiensemble Markov model (MEMM) with full thermodynamic and kinetic information at all ensembles. The approach combines the benefits of Markov state models—clustering of high-dimensional spaces and modeling of complex many-state systems—with those of the multistate Bennett acceptance ratio of exploiting biased or high-temperature ensembles to accelerate rare-event sampling. TRAM does not depend on any rate model in addition to the widely used Markov state model approximation, but uses only fundamental relations such as detailed balance and binless reweighting of configurations between ensembles. Previous methods, including the multistate Bennett acceptance ratio, discrete TRAM, and Markov state models are special cases and can be derived from the TRAM equations. TRAM is demonstrated by efficiently computing MEMMs in cases where other estimators break down, including the full thermodynamics and rare-event kinetics from high-dimensional simulation data of an all-atom protein–ligand binding model.

Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics

Inspired by the success of deep learning techniques in the physical and chemical sciences, we apply a modification of an autoencoder type deep neural network to the task of dimension reduction of molecular dynamics data. We can show that our time-lagged autoencoder reliably finds low-dimensional embeddings for high-dimensional feature spaces which capture the slow dynamics of the underlying stochastic processes-beyond the capabilities of linear dimension reduction techniques.

MHC class II complexes sample intermediate states along the peptide exchange pathway

The presentation of peptide-MHCII complexes (pMHCIIs) for surveillance by T cells is a well-known immunological concept in vertebrates, yet the conformational dynamics of antigen exchange remain elusive. By combining NMR-detected H/D exchange with Markov modelling analysis of an aggregate of 275 microseconds molecular dynamics simulations, we reveal that a stable pMHCII spontaneously samples intermediate conformations relevant for peptide exchange. More specifically, we observe two major peptide exchange pathways: the kinetic stability of a pMHCIIs ground state defines its propensity for intrinsic peptide exchange, while the population of a rare, intermediate conformation correlates with the propensity of the HLA-DM-catalysed pathway. Helix-destabilizing mutants designed based on our model shift the exchange behaviour towards the HLA-DM-catalysed pathway and further allow us to conceptualize how allelic variation can shape an individuals MHC restricted immune response.

Protein-peptide association kinetics beyond the seconds timescale from atomistic simulations

Understanding and control of structures and rates involved in protein ligand binding are essential for drug design. Unfortunately, atomistic molecular dynamics (MD) simulations cannot directly sample the excessively long residence and rearrangement times of tightly binding complexes. Here we exploit the recently developed multi-ensemble Markov model framework to compute full protein-peptide kinetics of the oncoprotein fragment 25–109Mdm2 and the nano-molar inhibitor peptide PMI. Using this system, we report, for the first time, direct estimates of kinetics beyond the seconds timescale using simulations of an all-atom MD model, with high accuracy and precision. These results only require explicit simulations on the sub-milliseconds timescale and are tested against existing mutagenesis data and our own experimental measurements of the dissociation and association rates. The full kinetic model reveals an overall downhill but rugged binding funnel with multiple pathways. The overall strong binding arises from a variety of conformations with different hydrophobic contact surfaces that interconvert on the milliseconds timescale.Binding and unbinding kinetics are important determinants of protein-protein or small molecule protein functional interactions that can guide drug development. Here the authors exploit the multi-ensemble Markov model framework to develop a computational approach that allows the estimation of binding kinetics reaching into the seconds timescale.

The structure and IR signatures of the arginine-glutamate salt bridge. Insights from the classical MD simulations

Salt bridges and ionic interactions play an important role in protein stability, protein-protein interactions, and protein folding. Here, we provide the classical MD simulations of the structure and IR signatures of the arginine (Arg)-glutamate (Glu) salt bridge. The Arg-Glu model is based on the infinite polyalanine antiparallel two-stranded β-sheet structure. The 1 μs NPT simulations show that it preferably exists as a salt bridge (a contact ion pair). Bidentate (the end-on and side-on structures) and monodentate (the backside structure) configurations are localized [Donald et al., Proteins 79, 898-915 (2011)]. These structures are stabilized by the short (+)N-H⋯O(-) bonds. Their relative stability depends on a force field used in the MD simulations. The side-on structure is the most stable in terms of the OPLS-AA force field. If AMBER ff99SB-ILDN is used, the backside structure is the most stable. Compared with experimental data, simulations using the OPLS all-atom (OPLS-AA) force field describe the stability of the salt bridge structures quite realistically. It decreases in the following order: side-on > end-on > backside. The most stable side-on structure lives several nanoseconds. The less stable backside structure exists a few tenth of a nanosecond. Several short-living species (solvent shared, completely separately solvated ionic groups ion pairs, etc.) are also localized. Their lifetime is a few tens of picoseconds or less. Conformational flexibility of amino acids forming the salt bridge is investigated. The spectral signature of the Arg-Glu salt bridge is the IR-intensive band around 2200 cm(-1). It is caused by the asymmetric stretching vibrations of the (+)N-H⋯O(-) fragment. Result of the present paper suggests that infrared spectroscopy in the 2000-2800 frequency region may be a rapid and quantitative method for the study of salt bridges in peptides and ionic interactions between proteins. This region is usually not considered in spectroscopic studies of peptides and proteins.

Foraging on the potential energy surface: A swarm intelligence-based optimizer for molecular geometry

We present a stochastic, swarm intelligence-based optimization algorithm for the prediction of global minima on potential energy surfaces of molecular cluster structures. Our optimization approach is a modification of the artificial bee colony (ABC) algorithm which is inspired by the foraging behavior of honey bees. We apply our modified ABC algorithm to the problem of global geometry optimization of molecular cluster structures and show its performance for clusters with 2-57 particles and different interatomic interaction potentials.

A Coupled Molecular Dynamics/Kinetic Monte Carlo Approach for Protonation Dynamics in Extended Systems.

We propose a multiscale simulation scheme that combines first-principles Molecular Dynamics (MD) and kinetic Monte Carlo (kMC) simulations to describe ion transport processes. On the one hand, the molecular dynamics trajectory provides an accurate atomistic structure and its temporal evolution, and on the other hand, the Monte Carlo part models the long-time motion of the acidic protons. Our hybrid approach defines a coupling scheme between the MD and kMC simulations that allows the kMC topology to adapt continuously to the propagating atomistic microstructure of the system. On the example of a fuel cell membrane material, we validate our model by comparing its results with those of the pure MD simulation. We show that the hybrid scheme with an evolving topology results in a better description of proton diffusion than a conventional approach with a static kMC transfer rate matrix. Furthermore, we show that our approach can incorporate additional dynamical features such as the coupling of the rotation of a side group in the molecular building blocks. In the present implementation, we focus on ion conduction, but it is straightforward to generalize our approach to other transport phenomena such as electronic conduction or spin diffusion.

Markov state models from short non-equilibrium simulations—Analysis and correction of estimation bias

Many state-of-the-art methods for the thermodynamic and kinetic characterization of large and complex biomolecular systems by simulation rely on ensemble approaches, where data from large numbers of relatively short trajectories are integrated. In this context, Markov state models (MSMs) are extremely popular because they can be used to compute stationary quantities and long-time kinetics from ensembles of short simulations, provided that these short simulations are in “local equilibrium” within the MSM states. However, over the last 15 years since the inception of MSMs, it has been controversially discussed and not yet been answered how deviations from local equilibrium can be detected, whether these deviations induce a practical bias in MSM estimation, and how to correct for them. In this paper, we address these issues: We systematically analyze the estimation of MSMs from short non-equilibrium simulations, and we provide an expression for the error between unbiased transition probabilities and the expected estimate from many short simulations. We show that the unbiased MSM estimate can be obtained even from relatively short non-equilibrium simulations in the limit of long lag times and good discretization. Further, we exploit observable operator model (OOM) theory to derive an unbiased estimator for the MSM transition matrix that corrects for the effect of starting out of equilibrium, even when short lag times are used. Finally, we show how the OOM framework can be used to estimate the exact eigenvalues or relaxation time scales of the system without estimating an MSM transition matrix, which allows us to practically assess the discretization quality of the MSM. Applications to model systems and molecular dynamics simulation data of alanine dipeptide are included for illustration. The improved MSM estimator is implemented in PyEMMA of version 2.3.

Efficient implementation and application of the artificial bee colony algorithm to low-dimensional optimization problems

Abstract We adapt a swarm-intelligence-based optimization method (the artificial bee colony algorithm, ABC) to enhance its parallel scaling properties and to improve the escaping behavior from deep local minima. Specifically, we apply the approach to the geometry optimization of Lennard-Jones clusters. We illustrate the performance and the scaling properties of the parallelization scheme for several system sizes (5–20 particles). Our main findings are specific recommendations for ranges of the parameters of the ABC algorithm which yield maximal performance for Lennard-Jones clusters and Morse clusters. The suggested parameter ranges for these different interaction potentials turn out to be very similar; thus, we believe that our reported values are fairly general for the ABC algorithm applied to chemical optimization problems.