Markov State Models in Molecular Dynamics Simulations: A Comprehensive Technical Review

Markov State Models (MSMs) provide a rigorous statistical framework for extracting long-timescale conformational dynamics from ensembles of short molecular dynamics (MD) simulations [1, 2]. By discretizing the continuous conformational space into a set of states (microstates) and estimating the probabilities of transitioning between these states over a discrete lag time, MSMs enable the construction of a kinetic model that approximates the true dynamics of the system [3, 4]. The resulting model can be used to identify metastable conformational states (macrostates), compute free energy landscapes, and simulate long-time behavior that would be computationally prohibitive to obtain directly from MD [5, 6].

MSMs are particularly valuable in structural bioinformatics for studying proteins relevant to veterinary medicine, such as viral glycoproteins, bacterial virulence factors, and host receptor complexes [7, 8]. The method bridges the gap between atomistic detail and biologically relevant timescales (microseconds to milliseconds) without requiring specialized hardware [9, 10]. This review details the core computational components of MSM construction: dimensionality reduction via time-lagged independent component analysis (TICA), clustering of trajectory frames into microstates, estimation of transition probability matrices, identification of macrostates, and visualization of structural transitions.

1. Dimensionality Reduction with Time-Lagged Independent Component Analysis (TICA)

MD trajectories typically contain many degrees of freedom (e.g., Cartesian coordinates of thousands of atoms). Direct clustering in this high-dimensional space is inefficient and prone to noise [2]. TICA is a linear dimensionality reduction method that identifies the slowest collective coordinates (independent components) in the data [11, 12]. Given a feature vector ( \mathbf{x}(t) ) (commonly dihedral angles, residue–residue distances, or root-mean-square deviations), TICA solves a generalized eigenvalue problem involving the time-lagged covariance matrix ( \mathbf{C}(\tau) = \langle \mathbf{x}(t+\tau) \mathbf{x}(t)^T \rangle ) and the instantaneous covariance matrix ( \mathbf{C}(0) = \langle \mathbf{x}(t) \mathbf{x}(t)^T \rangle ) [13, 14]. The eigenvectors corresponding to the largest eigenvalues define components that decorrelate most slowly with lag time ( \tau ) [15].

The number of TICA components retained is chosen based on the eigenvalue spectrum, typically keeping components with eigenvalues significantly above the noise floor [16]. TICA reduces the dimensionality to 5–50 dimensions, capturing the essential slow dynamics while discarding fast vibrational motions [17, 18]. This step is critical for the subsequent clustering to produce kinetically meaningful microstates [19].

2. Clustering Trajectory Frames into Microstates

After TICA projection, the reduced trajectory is discretized into a finite number of microstates using clustering algorithms [20]. The most common choices are k-means clustering and k-centers clustering (also called k-medoids) [21]. k-means partitions the data into ( N ) clusters by minimizing the sum of squared distances between points and their nearest cluster centroid [22]. The number of microstates ( N ) is typically on the order of hundreds to thousands, chosen to balance resolution (more states better approximate the continuous space) and statistical robustness (too many states lead to poorly sampled transition counts) [1].

Other clustering approaches include hierarchical clustering and density-based methods, but k-means remains popular due to its efficiency and straightforward interpretation [2]. The choice of distance metric (e.g., Euclidean distance in TICA space) and the number of microstates can be validated by examining the implied timescales as a function of ( N ) [3]. A good discretization yields consistent implied timescales across different state decompositions [4].

3. Estimating the Transition Probability Matrix

From the discretized trajectory, a count matrix ( \mathbf{C} ) is constructed, where ( C_{ij} ) is the number of transitions observed from microstate ( i ) to microstate ( j ) at a chosen lag time ( \tau ) [5]. The transition probability matrix ( \mathbf{T}(\tau) ) is then estimated by row-normalizing the count matrix: ( T_{ij}(\tau) = C_{ij} / \sum_k C_{ik} ) [6]. To avoid poor statistics from rarely visited states, a reversible transition matrix is often enforced using maximum likelihood estimation or the "reversible count" symmetrization [7, 8].

Validation of the Markovian assumption is performed by computing the implied timescales (ITS) as a function of lag time [9]. The ITS ( t_i ) are defined as ( t_i = -\tau / \ln \lambda_i ), where ( \lambda_i ) are the eigenvalues of ( \mathbf{T}(\tau) ) [10]. If the ITS become constant for increasing ( \tau ), the model is considered Markovian [11]. A typical analysis selects the smallest lag time for which the ITS plateau, balancing model accuracy with temporal resolution [12].

The transition matrix also yields the stationary (equilibrium) probability distribution ( \pi ), which is the eigenvector of ( \mathbf{T} ) corresponding to eigenvalue 1 [13]. From ( \pi ) and ( \mathbf{T} ), one can compute the free energy surface: ( F_i = -k_B T \ln \pi_i ) [14].

4. Identifying Macrostates: PCCA+ and Spectral Clustering

Although microstates provide a kinetic model, their large number (often >100) makes biological interpretation difficult. Macrostates are kinetically metastable sets of microstates that interconvert slowly [15]. The most widely used algorithm for coarse-graining microstates into macrostates is Perron Cluster Cluster Analysis (PCCA+) [16]. PCCA+ uses the eigenvectors of the transition matrix to partition microstates into a user-specified number of macrostates, maximizing the metastability of each set [17].

The number of macrostates is typically chosen by examining the gap in the eigenvalue spectrum of ( \mathbf{T} ); a distinct gap after ( m ) eigenvalues suggests ( m ) metastable states [18]. Other approaches include spectral clustering using the diffusion map or the robust Perron cluster analysis [19]. The resulting macrostates correspond to distinct conformational ensembles (e.g., open versus closed active site, ligand-bound versus unbound) and are often visualized as nodes in a free energy landscape [20].

5. Visualizing Structural Transitions in 3D

Once macrostates are defined, representative structures (e.g., centroid frames or medoid structures) for each macrostate can be extracted from the trajectory [21]. Animating transitions between these structures reveals the conformational pathways that the system follows under equilibrium conditions [22]. Common visualization pipelines use tools such as VMD, PyMOL, or modern web-based viewers (e.g., NGL Viewer) to load coordinate files for each macrostate and interpolate between them [1, 2].

To generate a transition animation, the user identifies pairs of macrostates with high transition probabilities, extracts the most probable transition path (e.g., using the transition path theory or the "reactive flux" decomposition) [3], and renders a smooth morph between the starting and ending structures. This process is often combined with free energy surface plots (e.g., as a function of the first two TICA components) to show the global landscape and the relative populations of states [4, 5].

6. Application in Veterinary Structural Biology

MSMs have been applied to a wide range of protein systems relevant to veterinary medicine and host–pathogen interactions [6]. For example, MSMs have elucidated the conformational dynamics of KRAS mutants, which are important in understanding oncogenic signaling in both human and animal cancers [7, 21, 22]. Studies on the μ-opioid receptor (MOR) have used MSMs to characterize partial agonism of compounds such as mitragynine and zerumbone, providing insights relevant to pain management in livestock and companion animals [3, 13]. Allosteric regulation of Hsp90, a chaperone targeted in cancer and parasitic infections, has been probed using MSMs combined with Gaussian accelerated MD [11]. In virology, MSMs have been used to study the ACE2–RBD interaction dynamics [9] and the SARS-CoV-2 papain-like protease [16], demonstrating the power of the method for understanding viral entry and drug inhibition. These examples underscore the versatility of MSMs in studying proteins from veterinary pathogens, including avian influenza, coronavirus, and bacterial virulence factors.

Workflow Overview

flowchart TD
    A[MD Simulation Trajectories], > B[Feature Extraction<br/>distances, dihedrals]
    B, > C[TICA Dimensionality Reduction<br/>Slow collective coordinates]
    C, > D[Clustering into Microstates<br/>k-means / k-centers]
    D, > E[Count Matrix C_ij]
    E, > F[Transition Probability Matrix T(τ)]
    F, > G[Validation: Implied Timescales]
    G, > H{Markovian?}
    H, >|No| I[Adjust lag time or clustering]
    H, >|Yes| J[PCCA+ Macrostate Identification]
    I, > D
    J, > K[Free Energy Landscape]
    J, > L[Transition Path Analysis]
    J, > M[3D Animation of Conformational Transitions]

Key Parameters in MSM Construction

Conclusion

Markov State Models provide a statistically rigorous and computationally efficient way to analyze the conformational dynamics of biological macromolecules from MD simulations. By combining TICA for dimensionality reduction, clustering for state discretization, and PCCA+ for macrostate identification, MSMs yield interpretable kinetic models and free energy landscapes. The ability to animate structural transitions between metastable states makes MSMs a powerful tool for studying mechanisms of protein function, ligand binding, and allostery. The method is broadly applicable to veterinary structural biology, enabling detailed analysis of host–pathogen interactions, receptor dynamics, and drug target conformational ensembles.

References

[1] Qaisrani MN, Kirsch C, Flötotto A et al. Bridging Atomistic and Mesoscale Lithium Transport via Machine-Learned Force Fields and Markov State Models. J Chem Theory Comput 2026. PubMed PMID: 42157717.

[2] Seki T, Ohnuki J, Okazaki KI. Conformational Dynamics of Na(+)-Pumping NADH-Quinone Oxidoreductase during Na(+) Translocation from AlphaFold-Facilitated Markov State Modeling. J Chem Inf Model 2026. PubMed PMID: 41930439.

[3] Bahari MNA, Azmi L, Fei LC et al. Mapping partial agonism of mitragynine at the µ-opioid receptor through molecular dynamics and Markov state modelling analysis. Sci Rep 2026. PubMed PMID: 41794923.

[4] Mahapatra S, Kar P. Allosteric activation of FGFR2 kinase in endometrial cancer: insights from Gaussian accelerated molecular dynamics and the Markov state model. Phys Chem Chem Phys 2026. PubMed PMID: 41736675.

[5] Mitra A, Paul S. Markov State Models Reveal How hLL-37(17-29) Tetramers Relax from a Disordered State to a Cross-α Amyloid Geometry. J Phys Chem B 2026. PubMed PMID: 41665928.

[6] Zhao L, Wang J, Ping H et al. Revealing conformational changes of GTP-bound NRAS mutants probed by GaMD and Markov state models. Phys Chem Chem Phys 2026. PubMed PMID: 41660764.

[7] Chen J, Wang J, Wang W et al. Allosteric Binding-Mediated Suppression on Activity of G12D KRAS Recognized via Markov State Model and Communication Pathway. J Phys Chem B 2026. PubMed PMID: 41610434.

[8] Zhao J, Yang Y, Zhang Z et al. Mapping the Allosteric Landscape of PPAR(γ): a Markov State Modeling and Energetic Analysis Approach. J Chem Inf Model 2026. PubMed PMID: 41498383.

[9] Zhou Y, Wang T. Dynamic insights into the structural evolution of ACE2-RBD interactions through molecular dynamics simulation, Markov state modeling, and large language model mutation prediction. J Chem Phys 2025. PubMed PMID: 41268951.

[10] Chakraborty G, Patra N. Mapping Conformational and Kinetic Landscapes of 14-3-3σ/α-Synuclein Assembly via Variational Autoencoders and Markov State Models. J Chem Inf Model 2025. PubMed PMID: 41225715.

[11] Bao H, Wang J, Zhao L et al. Ligand-mediated conformation diversity of Hsp90 revealed by GaMD simulations and Markov model. Mol Divers 2026. PubMed PMID: 41217569.

[12] Dutta S, Shukla D. Characterization of binding kinetics and intracellular signaling of new psychoactive substances targeting cannabinoid receptor using transition-based reweighting method. Elife 2025. PubMed PMID: 41181929.

[13] Ayub WMW, Marjohan NA, Khalid MH et al. Elucidating zerumbone's low-efficacy agonism at the μ-opioid receptor via molecular dynamics simulation and Markov state modeling. J Comput Aided Mol Des 2025. PubMed PMID: 41136794.

[14] Jiang Y, Duan X, Zheng J et al. Investigation of the Substrate Selection Mechanism of Poly (A) Polymerase Based on Molecular Dynamics Simulations and Markov State Model. Int J Mol Sci 2025. PubMed PMID: 41096778.

[15] Chen J, Wang J, Wang W et al. Conformational Transition and Recognition Mechanism of Eukaryotic Riboswitches Powered by Thiamine Pyrophosphate Analogues: An Elucidation through Multiple Short Molecular Dynamics Simulations and Markov Model. J Phys Chem B 2025. PubMed PMID: 40997929.

[16] Qiu Y, Hong WL, Wang X et al. Discovery of Inhibitors of SARS-CoV-2 Papain-like Protease Using Molecular Simulation and Markov State Model. J Phys Chem A 2025. PubMed PMID: 40970671.

[17] Ricci CG, Philpott JM, Torgrimson MR et al. Markovian state models uncover casein kinase 1 dynamics that govern circadian period. Biophys J 2025. PubMed PMID: 40968534.

[18] Shi D, Lv E, Li W et al. Conformational dynamics of peptide substrate recognition in polypeptide N-acetylgalactosaminyltransferase 2: Markov state models reveal the anchoring and stabilizing determinants. Int J Biol Macromol 2025. PubMed PMID: 40962091.

[19] O'Connor MS, Konovalov KA, Duvall JL et al. Role of Hsp70 chaperone in client-protein folding elucidated by Markov state modeling and NMR restraint-assisted molecular dynamics simulations. Biophys J 2025. PubMed PMID: 40873036.

[20] Hardie A, Powell FG, Lovera S et al. Elucidation of the mechanism of partial activation of EPAC1 allosteric modulators by Markov state modelling. Chem Sci 2025. PubMed PMID: 40687697.

[21] Xiao S, Alshahrani M, Hu G et al. Exploring binding and allosteric energy landscapes for the KRAS interactions with effector proteins using Markov state modeling of conformational ensembles and allosteric network modeling. Protein Sci 2025. PubMed PMID: 40671268.

[22] Zhang H, Ni D, Fan J et al. Markov State Models and Molecular Dynamics Simulations Reveal the Conformational Transition of the Intrinsically Disordered Hypervariable Region of K-Ras4B to the Ordered Conformation. J Chem Inf Model 2022. PubMed PMID: 35994329. *** Disclaimer: This article is for educational and informational purposes only. It is not intended to substitute for professional veterinary advice, diagnosis, treatment, or regulatory guidance. Always consult a licensed veterinarian or qualified specialist regarding animal health, disease diagnosis, and therapeutic decisions.