Accelerated estimation of long-timescale kinetics by combining weighted ensemble simulation with Markov model "microstates" using non-Markovian theory

TL;DR

This paper introduces a method combining weighted ensemble simulation with non-Markovian theory and Markov State Models to efficiently estimate long-timescale kinetics, such as protein folding, without requiring steady-state convergence.

Contribution

The authors develop a post-simulation clustering approach using haMSM to extract unbiased kinetics from WE data before steady state is reached, improving efficiency.

Findings

Validated on a 2D diffusive process showing unbiased MFPT estimation.

Applied to protein folding simulations, achieving progress in estimating folding times.

Demonstrated that the method bypasses convergence issues of traditional WE simulations.

Abstract

The weighted ensemble (WE) simulation strategy provides unbiased sampling of non-equilibrium processes, such as molecular folding or binding, but the extraction of rate constants relies on characterizing steady state behavior. Unfortunately, WE simulations of sufficiently complex systems will not relax to steady state on observed simulation times. Here we show that a post-simulation clustering of molecular configurations into "microbins" using methods developed in the Markov State Model (MSM) community, can yield unbiased kinetics from WE data before steady-state convergence of the WE simulation itself. Because WE trajectories are directional and not equilibrium-distributed, the history-augmented MSM (haMSM) formulation can be used, which yields the mean first-passage time (MFPT) without bias for arbitrarily small lag times. Accurate kinetics can be obtained while bypassing the often prohibitive convergence requirements of the non-equilibrium weighted ensemble. We validate the method in a simple diffusive process on a 2D random energy landscape, and then analyze atomistic protein folding simulations using WE molecular dynamics. We report significant progress towards the unbiased estimation of protein folding times and pathways, though key challenges remain.

Figures

Equations

$$k_{AB}=\mathrm{Flux}(A\to B;\mathrm{SS})$$

$$T_{ij}^{\alpha }=P\{X_{t+\tau }^{\alpha }=j|X_t^{\alpha }=i\},$$

$$k_{AB}=\mathrm{Flux}(A\to B;\mathrm{SS})=\sum _{i\notin B,j\in B}{p}_i^{\alpha }{T}_{ij}^{\alpha }\text{mbox}(haMSM)$$

$$\begin{array}[]{lr}T_{ba}=1,&T_{b\{i\neq a\}}=0,\end{array}$$

$$J_i=\frac{1}{{\tau }_{WE}}\sum _{j>i,k\le i}⟨w_{jk}-w_{kj}⟩$$

$$\begin{array}[]{lr}\mathrm{m}\frac{\partial^{2}\vec{x}}{\partial t^{2}}=-\gamma\frac{\partial\vec{x}}{\partial t}-\vec{\nabla}U(\vec{x})+\vec{f}(t),&\langle f_{i}(t)\cdot f_{j}(t^{\prime})\rangle=2\gamma k_{B}T\delta(t-t^{\prime})\delta{ij}\end{array}$$

$$\begin{array}[]{lr}\frac{\partial p}{\partial t}=-\vec{\nabla}\cdot\vec{J}+\gamma_{A}(\vec{x})\int_{d\Omega_{B}}\vec{J}\cdot\hat{n}\,d\vec{x}\\ \vec{J}=\beta D\vec{f}(\vec{x})p-D\vec{\nabla}p\end{array}$$

Full text

Non-Markovian, Markov State Models, Adaptive Sampling, Non-Equilibrium Dynamics

Introduction

The weighted ensemble (WE) is a parallel path sampling strategy for rare events, proposed by Huber and Kim in 1996Huber and Kim (1996), echoing earlier computational “splitting” strategiesKahn and Theodore (1951). WE entails running a set of simulations, accompanied by replication (“splitting”) of promising trajectories and pruning (via “merge” events) of less important trajectories. As long as simple statistical rules regarding the splitting and merging of trajectories are followed, the evolution of the trajectory ensemble is unbiasedZhang, Jasnow, and Zuckerman (2010); Aristoff (2018). By performing this splitting and merging based on a set of bins spanning the transition of interest, sampling can be greatly enhanced in rare or rate-limiting regionsAristoff and Zuckerman (2018). Recent efforts have utilized WE simulations to determine the pathways and effective kinetic rates of increasingly challenging long-timescale biological processes such as cell and genetic switchesDonovan et al. (2013, 2016); Tse et al. (2018), ion transportAdelman and Grabe (2015), protein-peptide bindingZwier et al. (2016), protein-ligand unbindingDixon, Lotz, and Dickson (2018); Dickson (2018), protein-protein associationSaglam and Chong (2019), and protein foldingAdhikari et al. (2019).

The Hill relation Hill (2005); Suarez et al. (2014) permits unbiased estimation of the mean first-passage time (MFPT), which is a proxy for the inverse rate constant, from steady-state probability flow. However, all applications of WE simulation for unbiased MFPT estimation are limited by the capability to obtain convergence to steady state (SS). In this study, we therefore address a question of considerable potential impact: Can SS kinetic information (i.e., the rate constant or MFPT) be extracted from unbiased transient data obtained via WE simulation prior to reaching steady state?

The data examined below will show that this is indeed the case, and that rate-constant estimates are most reliable when extracted using a non-Markov analysis Suarez et al. (2014); Suárez et al. (2016); Suarez, Adelman, and Zuckerman (2016) of a very fine discretization of configuration space – typically much finer than was used to run the original WE simulations themselves. WE simulations typically use a discretization of configuration space into bins, each of which may contain several trajectories running in parallel Huber and Kim (1996); Zuckerman and Chong (2017) – thus limiting the number of bins which can be used to tens or hundreds in practical cases.

We are considering the rate constant defined by $k_{AB}=1/\mathrm{MFPT}(A\to B)$, from one macrostate (A) to another (B), where A and B could be two non-overlapping conformational states, folded and unfolded states of a protein, or unbound and bound states of a complex. It has been established that $k_{AB}$ can be obtained from WE simulations which have reached SS based on the Hill relationHill (2005); Suarez et al. (2014); Suárez et al. (2016),

[TABLE]

which uses the probability flux into state B – i.e., the probability arriving per unit time – based on trajectories initiated in A. Analytical and numerical investigation of the transient relaxation to SS in the Smoluchowski equation for the probability flux has shown the conditions in which to expect monotonic relaxation to SS, and generically supports the merit of sampling the one-way transition ensemble necessary to haMSM constructionCopperman et al. (2019). In simple 1D systems with a single energy barrier, the SS relaxation time can be arbitrarily faster than the MFPT. However, complex systems inevitably will require long relaxation times to reach SS, and during the transient relaxation period “direct” WE estimates obtained from the probability flux will typically underestimate the true SS fluxCopperman et al. (2019). A complementary approach utilizes the event duration distribution to boost the measured flux during the transient regimeDeGrave and Chong (2018). Our concern is to extract accurate estimates of the SS flux and the entire SS distribution characterizing the transition, based on transient WE data – i.e., before steady state has been reached.

In prior work, we showed that a history-dependent non-Markov analysis of WE bins could be used to estimate the SS rate based on stationary solution of an appropriate transition matrix Warmflash, Bhimalapuram, and Dinner (2007); Dickson, Warmflash, and Dinner (2009); Vanden-Eijnden and Venturoli (2009); Suarez et al. (2014). Here, we term the non-Markov formulation a “history-augmented Markov state model” (haMSM) to emphasize both its relation to, and difference from, standard MSMs. In a haMSM, one constructs a separate transition matrix for each direction (A-to-B or B-to-A) based on the subset of trajectories which were most recently in macrostate A for the A-to-B direction (or B, for B-to-A). The history is equivalent to the directionality – i.e., the labeling of which macrostate was visited more recently. Once the history has been used to select the trajectory subset, rate estimation proceeds much as it would in a standard MSM calculation. Importantly, however, haMSMs provide unbiased estimates of the rates (inverse $\mathrm{MFPT}$s), at arbitrary lag timesSuarez, Adelman, and Zuckerman (2016). In this study we perform WE sampling of the one-way A-to-B ensemble, and the haMSM is itself just an MSM of the one-way WE trajectory ensemble; thus there is no direct comparison of an haMSM to a standard MSM in this context. We also emphasize that we are not computing “implied timescales” based on eigenvalues but the MFPT corresponding to a well-defined physical process.

Here we show that using finer (smaller) bins in a haMSM gives better performance for rate estimation from WE data, and we employ the clustering/binning processes which have been extensively developed in the MSM community Chodera and Noé (2014); Scherer et al. (2015); Paul et al. (2019). The motivation for using smaller bins is that larger bins are likely to possess internal free energy barriers leading to slower internal relaxation, which in turn will bias the transition probabilities (transition matrix elements) of the haMSM and ultimately the macroscopic transition rate estimate. Put another way, the distributions within smaller bins are more likely to be similar to the SS distribution, and the intra-bin distributions determine the transition probabilities.

Procedurally as shown in Fig. 1, we use fine “microbins” generated by analyzing uni-directional (A-to-B) data with the pyEMMMA software Scherer et al. (2015). These microbins are used to generate a haMSM, whose stationary solution provides the desired rate constant $k_{AB}$ Suarez et al. (2014). The stationary solution can also be used to initialize new WE simulations in the haMSM-estimated SS for validation. Our primary focus is showing that finer/smaller bins yield more accurate rate estimates, especially when compared to the relatively large bins used to run WE simulations. We note that practical WE simulations are limited in the number of bins which can be used because computing cost scales linearly with the number of bins Zuckerman and Chong (2017).

Methods

WE simulation of the A $\to$ B ensemble

We wish to describe the kinetics of macrostate (“state”) transitions in a molecular system using the Hill relation (1). We denote the initial/source state A and the target/sink state B, which are arbitrary non-overlapping regions of phase space. For the $A\to B$ transition we only need the “$\alpha$” (A-to-B directed) subset of trajectories which were most recently in A; those most recently in B (B-to-A) are denoted $\beta$. We employ WE simulation to enable parallel sampling of the $\alpha$ trajectory ensemble: trajectories are initiated in A, and those which subsequently arrive at the absorbing sink at B, are restarted in A. The starting distribution in A is noted below for each system. WE provides an unbiased representation of the $\alpha$ reactive trajectory ensemble Bhatt, Zhang, and Zuckerman (2010); Aristoff (2018), and the time-resolved probability flux into the sink state.

We report the variability between individual WE simulations by calculating nominal 95% credibility regions (CR) utilizing a Bayesian bootstrapping procedureMostofian and Zuckerman (2019); additional details about estimating the CR from haMSM estimates can be found in Adhikari et al.Adhikari et al. (2019)

haMSM formulation

We wish to post-analyze molecular dynamics data harvested from WE simulation in a discretized configurational space using a transition matrix $T$. The matrix encodes (conditional) transition probabilities $T_{ij}=T_{i\to j}$ among bins or “microbins” $i$ and $j$. Markov state models (MSMs) have been used to stitch together many independent simulations to approximate long-timescale processes Singhal, Snow, and Pande (2004); Noé et al. (2007); Voelz et al. (2010a) but do not distinguish the $\alpha$ and $\beta$ trajectory subsets Suárez et al. (2016).

The haMSM is a transition matrix formulation containing history labels, namely $\alpha$ or $\beta$, so that transition matrix elements are calculated solely from the corresponding trajectory subsets Suarez et al. (2014). Compared to a standard MSM, a haMSM expands the transition matrix formulation for an $N$-state system from $N\times N$ elements into a $2N\times 2N$ labeled rate matrix. Here, we are concerned solely with the $\alpha$ (last in A) ensemble of the source/sink system, which limits our attention to the $N\times N$ transition (sub)matrix

[TABLE]

where $X^{\alpha}_{t}$ is a trajectory in the $\alpha$ subset and $\tau$ is the lag time of the transition matrix.

The discretized version of the Hill relation (1) then becomes Suarez et al. (2014)

[TABLE]

with $p^{\alpha}_{i}$ the SS probability of microbin $i$ based on SS solution of the transition matrix $T^{\alpha}$. This non-equilibrium $\alpha$ SS does not obey detailed balance because of the net flux into the B state originating in the A state. Remarkably, the haMSM formulation (3) yields the correct MFPT independent of the microbin set or lag-time used to construct the transition matrix Suarez et al. (2014); Suárez et al. (2016) so long as the matrix elements are calculated based on trajectories launched in every bin $i$ according to the SS distribution. The lag-time independence is a powerful distinction from the traditional MSM because it means that all transitions collected in the $\alpha$ trajectory ensemble can be used to train the haMSM. In the atomistic protein folding example we use a $10\mathrm{ps}$ lag time, which should be contrasted with the $\sim$10-100ns lag times needed for accurate MSMs of molecular systemsVoelz et al. (2010a); Chodera and Noé (2014); Plattner and Noé (2015). In practice, then, training the haMSM should require significantly less trajectory data than for a standard MSM.

Because the haMSM transition matrix is built from the conditional transition probabilities between states, the SS distribution need only be reached locally within the defined bins, not globally over all configuration space. As noted above, we expect faster relaxation for smaller bins – that is, for haMSMs with more microbins. This expectation is borne out by the data shown below.

haMSM construction

Construction of haMSMs proceeds in a highly similar manner to constructing a standard MSM, except via WE weights instead of simple transition counts. All weights from WE simulation are tracked Zuckerman and Chong (2017) and available for this analysis. Microbins which are not fully connected (e.g., having transitions in but not out) were removed from the analysis. To build the transition matrix, the transitioning weight $w_{ij}$ from microbin $i$ to $j$ after lag $\tau$ was averaged over WE iterations and runs to construct the haMSM transition matrix using Suarez et al. (2014) $T_{ij}=\frac{\langle w_{ij}\rangle}{\langle w_{i}\rangle}$ with $w_{i}=\sum_{j}w_{ij}$. Here all transitions are used – that is, transitions only observed a single time were included in the analysis. We used the same lag time $\tau$ for the transition matrix as the WE lag time, 10.0 ps in the protein folding systems and 50.0 ps in the 2D test system. Extracting the structural transition information at each step requires a large amount of file input/output, and can be slow depending on file access capability/concurrency, but is not computationally demanding.

In using the transition matrix to extract non-equilibrium SS kinetics via Eq. (3), we must apply suitable boundary conditions. The source/sink boundary conditions were enforced by employing the exact target state (sink, or macrostate B) and source state (A) definitions used in the WE simulation when constructing the haMSM, with the transition matrix row for the target state enforcing all probability transfers directly back to the source state. That is, with $b$ the index of the sink state, and $a$ the index of the source state,

[TABLE]

The SS flux and mean first-passage time were estimated from the SS of the haMSM transition matrix via Eq. (3).

Microbin construction via clustering

History augmented Markov State Models (haMSMs) are constructed directly from the WE simulation trajectories. We employ “MSM-style” microbins, which is a new development in this work. Specifically, we apply unsupervised clustering methods to extract a set of microbins (Voronoi centers) representing the configuration space visited by the WE trajectories. For each haMSM based on clustering, the latest-occurring 100,000 structures from the training window are used as input to the clustering. The training window is chosen in different ways to address different questions (see below), and the windows used are always indicated in the results section. A new clustering was performed for each training window considered, so that only the information available inside the training window is ever used for the analysis. K-means clustering with a minimum RMSD metric based on all protein atom Cartesian coordinates (rotationally and translationally minimized) was performed using the pyEMMA software package Scherer et al. (2015), requiring $\sim$24 hours of computation over 12 CPUs for each clustering calculation. Given a set number of desired microbins and an initial choice of bin centers, a k-means clustering is deterministic, but accordingly varies given a change in the desired number of microbins or initialization.

In the NTL9 low-friction protein folding system, we also explored the use of dimensionality-reduction methods, to extract a lower dimensional subspace spanning the slow conformational degrees of freedom. We first processed the atomic coordinates to extract the matrix of intra $\alpha$-carbon distances. We tested Principal Component Analysis (PCA), as well as the Variational Approach for Markov Processes (VAMP)Paul et al. (2019) dimensionality reduction method which leverages the time-lagged covariance matrix to extract a representation of the generator and the subspace of slow degrees of freedom, appropriate for non-equilibrium ensembles. Trajectories of every segment from the last iteration of the training window were traced back to the first step (unfolded structure). These trajectories were used as input to the VAMP clustering algorithm implemented in pyemmaScherer et al. (2015) with an input lag time of half the training window.

Initializing WE simulations in estimated SS

A direct approach to validating the stationary solution of the haMSM is to initialize new WE simulations in the haMSM estimated SS and see if SS properties are observed. Re-initializing WE simulations from the steady-state distribution estimated from prior WE runs requires extracting structures and weights representative of the estimated SS, presented in pseudocode below.

Flux profiles

A stringent test of SS is to monitor the flux across collective iso-surfaces which separate source (A) and sink (B) states. At SS, the flux through any surface completely separating A from B will be identical, and constant– the SS flux $J_{\mathrm{SS}}$. To calculate the flux profile from a WE simulation, we calculate values of the collective variable (specified below for each system) for WE trajectories as follows. The collective variable is split into a set of monotonically ordered bins, and we average the weight transitioning between these bins at each WE iteration. The flux passing bin $i$, $J_{i}$ is then

[TABLE]

where $\tau_{WE}$ is the WE lag time (i.e., interval between iterations) and bin ordering is set to define positive flux towards the sink (B) state.

Systems and Results

Our goal is to validate the use of the haMSM approach to extract SS kinetics from transient trajectory data in complex systems. To that end, we apply the haMSM analysis to WE simulations of diffusion in a 2D random energy landscape, and to atomistic protein folding, described in detail in Ref. 14. We are interested in the extent to which haMSMs with many (thousands of) microbins can enable the estimation of the SS flux at molecular times less than the SS relaxation time. We will consider the ratio of the SS relaxation time ($\tau_{\mathrm{\mathrm{SS}}}$) to the latest molecular time (WE simulation time) $t_{\mathrm{mol}}$ used in the haMSM training window where the haMSM can predict the SS flux, to be the computational acceleration. This is not a computational acceleration as compared to any other method but rather assesses our ability to exploit transient information.

Our validation has three stages. (i) When the haMSM is trained with a full set of data which approaches SS, then the predicted $\mathrm{MFPT}$ should be independent of the clustering. Even coarse clusterings produce exact resultsSuarez et al. (2014); Suarez, Adelman, and Zuckerman (2016). (ii) When the haMSM is trained solely on transient data, we seek haMSMs which reliably predict the $\mathrm{MFPT}$ found from the more complete training – i.e., the SS value. (iii) The stationary solution of the haMSMs are used to re-initialize new WE simulations in the haMSM estimated SS, and these WE trajectories are analyzed for SS convergence.

Diffusion in a 2D Random Energy Landscape

We first tested our approach by simulating a particle diffusing in the 2D random energy landscape shown in Fig. 2. Parameters were chosen to emulate an amino acid in water at 300K ($m=100$ daltons, $D=\frac{k_{B}T}{\gamma}=6.08610^{-10}\frac{\mathrm{m}^{2}}{\mathrm{s}}$). Two stable states ($6k_{B}T$ deep Gaussian wells of $1.0\mathrm{nm}$ width) were separated by 10$\mathrm{nm}$ with the addition of 40 randomly placed Gaussians (python code to generate the energy surface provided in the supplemental information). A confining radial potential was also placed at $9\mathrm{nm}$ from the domain center. Particle trajectories were evolved according to the Langevin equation,

[TABLE]

with $\vec{x}(t)$ the particle position at time $t$, $\mathrm{m}$ the mass, $\gamma$ the friction, $U(\vec{x})$ the potential energy surface, $T$ the temperature ($k_{B}$ Boltzmann’s constant). Additionally, “recycling” was performed in the weighted ensemble implementation: trajectories which are found in the defined target state (B) are “recycled” back to the source state (A). State definitions were defined as all points within $0.1\mathrm{nm}$ of the metastable states at A $(x_{A},y_{A})=(4,-3)\mathrm{nm}$ and B $(x_{B},y_{B})=(-3,4)\mathrm{nm}$, shown in Fig. 2.

Weighted ensemble requires several implementation choices. Bins were based upon the radial distance to the target state and placed every $0.2\mathrm{nm}$ up to $10.0\mathrm{nm}$ and a final bin edge at $12.0\mathrm{nm}$ encompassing all trajectories with distances greater than $12.0\mathrm{nm}$. See Fig. 2 for state definitions and weighted ensemble binning. The system was evolved according to Eq. (6) using the OpenMM molecular dynamics Langevin integratorEastman et al. (2017). The weighted ensemble sampling was implemented with a $50.0\mathrm{ps}$ lag time in WESTPAZwier et al. (2015) with 4 trajectories per WE bin (200 trajectories for each WE simulation once full occupancy is reached). WE simulations, and numerical solution of the probability evolution were initialized from the source state.

We also computed the solution numerically in the Smoluchowski picture. Though inertial Langevin simulations were performed, the parameters chosen emulate an amino acid in water evolving on a potential with nanometer scale features and the system is well in the overdamped regime; inertial effects can be safely neglected as is borne out by the agreement with simulation data. Evolution of the probability distribution is thus approximated by the Smoluchowski equation with the addition of the recycling from target (sink, B) and initial (source, A) states,

[TABLE]

with $p(\vec{x},t)$ the probability distribution in the domain at time $t$ with absorbing boundary at the sink $p(\vec{x}\in B)=0$, $\vec{J}$ the current. The source distribution $\gamma_{A}(\vec{x})$ is a Dirac delta function: trajectories are recycled to the center of the gaussian well at $[4.0,-3.0]\mathrm{nm}$. Lastly, $\vec{f}=-\vec{\nabla}U$, $D=\frac{k_{B}T}{\gamma}$, and $\Omega_{B}$ is the boundary of B with inward-facing normal $\hat{n}$. Eq. (7) is the standard Smoluchowski equation with the addition of the source/sink boundary conditions for recyclingCopperman et al. (2019), and was solved numerically using the fipyGuyer, Wheeler, and Warren (2009) package. Slight variation was observed between different choices of grid sizes (400x400 - 800x800) and timestep ($10.0-100.0\mathrm{ps}$) and this variation set the minimum and maximum range of the SS flux into the target state shown in Fig. 4.

Transient evolution between the WE trajectory ensemble and the determistic probability evolution given from numerical solution of Eq. (7) are consistent (see supplementary Fig. S1), which is not surprising since WE is an unbiased path sampling procedureZhang, Jasnow, and Zuckerman (2010). A validation set of 50 long WE simulations of length $100\mathrm{ns}$ were performed, showing relaxation at around $0.1\mu\mathrm{s}$ to a SS indicating an $\mathrm{MFPT}$ of $\sim 1.0\mu\mathrm{s}$ via Eq. (1). A second test set of 50 WE simulations of length 2.5ns were performed, and were used as input to build haMSM models and estimate the steady state distribution and flux into the target. Fig. 2C shows the transient probability distribution from the WE test set (no haMSM) at the end of the training window at $2.5\mathrm{ns}$, about 1/40 of the SS relaxation time and very far from the SS.

The capability to estimate the correct SS flux value from transient data before SS was dependent upon the size of the microbins, confirming our expectations. Fig. 3 shows that while the coarsest bins yielded haMSM SS flux estimation which were the same as the direct transient flux, the haMSM prediction increased monotonically up to the correct SS flux value, at about 1024 bins, and remained consistent for all finer microbins calculated (1024-16,384 microbins). In the limit of very many microbins, the haMSM microbins evidently become Markovian in the sense that the microbins are small compared to the size of the features on the 2D landscape. Since the very first trajectories in the WE simulation only reach the target state by $2.0\mathrm{ns}$, this represents about the earliest possible estimation of the SS flux possible using roughly 1/40 of the SS relaxation time. The aggregate simulation time in the training set was about $1.3\mu\mathrm{s}$ which is roughly the $\mathrm{MFPT}$ itself. We note that it is likely that fewer WE simulations with fewer trajectories per WE simulation could have been used, although we do not explore this limit here.

The haMSM-predicted SS distribution does capture the overall scale and many important features of the SS distribution, shown in Fig. 2D. It is apparent that only approximate estimation of the entire SS distribution is necessary for accurate estimation of the SS target flux, and hence the $\mathrm{MFPT}$. With more training data and a longer training window, the predicted haMSM distribution would converge upon the SS distribution (as would the direct WE distribution). With the limited training set at $t=2.5\mathrm{ns}$, the haMSM estimate is clearly not an exact reproduction of the SS distribution (from numerical solution of Eq. (7) shown in Fig. 2B).

WE simulations re-initialized in the haMSM-estimated SS demonstrated clear convergence to SS (Fig. 4). We used the haMSM SS distribution to initialize a set of fully independent 100 WE simulations in the estimated SS (see section “Initializing WE simulations in estimated SS” for details). The target flux remains steady, consistent with the numerical solution of the Smoluchowski equation, haMSM SS estimations, and the direct flux from WE simulations longer than the SS relaxation time, shown in Fig. 4.

For a more granular view to confirm SS, in Fig. 4B we plot the flux profile along the 1D reaction coordinate (distance to the target state). This flux profile should become flat at SSGardiner (2009); Risken and Frank (1996). While indeed the flux profile becomes much flatter after reweighting/restarting, wrong-way fluxes away from the target state are transiently observed after reweighting, indicative of the errors in the estimated SS distribution at large distances to the target. Within an additional $5\mathrm{ns}$ these wrong-way fluxes relax to a nearly flat flux profile.

When metastable intermediates exist along the transition pathway between macrostates of interest, we expect SS convergence to be slow and to approach the MFPT itself. and in this situation we expect MSM microbin clustering will have the most utility to accelerate the estimation of SS properties. The 2D diffusion of a particle on a random energy landscape explored in this section is a simple toy model, although it does capture some of the complexity expected for more challenging systems such as the atomistic protein folding systems. To investigate the dependence of the haMSM performance with energy landscape ruggedness, we also examined the A to B transition in the 2D random energy landscape at 150K, or T/2. As expected, the mean first-passage time and the steady-state relaxation time increase exponentially with an Arrhenius-like exponential dependence, and are  300x slower. Meanwhile, the haMSM with thousands of microstates can predict the SS flux with only 10x more WE simulation time than at 300K, so the acceleration is increased by a factor of 30, see Fig. S2 of the supplementary information. Generically then, we should expect a greater speedup when processes are slower, and more activated (energy-barrier dependent).

Atomistic Protein Folding

In protein folding the configurational space is very high dimensional, and dynamical motion spans timescales across many orders of magnitude. We show here that in these very challenging systems, haMSMs with thousands of microstates are able to accelerate estimation of SS and protein folding times. In the toy model system, it is possible to reach a point where microbins are truly Markovian in the sense that the microbins are so small there are no important intrabin features of the energy landscape. Here we construct protein-folding haMSMs with $10.0\mathrm{ps}$ lag times, and do not expect this Markovian property to be satisfied for any construction method or number of microbinsSuarez, Adelman, and Zuckerman (2016). MSM models of NTL9 folding utilizing 100,000 microbins indicated that $10\mathrm{ns}$ lag times were necessary before implied timescales leveled off, suggesting Markovian behaviorVoelz et al. (2010a). Variationally optimized MSM construction of protein folding trajectories utilized lag times of $50\mathrm{ns}$Husic et al. (2016); Husic and Pande (2017) and $100\mathrm{ns}$Scherer et al. (2019). Despite the lack of truly Markovian microbins, we do find that haMSMs with many (thousands of) microbins accelerates SS estimation compared to the brute-force relaxation time to SS. While a full theoretical discussion of the SS relaxation time is beyond the scope of this work, we note that it should be very sensitive to the existence of metastable intermediatesBolhuis et al. (2002).

We study the folding of two proteins, the N-terminal domain of the ribosomal protein L9 (NTL9) and the IgG binding domain of streptococcal protein G. We analyze both weighted ensemble trajectories from Ref. 14 and from additional new WE simulations, following the same protocol. Computational wall time in these systems is discussed in more detail in Ref. 14 but we note that a single WE simulation with $\sim 1000$ trajectories accumulates roughly 0.5 ns/day of molecular time utilizing a single GPU/CPU (NVIDIA V100 gpu, Intel Xeon processor). The fast-relaxing low-friction NTL9 folding system serves as a well-converged system in which to validate the haMSM capability to accelerate SS estimation and explore the use of dimensionality reduction methods (PCA and VAMP) in the haMSM microbinning process. Utilization of dimensionality reduction methods significantly reduces the computational cost of the clustering but has a more subtle impact on the haMSM estimated SS, and in this work we explore their use only in the NTL9 folding system; we analyze the protein G folding system using haMSMs with microbins constructed using the all-atom based clustering with a minimum RMSD distance metric. We further attempt to validate the haMSM-estimated steady-state of the protein G folding system by reweighting and restarting a set of WE simulations in the haMSM estimated SS.

NTL9

NTL9 is a fast-folding globular proteinVoelz et al. (2010a); Lindorff-Larsen et al. (2011); Schwantes and Pande (2013); Nguyen et al. (2014) without experimentally measurable folding intermediatesHorng, Moroz, and Raleigh (2003). The previously performed weighted ensemble simulations of the protein folding process, Adhikari et al. (2019) determined a protein folding time utilizing Eq. (1) consistent with the experimentally measured 1.2 - 1.4ms folding timesHorng, Moroz, and Raleigh (2003). In those implicit solvent atomistic simulations, flux profiles along the RMSD to the folded state indicated that the WE simulations effectively approached SS in both the low-friction ($\gamma=1/5\mathrm{ps}$) and high-friction ($\gamma=\gamma_{\mathrm{water}}=1/80\mathrm{ps}$) systems within nanoseconds. Adhikari et al. (2019) Applying the haMSM analysis to both NTL9 systems (with a training window utilizing the final portion of the WE simulation), the haMSM estimated SS was independent of the method and number of microbins, confirming that the SS relaxation time was indeed short enough to approach the SS values via the brute force WE simulation. Adhikari et al. (2019)

We now take advantage of this well validated system to test the capability of haMSMs to estimate SS from pre-SS training data. The previously reported WE simulationsAdhikari et al. (2019) of 10 independent low-friction NTL9 folding simulations with a 2D progress coordinate (2D-WE) were run for $t_{\mathrm{mol}}=12\mathrm{ns}$ and are here considered the validation set. Five newly run independent 2D-WE simulations are used as the test set for building haMSMs. We estimate the SS relaxation time to be $\sim 5\mathrm{ns}$ in this system, which is the time needed for the average flux of the validation set to reach the 95% confidence region (CR) for the SS flux. This is also the time needed for the upper bound of the test set direct flux to reach the lower bound of the validation set SS flux.

Fig. 5 shows that haMSMs can indeed accelerate the estimate of SS, dependent upon the microbin clustering and number of microbins in the haMSM. The haMSM estimated SS flux, from haMSM models utilizing WE simulation data from the test set only up to the training window indicated, increases with the number of microbins in the haMSM when the training window is significantly less than the SS relaxation time. As the SS relaxation time of $\tau_{\mathrm{SS}}\sim 5\mathrm{ns}$ is approached, the haMSM estimate becomes independent of the number of microbins, for all clustering methods studied here, shown in Fig. 5 for both PCA and VAMP dimensionality reduction. Results for all-atom clustering using a minimum RMSD distance metric are quantitatively similar to the results utilizing PCA dimensionality reduction: see Fig. 6.

In principle, with sufficient training data and sufficiently small microbins, it should be possible to extract SS estimates from arbitrarily short (small $t_{\mathrm{mol}}$) trajectory data. In practice, the amount of trajectory and the ability to cluster configurational space into effective microbins will limit one’s ability to leverage transient information. As shown in Fig. 6A, all clustering methods reach order-of-magnitude estimates of the SS flux within $1\mathrm{ns}$ of molecular time ($\sim$1/5 of the transient, using aggregate simulation time $3.7\mu\mathrm{s}$) and a flux estimate within the validation set CR within 2.5ns of molecular time ($\sim$1/2 of the transient, using aggregate simulation time $11.8\mu\mathrm{s}$).

Flux profiles from the training set and the validation set (Fig. 6B) indicate continued relaxation towards SS even when the target flux and MFPT estimation have stabilized. The flux profile does not become completely flat, indicating the presence of a robust regime where flux into the sink (folded) state approaches SS before global convergence to SS. Steady-state estimation is most sensitive to the clustering procedure at early times in the transient, and becomes independent of the clustering method and number of microbins as SS is approached in the trajectory ensemble: see Figs. 5 and 6.

The haMSM steady-state probability distribution systematically redistributes weight to configurations similar to the folded structure (low RMSD to folded state) from unfolded configurations, as shown in Fig. 7, while a comparison of RMSD values and the fraction of native contacts on the landscapes are shown in the supplementary information Fig. S3. The PCA dimensionally-reduced landscapes show a landscape geometry where RMSD increases mostly monotonically. Meanwhile the VAMP landscape, which attempts to separate the subspace of “slow” degrees of freedom, shows a landscape where some unfolded large-RMSD structures are geometrically near small RMSD folded structures, with the largest geometrical distances between protein configurations of intermediate RMSD. It is not clear here which kinetic landscape is more correct.

For flux (Fig. 6A), VAMP performed similarly to PCA when many independent components (ICs) were retained; when only the leading ICs were retained, accelerated order of magnitude estimation of SS rates were obtained but overshot the true SS rate, and approached the correct SS rate from above. This counter-intuitive non-monotonic behavior illustrates that different dimensionality reduction and clustering procedures can lead to systematic variation in the estimated steady-state when applied to transient training data, supporting testing multiple approaches. A related consideration is to understand how haMSM models fail under training data sparsity. We took representative haMSM models with 10,000 microbins at different training windows and systematically reduced the amount of training data (while maintaining the same final iteration of the training window). Generically, models are well-behaved and predict similar SS target flux until they fail through a loss of connectivity in the matrix, with some dependence on the dimensionality reduction method, see Fig. S4 in the supplementary information.

A concern in building any complex model from limited training data is the possibility of overfitting, and choosing optimal hyperparameters. The choice of dimensionality reduction method, and the optimal number of microbins in a haMSM model given a finite set of trajectory (training) data, requires attention in future work, and could perhaps be guided by cross-validation procedures inspired by those developed in the MSM communityKellogg, Lange, and Baker (2012a); McGibbon and Pande (2015); Scherer et al. (2019), but this will require further development for application in this context. In the NTL9 system, we have an independent steady-state WE validation set, and beyond the validation of the probability flux into the folded state, we present additional validation of the steady-state distribution Fig. S5 and the likelihood of the validation set trajectories in the haMSM models, see Fig. S6 of the supplementary information.

Protein G

We expect that in larger, slower-to-relax, and more computationally expensive systems the haMSM method will have the most utility. We know that protein G, sampled out to $t_{\mathrm{mol}}=15\mathrm{ns}$, was not yet in SS by observing the flux profile along the WE progress coordinate shown in Fig. 8B. Experimental stopped flow kinetic measurementsPark, O’Neil, and Roder (1997); Park, Shastry, and Roder (1999) and coarse-grained structure-based simulationsShimada and Shakhnovich (2002) suggest that protein G has long-lived metastable on-pathway folding intermediates. Here we extend upon the protein G WE simulations and haMSM analysis we reported in Ref. 14, finding that in our atomistic WE folding simulations we observe slow relaxation to steady-state indicative of the presence of long-lived intermediates on the folding pathway. Using the haMSM estimated steady-state distribution projected along the fraction of native contacts, we find evidence for multiple metastable intermediates, see supplementary information Fig. S7.

Applying haMSM analysis to the original WE simulations (2D progress coordinate, friction $\gamma=\frac{5}{ps}$, 15 simulations of $15\mathrm{ns}$, aggregate simulation time $225\mu\mathrm{s}$) haMSMs with $10^{4}$ microbins constructed from the last $5\mathrm{ns}$ of WE simulation predict a SS flux about 3 orders of magnitude higher than the direct flux, shown in Fig. 8. Constructing haMSMs from this trajectory ensemble, we observe that the predicted $\mathrm{MFPT}$ depends strongly upon the number of microbins, consistent with the use of training data in the transient regime seen above.

To validate haMSM estimates, WE simulations were re-initialized in the haMSM-estimated SS and seen to maintain a consistent probability flux into the folded state consistent with the experimentally derived folding rate: see Fig. 8A. We initialized a set of 15 independent WE simulations (1D progress coordinate 1D-WE) in the haMSM-estimated SS; see section “Initializing WE simulations in estimated SS” for details. These reweighted simulations (rw1) were run for an additional $15\mathrm{ns}$. Aggregate simulation time of the reweighted simulations was $76\mu\mathrm{s}$. The lower bound of the direct flux CR remained inside the initial haMSM prediction CR, while the upper bound increased by a factor of $\sim 10$. These flux values remained consistent for a second round of restarted WE simulations (rw2), as well as haMSM analysis of the reweighted data. Figure 8 shows that the haMSM estimated SS flux did not vary systematically with the number of microbins (10-10,000) and were within the direct flux CR, indicative of convergence to SS. The protein folding time estimated via the Hill relation Eq. (1) of $0.3-4.3\mathrm{ms}$ is consistent with the experimentally derived folding time, which we take to be somewhere between the reported value of $3.1\mathrm{ms}$ at pH 4.0Park, O’Neil, and Roder (1997) and $39\mathrm{ms}$ at pH 11.2Alexander, Orban, and Bryan (1992). Note that the lowered friction (1/16 of water) of the Langevin simulation should decrease the observed $\mathrm{MFPT}$ of the WE simulation. Also note that the experimental protein folding time stated in our prior work Adhikari et al. (2019) was erroneous.

Fig. 8B indicates that the convergence to SS improves with each reweighting, but even after an additional $15\mathrm{ns}$ of 1D-WE simulation in the haMSM estimated SS (rw1), and a second round of reweighting and $2\mathrm{ns}$ of 1D-WE simulation, the flux profile is not yet flat: true SS is not reached – see Discussion. We speculate that relaxation processes far from the folded state may be very slow but not strongly contributing to the folding process. However, a weaker property of SS is that the haMSM-estimated SS flux becomes independent of the binning, which is observed when analyzing the WE simulations (rw1) initialized in the haMSM estimated SS (see Fig. 8A).

The protein G folding system, where long-lived metastable intermediates make the relaxation to SS very slow, is a system where there is a large benefit of using a haMSM to accelerate SS estimation. The SS relaxation time in the in silico system is not known, but ultrarapid stopped flow kinetic measurements indicated the presence of an on-pathway metastable intermediate with a lifetime of $600-700\mu\mathrm{s}$Park, Shastry, and Roder (1999) experimentally, which we take to be an estimate for the SS relaxation time. This $\tau_{SS}>10\mu$s would make WE simulation (without haMSM analysis) prohibitively costly, because hundreds of trajectories are integrated in parallel.

Discussion

The preceding data shows that the use of short-lag-time haMSMs with thousands of “microbins” can greatly accelerate the estimation of steady-state target fluxes, which translate directly into MFPT estimates. Not only does the haMSM analysis approach remove lag-time uncertainty, but we expect short lag times to be critical to analyzing mechanistic details on the ns scale and below. However, even though MFPT values evidently can be obtained prior to full SS relaxation, key challenges remain to determine the absolute convergence to the unique SS distribution.

Our investigation of flux profiles suggest that the determination of robust metrics to measure SS convergence remains an open problem which limits the ability to estimate the systematic error of rate estimation using the Hill relation. A flat flux profile along a reaction coordinate which separates the macrostates (folded and unfolded) of interest is a requirement for SS convergenceGardiner (2009); Risken and Frank (1996); Copperman et al. (2019). However, this requires global convergence even in regions of configurational space which may not be important to the transition of interest, and our data indicates a robust temporal regime where accurate kinetic rates can be obtained before this global convergence is reached: see Figs. 6 and Fig. 8. Nevertheless, we believe it will be valuable to continue to examine the flux profile to understand whether steady state has truly been achieved.

We expect that iterative procedures which alternate between WE path sampling simulation and haMSM SS estimation steps (schematic shown in Fig. 1) will be the most efficient in determining converged haMSM transition rates. WE trajectory analysis using haMSMs at many levels of resolution (e.g. numbers of microbins) can be used as a practical measure of convergence, where rate estimates become insensitive to the details underlying the haMSM construction as SS is effectively approached, see Figs. 5 and 8. This property can be used as a less stringent indicator of effective SS convergence, similar to how the leveling off of implied timescales is an indicator of Markovian behavior in MSM model building. In contrast, haMSM estimation of SS during the initial transient regime is very sensitive to the number of microbins and the dimensionality reduction method performed preceding the clustering, as suggested by our data (Figs. 4, 5, 6, 8). The efficiency of convergence to SS, beyond the system-specific presence or lack of metastable intermediates, will depend upon how rapidly and effectively the configurational space can be mapped and sampled.

The procedure of launching new WE simulations in the haMSM estimated steady-state and probing for steady-state relaxation, as described and performed in the 2D random energy model and the protein G folding model, is a rigorous validation of haMSM models and steady-state probability flux prediction. However, since it is not always going to be possible to achieve a completely relaxed steady-state, and moreover since additional simulations are computationally expensive and probably not practical to guide haMSM hyperparameter selection, it will be important to develop and apply cross-validation proceduresKellogg, Lange, and Baker (2012b); McGibbon and Pande (2015); Scherer et al. (2019); Wu and Noé (2020) which rely on only the haMSM model and existing WE trajectories in the transient regime, without a steady-state validation set.

Although the force field does not appear to have biased our results to a significant degree, a few issues deserve mention. Here we have studied implicit solvent low-friction models of protein folding which have allowed for the proof of principle determination of folding times in qualitative agreement with experimental rates, while maintaining computational tractability. The validity of the rates and mechanisms observed in any biomolecular simulation are dependent on the quality of the system description (force-field, solvent, etc.). The protein folding times estimated in this work, rescaled by the ratio of the friction coefficient of water at room temperature to the friction used in the implicit solvent Langevin simulations (a factor of 16), are encouragingly within an order of magnitude of experimental folding times. There are many reasons why this comparison is qualitative only. For Brownian dynamics in simple systems, the friction simply rescales time, but the behavior is less straightforward in complex biomolecular systemsShen and Freed (2002); Anandakrishnan et al. (2015). It is certainly physically reasonable to expect that the low-friction systems should have a shorter $\mathrm{MFPT}$ compared to high friction, but it is beyond the scope of this work to validate a time rescaling for the implicit solvent simulations. Moreover, under any simulation protocol, the protein folding $\mathrm{MFPT}$ is sensitive to the exact definition of the folded/unfolded states.

While the original WE simulationsAdhikari et al. (2019) were initialized from a single unfolded state, the ensemble of unfolded protein structures comprising the source state evolves toward a SS distribution in the WE bin encompassing the initial unfolded structure. Though trajectories which reach the folded state are fed back into the single initial unfolded state, the WE procedure merges trajectories randomly by weight, and trajectories near the folded state (sink) invariably have low weight due to the sink boundary condition (see supplementary Fig. S8). Hence, restarted trajectories in the single unfolded state immediately get pruned – i.e., merged with trajectories which have been allowed to evolve within the WE bin of the unfolded state; correspondingly, the nominal (and artificial) unfolded starting configuration is never observed in the SS re-initialized simulations. Thus it is more accurate to consider the source state of the folding simulations as being defined by the sub-ensemble of structures which map to the WE bin from which the single unfolded state belonged. Supplementary Fig. S9 shows the difference between the initial unfolded structure and the ensemble of unfolded structures from the same WE bin for the SS re-initialized protein G simulations.

In this work, we have focused on kinetics and the $\mathrm{MFPT}$, but haMSMs can provide mechanistic insight in the same way that MSMs are often used to elucidate pathways and mechanisms in biomolecular processesNoé and Fischer (2008); Voelz et al. (2010b); Schwantes and Pande (2013); Plattner and Noé (2015). We expect that the short lag times and granular structural detail of the haMSMs developed here will allow more realistic determination of kinetic pathways and mechanismsSuárez and Zuckerman (2018) complementing the unbiased rate determination. We speculate that the identification of structural states/pathways which are low direct weight during the transient regime, but important for the SS folding process, are key to the capability of haMSMs with fine structural detail to estimate the SS distribution– see comparison of transient and SS distributions for the 2D system in Fig. 2, NTL9 folding in 7, and protein G folding in supplementary Fig. S8.

As a proof of principle, we focused on the application of the WE sampling and haMSM modeling method to two protein folding conformational transitions we had previously investigated. However, the framework applies to any dynamical process which can be framed as an A-to-B transition, notably protein-ligand binding/unbinding. We are actively developing well-validated protein+ligand WE simulations to develop the application of the haMSM modeling for the determination of protein-ligand binding/unbinding rates and mechanism.

Conclusions

By using fine bins or “microbins” generated from pyEMMA clustering, we show that WE trajectory data from the transient (pre-steady-state) regime can be used to construct unbiased haMSMs which reliably estimate steady state kinetics. The fine bins can sidestep potentially long relaxation times that will occur if there are internal barriers within large WE bins. If a haMSM is not used, standard or “direct” rate calculation from WE would otherwise require reaching the steady state which likely is impractical for many systems of interest. Therefore, the approach developed here could be of considerable practical importance.

We validated the fine-bin haMSM approach in a 2D toy model of a particle in a random energy landscape, obtaining accurate estimation of the $\mathrm{MFPT}$ using a training window only up to 1/40 of the SS relaxation time, and in atomistic WE simulation of NTL9 folding, using a training window 1/10 to 1/2 of the SS relaxation time. In the case of the more challenging protein G folding system, we initialized new WE simulations in the haMSM estimated steady-state which determined a folding $\mathrm{MFPT}$ consistent with experimentally measured values, in a tiny fraction of the experimentally measured lifetime of metastable intermediates. Nevertheless, key challenges remain to obtain global SS convergence as we have detailed.

Overall, the accelerated estimation of long-timescale kinetics using WE simulation combined with haMSM analysis, demonstrated here, supports the ongoing application to more accurate but computationally expensive explicit solvent models, and larger systems.

Acknowledgements

We are very appreciative of helpful discussions with David Aristoff, Gideon Simpson, Barmak Mostofian, Ernesto Suarez and Lillian Chong. Computational resources were provided by the University of Pittsburgh Center for Research Computing, and the Advanced Computing Center and Exacloud cluster at Oregon Health & Science University. This work was supported by NIH Grant R01GM115805.

Supporting information containing supplementary figures, algorithms, and cross-validation is available.

Supplementary Information: Accelerated estimation of long-timescale kinetics from weighted ensemble simulation via non-Markovian “microbin” analysis

2D random energy landscape

NTL9

In the manuscript, we performed cross-validation of the observable of interest, namely the probability flux into the folded state. As a second cross-validation of this result, we show the Kullback-Leibler (KL) divergence $D(p|q)=\sum p\log{p/q}$ between the haMSM estimated steady-state distribution and the validation set steady-state distribution, as well as between the direct WE transient distribution, and the validation set, for various training windows and haMSM models with varying microbins and clustering methods shown in Fig. S6 (left). The KL-divergence is a measure of similarity between distributions, or relative entropy.

The KL-divergence decreases monotonically with longer training windows, as it should. Deep in the transient regime at 0.5ns, all haMSM models monotonically improve their ability to capture the steady-state as the number of microbins increases. As the training set data becomes closer to steady-state itself, overfitting can be observed for models somewhere between 1000 and 10,000 microbins. However, the overfit contribution is small compared to the vast gains in model prediction when the training data itself is deep in the transient regime. At intermediate times at 1.0ns training window, pca clustering shows monotonic improvement with increased number of microbins. However, vamp clustering shows subtle signs of overfitting, namely the KL-divergence begins to increase for models exceeding 100 microbins. Most tellingly, the direct and haMSM estimated steady-state distributions are equally good models of the steady-state for these VAMP clustered haMSMs. These happen to be the VAMP haMSM results where the steady-state flux is over-estimated. Further exploration is required to show if this phenomenom is due to overfitting of statistical noise, or a consequence of the VAMP dimensionality reduction at this lag time.

As another measure of model quality, and in the spirit of trajectory likelihood estimationKellogg, Lange, and Baker (2012b), we examine the likelihood of the validation trajectories in the haMSM models, relative to the likelihood of the training set trajectories. Because the haMSM models are only required to be a model of single-step probability transfer, we take the set of 2-step trajectories and calculate the log-likelihood as $\sum_{ij}\log{p_{i}}+\log{T_{ij}}$ filtered upon $p_{i},T_{ij}\neq 0$ where $p_{i}$ is the trajectory weight distribution or haMSM estimated steady-state weight in the microbin, and $T_{ij}$ is the haMSM transition matrix. Note that in the decomposition of the trajectory likelihood, we can ignore the contribution of the likelihood of a conformation in a state because those quantities are identical in the validation and training sets. Any model which is truly overfit should demonstrate such overfitting by making the training trajectories more likely than the validation trajectories. When applied to NTL9 folding, shown in Fig. S6 (right), we see that deep in the transient regime, there is no sign of this type of extreme overfitting. Only when the training data becomes close to steady-state itself at does this type of overfitting become apparent. We note that even by this definition of an overfit model, the overfit model better reproduces the validation set flux.

These analysis demonstrate that deep in the transient regime, increasing the resolution (number of microbins) in the haMSM model results in quantifiable model improvement without overfitting, beyond the improvement in the estimated steady-state probability flux.

Protein G