Non-equilibrium Markov state modeling of periodically driven biomolecules

TL;DR

This paper introduces a method for constructing Markov state models for biomolecules under periodic external driving, enabling the study of non-equilibrium dynamics beyond traditional equilibrium-based models.

Contribution

The paper presents a novel approach to build Markov state models for periodically driven biomolecular systems, extending modeling capabilities to non-equilibrium conditions.

Findings

Successfully applied to alanine dipeptide under electric field

Demonstrates modeling of non-equilibrium biomolecular dynamics

Extends Markov state modeling to driven systems

Abstract

Molecular dynamics simulations allow to study the structure and dynamics of single biomolecules in microscopic detail. However, many processes occur on time scales beyond the reach of fully atomistic simulations and require coarse-grained multiscale models. While systematic approaches to construct such models have become available, these typically rely on microscopic dynamics that obey detailed balance. In vivo, however, biomolecules are constantly driven away from equilibrium in order to perform specific functions and thus break detailed balance. Here we introduce a method to construct Markov state models for systems that are driven through periodically changing one (or several) external parameter. We illustrate the method for alanine dipeptide, a widely used benchmark molecule for computational methods, exposed to a time-dependent electric field.

Figures

Equations

$${\dot{p}}_j(t)=\sum _i{p}_i(t)W_{ij}(t),$$

$$\dot{\mathbf{K}}(t)=\mathbf{K}(t)\cdot \mathbf{W}(t)$$

$${\bar{p}}_i\equiv \frac{1}{T}{\int }_t^{t+T}\mathrm{d}s{p}_i(s),$$

$${\bar{\Phi }}_{ij}\equiv \frac{1}{T}{\int }_t^{t+T}\mathrm{d}s{p}_i(s)W_{ij}(s),$$

$$\mathbf{K}(t+nT)=\mathbf{K}(t)\cdot [\mathbf{K}(T){]}^n$$

$$\bar{\mathbf{p}}\cdot \mathbf{P}=\frac{1}{T}{\int }_t^{t+T}\mathrm{d}s\mathbf{p}(s+T)=\bar{\mathbf{p}}$$

$${\Phi }_{ij}\equiv p_i^{\text{s}}{\tilde{W}}_{ij},$$

$$\mathbf{P}=\mathbf{K}(T)\approx 1+\sum _{k=0}^{N-1}{\int }_{t_k}^{t_{k+1}}\mathrm{d}s\mathbf{K}(s)\cdot {\mathbf{W}}^{(k)}(s)$$

$$\tilde{\mathbf{W}}$$

$$\approx \frac{1}{T}\left[(\mathbf{P}-1)-\frac{(\mathbf{P}-1{)}^2}{2}+\frac{(\mathbf{P}-1{)}^3}{3}-...\right].$$

$$\dot{S}=\frac{1}{2}\sum _{ij}{J}_{ij}{A}_{ij}⩾0,$$

$${\Phi }_{ij}={\Phi }_{ij}^0+J_{ij}={\Phi }_{ij}^0+\sum _{\alpha \ni (i\to j)}{\phi }_{\alpha }$$

$${\dot{S}}_l=\sum _{\alpha \in C_l}{A}_{\alpha }{\phi }_{\alpha },$$

$${\dot{S}}_l\stackrel{!}{=}{A}_l{\phi }_l^{\text{cg}}\Rightarrow {\phi }_l^{\text{cg}}=\frac{{\dot{S}}_l}{A_l}$$

$${\Phi }_{ij}^{0,\text{cg}}=\frac{{\sum }_{l\ni (i\to j)}{\phi }_l^{\text{cg}}}{e^{A_{ij}}-1}.$$

$$W_{ij}^{\text{cg}}=\frac{{\Phi }_{ij}^{\text{cg}}}{p_i^{\text{cg}}}.$$

$$\left(\begin{array}{c}\varphi \\ \psi \end{array}\right)⟼\left(\begin{array}{c}\cos \varphi \\ \sin \varphi \\ \cos \psi \\ \sin \psi \end{array}\right).$$

$${\chi }_{\alpha }^i=\{\begin{array}{ll}1& \text{if state }i\text{ is in cycle }\alpha \\ 0& \text{otherwise}\end{array}$$

$$c_{\alpha }^{(i)}\equiv \left(\frac{1}{|\alpha |}\sum _k{R}_k^{(i)}{\chi }_{\alpha }^k\right),$$

$$d^{(i)}_{\alpha}\equiv\max_{k,l}\Bigg{\{}\chi^{k}_{\alpha}\chi^{l}_{\alpha}~{}\Big{|}R^{(i)}_{k}-R^{(i)}_{l}\Big{|}\Bigg{\}}.$$

$$P_X\equiv \sum _{i\in X}{p}_i^{\text{s}}$$

$$p_i(0)=\{\begin{array}{ll}1/{\sum }_{j\in {\alpha }_R},& i\in {\alpha }_R\\ 0.& \text{otherwise}.\end{array}$$

$${\phi }_{\alpha }\equiv \underset{(i\to j)\ni \alpha }{\min }\{J_{ij}\}.$$

$$\text{FPC}=\frac{1}{n}\sum _{l=1}^k\sum _{\alpha =1}^n{u}_{\alpha l}^2.$$

$$\Bigg{|}\frac{\text{MFPTs}^{\text{cg}}-\text{MFPTs}^{\text{full}}}{\text{MFPTs}^{\text{full}}}\Bigg{|}\leq\text{threshold},$$

Full text

I Introduction

The accurate modeling of large-scale complex systems such as proteins is a central challenge of computational (bio)chemistry. Molecular dynamics (MD) simulations have become an invaluable tool to gain microscopic insights into these pathways, and to calculate free energies and rate constant Frenkel and Smit (2002). However, fully atomistic unbiased simulations are still limited to hundreds of nanoseconds, while processes such as protein folding occur on much longer time scales. Markov state models (MSMs) are able to bridge time scales from nanoseconds up to milliseconds and beyond, and have been successfully applied to investigate folding pathways of large proteins Bowman and Pande (2010); Voelz et al. (2010), but also ligand-binding kinetics Plattner and Noé (2015). The dynamics of MSMs obeys a discrete-time master equation on a discrete space of long-lived (metastable) mesostates. MSMs are constructed from an ensemble of MD trajectories, but can also be augmented with experimental data Rudzinski et al. (2016). A key ingredient is the detailed balance condition linking forward and backward transition probabilities to energy changes, which drastically facilitates the construction of MSMs Prinz et al. (2011). Recently, some efforts have been undertaken to extend MSMs towards driven systems such as modelling sliding friction Pellegrini et al. (2016); Teruzzi et al. (2017), time-dependent external fields Wang and Schütte (2015), and polymer collapse in shear flow Knoch and Speck (2017).

In living cells, large proteins also perform specific functions such as assembling the ribosome and the generation of forces, e.g. the motor protein kinesin Kolomeisky and Fisher (2007). These molecular machines operate cyclically through converting chemical energy (typically through hydrolyzing adenosine triphosphate), which drives conformational changes. For example, in F1-ATPase these changes lead to the rotation of a central stalk through which the motor can exert mechanical forces and thus perform work on the environment Noji et al. (1997). The typical approach is to model such motors as operating in a non-equilibrium steady state (NESS) with stochastic dynamics that break detailed balance Liepelt and Lipowsky (2007); Zimmermann and Seifert (2012). Much progress has been made in formulating a coherent thermodynamic framework including the ever-present thermal fluctuations, which goes by the name stochastic thermodynamics Seifert (2012). In contrast, traditional designs of (macroscopic) engines are based on periodic variations of, e.g. the temperature, with the system returning to thermal equilibrium once the external stimulus is turned off. Stochastic thermodynamics has been extended to such periodic processes Brandner et al. (2015); Ray and Barato (2017). In the following, we will exploit that periodic driving can be represented as a NESS Raz et al. (2016); Rotskoff (2017), which greatly facilitates the theoretical analysis.

Here we develop a method to systematically construct periodically driven coarse-grained MSMs from atomistic data. To this end, we combine previous work on non-equilibrium Markov state modeling (NE-MSM) for time-dependent protocols Knoch et al. (2018) and non-equilibrium steady states Knoch and Speck (2015, 2017). Specifically, we construct a series of equilibrium fine-grained MSMs at different values of the external parameter. Exploiting the mapping to a NESS, coarse-graining is then performed in the space of cycles, thereby respecting the microscopic pathways along which entropy is produced. In Sec. II we provide some theoretical background on the mapping and stochastic thermodynamics Seifert (2012). In Sec. III, we describe the details of the implementation before discussing results in Sec. IV for a specific illustration of our method: alanine dipeptide in a time-dependent, periodic electric field. A possible application of this method is the efficient calculation of infrared spectra of biomolecules Krimm and Bandekar (1986).

II Theory

II.1 Mapping to NESS

We are concerned with classical systems that are well described by stochastic transitions between a (possibly very large) set of discrete microstates. These states represent small volumes in phase space. The time evolution is assumed to be governed by the master equation

[TABLE]

where $p_{i}(t)$ is the probability to occupy microstate $i$, and $W_{ij}(t)$ with $i\neq j$ is the transition rate to jump from state $i$ to state $j$. The diagonal entries $W_{ii}=-\sum_{j\neq i}W_{ij}$ ensure probability conservation. Eq. (1) is solved formally by the fundamental matrix $\mathbf{K}(t)$, which, in matrix notation, obeys

[TABLE]

with initial condition $\mathbf{K}(0)=\mathbf{1}$ and $(\mathbf{W})_{ij}=W_{ij}$. For a time-dependent rate matrix $\mathbf{W}(t)$, no general closed solution of Eq. (2) exists.

Here we are interested in rates with a periodic time dependence, $\mathbf{W}(t+T)=\mathbf{W}(t)$, where $T$ is the period. After a transient period, the actual probabilities acquire the same periodicity as the driving, $\mathbf{p}(t+T)=\mathbf{p}(t)$. We seek a stationary Markov processes described by a time-independent rate matrix $\tilde{\mathbf{W}}$ that mimics this (quasi)-stationary periodic process. The problem is typically formulated as follows: Given the time-averaged occupation probabilities and fluxes

[TABLE]

respectively, what are the rates $\tilde{W}_{ij}$ that yield the same probabilities $\bar{p}_{i}$ and currents $\bar{J}_{ij}\equiv\bar{\Phi}_{ij}-\bar{\Phi}_{ji}$? As realized by Zia and Schmittmann Zia and Schmittmann (2007), while the currents fix the antisymmetric part, there is some freedom in choosing the symmetric part of the fluxes. Two choices have been discussed recently, one that preserves the entropy production along edges Raz et al. (2016) and one that preserves the fluctuations of the currents Rotskoff (2017).

Here we follow a different route and make use of Floquet’s theorem Kuchment (2013), which states that the solution of Eq. (2) can be split into a periodic and a nonperiodic part leading to

[TABLE]

with integer $n$. Clearly, $\mathbf{P}=\mathbf{K}(T)$ is the transition matrix of a time-discrete Markov chain. We assume that the Markov chain can be reproduced by a time-continuous Markov process with time-independent rate matrix $\tilde{\mathbf{W}}$ so that $\mathbf{P}\equiv\exp(\tilde{\mathbf{W}}T)$. By construction, both $\mathbf{P}$ and $\tilde{\mathbf{W}}$ have the same eigenvectors and thus the same stationary probabilities $\mathbf{p^{\text{s}}}$. It is easy to check that

[TABLE]

independent of $t$ with $\mathbf{p}(t)\cdot\mathbf{P}=\mathbf{p}(t+T)$. Hence, we conclude that $\mathbf{p^{\text{s}}}=\bar{\mathbf{p}}$ due to the uniqueness of the stationary solution of $\mathbf{P}$. However, the fluxes

[TABLE]

and thus the currents $J_{ij}\equiv\Phi_{ij}-\Phi_{ji}$, are not necessarily the same as the time-averaged fluxes and currents of the periodic process. The computational advantage of our approach is that we do not have to measure currents in the periodic process, but it is sufficient to estimate the transition probabilities $P_{ij}$, a task that is well suited for Markov state modeling.

II.2 Determining the effective dynamics

The challenge is to determine the effective transition probabilities $\mathbf{P}$ in non-trivial systems with many microstates. One approach described by Wang and Schütte Wang and Schütte (2015) is to explicitly incorporate some time-dependent, periodic signal into molecular dynamics simulations and to build a non-equilibrium Markov chain from these simulations. This fine-grained model is described by the transition matrix $\mathbf{P}$ and is then further coarse-grained into metastable states identified from milestoning Schütte et al. (2011); Sarich and Schütte (2014). However, this approach is somewhat limited because changing the protocol, or even the period, requires to rerun the atomistic MD simulations.

Our strategy is different. From stationary atomistic MD simulations, and using standard tools, we construct a series of $N$ equilibrium Markov state models parametrized by the value of the external parameter $\lambda$ Knoch et al. (2018). The procedure is sketched in Fig. 1(a). For each of these MSMs indexed by $k=1\dots N$ with $\lambda_{k}$, we determine the rate matrix $\mathbf{W}^{(k)}$. The time-dependent matrix $\mathbf{W}(t)$ is then approximated step-wise by these matrices through discretizing the protocol $\lambda(t)\approx\lambda_{k}$ for $t_{k}<t<t_{k+1}$. This allows us to implement different protocols $\lambda(t)$ and, in particular, different periods $T$ reusing the same equilibrium MSMs. The transition probabilities

[TABLE]

are obtained through integrating Eq. (2).

Formally, the rate matrix $\tilde{\mathbf{W}}$ can be obtained by employing the matrix logarithm, $\tilde{\mathbf{W}}=\frac{1}{T}\ln(\mathbf{P})$. However, this expression is not guaranteed to return a correct rate matrix $\tilde{\mathbf{W}}$, i.e., every non-diagonal element is non-negative, which is known as the embedding problem. It is not yet known what conditions $\mathbf{P}$ must fulfill to be embeddable. When, however, a given matrix $\mathbf{P}$ is not embeddable, it is still useful to determine an auxiliary rate matrix approximating the time-continuous dynamics. One approach described in Ref. Israel et al. (2001) is to approximate the matrix logarithm by its series expansion

[TABLE]

The expansion can be computed recursively, while we stop either if any non-diagonal entry of the auxiliary matrix becomes negative, or the change in the next order approximation is small enough. Note that the linear approximation always returns a correct rate matrix. More powerful methods to calculate matrix logarithms are available and will be explored in the future.

II.3 Detailed balance and entropy production

The original dynamics obeys detailed balance at each point in time in the sense that for the vectors $\mathbf{p^{\text{eq}}}$ obeying $\mathbf{p^{\text{eq}}}\cdot\mathbf{W}(t)=0$ we have $p^{\text{eq}}_{i}W_{ij}(t)=p^{\text{eq}}_{j}W_{ji}(t)$ for any value of $t$ (which merely parametrizes the rates). Hence, freezing the protocol, the system would return to thermal equilibrium. For periodic driving, entropy is produced with rate $\dot{S}_{\text{tot}}$ due to the actual probabilities $\mathbf{p}(t)$ being different from $\mathbf{p^{\text{eq}}}$.

Mathematically, the network of states $\{i\}$ can be represented as an undirected graph, where the edges correspond to possible transitions. A single system jumps randomly from state to state so that the ensemble average obeys the master equation (1). In this graph we can identify cycles $\alpha\equiv(i\to j\to\dots\to i)$, which are defined as ordered sets of states, at the end of which the starting state is reached again and all other states are visited exactly once. Stationary Markov processes that break detailed balance, such as the effective dynamics obtained in Sec. II.2, imply non-zero currents $J_{ij}$ which have to flow along these cycles (which is a consequence of Kirchhoff’s current law $\sum_{j}J_{ij}=0$). These currents imply an entropy production with average rate Seifert (2012)

[TABLE]

where we have introduced the microscopic affinities $A_{ij}\equiv\ln(\Phi_{ij}/\Phi_{ji})$. In contrast, in the original time-dependent process probability can be temporarily accumulated within a period, see Ref. Rotskoff (2017) for a discussion. In the following, we ignore this extra contribution and focus on the entropy production $\dot{S}$ defined in Eq. (10). In particular, our coarse-graining scheme will preserve this entropy production of the effective dynamics.

II.4 Thermodynamically consistent coarse-graining

As already mentioned, in the effective stationary process all entropy has to be produced in cycles. Coarse-graining, i.e. the removal of states, most likely will also remove cycles and thus reduce the entropy production Puglisi et al. (2010). We have devised a coarse-graining scheme that circumvents this problem through performing the coarse-graining in the space of cycles Knoch and Speck (2015, 2017).

Our main tool is the cycle-flux decomposition Kalpazidou (2007); Altaner et al. (2012). We decompose fluxes as

[TABLE]

with cycle weights $\varphi_{\alpha}$, where the sum is over all cycles $\alpha$ that contain the oriented edge $i\to j$. Here, $\Phi^{0}_{ij}=\Phi^{0}_{ji}$ is the symmetric part of the fluxes. Details on the algorithm and the efficient determination of the cycle weights are provided in appendix A. Inserting Eq. (11) into Eq. (10), we obtain $\dot{S}=\sum_{\alpha}\varphi_{\alpha}A_{\alpha}$ with cycle affinities $A_{\alpha}$ (summing the edge affinities $A_{ij}$ along cycle $\alpha$).

The second step is to identify similar cycles and group them into communities, each of which is represented by a single cycle (details on how to do this are given in Sec. III.3). We keep only the set of these cycle representatives and delete all remaining cycles. The surviving states (being part of at least one cycle representative) then form the coarse-grained model. The new transition rates are obtained under two restrictions: (i) the mean entropy production rate $\dot{S}$ is preserved and (ii) all edge affinities $A_{ij}$ are preserved. The entropy production rate of the $l$th community follows as

[TABLE]

where $C_{l}$ contains all cycles belonging to the $l$th cycle community. The cycle affinity of the cycle representative (which is an element of $C_{l}$) is indicated as $A_{l}$. The only quantity that can be modified then is the cycle weight $\varphi_{l}\to\varphi^{\text{cg}}_{l}$ of the cycle representative, which is chosen as

[TABLE]

such that $\dot{S}_{l}$ is preserved. The symmetric contribution to the fluxes can be shown to be given by Knoch and Speck (2015)

[TABLE]

Finally, preserving affinities $A_{ij}$ implies that the logarithm of the ratio of conjugated nonzero rates $\ln W^{\text{cg}}_{ij}/W^{\text{cg}}_{ji}$ is conserved, which is fulfilled if $p^{\text{s}}_{i}/p^{\text{s}}_{j}=p^{\text{cg}}_{i}/p^{\text{cg}}_{j}$. We rescale all probabilities by the same factor so that $\sum_{i}p_{i}^{\text{cg}}=1$. Having obtained $\Phi_{ij}^{\text{cg}}$ and $p_{i}^{\text{cg}}$, the final elements of the coarse-grained rate matrix are determined through

[TABLE]

Further details can be found in Ref. Knoch and Speck (2015).

III Implementation

III.1 Molecular dynamics simulations

As a concrete example, we study alanine dipeptide in aqueous solution, which serves as a paradigm for computational studies of free energy calculations and rate estimations of protein dynamics Esque and Cecchini (2015); Chodera et al. (2007); de Oliveira et al. (2007); Velez-Vega et al. (2009). The molecular conformations are well characterized by its two dihedral angles $\phi$ and $\psi$ illustrated in Fig. 2(a).

All molecular dynamics simulations of alanine dipeptide have been performed employing the molecular dynamics software Gromacs 5.1.2 Abraham et al. (2015). The alanine dipeptide molecule is modeled with the CHARMM27 MacKerell et al. (2000) force field, dissolved in 362 TIP3P Jorgensen et al. (1983) water molecules and set up in a 2.25 nm $\times$ 2.25 nm$~{}\times~{}$2.25 nm box with periodic boundary conditions. Long range electrostatics were treated using particle mesh ewald summation with cubic interpolation and Fourier grid spacing of 0.16 nm. The cutoff for all short-ranged interactions was set to 1.0 nm. All hydrogen-involving covalent bonds were constraint by the LINCS algorithm Hess et al. (1997). For the time step we chose 2 fs. The temperature was set to 300 K using velocity-rescaling Bussi et al. (2007) thermostat with $\tau_{T}=1~{}$ps, while for the isotropic pressure coupling we used the Parrinello-Rahman Parrinello and Rahman (1981) barostat with $\tau_{p}=2~{}$ps.

III.2 Constructing the equilibrium MSM

MD simulations are performed at a constant electric field. We harvest 500 independent equilibrium trajectories ($\approx 10~{}$ns per trajectory) for each of the ten considered field strengths [see Fig. 2(b)] and record configurations every 1 ps. All simulations are started from the same initial configuration with the first 2 ns of each trajectory being discarded. From these trajectories, we construct a series of equilibrium MSMs based on both dihedral angles ($\psi,\phi$). Instead of directly discretizing the dihedral angles, we employ the cos/sin projection for both angles, returning a four-dimensional space, i.e.,

[TABLE]

The benefit of doubling the dimensionality is that, in contrast to the dihedral space, the cos/sin space offers a distance metric which allows the application of distance-based clustering algorithms and is required for our definition of cycle centers and diameters.

To further improve the configuration space representation, we employ the time-lagged independent components analysis (TICA) Pérez-Hernández et al. (2013) with lag time $\tau_{\text{tica}}=1$ ps. TICA is an orthogonal linear transformation similar to a principal component analysis (PCA). However, its orthogonal components point in the direction in which the slowest conformational changes occur for a given input lag time $\tau_{\text{tica}}$. Although the kinetic variance ($\approx 96\%$) of the TICA transformation suggests a reduced two-dimensional configuration space, we keep all four TICA dimensions as it allows us to back-transform from TICA space into dihedral space. Note that our analysis can be equally applied for a reduced configuration space. However, for illustrative purposes it is instructive to keep all dimensions.

The transformed configuration space is then clustered (using the $k$-means algorithm) into 250 states represented by centroids, which are positions in the underlying configuration space (here TICA space). We have chosen 250 states as a compromise between a too fine-grained rate matrix (with possibly badly estimated entries due to insufficient statistic for states that are almost never visited) and a too-coarse rate matrix. Because the output of the $k$-means algorithm depends on the starting configuration of centroids (zeroth iteration), we perform the clustering with 50 different randomly selected configurations and select the discretization exhibiting the smallest error. We verified that the best three discretizations do not show any qualitative difference. The final centroids are labeled as $R^{(i)}_{k}$ indicating the $i$th component in TICA space of the $k$th centroid, with $k=1\dots 250$ and $i=1\dots 4$.

For constructing the equilibrium MSMs, we identify along every recorded trajectory a sequence of states, i.e., every position in TICA space is associated with its nearest (Euclidean distance) centroid. To estimate the rate matrix $\mathbf{W}^{(k)}$, we employ the maximum likelihood estimator McGibbon and Pande (2015) provided in the software package MSMBuilder Beauchamp et al. (2011), for which use a lag time $\tau=2~{}$ps.

III.3 Cycle space

To identify similarities among cycles, we need to quantify their properties. Here we follow the approach of Ref. Knoch and Speck (2017) and focus on “geometric properties”, but other measures are possible. Before defining these cycle space coordinates, it is convenient to introduce the indicator function

[TABLE]

whether or not a given state is part of a specific cycle.

We define two geometric collective variables to characterize cycles: the cycle centers (or rather their $i$th component)

[TABLE]

where $|\alpha|=\sum_{k}\chi^{k}_{\alpha}$ denotes the length (number of states) of cycle $\alpha$ and $R^{(i)}_{k}$ the $i$th component of the $k$th centroid, and the cycle diameters

[TABLE]

Together with the cycle affinities $A_{\alpha}$, the cycle centers and cycle diameters of all cycles form the cycle space. Note that, although strictly speaking the cycle centers and diameters are four-dimensional, we only make use of the first two dimensions in accordance with the two largest TICA components. Based on these features, we identify similiar cycles and group them into communities, each of which is represented by a single cycle representative. The details of the algorithm are given in Appendix B, which yields the final coarse-grained model with rates $\mathbf{W}^{\text{cg}}$.

IV Results

IV.1 Static electric field

Applying a static electric field leads to a new equilibrium state. To show the influence of a static electric field, we compute the free energy landscape projected onto the dihedral angles for three different field strengths shown in Fig. 2(c). Without an applied electric field, we identify three different regions known as the extended conformation of the backbone ($C7_{\text{eq}}$), left-handed ($\alpha_{L}$) and right-handed ($\alpha_{R}$) $\alpha$-helical conformers. Note that other studies split the $C7_{\text{eq}}$ Wang and Schütte (2015) or $\alpha_{R}$ Esque and Cecchini (2015) regions, or both Chodera et al. (2007), into further subregions. The classification employed here is the coarsest that still allows to identify the dominant conformations.

Already when applying a static electric field, a significant shift in the depths of the free energy basins can be observed. For $E=0~{}$V/nm the $C7_{\text{eq}}$ and $\alpha_{L}$ conformations are equally dominant, while the $\alpha_{R}$ conformations are comparatively rare. However, for $E=1~{}$V/nm the $\alpha_{R}$ conformations become more dominant while $C7_{\text{eq}}$ and $\alpha_{L}$ are less populated.

IV.2 Oscillating electric field

The periodic protocol that drives the system out of equilibrium has the physical meaning of an oscillating electric field caused, for example, by a laser beam. For concreteness, we assume a polarized electric field $\mathbf{E}(t)=E(t)\mathbf{e}_{x}$ with periodic $E(t)=E_{0}[\cos(\omega t)+1]$, which is characterized by the amplitude $E_{0}$ and frequency $\omega=2\pi/T$ with $T$ the oscillation period 111We choose the electric field to be positive because the $(\phi,\psi$) configuration space would otherwise be degenerated with respect to positive/negative field directions, making it impossible to compare both approaches within the $(\phi,\psi)$ representation..

We start our modeling approach by discretizing the time-continuous electric field $E(t)$ into 20 steps. Both continuous and discrete electric fields are shown in Fig. 2(b). Because of the symmetry of the cosine function, however, the number of effective discretized field strengths are reduced by a factor of two. As described in Sec. III.2, we conduct independent MD simulations for these 10 different field strengths, each yielding an equilibrium MSM defined through the rate matrix $\mathbf{W}^{(k)}=\mathbf{W}(E_{k})$.

IV.3 Stochastic pumping

To study the effect of an oscillating field on the dynamics of alanine dipeptide, we solve Eq. (8) for different periods and extract the corresponding effective rate matrix $\tilde{\mathbf{W}}$ from which we compute the respective steady-state probability distributions. In Fig. 3 the cumulated stationary probability distribution for all three metastable sets ($C7_{\text{eq}}$, $\alpha_{L}$, $\alpha_{R}$) is shown for different driving periods $T$, whereas the solid/dashed lines illustrate the equilibrium distributions for no ($E=0~{}$V/nm) and the strongest ($E=1~{}$V/nm) statically applied electric field. The cumulated probabilities are defined by summing over the corresponding steady-state probabilities,

[TABLE]

with $X=C7_{\text{eq}},\alpha_{L},\alpha_{R}$. All three probabilities seem to saturate for longer periods. The probabilities $P_{\alpha_{L}}$ and $P_{\alpha_{R}}$ increase with longer periods $T$, while $P_{C7_{\text{eq}}}$ decreases. The absolute change, however, is moderate with about 10-15%.

Interestingly, the stationary distributions compared with the equilibrium distributions do not follow any general trend, e.g. $P_{\alpha_{L}}$ is located between both static distributions with $P_{\alpha_{L}}(E=0\,\text{V/nm})<P_{\alpha_{L}}<P_{\alpha_{L}}(E=1\,\text{V/nm})$. For $P_{C7_{\text{eq}}}$ we observe the opposite trend, i.e., $P_{C7_{\text{eq}}}(E=0\,\text{V/nm})>P_{C7_{\text{eq}}}>P_{C7_{\text{eq}}}(E=1\,\text{V/nm})$, and finally for $P_{\alpha_{R}}$ none is true as both values are exceeded. In particular, the latter is of importance as it demonstrates the concept of stochastic pumping Astumian (2011), where an oscillating protocol is used to increase the occupation of certain molecular conformations above their equilibrium values.

To compare the approach introduced in this study with the complementary approach recently introduced by Wang and Schütte Wang and Schütte (2015) (cf. Sec. II.2), we conduct 1000 MD simulations (10 ns each, removing the first 2 ns for equilibration) explicitly including the oscillating electric field with oscillation period $T=5~{}$ps. From the recorded trajectories, we directly estimate the matrix $\mathbf{K}(T)$ (lag time 5 ps) based on the same configuration space discretization as described in Sec. III.2, and convert it to $\tilde{\mathbf{W}}$ following Eq. (9). To compare the quality of both approaches, we start in the initial state

[TABLE]

and monitor the time evolution of the cumulative probability distribution of all three metastable sets, which are now time-dependent. Fig. 4 depicts the time evolution of all three probability distributions for both approaches (gray lines represent the approach introduced by ref. Wang and Schütte (2015) and colored lines represent the approach introduced in this study), illustrating that both approaches are in excellent agreement. Moreover, results are indistinguishable from the averaged MD trajectories.

IV.4 Cycle space

To infer information about the non-equilibrium dynamics from the cycles detected by the cycle-flux decomposition, we compute the cycle space introduced in Sec. III.3, which allows us to identify patterns in the distribution of cycles. To this end, we determine for every identified cycle its center ($c^{(1)}_{\alpha},c^{(2)}_{\alpha}$) and diameter ($d^{(1)}_{\alpha},d^{(2)}_{\alpha}$) in configuration space, and the cycle affinity $A_{\alpha}$.

In Fig. 5(a), we illustrate four selected projections of these five cycle space coordinates determined for $T=5~{}$ps, which clearly show that the distribution of cycles indeed exhibits an internal structure and is not random. While the left panels exhibit four distinguishable clusters, the right panels only offer two or three separated point clouds. Our next task is to group cycles (points in cycle space) together forming cycle communities. To this end we employ $c$-means clustering (cf. Sec. III.3). For the shown example we find $k=5$ communities to fit the data best (highlighted by different colors).

To understand their significance with respect to the change in molecular configurations, we highlight in dihedral space ($\phi,\psi$) all states that belong to cycles of a specific cycle community by the same color [see Fig. 5(b)]. Apparently, cycles and thus states belonging to the red and blue community coincide with the metastable set $C7_{\text{eq}}$ and $\alpha_{L}$, respectively. The purple and green community, on the other hand, connect each two metastable sets, with green cycles linking $C7_{\text{eq}}$ and $\alpha_{L}$ configurations and purple cycles linking $\alpha_{L}$ and $\alpha_{R}$ configurations. Since the red and blue cycle community only contain states belonging to a single metastable set, we refer to them as local cycle communities. On the contrary, cycles that connect two or more metastable sets, here green and purple, are referred to as global cycle communities.

IV.5 Coarse-grained model

At this stage, we understand that cycle communities emerge through the interplay between the underlying potential energy landscape and the oscillating electric field that drives the system out of equilibrium. To be able to effectively examine the non-equilibrium dynamics when changing the oscillation period, we coarse-grain the effective rate matrix $\tilde{\mathbf{W}}$, returning the coarse-grained rate matrix $\mathbf{W}^{\text{cg}}$ with only a handful of states that represent the dominant (cyclic) dynamics. The basic idea is that all cycles forming a specific community are lumped together into a single cycle, the community representative, and all remaining cycles are deleted. The coarse graining is completed when all nonvanishing transition rates are dynamically consistently rescaled while preserving dominant time scales of the system.

To examine how the dynamics is altered when changing the driving period $T$, we construct coarse-grained NE-MSMs for three different periods: $T=5~{}$ps, $50~{}$ps and $100~{}$ps. In Fig. 6(a), the final coarse-grained rate matrices are visualized where colored cycles are representatives of their respective community. For $T=5~{}$ps the coarse-grained model quantitatively confirms the cycle distribution shown in Fig. 5(b). Interestingly, when the oscillation period is changed, the communities and therefore their representatives change too. For example, for $T=50~{}$ps the former ($T=5$ ps) local blue and red community vanish and two new communities, also labeled in red and blue, appear. The blue community is of importance as it encloses all three metastable sets $\alpha_{R}\to\alpha_{L}\to C7_{\text{eq}}\to\alpha_{R}$. For $T=100~{}$ps no new communities emerge. However, the former ($T=5$ ps) purple and red community disappear, leaving two effective cyclic modes. The green cycle represents a clockwise rotation of the $\psi$ angle, whereas the blue cycle connects all three metastable sets in clockwise direction. Both cycles are illustrated in Fig. 6(b).

V Conclusions

To conclude, we have introduced a computational method to construct coarse-grained Markov state models with a few states for systems that are periodically driven. This is achieved through exploiting Floquet’s theorem and mapping the time-periodic dynamics onto a stationary effective dynamics with the same occupation probabilities. For such a non-equilibrium steady state, we have previously developed a method that is based on coarse-graining in cycle space through clustering, i.e., the identification of similar clusters as measured through suitable collective variables. This approach has the advantage that the entropy production of the effective dynamics is preserved in the coarse-grained model. Through basing the selection of cycle representatives on dynamical data (in particular mean-first passage times) also the dynamics of the coarse-grained system is preserved, a feature that is very hard to achieve in structural coarse-graining (i.e., the reduction of continuous degrees of freedom through combining atoms into beads) Müller-Plathe (2002).

As a specific system, we have studied alanine-dipeptide, a paradigm for computational studies of protein dynamics. In a time-dependent electric field, the dynamics averaged over one period breaks detailed balance and becomes cyclic, i.e., transitions have a preferred direction with non-zero probability currents. These cycles can be classified into local cycles residing in a metastable basin, and global cycles connecting these basins. For longer periods, the peptide effectively equilibrates and consequently the local cycles vanish. Our method is not restricted to biomolecules, but could be applied to the stochastic dynamics in other networks.

The major challenge in constructing Markov state models is to identify a suitable low-dimensional space of collective variables (like the cycle center and diameters employed in this work). Major efforts in this field are currently directed towards this dimensionality reduction. In particular machine learning approaches seem to be a promising route to reduce the complexity involved in determining suitable collective variables Mardt et al. (2018). However, such an approach for driven systems will again have to respect the underlying cycle structure, and it will be interesting to see how our results can be incorporated into learning strategies.

Acknowledgements.

We acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through the collaborative research center TRR 146 (project A7).

Appendix A Cycle-flux decomposition

To decompose the current matrix $\mathbf{J}$, we introduce an algorithm that can be split into two parts. First, we identify nontrivial cycles, i.e., cycles with more than two different states, and secondly we determine their cycle weights. Theoretically, the number of possible cycles in a graph grows exponentially with the number of vertices. However, if both algorithmic steps – searching for a nontrivial cycle and determine its cycle weight – are combined by running them alternately, the decomposition becomes computationally affordable even for a large number of states.

To detect a nontrivial cycle, we perform the following steps:

Find the oriented edge $i\to j$ of the largest element of $\mathbf{J}$: $\operatorname{arg\,max}\,(J_{ij})$. 2. 2.

Identify the shortest path (smallest number of transitions) from state $j$ leading back to state $i$ by only following the nonzero transitions of $\mathbf{J}$. This step can be efficiently realized by applying a breadth-first search Newman (2010). 3. 3.

Return the cycle, $\alpha=\{i\to j\to\,$found path$\}$.

To determine the corresponding cycle weight $\varphi_{\alpha}$, we consider all fluxes along cycle $\alpha$ and determine their smallest value,

[TABLE]

which becomes the cycle weight $\varphi_{\alpha}$.

Summing up both steps, the final cycle-flux decomposition algorithm reads

Find a nontrivial cycle in $\mathbf{J}$. 2. 2.

Compute its cycle weight $\varphi_{\alpha}$. 3. 3.

Update all $J_{ij}$ by subtracting $\varphi_{\alpha}$ along cycle $\alpha$, $J_{ij}\leftarrow J_{ij}-\varphi_{\alpha},~{}\forall(i\to j)\in\alpha$ and continue with step 2. 4. 4.

The algorithm stops when the residuum $||\mathbf{J}_{\text{max}}||$ has become smaller than a threshold.

An implementation of this algorithm is provided with the Supplementary Information of Ref. Knoch and Speck (2017).

Appendix B Cycle representatives

All input quantities (i.e., the cycle features affinity, diameter, etc.) are normalized by their variance to make them comparable. These features are then used as input for the fuzzy $c$-means clustering algorithm returning membership degrees $u_{\alpha l}$ that express the probability that observation $\alpha$ belongs to community $l$. To obtain an indicator of how good the clustering results are we compute the fuzzy partition coefficient (FPC) that is defined as the Frobenius norm of the membership matrix

[TABLE]

Here, $k$ is the number of chosen communities and $n$ the number of observations (i.e., number of cycles in our case). The best fitting number of cycle communities is chosen such that the FPC becomes maximal/large. Note that it is not necessarily the best choice to always pick $k$ corresponding to the maximal FPC. Sometimes the determined $k$-value is too low, causing a too coarse description of the cyclic structure. In such a case it is worth checking the second/third best choice.

Once we have determined the number of cycle communities, we require a suitable set of representatives, each representing one community. In order to preserve the large-scale dynamics, we propose the following stochastic algorithm: The general idea is that any set of cycle representatives is thought to be appropriate if the graph spanned by $\mathbf{W}^{\text{cg}}$ is ergodic and the mean first passage times (MFPT) between metastable configurations are preserved. Metastable configurations, if not known beforehand, can be determined for instance by considering all configurations / states found to belong to “local” cycle communities, i.e., communities with cycles which states do not overlap with states of other communities [see for example the red and blue colored states in Fig. 5(b)]. Another approach for NE-MSMs is to employ a kinetic clustering scheme using hitting times Sarich and Schütte (2014); Wang and Schütte (2015), which is in analogy to the often used PCCA+ algorithm Röblitz and Weber (2013) but also valid for transition matrices breaking detailed balance. For alanine dipeptide, however, we keep the manual selected metastable configurations as shown in Fig. 2(c).

Summing up both parts, the full coarse-graining algorithm is outlined as follows:

Pick one representative per cycle community by drawing a random number. 2. 2.

Check if the set of cycle representatives span an ergodic (connected) transition network. If yes, determine the coarse-grained transition rates (cf. Sec. II.4), else go back to step (1). 3. 3.

Compute MFPTs of the coarse-grained model and compare to MFPTs of the full model. If

[TABLE]

return $\mathbf{W}^{\text{cg}}$, else go back to step (1).