Multiensemble Markov models of molecular thermodynamics and kinetics - PubMed (original) (raw)

Multiensemble Markov models of molecular thermodynamics and kinetics

Hao Wu et al. Proc Natl Acad Sci U S A. 2016.

Erratum in

Abstract

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.

Keywords: Bennett acceptance ratio; Markov state models; enhanced sampling; molecular dynamics; transition-based reweighting analysis method.

PubMed Disclaimer

Conflict of interest statement

The authors declare no conflict of interest.

Figures

Fig. 1.

Fig. 1.

Multiensemble Markov models (MEMMs) contain thermodynamic (e.g., free energies) and kinetic (e.g., transition probabilities) information among all configuration states (subscript index) and ensembles (superscript index). TRAM estimates MEMMs by combining transitions observed between configuration states, and the statistical weights/reduced energies of samples in all ensembles (reweighting).

Fig. 2.

Fig. 2.

Relationship of different statistically optimal reweighting estimators. TRAM is the most general method considered here, it can be specialized to MBAR, discrete TRAM, and WHAM by adding the assumption of global equilibrium or performing a binning of sampled configurations.

Fig. 3.

Fig. 3.

MBAR versus TRAM using REMD simulations of alanine dipeptide; mean and SD over seven independent simulations are shown. (A) Histogram-based free-energy landscape in the backbone torsions ϕ, ψ at a temperature of 300 K. Thin lines: borders of 20 clusters used to partition the space {cos⁡ϕ,sin⁡ϕ,cos⁡ψ,sin⁡ψ} (hence the nonlinear boundaries). Thick lines: borders of four metastable sets analyzed in B. (B) Convergence of equilibrium probability estimates of sets I–IV at 300 K, as a function of simulation length. (C–F) TRAM estimates. (C) Equilibrium probabilities as a function of the lag time τ. (D) Slowest relaxation timescale at 366 K as a function of the lag time. (E) Equilibrium probabilities of sets I–IV as a function of temperature. (F) Mean first passage times at temperatures where transitions were found.

Fig. 4.

Fig. 4.

Comparison of MBAR, xTRAM, and TRAM for estimating the equilibrium distribution of a three-well potential from US data. (A) Potential function u(x,y), where thin white lines represent the borders of 20 discrete states, thick white lines represent the borders of the three potential wells, and the dashed black lines indicate the umbrella centers. (B) A simulation trajectory with the bias potential centered at x=8.33. (C) Autocorrelation times τess(x) and τess(y) with respect to x axis and y axis for different umbrella centers. (D) Average estimation errors and their SDs of MBAR, xTRAM, and TRAM for different simulation trajectory lengths over 30 independent realizations of US.

Fig. 5.

Fig. 5.

Thermodynamics and kinetics of all-atom protein–ligand binding model for trypsin–benzamidine. (A–C) US simulations. (D and E) MEMM using both unbiased and US simulations. (A) Trajectories projected on the space of the umbrella sampling (US) coordinate and the second independent component (IC). The US coordinate describes a transition from benzamidine bound to Asp-189 to benzamidine located outside the binding pocket on the surface of trypsin. The second IC corresponds to concerted opening of loop (Trp-215-Gln-221) and flipping of Trp-215. The Voronoi centers of the Markov states are shows as disks. Markov states that are irreversibly connected to the data set are shown as red disks and are excluded from the MEMM. (B) Potential of mean force (PMF) in the same coordinate space computed with MBAR; (C) PMF computed with TRAM. Besides a higher barrier along the US coordinate, the TRAM-PMF gives the Trp-occluded conformation a lower free energy compared with the MBAR result. (D) Coarse-grained kinetic network of the MEMM. Structures (i, ii, iii, and iv) are found in the four quadrants of A. The largest transition rates (where at least one direction exceeds 1/ms) between these macrostates, the unbound state and two alternatively bound states are shown as arrows. Units are events per millisecond. (E) Efficiency of TRAM in the estimation of unbinding kinetics compared with an MSM built from the same unbiased data. Shown is the probability that log⁡koff calculated from a bootstrap sample falls into the interval [0.5⁡log⁡koffall,2⁡log⁡koffall] where koffall is the TRAM estimate calculated using all data.

References

    1. Jensen MO, et al. Mechanism of voltage gating in potassium channels. Science. 2012;336(6078):229–233. -PubMed
    1. Zhu F, Hummer G. Pore opening and closing of a pentameric ligand-gated ion channel. Proc Natl Acad Sci USA. 2010;107(46):19814–19819. -PMC -PubMed
    1. Bernèche S, Roux B. Energetics of ion conduction through the K+ channel. Nature. 2001;414(6859):73–77. -PubMed
    1. Köpfer DA, et al. Ion permeation in K+ channels occurs by direct Coulomb knock-on. Science. 2014;346(6207):352–355. -PubMed
    1. Kohlhoff KJ, et al. Cloud-based simulations on Google Exacycle reveal ligand modulation of GPCR activation pathways. Nat Chem. 2014;6(1):15–21. -PMC -PubMed

LinkOut - more resources