---
author:
- Frank Noé
- Edina Rosta
bibliography:
- 'all.bib'
- 'own.bib'
- 'special\_issue.bib'
title: 'Special Topic: Markov Models of Molecular Kinetics'
---

#### Introduction

The *Journal of Chemical Physics* article collection on [Markov Models of Molecular Kinetics (MMMK)](https://aip.scitation.org/toc/jcp/collection/10.1063/jcp.2019.MMMK.issue-1) features recent advances developing and using Markov State Models (MSMs) [@SchuetteFischerHuisingaDeuflhard_JCompPhys151_146; @SwopePiteraSuits_JPCB108_6571; @NoeHorenkeSchutteSmith_JCP07_Metastability; @ChoderaEtAl_JCP07; @BucheteHummer_JPCB08; @PrinzEtAl_JCP10_MSM1] in atomistic molecular simulations and related applications – see [@BowmanPandeNoe_MSMBook; @SarichSchuette_MSMBook13; @ChoderaNoe_COSB14_MSMs; @HusicPande_JACS18_MSMReview] for recent MSM reviews.  MSMs have been an important driving force in molecular dynamics (MD), as they facilitate divide-and-conquer integration of short, distributed MD simulations into long-timescale predictions, they are conceptually simple and provide readily-interpretable models of kinetics and thermodynamics. 

Most MSM estimation approaches proceed by a sequence, or pipeline, of data processing steps that is also represented by MSM software packages [@SchererEtAl_JCTC15_EMMA2; @DoerrEtAl_JCTC16_HTMD; @HarriganEtAl_BJ17_MSMbuilder], and typically includes:

1.   **Featurization**: The MD coordinates are transformed into features, such as residue distances, contact maps or torsion angles [@HumphreyDalkeSchulten_JMG96_VMD; @McGibbon_BJ15_MDTraj; @SchererEtAl_JCTC15_EMMA2; @DoerrEtAl_JCTC16_HTMD], that form the input of the MSM analysis. 

2.   **Dimension reduction**: The dimension is reduced to much fewer (typically 2-100) slow collective variables (CVs), [@NoeNueske_MMS13_VariationalApproach; @NueskeEtAl_JCTC14_Variational; @PerezEtAl_JCP13_TICA; @SchwantesPande_JCTC13_TICA; @Molgedey_94; @ZieheMueller_ICANN98_TDSEP; @Mezic_NonlinDyn05_Koopman; @SchmidSesterhenn_APS08_DMD; @WilliamsKevrekidisRowley_JNS15_EDMD; @TuEtAl_JCD14_ExactDMD; @NoeClementi_COSB17_SlowCVs].  The resulting coordinates may be scaled, in order to embed them in a metric space whose distances correspond to some form of dynamical distance [@NoeClementi_JCTC15_KineticMap; @NoeClementi_JCTC16_KineticMap2]. 

3.   **Discretization**: The space may be discretized by clustering the projected data [@BowmanPandeNoe_MSMBook; @SchererEtAl_JCTC15_EMMA2; @HusicPande_JCTC17_Ward; @SheongEtAl_JCTC15_APM; @ChoderaEtAl_JCP07; @WuNoe_JCP15_GMTM; @HarriganPande_bioRxiv17_LandmarkTICA; @WeberFackeldeySchuette_JCP17_SetfreeMSM], typically resulting in 100-1000 discrete “microstates”. 

4.   **MSM estimation**: A transition matrix or rate matrix describing the transition probabilities or rate between the discrete states at some lag time $\tau$ is estimated [@Bowman_JCP09_Villin; @PrinzEtAl_JCP10_MSM1; @BucheteHummer_JPCB08; @TrendelkampSchroerEtAl_InPrep_revMSM]. 

5.   **Coarse-graining**: In order to get an easier interpretable kinetic model, the MSM from step 5 is often coarse-grained to a few states [@DeuflhardWeber_LAA05_PCCA+; @KubeWeber_JCP07_CoarseGraining; @YaoHuang_JCP13_Nystrom; @FackeldeyWeber_WIAS17_GenPCCA; @GerberHorenko_PNAS17_Categorial; @HummerSzabo_JPCB15_CoarseGraining; @OrioliFaccioli_JCP16_CoarseMSM; @NoeEtAl_PMMHMM_JCP13; @Schuette_HandbookNumAnal_2003]. 

Some method skip or combine some of these steps, novel machine learning methods attempt to integrate most or all of them in an end-to-end learning framework. 

Key in much of the methodological progress in Markov modeling has been the mathematical theory of conformation dynamics pioneered by Schütte [@SchuetteFischerHuisingaDeuflhard_JCompPhys151_146] and further developed by many contributors.  This theory models the dynamics of molecules by a Markov propagator $\mathcal{T}(\tau)$ that describes how an ensemble of molecules $\rho_{0}$ evolves in a time step $\tau$. $$\rho_{\tau}=\mathcal{T}(\tau)\rho_{0}.$$ The MSM transition matrix is a discrete version of $\mathcal{T}(\tau)$. Discretization in high-dimensional spaces is difficult to impossible, so it is important to understand the structure underlying these dynamics in order to make the MSM estimation problem feasible.  If the dynamics are in equilibrium, this propagation can further be approximated by a sum of processes $\psi_{i}$ that relax the initial distribution towards the equilibrium distribution with characteristic time scales $t_{i}$. $$\rho_{\tau}(\mathbf{x})\approx\sum_{i}\mathrm{e}^{-\tau/t_{i}}\langle\psi_{i},\rho_{0}\rangle\psi_{i}(\mathbf{x})\label{eq:EVD}$$ As $\mathrm{e}^{-\tau/t_{i}}$ decays exponentially fast in the time step $\tau$, only few terms are needed for Eq. ( \[eq:EVD\]) to be an accurate description if we focus on the long-time dynamics, i.e., the kinetics.  The key insight from this theory is no less that Markov modeling is possible even for complicated and very high-dimensional molecular systems: We cannot sample or discretize truly high-dimensional spaces, but we can do that for metastable molecular systems because we are ultimately only interested in the low-dimensional manifold spanned by a few eigenfunctions $\psi_{i}$ of the Markov operator.  Characterizing this manifold more compactly than by modeling all relevant eigenfunctions $\psi_{i}$ explicitly is subject of current research [@BittracherEtAl_arXiv2017_TransitionManifols]. 

An important cornerstone for improving MSM estimators, developing new ones and for turning MSM estimation into generic machine learning problem that can be combined with kernel machines or neural networks, is the development of variational optimization methods.  The variational approach for conformation dynamics (VAC) [@NoeNueske_MMS13_VariationalApproach; @NueskeEtAl_JCTC14_Variational] shows that eigenvalues of MSMs (the same is true for other linear Markovian models such as TICA) systematically underestimates timescales $t_{i}$ and eigenvalues $\mathrm{e}^{-\tau/t_{i}}$, and defines a variational score – essentially the sum of eigenvalues of an estimated Markov model – that ought to be maximized to optimally approximate the unknown eigenvalues and eigenfunctions in (\[eq:EVD\]).  The variational formulation is key to many contributions in the MMMK collection, and remains to be the subject of further development and application, e.g. in the context of MSM hyperparameter optimization [@McGibbonPande_JCP15_CrossValidation; @HusicPande_JCTC17_Ward]. 

While VAC only describes the scenario of equilibrium dynamics, i.e. where our dynamics have a unique equilibrium distribution and obey detailed balance, much recent research has focused on non-equilibrium Markov models [@KoltaiEtAl_Computation18_NonrevMSM; @ReuterEtAl_JCPMSM19_Schur; @KnochSpeck_JCPMSM19_Noneq; @KnochSpeck_NJP15_NoneqMSM; @WuNoe_VAMP; @MardtEtAl_VAMPnets].  While nonequilibrium MSM studies are still in their infancy, several theoretical principles are known in oder to make progress here.  An important framework for the description of such processes is that of nonequilibrium work as covered by the Jarzynski fluctuation theorem [@Jarzynski_PRL97].  From a machine learning and optimization perspective, we can replace the eigenvalue decomposition in (\[eq:EVD\]) with a singular value decomposition of the operator whose components can be approximated with the variational approach of Markov processes (VAMP) [@WuNoe_VAMP].  VAMP is exploited in several contributions in the MMMK collection. 

#### Feature selection

One step in MSM estimation that had not yet undergone systematic analysis or optimization is the selection of features used as an input.  For solvated molecules, roto-translationally invariant features were usually chosen based on what works best for a given application – including intramolecular distances, angles, contact matrices or features implicitly defined by pairwise metrics, such as minimal root mean square distance.  The VAC approach has previously been invoked to define variationally optimal features for short peptides, in the spirit of defining optimized basis sets in quantum chemistry [@VitaliniNoeKeller_JCTC15_BasisSet].  In the MMMK collection, Scherer et al. [ @SchererEtAl_JCPMSM19_VAMPFeatures] propose to use the VAMP approach to variationally score different candidates of features for a given MD analysis task.  Considering a large list of candidate features and all 12 fast-folding protein simulations published by DESRES [@LindorffLarsenEtAl_Science11_AntonFolding], the authors of [@SchererEtAl_JCPMSM19_VAMPFeatures] conclude that a combination of residue-residue contact signals that decay exponentially in the distance and backbone torsions performs best for protein folding. 

#### Slow collective variables

A major leap forward in MSM construction was the finding that machine-learning methods that identify a manifold of slow collective variables (CVs) [@NoeClementi_COSB17_SlowCVs], such as the time-lagged independent analysis (TICA) method [@Molgedey_94], led to superior MSMs [@PerezEtAl_JCP13_TICA; @SchwantesPande_JCTC13_TICA; @NaritomiFuchigami_JCP11_TICA].  Intuitively, this success is due to the fact that MSMs aim a modeling the kinetics between metastable states, and first reducing the dimension to the manifold of slow (kinetic) processes makes subsequent geometric operations such as clustering much simpler and faster than directly working in a high-dimensional space.  Theoretically, these methods can indeed be derived from VAC [@NoeNueske_MMS13_VariationalApproach; @PerezEtAl_JCP13_TICA; @NueskeEtAl_JCTC14_Variational], and thus be showed to variationally approximate the eigenfunctions of the Markov propagator (\[eq:EVD\]). 

Several papers in the MMMK collection develop this approach further.  Karasawa et al. [ @KarasawaEtAl_JCPMSM19_RMA] propose and extension to relaxation mode analysis (RMA) [@TakanoMiyashita_JPSJ95_RelaxationModeAnalysis], a close sibling of TICA, in which one first solves an eigenvalue problem of the time correlation matrix of features to identify the manifold of slow CVs, and then finds the subspace in which the matrix is positive definite, promoting numerically stable estimate of relaxation rates. 

A close relative of VAC is the spectral gap optimization of order parameters (SGOOP) developed by Tiwary [@TiwaryBerne_PNAS16_SGOOP].  While both SGOOP and VAC find a manifold of slow CVs by maximizing the largest eigenvalues, SGOOP combines this principle with a maximum Caliber based estimation of the transition or rate matrix, and is thus applicable to enhanced-sampling simulations where dynamics are not readily available.  In the MMMK collection, Smith et al.  develop a multidimensional version of SGOOP by introducing conditional probability factorization, and demonstrate its usefulness on the rare-event dissociation pathway of benzene from Lysozyme [@SmithEtAl_JCPMSM19_SGOOP2]. 

Paul et al. [ @PaulEtAl_VAMP_dimred] generalize the idea of VAC and TICA to nonequilibrium processes: They show that VAMP can be used to variationally find slow CVs in systems that are driven by external forces, such as an ion channel in an electrical field.  Operationally, the approach is as easy as TICA: time-correlation matrices between features are computed and a singular value decomposition yields the slow CVs.  The MSMs estimated in this manifold reveal the circular fluxes between long-lived states driven by the external potential [@PaulEtAl_VAMP_dimred]. 

#### Estimating transition matrices and other quantities

It is easy to show that the estimation of MSM transition matrices by maximizing the Markov chain likelihood is statistically unbiased if all simulations are in global equilibrium.  However, MSMs are usually estimated from short simulation trajectories that may be simulated under equilibrium dynamics, but whose starting points do not start from a global equilibrium distribution.  Nüske et al. [ @NueskeEtAl_JCP17_OOMMSM] have derived the mathematical form of the MSM estimation error for such data, and have proposed a reweighting method that allows the user to estimate MSMs without this initial state bias.  In the MMMK collection, [@BacciEtAl_JCPMSM19_Initial] provide a new estimation method for the same aim which is based on statistical resampling.  As always in machine learning, there is a trade-off between bias and variance of estimators [@Vapnik_IEEE99_StatisticalLearningTheory] that all these bias-reducing estimators must face.  While most MSM estimators have been developed with the aim of reducing the bias, a systematic account for MSM estimators with an optimal bias-variance trade-off is still an open issue for the future. 

As described above, slow CVs, transition rates and MSM transition matrices can be viewed as the result of a variational optimization process, e.g. using VAC, VAMP or SGOOP.  From a mathematical point of view, all methods which do this via some form of linear combination of basis functions – and this includes TICA and standard MSM transition matrix estimators – can also be described by the Galerkin approximation framework [@SarichSchuette_MSMBook13].  The idea of the Galerkin approach is as follows: we define basis functions – the mean-free feature functions in TICA or indicator functions denoting where Markov states are in discrete MSMs – and consider the projection of the dynamics onto this basis set.  The Galerkin approach then gives us expression for dynamical quantities, such as the Markov propagator eigenfunctions and its eigenvalues / relaxation rates, based on linear combinations of these basis functions.  While the Galerkin approach is primarily a mathematical explanation of what happens in TICA or MSM estimation algorithms, Thiede et al show in the MMMK collection how it can be exploited and expanded to develop better MSM estimators, and also obtain direct estimators for quantities that are usually estimated via MSMs, such as committor functions [@ThiedeEtAl_JCPMSM19_Galerkin]. 

As mentioned in the context of core-based MSMs, committors, mean first-passage times (MFPTs), milestones and MSMs are deeply connected.  In the MMMK collection, Berezhkovskii and Szabo [@BerezhkovskiiSzabo_JCPMSM19_MSMcommitttor] further our theoretical understanding of these relationships by showing why exact MFPTs can be computed via a milestoning MSM and provide a relationship between the equilibrium population of milestoning MSM states and the committor functions. 

Building upon previous work done on a variational framework for the identification of Markovian transition states, [@MartiniEtAl_PRX17_VariationalTS], Kells et al. [ @KellsEtAl_JCPMSM19_MFPT] develop a variationally optimal coarse-graining framework for MSM transition matrices that has broad applicability and for time series analysis of large datasets in general.  They demonstrate that coarse-graining an MSM into two or three states with this method has a simple physical interpretation in terms of mean first passage times and fluxes between the coarse grained states.  Results are presented using both analytic test potentials and MD simulations of pentalanine. 

#### Markov model estimation with rare-events

While MSMs effectively turn the problem of estimating molecular kinetics and thermodynamics into an embarrassingly parallel process, estimating a statistically precise or even connected MSM is still hampered by sampling the rare transition events sufficiently often.  A manifold of MSM-based approaches have been proposed to address the sampling problem, most prominently: 1) Adaptive sampling approaches, where an MSM is used which starting points for new simulations are most promising to discover new states or reduce statistical error [@Singhal_JCP07; @HuangBowmanPande_PNAS09_AdaptiveSeeding; @BowmanEnsignPande_JCTC2010_AdaptiveSampling; @PretoClementi_PCCP14_AdaptiveSampling; @DoerrDeFabritiis_JCTC14_OnTheFly; @DoerrEtAl_JCTC16_HTMD; @PlattnerEtAl_NatChem17_BarBar; @HruskaEtAl_JCPMSM19_Adaptive], and 2) multi-ensemble Markov modeling approaches, which estimate MSMs with the aid of generalized ensemble simulations (multiples temperatures or biases), in order to exploit expedited rare-event sampling [@WuMeyRostaNoe_JCP14_dTRAM; @RostaHummer_DHAM; @MeyWuNoe_xTRAM; @WuEtAL_PNAS16_TRAM; @StelzlHummer_JCTC17_KineticsREMD; @PaulEtAl_PNAS17_Mdm2PMI].  In the MMMK collection, several new methods are proposed to construct MSMs without sampling the rare events by brute force. 

Adaptive sampling is considered in Hruska et al. [ @HruskaEtAl_JCPMSM19_Adaptive].  While a variety of adaptive sampling methods have been developed before, the authors conduct a systematic study of the effectiveness of different adaptive sampling strategies on several fast folding proteins. [ @HruskaEtAl_JCPMSM19_Adaptive] provides theoretical limits for the adaptive sampling speed-up and shows that different adaptive sampling strategies are optimal, depending on sampling starts without prior knowledge of the metastable states, or whether some states are already known and finding new ones is the aim. 

By combining the maximum caliber approach [@Jaynes_ARPC80_MaxCal; @PresseEtAl_RMP13_MaxCal] with optimal transport theory, Dixit and Dill develop an approach to approximate MSM rate matrices from short non-equilibrium simulations [@DixitDill_JCPMSM19_MSMtransport].  Maximum caliber-based estimation of MSMs is used in Meral et al. [ @MeralEtAl_JCPMSM19_DNA] in combination with enhanced sampling using well-tempered Metadynamics [@LaioParrinello_PNAS99_12562; @BarducciBussiParrinello_PRL08_WellTemperedMetadynamics].  The authors apply their framework to the challenging problem of studying the activation of a G Protein-Coupled Receptor (GPCR), here the $\mu$-opioid receptor.  They demonstrate that the caliber-derived transition rates are in agreement with those obtained from adaptive sampling, suggesting that the framework is of general usefulness. 

Another approach to avoid waiting for the rare events to happen is to speed up sampling between known metastable states using transition path methods, such as transition path sampling [@BolhuisChandlerDellagoGeissler_AnnuRevPhysChem02_TPS], transition interface sampling (TIS) [@Du2011], or Forward Flux Sampling (FFS) [@AllenFrenkelWolde_JCP06_FFS], and subsequently constructing coarse-grained MSMs from the transition path statistics.  Recently developed software for transition path simulations facilitates this task [@Swenson_JCTC19_OPS1].  In the MMMK collection, Qin et al. [ @QinEtAl_JCPMSM19_TIS] develop the reweighted partial path (RPP) method approach which can efficiently reweight TIS or FFS simulations in order to derive equilibrium distributions of states or free energy profiles. 

Path-based sampling is also considered in Zhu et al. [ @ZhuEtAl_JCPMSM19_TAPS].  The authors develop a new path-searching method for connecting different metastable states of biomolecules that employs ideas from the traveling-salesman problem.  Their TAPS algorithm outperforms the string method by 5 to 8 times for peptides in vacuum and solution, suggesting that it is an efficient method to obtain initial pathways and intermediates that facilitate the construction of MSMs and thereby full kinetics of complex conformational changes. 

#### Clustering and coarse-graining

A successful class of kinetic models are core-based MSMs, originally proposed in [@BucheteHummer_JPCB08].  Core-based MSMs directly go from a low-dimensional manifold of feature space to an MSM of few metastable states, skipping over the traditional approach of clustering that space into microstates and coarse-graining the microstate transition matrix.  The basic idea of core-based MSMs is to identify dynamical cores – the most densely populated regions of state space which are parts of metastable states – and estimate an MSM from the rare transition paths between cores.  Theoretically, core MSMs are closely related to milestoning [@FaradjianElber_JCP04_Milestoning], can be shown to approximate committor functions between metastable states, which are in turn approximating the eigenfunctions of the Markov propagator in metastable systems [@SchuetteEtAl_JCP11_Milestoning]. 

A natural approach to identify cores are density-based clustering algorithms [@LemkeKeller_JCP16_DensityBasedClustering; @SittelStock_JCP16_PathBasedClustering].  In the MMMK collection, Nagel et al. [ @NagelEtAl_JCPMSM19_DynamicalCoring] propose an extension of their previous density-based coring algorithm [@SittelStock_JCP16_PathBasedClustering], which avoids misclassification of MD simulation frames to cores by requiring a minimum time spend in a new core to qualify as a core transition.  They demonstrate that dynamical coring obtains better MSMs using alanine dipeptide dynamics and Villin headpiece folding as examples. 

#### Nonequilibrium Markov models

Deviations from equilibrium can come in different forms: An ion channel in an electric field may be in steady-state, i.e. it has an unique stationary distribution, but does not obey detailed balance.  A spectroscopically probed molecule may be subject to a period external force.  When a molecular system is expanded by pulling it with a nonequilibrium optical tweezer experiment, even the dynamics themselves become time-dependent and neither a stationary distribution exists nor detailed balance is obeyed.  These different degrees of nonequilibrium call for different analysis methods that are only beginning to unfold now.  VAMP-based identification of the slow kinetics manifold for nonequilibrium has been discussed above [@PaulEtAl_VAMP_dimred]. 

In [@ReuterEtAl_JCPMSM19_Schur], Reuter et al.  generalize the popular robust Perron Cluster Cluster Analysis (PCCA+) method for coarse-graining transition matrices to obtain metastable states.  Their generalized method (G-PCCA) decomposes the MSM transition matrix with a Schur decomposition instead of an eigenvalue decomposition, and can obtain metastable states as well as slow cyclical processes from transition matrices that do not obey detailed balance. 

Knoch and Speck [@KnochSpeck_JCPMSM19_Noneq] develop a method to construct MSMs for systems that are periodically driven, and illustrate the method using a alanine dipeptide molecule that is exposed to a periodic electric field. 

#### Hidden Markov Models and experimental data

Hidden Markov models (HMMs) are an alternative to MSMs and have been used to obtain few-state kinetic models from MD data [@NoeEtAl_PMMHMM_JCP13; @McGibbon_ICML14_HMM].  They have also been used to extract kinetics directly from experimental trajectories, such as single-molecule FRET or optical tweezer measurements, which only track one or a few experimental observables over time instead of the full configuration vector [@GebhardRief_PNAS10_AFMEnergyLandscapeProtein; @PirchiHaran_NatureComms11_SingleMoleculeFRET; @KellerNoe_JACS14_HMM-FRET].  Similar methods have been used in order to analyze the kinetics from short FRET trajectories of single molecules diffusing through a confocal volume [@GopichSzabo_FRETCorrelation_JCP09; @GopichSzabe_HMMFRET_JCP09].  Much of MSM theory can be reused when dealing with HMMs, for example in order to compute relaxation times and a hierarchical decomposition of the system kinetics into metastable states with different lifetimes [@KellerNoe_JACS14_HMM-FRET]. 

In the MMMK collection, Jazani et al.  develop a HMM which analyzes fluorescence data from molecules diffusing through a single confocal volume [@JazaniEtAl_JCPMSM19_Tracking].  Although the fluorescence data stems from a single sensor – in contrast to wide-field optical microscopy – [@JazaniEtAl_JCPMSM19_Tracking] shows that the intensity fluctuations resulting from the fact that the non-homogeneity of the confocal excitation volume bears information about the spatial location of the molecule that can be exploited to reconstruct molecular diffusion paths. 

#### Machine Learning

Recently, the classical approach of constructing MSMs has been disrupted by replacing the traditional estimation pipeline (see above) by VAMPnets, where a deep neural network is trained with VAMP to map from high-dimensional coordinate or feature space to a few-state MSM [@MardtEtAl_VAMPnets].  Deep neural networks are now routinely used for several MSM-related tasks, such as learning slow CVs or aiding rare event sampling [@RibeiroTiwary_JCP18_RAVE; @ZhangYanNoe_ChemRxiv19_TALOS; @Hernandez_PRE18_VTE; @Jung2019; @WangEtAl_NatComm19_PastFutureBottleneck]. 

Another important machine-learning framework are kernel methods.  Kernel methods have been previously used for TICA [@HarmelingEtAl_NeurComput03_KernelTDSEP; @SchwantesPande_JCTC15_kTICA], and are also underlying the diffusion map approach that is a popular MD analysis framework [@CoifmanLafon_PNAS05_DiffusionMaps; @RohrdanzClementi_JCP134_DiffMaps].  Following a similar approach as VAMPnets.  Klus et al.  develop a VAMP-optimal kernel method [@KlusEtAl_JCPMSM19_Kernel] to estimate conformation dynamics directly using a kernel function acting on molecular feature space.  They show that established linear models such as TICA and MSMs are special cases of their kernel model and demonstrate the computation of metastable states and kinetics for alanine dipeptide dynamics and NTL9 protein folding. 

#### Software

An important technology driving MSM method development, dissemination and application is publicly available and software.  Luckily, two large-scale and widely used open-source packages, PyEMMA [@SchererEtAl_JCTC15_EMMA2] and MSMbuilder [@HarriganEtAl_BJ17_MSMbuilder], exist that implement a wide range of MSM methods and welcome contributions from the community. 

Recently, new software packages have added that publish additional MSM methods and techniques, e.g. [@Swenson_JCTC19_OPS1; @DeSancho_JCIM19_MasterMSM].  In the MMMK collection, Porter et al. [ @PorterEtAl_JCPMSM19_Enspara] present the Enspara library which is geared towards scalability to large data or models, i.e. MSMs with many states or from very large datasets.  Enspara includes parallelized implementations of computationally intensive operations, and represents a flexible framework for MSM construction and analysis. 

#### Applications

While MSMs and related techniques in the molecular sciences have been primarily developed to study peptide and protein folding, they are now used for a wide range of dynamical processes, including the study of liquids, aggregation, and structural transitions in materials such as alloys.  The MMMK article collection is no exception and contains MSM applications to various interesting molecular processes.  In all of these examples, the MSM framework reveals new physical or biological insight by revealing structures and transition processes at an unprecedented degree of detail. 

Liquid water is a surprisingly rich and complex dynamical system.  Long-standing questions, for example, include which dynamical rearrangements lead to the picosecond dynamics observed in spectroscopic data of liquid water.  The difficulty in answering such questions – even with accurate molecular models at hand – lies in data analysis: how can we define “state space” in a practical way and which molecular features are suitable in a liquid of molecules that are diffusing around, constantly switching between states that are identical up to the exchange of labels.  In the MMMK collection, Schulz et al.  use MSM methodology to pursue a detailed analysis of liquid water [@SchulzEtAl_JCPMSM19_DNA].  They solve the permutation-invariance problem by considering each water trimer as a subsystem in a 12-dimensional space defined by aligning the coordinate system to one of the water molecules, and then perform an MSM analysis using all water trimer trajectories of a solvent box simulation in this space.  The analysis suggests which exact transition processes are observed by experiment and how elementary dynamical processes, such as hydrogen-bond exchange in liquid water occur in detail. 

Gopich and Szabo [@GopichSzabo_JCPMSM19_TwoSite] work out a detailed analysis of diffusion-limited kinetics of a ligand to a macromolecule with two competing binding sites.  Their results indicate that the kinetics of such a system are surprisingly rich, the presence of the second empty binding site can slow down binding to the first as a result of competition, or it can speed up binding when populated as direct transitions of ligands between the two binding sites are possible. 

Shin and Kolomeisky employ MSM methods in order to model the kinetics of a one-dimensional walker with conformational changes that affect its transition probabilities [@ShinKolomeisky_JCPMSM19_Walker].  Biological systems have many examples of such processes, such as the dynamics of molecular motors along filaments, whose motion depends on the current conformation of the motor protein.  The authors derive a phase diagram of such systems exhibiting several dynamical regimes of the one-dimensional search process that are determined by the ratios of the relevant length scales. 

Pinamonti et al.  combine an advanced clustering technique with core-based MSMs in order to analyze the process of RNA base fraying in detail [@PinamontiEtAl_JCPMSM19_DNA].  The dynamics of four different RNA duplexes are analyzed and an interesting interplay between the equilibrium probability of intermediate states and the overall fraying kinetics is described. 

Chakraborty and Wales [@ChakrabortyWales_JCPMSM19_DPS] obtain an MSM of the adenine-adenine RNA conformational switch using the discrete path sampling technique (DPS) [@Wales_MolPhys100_3285].  DPS allows the authors to probe very rare events, with interconversion time scale here predicted to be in the range of minutes.  Several competing structures, separated by high barriers are found but the two main energy funnels lead to the major and minor conformations known from NMR experiments. 

Similar issues with permutation invariances exist when studying aggregation and self-assembly of many identical molecules.  To this end, Sengupta et al. [ @SenguptaEtAl_JCPMSM19_MFPT] construct CVs from descriptors that are invariant with respect to permutation of identical molecules.  Using these CVs, the authors construct MSMs to describe the aggregation of a subsequence of Alzheimers amyloid-$\beta$ peptide.  The results suggest that disordered and $\beta$-sheet oligomers do not interconvert, and thus amyloid formation relies on having formed ordered aggregates from the very beginning. 

#### Conclusion

Markov modeling has come of age.  In the molecular sciences, it has grown from an activity practiced by a handful of groups to a technique used by a large fraction – if not the majority – of MD simulation groups.  Markov modeling has also gone beyond molecular sciences and found applications in other areas of dynamical systems.  Recently it has been found that key MSM techniques have evolved in parallel in other fields under different names [@ZieheEtAl_JMLR04_FastTDSEP; @SchmidSesterhenn_APS08_DMD; @TuEtAl_JCD14_ExactDMD; @WilliamsKevrekidisRowley_JNS15_EDMD], and consolidating these efforts is fruitful for all these fields. 

While basic aspects of Markov modeling, such as steps of the data processing pipeline outlined in the beginning of this editorial are now well established and largely under control, new research questions have emerged, such as how to treat nonequilibrium processes, how to deal with systems with permutation invariance such as liquid and membrane systems, and how to exploit modern machine learning methods for molecular thermodynamics and kinetics.  The contributions of the MMMK collection are a cross-section of this change. 

An important driving force for the development of the field was, and is, the availability of open-source well-maintained software.  Currently MSM softwares are primarily developed and maintained by individual groups.  We believe that a key to make this development sustainable and maintainable – and therefore to preserve the accumulated methodological knowledge for the community – is to move these softwares from single groups to communities, or in other words dissociate them from individual principle investigators, e.g. by merging packages from different groups.  The next years will be crucial in order to show whether this step succeeds, and therefore whether MSM research can proceed at full steam. 

#### Acknowledgements {#acknowledgements .unnumbered}

We are grateful to the staff and editors of Journal of Chemical Physics who proposed the MMMK collection and did the heavy lifting in collecting and editing papers, especially Erinn Brigham, John Straub and Peter Hamm.  F.N. acknowledges funding from Deutsche Forschungsgemeinschaft (CRC 1114, projects C03 and A04) and European Commission (ERC CoG 772230 “ScaleCell”).  E.R. acknowledges funding from EPSRC (EP/R013012/1, EP/L027151/1, and EP/N020669/1) and European Commission (ERC StG 757850 “BioNet”).