MSMBuilder2: Modeling Conformational Dynamics at the Picosecond to Millisecond Scale

Abstract

Markov State Models provide a framework for understanding the fundamental states and rates in the conformational dynamics of biomolecules. We describe an improved protocol for constructing Markov State Models from molecular dynamics simulations. The new protocol includes advances in clustering, data preparation, and model estimation; these improvements lead to significant increases in model accuracy, as assessed by the ability to recapitulate equilibrium and kinetic properties of reference systems. A high-performance implementation of this protocol, provided in MSMBuilder2, is validated on dynamics ranging from picoseconds to milliseconds.

Figure 1

(a). Relaxation timescales of models constructed with k-centers and hybrid clustering. (b). RMSD correlation functions as calculated by different clusterings. MSMs in (b) constructed with 90 ns lagtime. MSMs constructed from simulations of the WW protein; see Appendix 1.

Figure 2

Cluster centers (opaque) and randomly sampled conformations (transparent) are displayed for the most populated state from models based on the k-centers and hybrid clustering algorithms. Both models are based on simulations of the WW domain. The hybrid clusters (b) were constructed by improving the initial k-centers clustering in (a). Both clusterings have 806 states (f_max = 6.5Å).

Figure 3

Simulations of the WW protein were used to compare the performance of the symmetrized and MLE protocols. Folding free energies calculated using a two-state approximation $(-RT\text{log}(\frac{{\pi }_{\mathit{\text{folded}}}}{{\pi }_{\mathit{\text{unfolded}}}}))$, show good agreement (Δ ≤ 0.03 kcal / mol) between models constructed using the symmetrized and MLE protocols, as expected for long trajectories. The near-zero folding free energy is expected, as the simulations were performed near the melting temperature; the exact free energy depends weakly on how one defines the folded state. Here, the folded state is defined as all states with an RMSD (to crystal structure) below some cutoff value; the unfolded state is defined as the remaining states. The large RMSD values observed are due to the large conformational fluctuations observed in the high temperature (393 K) simulations.

Figure 4

Simulated two-state folding simulations generated from a reference transition matrix (a) were used to estimate transition matrices. The MLE reversible procedure (b) shows good agreement with the reference transition matrix, while the symmetrized procedure (c) shows poor agreement with the reference. Furthermore, as compared to the symmetrized estimate, the MLE estimate better recapitulates the reference equilibrium properties (d).

Figure 5

Ensemble average RMSD to native is calculated for a sequence of models constructed from WW simulations.

Figure 6

The two slowest relaxation timescales for alanine dipeptide are plotted as a function of lagtime.