Markov State Models for Rare Events in Molecular Dynamics

Abstract

Rare, but important, transition events between long-lived states are a key feature of many molecular systems. In many cases, the computation of rare event statistics by direct molecular dynamics (MD) simulations is infeasible, even on the most powerful computers, because of the immensely long simulation timescales needed. Recently, a technique for spatial discretization of the molecular state space designed to help overcome such problems, so-called Markov State Models (MSMs), has attracted a lot of attention. We review the theoretical background and algorithmic realization of MSMs and illustrate their use by some numerical examples. Furthermore, we introduce a novel approach to using MSMs for the efficient solution of optimal control problems that appear in applications where one desires to optimize molecular properties by means of external controls.

1. Introduction

Stochastic processes are widely used to model physical, chemical or biological systems. The goal is to approximately compute interesting properties of the system by analyzing the stochastic model. As soon as randomness is involved, there are mainly two options for performing this analysis: (1) Direct sampling and (2) the construction of a discrete coarse-grained model of the system. In a direct sampling approach, one tries to generate a statistically significant amount of events that characterize the property of the system one in which is interested. For this purpose, computer simulations of the model are a powerful tool. For example, an event could refer to the transition between two well-defined macroscopic states of the system. In chemical applications, such transitions can often be interpreted as reactions or, in the context of a molecular system, as conformational changes. Interesting properties are, e.g., average waiting times for such reactions or conformational changes and along which pathways the transitions typically occur. The problem with a direct sampling approach is that many interesting events are so-called rare events. Therefore, the computational effort for generating sufficient statistics for reliable estimates is very high, and particularly if the state space is continuous and high dimensional, estimation by direct numerical simulation is infeasible.

Available techniques for rare event simulations in continuous state space are discussed in [1]. In this article, we will discuss approach (2) to the estimation of rare event statistics via discretization of the state space of the system under consideration. That is, instead of dealing with the computation of rare events for the original, continuous process, we will approximate them by a so-called Markov State Model (MSM) with discrete finite state space. The reason is that for such a discrete model, one can numerically compute many interesting properties without simulation, mostly by solving linear systems of equations as in discrete transition path theory (TPT) [2]. We will see that this approach, called Markov State Modeling, avoids the combinatorial explosion of the number of discretization elements with the increasing size of the molecular system in contrast to other methods for spatial discretization.

The actual construction of an MSM requires one to sample certain transition probabilities of the underlying dynamics between sets. The idea is: (1) to choose the sets such that the sampling effort is much lower than the direct estimation of the rare events under consideration; and (2) to compute all interesting quantities for the MSM from its transition matrix, cf. [2,3]. There are many examples for the successful application of this strategy. In [4], for example, it was used to compute dominant folding pathways for the PinWW domain in explicit solvent. However, we have to make sure that the Markov State Model approximates the original dynamics well enough. For example, the MSM should correctly reproduce the timescales of the processes of interest. These approximation issues have been discussed since more than a decade now [5,6]; in this article, we will review the present state of research on this topic. In the algorithmic realization of Markov State Modeling for realistic molecular systems, the transition probabilities and the respective statistical uncertainties are estimated from short molecular dynamics (MD) trajectories only, cf.[7]. This makes Markov State Modeling applicable to many different molecular systems and processes, cf.[8–13].

In the first part of this article, we will discuss the approximation quality of two different types of Markov State Models that are defined with respect to a full partition of state space or with respect to so-called core sets. We will also discuss the algorithmic realization of MSMs and provide references to the manifold of realistic applications to molecular systems in equilibrium that are available in the literature today.

The second part will show how to use MSMs for optimizing particular molecular properties. In this type of application, one wants to steer the molecular system at hand by external controls in a way such that a pre-selected molecular property is optimized (minimized or maximized). That is, one wants to compute a specific external control from a family of admissible controls that optimizes the property of interest under certain side conditions. The property to be optimized can be quite diverse: For example, it can be (1) the population of a certain conformation that one wants to maximize under a side condition that limits the total work done by the external control or (2) the mean first passage time to a certain conformation that one wants to minimize (in order to speed up a rare event), but under the condition that one can still safely estimate the mean first passage time of the uncontrolled system. The theoretical background of case (1) has been considered in [14], for example, and of case (2) in [1,15]. There, one finds the mathematical problem that has to be solved in order to compute the optimal control. Here, we will demonstrate that one can use MSMs for the efficient solution of such a mathematical problem (for both cases). We will see that the spatial discretization underlying an MSM turns the high-dimensional continuous optimal control problem into a rather low-dimensional discrete optimal control problem of the same form that can be solved efficiently. Based on these insights, MSM discretization yields an efficient algorithm for solving the optimal control problem, whose performance we will outline in some numerical examples, including an application to alanine dipeptide.

2. MSM Construction

Let (X_t)_t≥0 be a time-continuous Markov process on a continuous state space, E, e.g., E ⊂ R^d. That is, X_t is the state of the molecular system at time t resulting from any usually used form of molecular dynamics simulation, be it based on Newtonian dynamics with thermostats or resulting from Langevin dynamics or other diffusion molecular dynamics models. The idea of Markov State Modeling is to derive a Markov chain, (X̂_k)_k∈ℕ, on a finite and preferably small state space Ê = {1, ..., n} that models characteristic dynamics of the continuous process, (X_t). For example, in molecular dynamics applications, such characteristic dynamics could refer to protein folding processes [16,17], conformational rearrangements between native protein substates [18,19], or ligand binding processes [20]. Since the approximating Markov chain, (X_k)_k∈ℕ, lives on a finite state space, the construction of an MSM boils down to the computation of its transition matrix, P:

$$P_{\mathit{ij}}=\mathbb{P}\left[{\hat{X}}_{k+1}=j|{\hat{X}}_k=i\right]$$

(1)

The main benefit is that for a finite Markov chain, one can compute many interesting dynamical properties directly from its transition matrix, e.g., timescales and the metastability in the system [5,21,22], a hierarchy of important transition pathways [2] or mean first passage times between selected states. With respect to an MSM, these computations should be used afterwards to answer related questions for the original continuous process. To do this, we must be able to link the states of the Markov chain back to the spatial information of the original process, and the approximation of the process (X_t) by the MSM must be valid in some sense.

Having this in mind, the first natural idea is to let the states of an MSM correspond to sets A₁, ..., A_n ⊂ E in continuous state space that form a full partition, i.e.:

$$A_i\cap A_j=\varnothing \hspace{0.17em}\text{for}\hspace{0.17em}i\ne j,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\underset{i=1}{\overset{n}{\cup }}{A}_i=E$$

(2)

Typical choices for such sets are box discretizations or Voronoi tessellations [23]. For such a full partition, it is trivial to also define a corresponding discretized process by the original switching dynamics between the sets. For a given lag time, τ > 0, we can define the index process:

$${\tilde{X}}_k=i\iff X_{k\tau }\in A_i$$

(3)

It is well known that this process is not Markovian, mainly due to the so-called recrossing problem. This refers to the fact that the original process typically crosses the boundary between two sets, A_i and A_j, several times when transitions take place, as illustrated in Figure 1. This results in cumulative transitions between indices i and j for the index process, that is, a not memoryless transition behavior.

The non-Markovianity of the index process is often seen as a problem in Markov State Modeling, because many arguments assume that X̃_k is a Markov process. In this article, we will not make this assumption. We interpret the process (X̃_k) as a tool to construct the following transition matrix, P^τ:

$$P_{\mathit{ij}}^{\tau }=\mathbb{P}\left[{\tilde{X}}_{k+1}=j|{\tilde{X}}_k=i\right]=\mathbb{P}\left[X_{(k+1)\tau }\in A_j|X_{k\tau }\in A_i\right]$$

(4)

and, hence, the MSM as the Markov chain, (X̂_k)_k∈ℕ, associated with this transition matrix. From above, it is clear that, in general, we have X̂_k ≠ X̃_k, and in [24] it was analyzed how these two processes relate in terms of density propagation. In the following, we will show under which assumptions and in which sense the MSM (X̂_k) will be a good approximation of the original dynamics given by (X_t). For convenience, we will usually write P^τ ≡ P and leave the τ-dependence implicit.

3. Analytical Results

In order to compare the MSM to the continuous process, we introduce one of the key objects for our analysis, the transfer operator of a Markov process. We assume that the Markov process (X_t) has a unique, positive invariant probability measure, μ, and that it is time-reversible. Then, for any time-step, t ≥ 0, we define the transfer operator, T_t, via the property:

$${\int }_A{T}_{t}v(y)\mu (\mathit{dy})={\int }_{E}v(x)p(t,x,A)\mu (\mathit{dx})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\text{all}\hspace{0.17em}\text{measurable}\hspace{0.17em}A$$

(5)

as an operator T_t : L²(μ) → L²(μ). Here, p(t, x, A) = ℙ[X_t ∈ A|X₀ = x] defines the transition probability measure and L²(μ) denotes the Hilbert space of functions v with:

$${\int }_{E}v{(y)}^2\mu (\mathit{dy})\ge \infty$$

(6)

and the scalar product:

$$\langle v,w\rangle ={\int }_{E}v(y)w(x)\mu (\mathit{dy})$$

(7)

Note that T_t is nothing else other than the propagator of densities under the dynamics, but the densities are understood as densities with respect to the measure, μ. That is, if the Markov process is initially distributed according to:

$$\mathbb{P}\left[X_0\in A\right]={\int }_A{v}_0(x)\mu (\mathit{dx})$$

(8)

its probability distribution at time t is given by:

$$\mathbb{P}\left[X_t\in B\right]={\int }_B{v}_t(x)\mu (\mathit{dx}),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{v}_t=T_t{v}_0$$

(9)

The benefit of working with μ-weighted densities is that the transfer operator, T_t, becomes essentially self-adjoint on L²(μ) for all cases of molecular dynamics satisfying some form of detailed balance condition. Hence, it has real eigenvalues and orthogonal eigenvectors with respect to Equation (7) (or, at least, the dominant spectral elements are real-valued). Moreover, the construction of an MSM can be seen as a projection of the transfer operator [25]. Assume Q is an orthogonal projection in L²(μ) onto an n-dimensional subspace, D ⊂ L²(μ), with 𝟙 ∈ D, and χ₁, ..., χ_n is a basis of D. Then, the so-called projected transfer operator, QT_τQ: D → D, has the matrix representation:

$$P_Q=P{M}^{-1}$$

(10)

with the non-negative, invertible mass matrix, M ∈ ℝ^n,n, with entries:

$$M_{\mathit{ij}}=\frac{\langle {\chi }_i,{\chi }_j\rangle }{\langle {\chi }_i,\mathrm{𝟙}\rangle }$$

(11)

The matrix, P ∈ ℝ^n,n, is also non-negative and has entries:

$$P_{\mathit{ij}}=\frac{\langle {\chi }_i,T_{\tau }{\chi }_j\rangle }{\langle {\chi }_i,\mathrm{𝟙}\rangle }$$

(12)

Full Partition MSM. If we choose χ_i = 𝟙_Ai to be the characteristic function of set A_i for i = 1, ..., n, one can easily check that we get M = I to be the identity matrix and:

$$P_{\mathit{ij}}={\mathbb{P}}_{\mu }\left[X_{\tau }\in A_j|X_0\in A_i\right]$$

(13)

as in Equation (4). The subscript, μ, shall indicate that X₀ ∼ μ. Therefore, the transition probabilities are evaluated along equilibrium paths.

The previously constructed transition matrix of the MSM based on a full partition can be interpreted as a projection onto a space of densities that are constant on the partitioning sets. This interpretation of an MSM is useful, since it allows one to analyze its approximation quality. For example, in [25,26], it is proven that we can reproduce an eigenvalue, λ, of a self-adjoint transfer operator, T_t, by the MSM by choosing the subspace appropriately. That is, if u is a corresponding normalized eigenvector, Q the orthogonal projection to a subspace, D, with 𝟙 ∈ D, then there exists an eigenvalue, λ̂, of the projected transfer operator, QT_tQ, with:

$$|\lambda -\hat{\lambda }|\le {\lambda }_1\delta {(1-{\delta }^2)}^{-\frac{1}{2}}$$

where λ₁ < 1 is the largest non-trivial eigenvalue of T_t and δ = ‖u − Qu‖.

In particular, for $\delta \le \frac{3}{4}$, one can simplify the equation to:

$$|\lambda -\hat{\lambda }|\le 2{\lambda }_1\delta$$

(14)

An eigenvalue, λ_i, of the transfer operator directly relates to an implied timescale, $𝒯$_i, of the system via:

$${𝒯}_i=-\frac{\tau }{\text{log}\left({\lambda }_i\right)}$$

(15)

Therefore, the transition matrix Equation (4) that we construct from transitions between the sets, A₁, ..., A_n, will generate a Markov chain that will reproduce the original timescales well if the partitioning sets are chosen such that the corresponding eigenvectors are almost constant on these sets. In this case, δ = ‖u − Qu‖; that is, the approximation error of the eigenvector by a piecewise constant function on the sets will be small.

The projection error, δ, depends on our choice of the discretizing sets. As an example, let us consider a diffusion in the potential that is illustrated in Figure 2, that is, the reversible Markov process given by the stochastic differential equation:

$$d{X}_t=-\nabla V\left(X_t\right)\mathit{dt}+\sqrt{2\epsilon }d{B}_t$$

(16)

where V is the potential, B_t denotes a Brownian motion and ε > 0.

The figure also shows a choice of three sets that form a full partition of state space. The computation of the transition matrix Equation (4) for σ = 0.7 and a lag time τ = 1 yields:

$$P_Q=P=\left(\begin{array}{lll}0.9877\hfill & 0.0123\hfill & 0.0000\hfill \\ 0.0420\hfill & 0.9160\hfill & 0.0419\hfill \\ 0.0000\hfill & 0.0123\hfill & 0.9877\hfill \end{array}\right)$$

that has three eigenvalues λ₀ = 1, λ₁ = 0.9877, λ₂ = 0.9037.

Table 1. Comparison of implied timescales

	T1	T2
Original	103.7608	11.9566
Full partition 3 sets	80.6548	9.8784

shows the two resulting implied timescales Equation (15) in comparison to the timescales of the original system.

As one can see, the timescales are strongly underestimated. This is a typical phenomenon. From a statistical point of view, the recrossing problem will lead to cumulatively appearing transition counts when one computes the transition probabilities, ℙ_μ[X_τ ∈ A_j|X₀ ∈ A_i], from a trajectory (X_t), as discussed above. Therefore, on average, transitions between sets seem to become too likely, and hence, the processes in the coarse-grained system get accelerated. We have seen in Equation (14) that this cannot happen if the associated eigenvectors can be approximated well by the subspace that corresponds to the MSM. Figure 3 shows the first non-trivial eigenvector, u₁, belonging to the timescale T₁ = 103.7608 and its best-approximation by a step function.

The eigenvector is indeed almost constant in the vicinity of the wells, but within the transition region between the wells, the eigenvector is varying and the approximation by a step function is not accurate. Therefore, we have two explanations of why the main error is introduced in the region close to shared boundaries of neighboring sets: (1) because of recrossing issues; and (2) because of the main projection error of the associated eigenvector. Of course, one solution would be an adaptive refinement of the discretization, that is, one could choose a larger number of smaller sets, such that the eigenvector is better approximated by a step function on these sets. In the following section, we will present an alternative solution for overcoming the recrossing problem and reducing the projection error without refining the discretization.

4. The Core Set Approach

From Equation (10), we know how to compute a matrix representation for a projected transfer operator for an arbitrary subspace, D ⊂ L²(μ). For a given basis, χ₁, ..., χ_n, we have to compute Equations (11) and (12), so:

$$M_{\mathit{ij}}=\frac{\langle {\chi }_i,{\chi }_j\rangle }{\langle {\chi }_i,\mathrm{𝟙}\rangle },\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{P}_{\mathit{ij}}=\frac{\langle {\chi }_i,T_{\tau }{\chi }_j\rangle }{\langle {\chi }_i,\mathrm{𝟙}\rangle }$$

(17)

In general, the evaluation of these scalar products for arbitrary basis functions is a non-trivial task. On the other hand, we have seen that for characteristic functions χ_i = 𝟙_Ai on a full partition, we do not have to compute the scalar products numerically, since the matrix entries have a stochastic interpretation in terms of transition probabilities between set Equation (13). This means they can be directly estimated from a trajectory of the process, which is a strong computational advantage, particularly in high-dimensional state spaces.

Now, the question is if there is another basis other than characteristic functions that: (a) is more adapted to the eigenvectors of the transfer operator; and (b) still leads to a probabilistic interpretation of the matrix entries Equation (17), such that scalar products never have to be computed. The basic idea is to stick to a set-oriented definition of the basis, but to relax the full partition constraint. We will define our basis with respect to so-called core sets, C₁, ..., C_n ⊂ E, that are still disjoint, so C_i ∩ C_j = ∅, but they do not have to form a full partition. Figure 4 suggests that this could lead to a reduction of the recrossing phenomenon, since the sets do not share boundaries anymore.

Now, we use the core sets to define our basis functions, χ₁, ..., χ_n. Assume T_τ is, again, a self-adjoint transfer operator and consider n core sets C₁, ..., C_n. For every i, take the committor function, χ_i, of the process with respect to core set C_i; that is, χ_i(x) denotes the probability to hit the core set, C_i, next, rather than the other core sets, when starting the process in x. If we now study the projection, Q, onto the space spanned by these committor functions, the two following properties hold [25,27].

(P1)

The matrices, M and P, in Equation (10) can be written as:

$$M_{\mathit{ij}}={\mathbb{P}}_{\mu }\left[{\tilde{X}}_k^{+}=j|{\tilde{X}}_k^{-}=i\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{P}_{\mathit{ij}}={\mathbb{P}}_{\mu }\left[{\tilde{X}}_{k+1}^{+}=j|{\tilde{X}}_k^{-}=i\right]$$

(18)

where ( ${\tilde{X}}_k^{+}$) and ( ${\tilde{X}}_k^{-}$) are forward and backward milestoning processes [25,28]; that is, ${\tilde{X}}_k^{-}=i$ if the process came at time t = kτ, last from core set C_i and ${\tilde{X}}_k^{+}=j$ if the process went next to core set C_j after time t = kτ.

(P2)

Let u_i be an eigenvector of T_τ that is almost constant on the core sets. Let the region C = E\∪_iC_i that is not assigned to a core set be left quickly enough, so 𝔼_x[τ(C^c)] ≪ $𝒯$_i for all x ∈ C, where $𝒯$_i is the timescale associated with u_i and 𝔼_x[τ(C^c)] is the expected hitting time of C^c = ∪_iC_i when starting in x ∈ C. Then, ‖u_i − Qu_i‖ is small; so, the committor approximation to the eigenvector is accurate.

The message behind (P1) is that it is possible to relax the full partition constraint and use a core set discretization that does not cover the whole state space. We can still define a basis for a projection of the transfer operator that leads to a matrix representation that can be interpreted in terms of transition probabilities.

Important Remark: The construction of the projection onto the committors is only necessary for theoretical purposes. In practice, neither the committor functions nor scalar products between the committors have to be computed numerically, since the matrix entries of M and P can be estimated from trajectories again.

Property (P2) yields that the relaxation of the full partition constraint should also lead to an improvement of the MSM if the region, C, between the core sets is typically left on a faster timescale than the processes of interest taking place. Let us get back to the example from above. We will see that we can achieve a strong improvement of the approximation by simply excluding a small part of state space from our discretization. In Figure 5, we have turned our initial full partition into a core set discretization by removing parts of the transition region between the wells.

The matrix P_Q = PM⁻¹ that represents the projection, QT_τ Q, of the transfer operator onto the committor space associated with the core sets is given by:

$$P_Q=\left(\begin{array}{lll}0.9897\hfill & 0.0103\hfill & 0.0000\hfill \\ 0.0352\hfill & 0.9298\hfill & 0.0351\hfill \\ 0.0000\hfill & 0.0103\hfill & 0.9897\hfill \end{array}\right)$$

Comparing to the MSM for the full partition one can see that transitions between indices i and j, i ≠ j are less likely.

Table 2. More accurate approximation if implied timescales

	T1	T2
Original	103.7608	11.9566
3 core sets	100.8066	11.9145
Full partition 3 sets	80.6548	9.8784

shows this leads to a far more accurate reproduction of the timescales in the system.

From the discussion above, this has to be expected, because the eigenvectors are almost constant in the vicinity of the wells, and we removed a part of state space from the discretization that is typically left quickly compared to the timescales, T₁ and T₂. Therefore, the committor functions should deliver a good approximation of the first two eigenvectors. Figure 6 underlines this theoretical result.

5. Practical Considerations and MD Applications

In the previous sections, we have interpreted the construction of an MSM as a projection of the dynamics onto some finite dimensional ansatz space. We have discussed two types of spaces that both have been defined on the basis of a set discretization. First, we chose a full partition of state space and the associated space of step functions, and second, we analyzed a discretization by core sets and the associated space spanned by committor functions. These two methods have the advantage that the resulting projections lead to transition matrices for the MSM with entries that are given in terms of transition probabilities between the sets. That is, one can compute estimates for the transition matrices from simulation data. This is an important property for practical applications, because it means that we never need to compute committor functions or scalar products between committors or step functions. We rather generate trajectories x₀,x₁, ...x_N of the process (X_t), let us say, for a time step h > 0, so x_i = X_hi. For example, we can then define for a full partition, A₁, ..., A_m, and a lag time τ = nh the discrete trajectory s_k = i ⇔ x_k ∈ A_i and compute the matrix, P̂:

$${\hat{P}}_{\mathit{ij}}=\frac{C_{\mathit{ij}}}{\sum _j{C}_{\mathit{ij}}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{C}_{\mathit{ij}}=\sum _{k=0}^{N-n}{\mathrm{𝟙}}_{\{s_k=i\}}{\mathrm{𝟙}}_{\{s_{k+n}=j\}}$$

(19)

It is well known [7] that P̂ is a maximum likelihood estimator for the full partition MSM transition matrix Equation (4). Similarly one can also compute estimates for a core set MSM by using the definition of milestoning processes [27,28]. That is, if we have core sets C₁, ..., C_m, a lag time τ = nh as before, and we define discrete milestoning trajectories by:

$$\begin{array}{l}{s}_k^{-}=i\iff x_k\in A_i\hspace{0.17em}\text{or}\hspace{0.17em}\text{came}\hspace{0.17em}\text{last}\hspace{0.17em}\text{from}\hspace{0.17em}{A}_i\hspace{0.17em}\text{before}\hspace{0.17em}\text{time}\hspace{0.17em}k\\ s_k^{+}=i\iff x_k\in A_i\hspace{0.17em}\text{or}\hspace{0.17em}\text{went}\hspace{0.17em}\text{next}\hspace{0.17em}\text{to}\hspace{0.17em}{A}_i\hspace{0.17em}\text{after}\hspace{0.17em}\text{time}\hspace{0.17em}k\end{array}$$

we can compute an estimator P̂_Q = P̂M̂⁻¹ of the core set MSM matrix Equation (10) by counting transitions:

$${\hat{P}}_{\mathit{ij}}=\frac{C_{\mathit{ij}}}{\sum _j{C}_{\mathit{ij}}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{C}_{\mathit{ij}}=\sum _{k=0}^{N-n}{\mathrm{𝟙}}_{\{s_k^{-}=i\}}{\mathrm{𝟙}}_{\{s_{k+n}^{+}=j\}}$$

(20)

$${\hat{M}}_{\mathit{ij}}=\frac{N_{\mathit{ij}}}{\sum _j{N}_{\mathit{ij}}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{N}_{\mathit{ij}}=\sum _{k=0}^N{\mathrm{𝟙}}_{\{s_k^{-}=i\}}{\mathrm{𝟙}}_{\{s_k^{+}=j\}}$$

(21)

Since, in practice, we will only have a finite amount of data available, we will have statistical errors when constructing an MSM. This is an additional error to the projection error related to the discretization that we have discussed above. On the other hand, one should note that these errors are not independent of each other. For example, it is clear that if we take a full partition of state space, and we let the partition become arbitrarily fine by letting the number of sets go to infinity, the discretization error will vanish. At the same time, for a fixed amount of statistics, the statistical error will become arbitrarily large, because we will need to compute more and more estimators for transition events between the increasing number of sets. For more information on statistical errors, we refer to the literature [7,29].

Besides the choice of discretization and the available statistics, the estimates above also depend on a lag time, τ. This dependence can be used to validate an MSM by a Chapman–Kolmogorov test [7]. This is based on the fact that the MSM matrices approximately form a semi-group for all large enough lag times τ > τ^*; although, for small lag times, this is typically not true, due to memory effects. These facts also motivate one to look at something, like an infinitesimal generator, that approximately generates these MSM transition matrices for large enough lag times. In [27], two types of generator constructions have been compared for a core set setting. The first generator, K, is simply constructed from the transition rates between the core sets in the milestoning sense, that is:

$$K_{\mathit{ij}}=\underset{T\to \infty }{\text{lim}}\frac{N_{\mathit{ij}}^T}{R_i^T},i\ne j\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{K}_{\mathit{ii}}=-\sum _{j\ne i}{K}_{\mathit{ij}}$$

(22)

where $N_{\mathit{ij}}^T$ is the amount of time in [0,T] the process has spent on its way from core set C_i to C_j and $R_i^T$ is the total time in [0,T] the process came last from C_i. On the other hand, one can see [27,30] that K^* = KM⁻¹ with the mass matrix, M, from above Equation (18), can be interpreted as a projection of the original generator of the process and, also, as a derivative of the core set MSM from above, i.e.:

$$K^{*}=\underset{\tau \to 0}{\text{lim}}\frac{P{M}^{-1}-I}{\tau }$$

(23)

where P depends on τ Equation (17).

Let us now analyze how the choice of core sets, particularly the size of the core sets, influences the resulting approximation. Therefore, we consider an MD example that was discussed in [27], namely one molecule of alanine dipeptide monitored via its ϕ and ψ backbone dihedral angles. Two core sets are defined as balls with radius r around the two points with angular coordinates x_α = (−80, −60) and x_β = (−80, 170). The stationary distribution of the process and the two centers of the core sets, x_α,x_β, in the angular space are shown in Figure 7.

For computing a reference timescale, several MSMs based on three different full partitions using 10, 15 and 250 sets have been constructed for increasing lag times. In [27], it is shown that in each setting, the estimate for the longest implied timescale of the process converged to ≈19 ps for large enough τ. Now, the implied timescales for the two different generators, K Equation (22) and K^* Equation (23), are computed. In Figure 8, the resulting timescales are plotted against the reference timescale ≈ 19 ps for varying size of the core sets.

One can see that the estimate by the milestoning generator, K, is rather sensitive to the size of core sets. It overestimates the timescales for small core sizes and underestimates it for larger core sizes. On the other hand, the projected generator, K^*, can never overestimate the timescale, due to its interpretation as projection. It is also rather robust against the choice of size of the core sets until the core sets become too large, e.g., r > 15. Then, the discretization becomes close to a full partition discretization using only two sets. In this case, the timescales have to be underestimated heavily, because of recrossing phenomena. On the other hand, the underestimation for very small core sets has to be explained by a lack of statistics. When the core sets are chosen as arbitrarily small, it is clearly more difficult for the process to hit the sets, and therefore, transition events become rare. Note that for the straightforward milestoning generator, K, the processes seem to become very slow, but for the projected generator K^* = KM⁻¹, this effect is theoretically corrected by the mass matrix, M. Nevertheless, in both cases, the generation of enough statistics will be problematic for too small core sets.

6. Further Applications in MD

Markov State Modeling has been show to apply successfully to many different molecular systems, like peptides, including time-resolved spectroscopic experiments [10–12], proteins and protein folding [4,9,13], DNA [31] and ligand-receptor interaction [32]. In most of the respective publications, full partition MSMs are used, and the underlying discretization is based on cluster finding methods (see [7] for a review), while the sampling issues are tackled by means of ideas from enhanced sampling [33] and based on ensembles of rather short trajectories instead of one long one, cf.[4]. Core set-based approaches have been used just recently [10,27]; related algorithms are less well developed. However, recent work has shown that and how every full partition MSM can be easily transformed into a core set-based MSM with significantly improved approximation quality [34], making core set MSMs the most promising next generation MSM tools.

Very Rare Transitions between Discretization Sets. When constructing a full partition or a core set MSM, we have to estimate transition probabilities between sets in state space, and it can happen that we cannot avoid that some of these transitions are very rare. That is, the transition probabilities for a lag time, τ, between some sets can be non-zero, but small even, if compared to the remaining transition probabilities that are small already. This is why it is important to note that neglecting these very rare transitions during the construction of an MSM does not harm its approximation quality. For example, we can define for a transition matrix, P, another transition matrix, P̃, by:

$${\tilde{P}}_{\mathit{ij}}=\{\begin{array}{ll}{P}_{\mathit{ij}},\hfill & i\ne j,(i,j)\notin R\hfill \\ 0,\hfill & i\ne j,(i,j)\in R\hfill \end{array}$$

(24)

where R denotes the set of pairs of indices for which the transition are very rare and for which we set the transition probability to zero. If the Markov chain is reversible and (i, j) ∈ R ⇔ (j, i) ∈ R, one can show that for all ordered eigenvalues, λ_k(P) and λ_k(P̃), it holds that:

$$|{\lambda }_k(P)-{\lambda }_k(\tilde{P})|\le \underset{i}{\text{max}}\sum _{j\ne i,(i,j)\in R}{P}_{\mathit{ij}}$$

(25)

That is, if we cannot estimate a very small transition probability, P_ij, for a very rare transition event between two sets, A_i and A_j, and even totally neglect this probability by setting it to zero, the timescales of the MSM remain almost unaffected. Thus, if we compute the set of the “first order” transition probability of a system correctly enough and ignore all “higher order” ones, then our accuracy will not be spoiled. This nicely illustrates the main advantage of MSM modeling compared to classical long-term simulation: since only neighboring core sets have to be connected by accurately estimated rates, the long residence time of long-term trajectories between and in core sets can be avoided, thus cutting down total simulation time.

Computation from Trajectories. Clearly, constructing and analyzing a core set MSM will only have a computational advantage compared to the direct sampling of a rare event if the transition events between neighboring core sets occur on a much shorter timescale than the rare event itself. One should note that from the purely theoretical point of view, it would be optimal to choose only very few core sets in the most metastable regions of state space, because this would minimize the projection error δ = ‖u − Qu‖ for each dominant eigenvector u, as discussed in Section 3. On the other hand, when estimating the MSM from trajectories, only a finite amount of statistics will be available, so there will also be a statistical error. In order to keep the total error small, additional core sets in less metastable parts of state space typically have to be introduced. In the end, this makes the estimation of a core set MSM possible without having to sample rare events. Note that the projection error is still under control, as long as there is a transition region between the core sets that is typically left very quickly (see Property (P2) in Section 4).

In practice, the statistics of the transition events between core sets will preferably be estimated from many short trajectories using milestoning techniques [27,28] and parallel computing. However, any algorithm for the construction of a core set MSM has to find a balance between sampling issues (not too many too long trajectories needed) and discretization issues (not too many core sets). Construction of such an algorithm still is ongoing research.

This article cannot give a detailed review on the algorithmic realization of MSMs for realistic molecular systems and on the findings resulting from such applications, since this is discussed to some extent elsewhere; see [7] for a recent review of the algorithmic aspects and [32,35] for ligand-receptor interaction.

7. MSM for Optimal Control Problems

In this section, we will borrow ideas from the previous section and explain how MSMs can be used to discretize optimal control problems that are linear-quadratic in the control variables and which appear in, e.g., the sampling of rare events. Specifically, we consider the case that (X_t)_t≥0 is the solution of:

$$d{X}_t=\left(\sqrt{2}{u}_t-\nabla V\left(X_t\right)\right)\mathit{dt}+\sqrt{2\epsilon }d{B}_t$$

(26)

with potential V, Brownian motion B_t and temperature ε > 0, as in Equation (16), and an unknown control variable, u: [0, ∞) → ℝ^d, that is chosen so as to minimize the cost function:

$$J(u;x)=𝔼\left[{\int }_0^{\tau }\left(f\left(X_s\right)+\frac{1}{2}|u_t{|}^2\right)\mathit{ds}|X_0=x\right]$$

(27)

(The factors of 1/2 and $\sqrt{2}$ in front of the control terms are for notational convenience.) Here, f ≥ 0 is a bounded continuous function called running cost and τ < ∞ (almost surely) is a random stopping time that is determined by X_t hitting a given target set, A ⊂ E, i.e., τ = inf{t > 0: X_t ∈ A}, in other words, we are interested in controlling $X_t=X_t^u$ until it reaches A. As an example, consider the case f = 1 with the potential considered in Figure 5 and the target region, A, around the left well. This situation is illustrated in Figure 9 and amounts to the situation that one seeks to minimize the time to reach A by tilting the potential towards A; tilting the potential too much is prevented by the quadratic penalization term in the cost functional that grows when too much force is applied.

Other choices of f in Equation (26) result in alternative applications. One obvious application would be to set τ = T to a fixed time and f to the characteristic function of the complement of a conformation set C, f = 𝟙_E\C. In this case, minimization of J wrt. the control u_t would mean maximization of the probability to find the system in the conformation, C, until time T under a penalty on the external work done to the system. See [14] for more details on such applications.

There are other types of cost functions, J, one might consider, e.g., control until a deterministic finite time τ = T is reached or, even, τ → ∞, and the construction would follow analogously. For compactness, we consider here only cost functions as in Equation (27).

Optimal Control and Equilibrium Expectation Values. It turns out that when minimizing J, it is sufficient to consider control strategies that are Markovian and depend only on X_t, i.e., we consider feedback laws of the form u_t = α(X_t) for some smooth function, α: E → ℝ^d. Moreover, only controls with finite energy are considered, for otherwise, J(u; x) = ∞. For control problems of the form Equations (26) and (27), the optimal feedback function can be shown to be α^*(x) = − $\sqrt{2}$∇W, where W is the value function or optimal-cost-to-go [1,15]:

$$W(x)=\underset{u}{\text{min}}J(u;x)$$

(28)

with the minimum running over all admissible Markovian feedback strategies. It can be shown that W satisfies the following dynamic programming equation of the Hamilton–Jacobi–Bellman type (see [36]):

$$\begin{array}{r}\hfill \mathit{LW}(x)-|\nabla W(x){|}^2+f=0\\ \hfill W{|}_A=0\end{array}$$

(29)

with the second-order differential operator:

$$L=\epsilon \Delta -\nabla V\cdot \nabla$$

that is the infinitesimal generator of the process, X_t, for u = 0. If the value function, W, is known, it can be plugged into the equation of motion, which then turns out to be of the form:

$$d{X}_t^{*}=-\nabla U\left(X_t^{*}\right)\mathit{dt}+\sqrt{2\epsilon }d{B}_t$$

(30)

with the new potential:

$$U(x)=V(x)+2W(x)$$

The difficulty is that Equation (29) is a nonlinear partial differential equation and for realistic high-dimensional systems, it is not at all obvious how to discretize it, employing any kind of state space partitioning. It has been demonstrated in [14,15] that Equation (29) can be transformed into a linear equation by a logarithmic transformation. Setting: W (x) = −ε log ϕ(x) it readily follows, using chain rule and Equation (29), that ϕ solves the linear equation:

$$\begin{array}{r}\hfill \left(L-{\epsilon }^{-1}f\right)\varphi =0\\ \hfill \varphi {|}_A=1\end{array}$$

(31)

The last equation is linear and can be solved by using MSMs, as we will show below. Moreover, by the Feynman–Kac theorem [37], the solution to Equation (31) can be expressed as:

$$\varphi (x)=\text{𝔼}\left[\text{exp}\left(-\frac{1}{\epsilon }{\int }_0^{\tau }f\left(X_t\right)\mathit{dt}\right)|X_0=x\right]$$

(32)

where X_t solves the control-free equation:

$$d{X}_t=-\nabla V\left(X_t\right)\mathit{dt}+\sqrt{2\epsilon }d{B}_t$$

That is, the optimal control for Equation (26) can be computed by solving Equation (31), which can be done in principle via Monte Carlo approximation of the expected value in Equation (32) if critical slowing down by rare events can be avoided.

Remark. The optimization problem Equation (28) admits an interpretation in terms of entropy minimization: let $Q=Q_x^u$ and $P=Q_x^0$ denote the path probability measures of controlled and uncontrolled trajectories starting at x at time t = 0, and set:

$$Z={\int }_0^{\tau }f\left(X_s\right)\mathit{ds}$$

Then, it follows that we can write:

$$W(x)=\underset{Q\ll P}{\text{min}}J(u;x),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}J(u;x)=\int \left\{Z+\epsilon \text{log}\left(\frac{dQ}{dP}\right)\right\}dQ$$

(33)

where the notation “Q ≪ P” means that Q has a density (That is, the density function, dQ/dP, exists and is almost everywhere positive and normalized) with respect to P. It turns out that for every such Q, there is exactly one control strategy, u, such that $Q=Q_x^u$ is generated by Equation (26); in this sense, the notation in Equation (33) is meaningful. The second term:

$$H(Q\Vert P)=\epsilon \int \text{log}\left(\frac{dQ}{dP}\right)dQ$$

is the relative entropy or Kullback–Leibler divergence between Q and P. For details on this matter that are based on Girsanov transformations for stochastic differential equations, we refer to [38]orthe article [1] in this special issue.

8. MSM Discretization of Optimal Control Problems

The basic idea is now to choose a subspace, D ⊂ L²(μ), with basis χ₁,...,χ_n as in Markov State Modeling and then discretize the dynamic programming Equation (29) of our optimal control problem by projecting the equivalent log transformed Equation (31) onto that subspace. As we will see, the resulting discrete matrix equation can be transformed back into an optimal control problem for a discrete Markov jump process (MJP).

We will do this construction for the full partition case χ_i = 𝟙_Ai and the core set case χ_i = q_i discussed earlier. We will see that in both cases, we arrive at a structure-preserving discretization of the original optimal control problem, where the states of the corresponding MJP will be related to the partition subsets, A_i. The first case will give us back a well-known lattice discretization for continuous control problems, the Markov chain approximation [39]. This is illustrated in the following diagram:

Subspace Projection. The key step for the discretization is that we pick a suitable subspace, D ⊂ L²(μ), that is adapted to the boundary value problem Equation (31). Specifically, we require that the subspace contains the constant function, 𝟙 ∈ D, and that it gives a good representation of the most dominant metastable sets. To this end, we choose basis functions χ₁,...,χ_n+1 with the following properties:

(S1)

The χ_i form a partition of unity, that is ${\sum }_{i=1}^{n+1}{\chi }_i=\mathrm{𝟙}$.

(S2)

The χ_i are adapted to the boundary conditions in Equation (31), that is χ_n+1|_A = 1 and χ_i|_A = 0 for i ∈ {1,...,n}.

Now, let Q be the orthogonal projection onto D, and define the matrices:

$$F_{\mathit{ij}}=\frac{\langle {\chi }_i,f{\chi }_j\rangle }{\langle {\chi }_i,\mathrm{𝟙}\rangle },\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{K}_{\mathit{ij}}=\frac{\langle {\chi }_i,L{\chi }_j\rangle }{\langle {\chi }_i,\mathrm{𝟙}\rangle }$$

Now, if ϕ solves the linear boundary value problem Equation (31), then the coefficients, ϕ̂₁,…,ϕ̂_n+1, of its finite-dimensional representation Qϕ = ∑_jϕ̂_jχ_j on the subspace, D, satisfy the constrained linear system:

$$\begin{array}{c}\sum _{j=1}^{n+1}\left(K_{\mathit{ij}}-{\epsilon }^{-1}{F}_{\mathit{ij}}\right){\hat{\varphi }}_j=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i\in \{1,\dots ,n\}\\ {\hat{\varphi }}_{n+1}=1\end{array}$$

(34)

that is the discrete analogue of Equation (31). The discrete solution ϕ̂ = Qϕ is optimal in the sense of being the best approximation of ϕ in the energy norm, i.e.:

$$\Vert \varphi -\hat{\varphi }{\Vert }_A=\underset{\psi \in D}{\text{inf}}\Vert \varphi -\psi {\Vert }_A$$

(35)

where

$$\Vert \varphi {\Vert }_A^2=\langle \varphi ,\left({\epsilon }^{-1}f-L\right)\varphi \rangle$$

is the energy norm on L²(μ), and the infimum runs over all functions, ψ ∈ L²(μ), that are of the form ψ(x) = ∑_jψ_jχ_j(x) with coefficients ψ_j ∈ ℝ. This is a standard result about projections of PDEs; see [40] for details. (By the same argument as in the previous sections, A = ε⁻¹f − L is symmetric and positive definite as an operator on the weighted Hilbert space, L²(μ). Moreover, ${\Vert \varphi \Vert }_A^2={\epsilon }^{-1}\langle \varphi ,f\varphi \rangle +\epsilon \langle \nabla \varphi ,\nabla \varphi \rangle .)$ In analogy with Equation (14), we can use the above result to get the error estimate:

$$\Vert \varphi -\hat{\varphi }{\Vert }_{\mu }^2\le \left(1+\frac{1}{{\delta }^2}\Vert {\mathit{QAQ}}^{\perp }{\Vert }^2\right)\underset{\psi \in D}{\text{inf}}\Vert \varphi -\psi {\Vert }_{\mu }^2$$

(36)

where A = ε⁻¹f − L is a shorthand for the operator appearing in Equation (31) and the constant δ > 0 is defined, such that ${\Vert v\Vert }_A^2\ge \delta {\Vert v\Vert }_{\mu }^2$ holds for all v ∈ L²(μ); see [41]. The bottom line of Equation (35) shows that discretizing Equation (31) via Equation (34) minimizes the projection error measured in the energy norm. Since all functions are μ-weighted, the approximation will be good in regions visited with high probability and less good in regions with lower probability. The error estimate Equation (36) is along the lines of the MSM approximation result: if we switch to the norm on L²(μ), the function ϕ̂ = Qϕ is still almost the best approximation of ϕ, provided that A leaves the subspace, D, almost invariant. As was pointed out earlier, this is exactly the case when the χ_i are close to the eigenfunctions of A (e.g., when the system is metastable).

The best approximation error ‖Q^⊥ϕ‖_μ = inf_ψ∈D‖ϕ − ψ‖_μ, which appears in Equation (36), will vanish if the χ_i form an arbitrarily fine full partition of E. If we follow the core set idea from Section 4 and choose the χ_i to be committor functions, we have good control over ‖Q^⊥ϕ‖_μ. Due to [41]:

$$\Vert Q^{\perp }\varphi {\Vert }_{\mu }\le \Vert P^{\perp }\varphi {\Vert }_{\mu }+\mu {(C)}^{1/2}\left[\kappa \Vert f{\Vert }_{\infty }+2\Vert P^{\perp }\varphi {\Vert }_{\infty }\right]$$

(37)

where C = E \ ∪_iC_i is the transition region, κ = sup_x∈C 𝔼_xτ_E\C is the maximum expected time of hitting the metastable set from outside (which is short) and P is the orthogonal projection onto the subspace V = {v ∈ L²(μ), v = const on every C_i} ⊂ L²(μ). Note that P^⊥ϕ = 0 on C. The errors, ‖P^⊥ϕ‖_μ and ‖P^⊥ϕ‖_∞, measure how constant the solution, ϕ, is on the core sets. Hence, Equation (37) together with Equation (36) gives us complete control over the approximation error of our projection method, even if we consider just a few core sets. In Section 9, we will investigate the full and core set partition cases further.

Properties of the Projected Problem. We introduce now the diagonal matrix, Λ, with entries Λ_ii = ∑_jF_ij (zero otherwise) and the full matrix G = K − ε⁻¹(F − Λ), and rearrange Equation (34) as follows:

$$\begin{array}{c}\sum _{j=1}^{n+1}\left(G_{\mathit{ij}}-{\epsilon }^{-1}{\Lambda }_{\mathit{ij}}\right){\hat{\varphi }}_j=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i\in \{1,\dots ,n\}\\ {\hat{\varphi }}_{n+1}=1\end{array}$$

(38)

This equation can be given a stochastic interpretation. To this end, let us introduce the vector, π ∈ ℝⁿ⁺¹, with nonnegative entries π_i = 〈χ_i, 𝟙〉 and notice that ∑_iπ_i = 1 follows immediately from the fact that the basis functions, χ_i, form a partition of unity, i.e., ∑_iχ_i =𝟙. This implies that π is a probability distribution on the discrete state space Ê = {1,...,n + 1}. We summarize properties of the matrices, K, F and G; see also [41]:

(M1)

K is a generator matrix of an MJP (X̂_t)_t≥0 (i.e., K is a real-valued square matrix with row sum zero and positive off-diagonal entries) with stationary distribution, π, that satisfies detailed balance:

$${\pi }_i{K}_{\mathit{ij}}={\pi }_j{K}_{\mathit{ji}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,j\in \hat{E}$$

(M2)

F ≥ 0 (entry-wise) with π_iF_ij = π_jF_ji for all i, j ∈ Ê.

(M3)

G has a row sum of zero and satisfies π^T G = 0 and π_iG_ij = π_jG_ji for all i, j ∈ Ê; furthermore, there exists a constant 0 < C < ∞, such that G_ij ≥ 0 for all i ≠ j if ‖f‖_∞ ≤ C. In this case, Equation (38) admits a unique and strictly positive solution ϕ̂ > 0.

It follows that if the running costs, f, are such that (M3) holds, then G is a generator matrix of an MJP that we shall denote by (X̂_t)_t≥0, and Equation (38) has a unique and positive solution. In this case, the logarithmic transformation Ŵ = −ε log ϕ̂ is well defined. It was shown in [42] that Ŵ can be interpreted as the value function of a Markov decision problem with cost functional (cf. also [36]):

$$\hat{J}(v;i)=𝔼\left[{\int }_0^{\tau }\left(\hat{f}\left({\hat{X}}_s\right)+k\left({\hat{X}}_s,v_s\right)\right)\mathit{ds}|{\hat{X}}_0=i\right]$$

(39)

that is minimized over the set of Markovian control strategies, v : Ê → (0, ∞), subject to the constraint that the controlled process ${\hat{X}}_t={\hat{X}}_t^v$ is generated by G^v, where:

$$G_{\mathit{ij}}^v=\{\begin{array}{cc}v{(i)}^{-1}{G}_{\mathit{ij}}v(j),& i\ne j\\ -{\sum }_{j\ne i}{G}_{\mathit{ij}}^v,& i=j\end{array}$$

(40)

with stopping time τ = inf{t > 0 : X̂_t = n +1} and running costs:

$$\hat{f}(i)={\Lambda }_{ii},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k(i,v)=\epsilon \sum _{j\ne i}{G}_{\mathit{ij}}\left\{\frac{v(j)}{v(i)}\left[\text{log}\frac{v(j)}{v(i)}-1\right]+1\right\}$$

(41)

Properties of the Projected Problem, Continued. From [42], we know that the optimal cost:

$$\hat{W}(i)=\underset{v}{\text{min}}\hat{J}(v;i)$$

is given by Ŵ = −∊log ϕ̂, where ϕ̂ solves Equation (38), with the optimal feedback strategy given by v^*(i) = ϕ̂_i (see [36]). We list additional properties:

(i)

The v-controlled system has the unique invariant distribution:

$${\pi }^v=\left({\pi }_1^v,\dots ,{\pi }_{n+1}^v\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_i^v=\frac{v{(i)}^2{\pi }_i}{Z_v}$$

with Z_v an appropriate normalization constant; in terms of the value function, π^* = π^v* reads:

$${\pi }^{*}=\left({\pi }_1^{*},\dots ,{\pi }_{n+1}^{*}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_i^{*}=\frac{1}{Z_{*}}{e}^{-2{\epsilon }^{-1}\hat{W}(i)}{\pi }_i$$

(ii)

G^v is reversible and stationary with respect to π^v, i.e., ${\pi }_i^v{G}_{\mathit{ij}}^v={\pi }_j^v{G}_{\mathit{ji}}^v$ for all i, j ∈ Ê.

(iii)

Ĵ admits the same interpretation as Equation (33) in terms of the relative entropy:

$$\hat{W}(i)=\underset{Q\ll P}{\text{min}}\hat{J}(v;i),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hat{J}(v;i)=\int \left\{\hat{Z}+\epsilon \text{log}\left(\frac{dQ}{dP}\right)\right\}dQ$$

where P denotes expectation with respect to the uncontrolled MJP, X̂_t, starting at X̂₀ = i, Q denotes the path measure of the corresponding controlled process with generator G^v and:

$$\hat{Z}={\int }_0^{\tau }\hat{f}\left({\hat{X}}_s\right)\mathit{ds}$$

A few remarks seem in order: Item (i) of the above list is in accordance with the continuous setting, in which the optimally controlled dynamics is governed by the new potential U = V +2W and has the stationary distribution, μ^* ∝ exp(−2∊⁻¹W)μ, with μ being the stationary distribution of the uncontrolled process. Hence, the effect of the control on the invariant distribution is the same in both cases. Further, note that optimal strategies change the jump rates according to:

$$G_{\mathit{ij}}^{v*}=G_{\mathit{ij}}{e}^{-{\epsilon }^{-1}\left(\hat{W}(j)-\hat{W}(i)\right)}$$

(42)

that is, Ŵ acts as an effective potential as in the continuous case, and the change in the jump rates can be interpreted in terms of Kramer’s law for this effective potential.

This completes our derivation of the discretized optimal control problem, and we now compare it with the continuous problem we started with for the case of a full partition of E and a core set partition of E.

9. Markov Chain Approximations and Beyond

Full Partitions. Let E be fully partitioned into disjoint sets, A₁,...,A_n+1, with centers x₁,...,x_n+1 and such that A_n+1 ≔ A, and define χ_i ≔ χ_Ai. These χ_i satisfy Assumptions (S1) and (S2) discussed in Section 8. Since they are not overlapping, F is diagonal, and:

$$\hat{f}(i)=\frac{1}{{\pi }_i}{\int }_{A_i}f(x)\mu (x)\mathit{dx}={𝔼}_{\mu }\left[f\left(X_t\right)|X_t\in A_i\right]$$

(43)

is just obtained by averaging f(x) over the cell, A_i. Equation (43) is also a sampling formula for f̂(i). It follows directly that G = K, and in particular, (M3) holds for any f. One can show that K has components:

$$K_{\mathit{ij}}\approx \frac{1}{{\Delta }_{\mathit{ij}}}{e}^{-\beta \left(V\left({\overline{x}}_{\mathit{ij}}\right)-V\left(x_i\right)\right)},{\Delta }_{\mathit{ij}}^{-1}={\beta }^{-1}\frac{m\left(S_{\mathit{ij}}\right)}{m\left(h_{\mathit{ij}}\right)m\left(A_i\right)}$$

(44)

if i and j are neighbors (K_ij = 0 otherwise). Here, m is the Lebesgue measure, and h_ij, S_ij and x̄_ij are defined as in Figure 10. K is the generator of an MJP on the cells, A_i, and coincides with the so-called finite volume approximation of L discussed in [43]. It is reversible with stationary distribution:

$${\pi }_i={\int }_{A_i}d\mu \approx m\left(A_i\right)e^{-\beta V\left(x_i\right)}$$

One can show that the approximation error vanishes for n → ∞. K and π can be computed from the potential, V, and the geometry of the mesh. By inspecting Equations (12) and (13), we see that K is connected to the transition matrix, P^τ, of a full partition MSM with lag time τ by

$$\underset{\tau \to 0}{\text{lim}}\frac{1}{\tau }\left(P_{\mathit{ij}}^{\tau }-M_{\mathit{ij}}\right)\underset{\tau \to 0}{\text{lim}}\frac{1}{{\pi }_i}\langle {\chi }_i,\frac{1}{\tau }\left(T_{\tau }-\mathrm{𝟙}\right){\chi }_j\rangle =\frac{1}{{\pi }_i}\langle {\chi }_i,L_{{\chi }_j}\rangle =K_{\mathit{ij}},$$

thus K is the generator of the semigroup of transition matrices, P^τ. Therefore we could obtain K by sampling in the same way we obtained P^τ through Equation (19) in Section 5. This is difficult, however, due to recrossing problems for small τ; see e.g., [44]. Finally, let us note in passing that we can drastically simplify k^v if the cells, A_i, are boxes of length h. Denote the elementary lattice vectors by e_n. Then,

$$k^v(i)=\frac{1}{2}{|u^v(i)|}^2+𝒪(h),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{u}_n^v(i):=\frac{1}{\sqrt{2}}\frac{∊}{2h}\left(\text{log}v\left(i+e_n\right)-\text{log}v\left(i-e_n\right)\right),$$

which establishes the connection to the continuous case. However, more is true: The whole discrete control problem reduces to first order in h to the well-known Markov chain approximation (MCA) [39], which allows us to use convergence theory for MCAs to conclude that, for n → ∞, the optimal control and value function of the discrete control problem converge to their continuous counterparts. More details can be found in [41].

Core Set Partition. Now, we choose core sets C₁,...,C_n+1 with C_n+1 = A, and we let χ_i = q_i be the committor function of the process with respect to C_i, as in Section 4. These χ_i satisfy Assumptions (S1) and (S2) discussed in Section 8. The projection onto the committor basis also allows for a stochastic interpretation. Recall the definition of the forward and backward milestoning process, ${\tilde{X}}_t^{\pm }$, from Equation (18). The discrete costs can be written as:

$$\hat{f}(i)=\frac{1}{{\pi }_i}\langle q_i,f\sum _j{q}_j\rangle =\int v_i(x)f(x)\mathit{dx}={𝔼}_{\mu }\left[f\left(X_t\right)|{\tilde{X}}_t^{-}=i\right]$$

(45)

where ${\nu }_i(x)=\frac{q_i(x)\mu (x)}{{\pi }_i}=\mathbb{P}\left(X_t=x|{\tilde{X}}_t^{-}=i\right)$ is the probability density of finding the system in state x given that it came last from i. Hence, f̂(i) is the average costs conditioned on the information ${\tilde{X}}_t^{-}=i$, i.e., X_t came last from A_i, which is the natural extension to the full partition case, where f̂(i) was the average costs conditioned on the information that X_t ∈ A_i.

The matrix $K={\pi }_i^{-1}\langle q_i,L{q}_j\rangle$ is reversible with stationary distribution

$${\pi }_i=\langle q_i,\mathrm{𝟙}\rangle ={\mathbb{P}}_{\mu }\left({\tilde{X}}_t^{-}=i\right)$$

and is related to core MSMs again:

$$K=\underset{\tau \to 0}{\text{lim}}\frac{1}{\tau }\left(P^{\tau }-M\right)$$

where P^τ and M are now the matrices for core MSMs, as in Equation (18). Formally, K is the generator of the P^τ, but these do not form a semigroup, since M ≠ 𝟙. Therefore, we cannot interpret K directly as, e.g., the generator of ${\tilde{X}}_t^{-}$. Nevertheless, the entries of K are the transition rates between the core sets, as defined in transition path theory [45]. We can sample P^τ and M using Equations (20) and (21); because we used an incomplete partition, the recrossing problem is removed, and there is no difficulty in sampling P^τ for all lag times, τ, and therefore, K directly. It is worth noting that F can also be sampled:

$$F_{\mathit{ij}}={𝔼}_{\mu }\left[f\left(X_t\right){\chi }_{\left\{{\tilde{X}}_t^{+}=j\right\}}|{\tilde{X}}_t^{-}=i\right]$$

Therefore, as in the construction of core MSMs, we do not need to compute committor functions explicitly. Note, however, that G ≠ L, there is a reweighting, due to the overlap of the q_i’s, which causes F to be non-diagonal. This reweighting is the surprising bit of this discretization. From properties (M1)–(M3) from Section 8, we see, however, that G and K are both reversible with stationary distribution, π. Finally, note that if the cost function, f(x), does not satisfy ‖f‖_∞ ≤ C from (M3), G will not even be a generator matrix. In this case, (34) still has a solution, ϕ̂, which is the best approximation to ϕ, but this solution may not be unique; it may not satisfy ϕ̂ > 0, and we have no interpretation as a discrete control problem.

10. Numerical Results

10.1. 1D Potential Revisited

Firstly, we study diffusion in the triple well potential, which is presented in Figure 2. This potential has three minima at approximately x_0/1 = ±3.4 and x₂ = 0. We choose the three core sets C_i = [x_i − δ, x_i + δ] around the minima with δ = 0.2. Take τ to be the first hitting time of C₀. We are interested in the moment generating function ϕ(x) = 𝔼[e^{−∊−1στ}] of passages into C₀ and the cumulant generating function W = ∊log ϕ. This is of the form Equation (32) for A = C_i and f = σ, a constant function.

In Figure 11a, the potential, V, and effective potential, U, are shown for β = 2 and σ = 0.08 (solid lines), cf. Equation (30). One can observe that the optimal control effectively lifts the second and third well up, which means that the optimal control will drive the system into C₀ very quickly. The reference computations here have been carried out using a full partition FEM (finite element method) discretization of Equation (31) with a lattice spacing of h = 0.01. Now, we study the MJP approximation constructed via the committor functions shown in Figure 11b. These span a three-dimensional subspace, but due to the boundary conditions, the subspace, D, of the method is actually two-dimensional. The dashed line in Figure 11a gives the approximation to U calculated by solving Equation (38). We can observe extremely good approximation quality, even in the transition region. In Figure 11c, the approximation to the optimal control, α^*(x) (solid line), and its approximation α̂^* = −$\sqrt{2}$▿Ŵ (dashed line) are shown. The core sets are shown in blue. We can observe jumps in α̂^* at the left boundaries of the core sets. This is to be expected and comes from the fact that the committor functions are not smooth at the boundaries of the core sets, but only continuous. Therefore, the approximation to U is continuous, but the approximation to α^* is not.

Next, we construct a core MSM to sample the matrices, K and F. One hundred trajectories of length T = 20, 000 were used to build the MSM. In Figure 11d, W and its estimate using the core MSM are shown for ∊ = 0.5 and different values of σ. Each of the 100 trajectories has seen about four transitions. For comparison, a direct sampling estimate of W using the same data is shown (green). The direct sampling estimate suffers from a large bias and variance and is practically useless. In contrast, the MSM estimator for W performs well for all considered values of σ, and always, its variance is significantly small. The constant, C, which ensures ϕ̂ > 0 when σ ≤ C, is approximately 0.2 in this case. This seems restrictive, but still allows one to capture all interesting information about ϕ and W.

10.2. Alanine Dipeptide

Lastly, we study α-β-transitions in alanine dipeptide, a well-studied test system for molecular dynamics applications. We use a 1μs long trajectory simulated with the CHARMM (Chemistry at HARvard Molecular Mechanics) 27 force field. The conformational dynamics is monitored as usual via the backbone dihedral angles, ϕ and ψ. The data was first presented in [27]. We construct a full partition MSM with 250 clusters using k-means clustering. We are interested in the MFPT (mean first passage time) t̂(i) = 𝔼_i[τ_α], where τ_α is the first hitting time of the α conformation, which we define as a circle with radius r = 45 around (ϕ_α,ψ_α) = (−80, −60). The MFPT vector, t̂, solves the boundary value problem

$$K\hat{t}=-1\hspace{0.17em}\text{outside}\hspace{0.17em}\text{of}\hspace{0.17em}\alpha ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hat{t}=0\hspace{0.17em}\text{in}\hspace{0.17em}\alpha$$

, but since K is not available directly via sampling, we have to consider the equation

$$\frac{1}{\tau }\left(P^{\tau }-1\right)\hat{t}=-1\hspace{0.17em}\text{outside}\hspace{0.17em}\text{of}\hspace{0.17em}\alpha ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hat{t}=0\hspace{0.17em}\text{in}\hspace{0.17em}\alpha$$

instead. The result will depend on the choice of lag time τ. In Figure 12a, the results are shown for τ = 5; we can identify the β-structure as the red cloud of clusters where t̂(i) is approximately constant. In Figure 12b, t̂_βα = 𝔼(t̂(i)|i ∈ β) is shown as a function of τ. We observe a linear behavior for large τ, which is due to the linear error introduced in the replacement of K with $\frac{1}{\tau }$ (P^τ − 1), and a nonlinear drop for small τ, which is due to non-Markovianity. Our best guess is, therefore, a linear interpolation to τ = 0, which is indicated by the solid line. The result is ${\hat{t}}_{\beta \alpha }^{(0)}$ = 35.5ps. As a comparison, the reference value ${\hat{t}}_{\beta \alpha }^{\mathit{ref}}$ = 36.1ps from [27] is shown as a dashed line. It was computed in [27] as an inverse rate, using the slowest ITS (implied time scale) and information about the equilibrium weights of the α and β structure. We see very good agreement. The result is, of course, dependent, though, on the assignment of clusters to the α and β structure. Some tests show that ${\hat{t}}_{\beta \alpha }^{(0)}$ as computed with the interpolation method is fairly insensitive to this choice.

In [14], it is demonstrated how to use the method presented herein for maximizing the population of the α-conformation of alanine dipeptide based on the MSM used here.

11. Conclusions

In this article, we have discussed an approach to overcome direct sampling issues of rare events in molecular dynamics based on spatial discretization of the molecular state space. The strategy is to define a discretization by subsets of state space, such that the sampling effort with respect to transitions between the sets is much lower than the direct estimation of the rare events under consideration. That is, without having to simulate rare events, we construct a so-called Markov State Model, a Markov chain approximation to the original dynamics. Since the state space of the MSM is finite, we can then calculate the properties of interest by simply solving linear systems of equations. Of course, it is crucial that these properties of the MSM can be related to the rare event properties of the original process that we have not been able to sample directly.

This is why we have analyzed the approximation quality of MSMs in the first part of the article. We have used the interpretation of MSMs as projections of the transfer operator to: (1) derive conditions that guarantee an accurate reproduction of the dynamics; and (2) show how to construct models based on a core set discretization by leaving the state space partly undiscretized.

In the second part of the article, we have used the concept of MSM discretization to solve MD optimal control problems in which one computes the optimal external force that drives the molecular system to show an optimized behavior (maximal possible population in a conformation; minimal mean first passage time to a certain conformation) under certain constraints. We have demonstrated that the spatial discretization underlying an MSM turns the high-dimensional continuous optimal control problem into a rather low-dimensional discrete optimal control problem of the same form that can be solved efficiently. This result allows two different types of applications: (1) if one can construct an MSM for a molecular system in equilibrium, then one can use it to compute optimal controls that extremize a given costs criterion; (2) if an MSM can be computed based on transition probabilities between neighboring core sets alone, then the rare event statistics for transitions between strongly separated metastable states of the system can be computed from an associated optimal control problem that can be solved after discretization using the pre-computed MSM.