Efficient Estimation of Rare-Event Kinetics

The efficient calculation of rare-event kinetics in complex dynamical systems, such as the rate and pathways of ligand dissociation from a protein, is a generally unsolved problem. Markov state models can systematically integrate ensembles of short simulations and thus effectively parallelize the computational effort, but the rare events of interest still need to be spontaneously sampled in the data. Enhanced sampling approaches, such as parallel tempering or umbrella sampling, can accelerate the computation of equilibrium expectations massively, but sacrifice the ability to compute dynamical expectations. In this work we establish a principle to combine knowledge of the equilibrium distribution with kinetics from fast “downhill” relaxation trajectories using reversible Markov models. This approach is general, as it does not invoke any specific dynamical model and can provide accurate estimates of the rare-event kinetics. Large gains in sampling efficiency can be achieved whenever one direction of the process occurs more rapidly than its reverse, making the approach especially attractive for downhill processes such as folding and binding in biomolecules. Our method is implemented in the PyEMMA software.

Popular Summary

The dynamics of proteins—the biomolecules central to the functioning of living organisms—are often characterized by infrequent, large-scale transitions between two or more overall geometric structures. Between these rare transition events, the protein only exhibits small fluctuations around a definite overall structure. Precise information about the probabilities of this stochastic switching behavior is very useful for understanding the role or the function of a specific protein. However, obtaining this information is unfortunately very difficult; the required length and time scales cannot be resolved simultaneously in an experiment. Here, we show that it is possible to combine standard simulations and a second class of simulations, from which the interesting switching probabilities cannot be directly extracted, to obtain reliable values with reasonable computational times.

Simulating the dynamics of a protein on a computer can resolve the system on atomistic length and time scales, but the rare-event nature of the interesting transitions requires impossibly long running times in order to obtain reliable information about the relevant quantities. Here, we conduct simulations of “downhill” relaxation paths. We focus on an alanine-dipeptide molecule, which is commonly used to study molecular dynamics, and we also simulate the motion of a Brownian particle. We are able to compute the kinetics of rare events using simulations that are between 10 and 1 million times faster than conventional simulations that wait for rare events to occur.

We anticipate that our approach will help with the computer-aided design of new drugs in which the rare unbinding of the drug molecule is invariably a central obstacle.

Figure 1

Mean and standard error of the largest implied time scale $t_2$ total simulation effort $N$ for metastable three-state system with barrier parameter $b=4$. (a) Convergence of the mean value, using either a single long trajectory starting in one of the metastable states or the equilibrium distribution together with an ensemble of short chains relaxing from the transition state. The latter approach allows us to obtain a reliable estimate before the average waiting time for a single rare event ${\tau }_{13}+{\tau }_{31}$ has elapsed. The comparison of the estimated standard error (b) indicates a 3 orders of magnitude speed-up when estimating the rare-event sensitive quantity $t_2$ using the equilibrium distribution in combination with short relaxation trajectories.

Figure 2

Mean and standard error of the largest relaxation time scale $t_2$ total simulation effort $N$ for metastable three-state system with barrier parameter $b=9$. (a) Convergence of the mean value using short trajectories relaxing from the transition state. A correct estimate can be obtained 6 orders of magnitude before a single rare event would have occurred on average. (b) Standard error of the estimate. The estimation using long trajectories is unfeasible.

Figure 3

Potential $V(x)$ and equilibrium distribution $\pi (x)$ for Brownian dynamics in double-well potential. The equilibrium distribution (shaded area) is scaled to fit the scale of the potential function. It can be seen that the equilibrium distribution is concentrated in the metastable regions around the two minima of the potential at $\pm \sigma$.

Figure 4

Mean and standard error of largest implied time scale $t_2$, given total simulation effort $N$, for Brownian dynamics in double-well potential. (a) Convergence of the mean value, using either a single long trajectory starting in one of the metastable states or the equilibrium distribution together with an ensemble of short chains relaxing from the transition state. The latter approach allows us to obtain a reliable estimate before the average waiting time for a single rare event ${\tau }_{AB}+{\tau }_{BA}$ has elapsed. A comparison of the standard error (b) indicates a more than 2 orders of magnitude speed-up when estimating the rare-event sensitive quantity $t_2$. By combining short trajectories with information about the equilibrium probabilities, reliable estimates of the slowest relaxation time scale can be obtained with a total amount of simulation data that is about 1 order of magnitude smaller than the expected waiting time for a forward and backward transition across the barrier.

Figure 5

Free-energy profile of alanine dipeptide as a function of the dihedral angles. Energies are given in $\text{kJ}/\mathrm{mol}$. The average thermal energy $k_{B}T$ at 300 K is $2.493\mathrm{kJ}/\mathrm{mol}$. One can identify five metastable sets on the dihedral angle torus, indicated here by black lines. There are three low-energy (high-probability) sets $C_5$, $P_{\mathrm{II}}$, and ${\alpha }_R$, with $\varphi <0$, and two high-energy (low-probability) sets ${\alpha }_R$ and $C_7^{ax}$, with $\varphi >0$.

Figure 6

Mean and standard error of mean first passage time (MFPT) ${\tau }_{AB}$, total simulation effort $N$, for alanine-dipeptide MSM on the $\varphi$, $\psi$ dihedral angles. The mean first-passage time ${\tau }_{AB}$ of the $C_5$ to ${\alpha }_L$ transition is used as an observable for a rare-event process. (a) Convergence of the mean value is shown for a small number of long chains starting in the $C_5$ region (blue), an ensemble of short chains starting in the ${\alpha }_L$ region combined with different amounts of umbrella sampling simulations (green, red, light blue). The correct value of the $C_5$ to ${\alpha }_L$ transition, ${\tau }_{AB}=43\mathrm{ns}$, can be obtained for a total simulation effort of $N=70\mathrm{ns}$ when short “down-hill” simulations are used in combination with umbrella sampling data. (b) The standard error shows almost 1 order of magnitude speed-up when estimating the kinetic characteristic of a rare event ${\tau }_{AB}$ using short trajectories in combination with umbrella sampling simulations compared to using long trajectories and no additional information about the equilibrium distribution.

Figure 7

Forward committor $q^{+}(x)$ for transition from $C_5$ to ${\alpha }_L$ region. (a) Nonreversible reference estimate for $N=10\mu \mathrm{s}$ of simulation data. Dark contour lines indicate the free-energy profile. (b) Difference between the reference estimate and a nonreversible estimate for $N=1\mu \mathrm{s}$ of simulation data. There is a large error in the transition region due to insufficient sampling in the short simulation. (c) Distance for an estimate using a combination of umbrella sampling and standard simulation data with $N=N_{\pi }+N_C=960\mathrm{ns}$. There is no significant error in the transition region; the small error close to the second saddle is probably due to insufficient sampling of this region by the reference simulation.

Figure 8

Free-energy profile for alanine dipeptide as a function of the $\varphi$ dihedral angle. One can identify three metastable sets: two low-energy (high-probability) sets with $\varphi <0$ and a single high-energy (low-probability) set with $\varphi >0$.

Figure 9

Mean and standard error of MFPT ${\tau }_{AB}$, total simulation effort $N$, for alanine-dipeptide MSM on the $\varphi$ dihedral angle alone. The mean first-passage time ${\tau }_{AB}$ of the transition from the low free-energy region, $A=\{\varphi |-162°<\varphi <-54°\}$, to the high free-energy region, $B=\{\varphi |36°<\varphi <72°\}$, is used as an observable for a rare-event process. (a) Convergence of the mean value is shown for a small number of long chains starting in the $A$ region (blue) and for an ensemble of short chains starting in the $B$ region combined with different amounts of umbrella sampling simulations (green, red, light blue). The correct value, ${\tau }_{AB}=43\mathrm{ns}$, for the $C_5$ to ${\alpha }_L$ transition can be obtained even if no information about the $\psi$ dihedral angle is used in the construction of the MSM. (b) The standard error shows almost 1 order of magnitude speed-up when estimating the kinetic characteristic of a rare event ${\tau }_{AB}$ using short trajectories in combination with umbrella sampling simulations compared to using long trajectories and no additional information about the equilibrium distribution.

Figure 10

Energy landscape for the different attachment modes, $m=0,1,\dots ,4$

Figure 11

Mean and standard error of mean first-passage time of vesicle dissociation on a coarse-grained non-Markovian state space. (a) Convergence of the mean value, (b) standard error. Estimates are obtained for an ensemble of association trajectories starting in the high-energy region and relaxing towards the low-energy region in combination with the coarse-grained equilibrium distribution. The dissociation time can be estimated orders of magnitude before a single dissociation event would have been observed.

Figure 12

(a) Implied time-scale test. Convergence of the largest relaxation time scale $t_2$ indicates a good Markov model fit; i.e., the slow eigenfunction of the associated dynamical operator is well approximated. (b) The Chapman-Kolmogorov test validates the Markov assumption by comparing the evolution of self-transition probabilities predicted by the MSM parametrized at lag time $\tau$ with direct estimates from the data at larger lag times $n\tau$.

Figure 13

(a) Implied time-scale test. Convergence of the largest relaxation time scale $t_2$ indicates a good Markov model fit; i.e., the slow eigenfunction of the associated dynamical operator is well approximated. (b) The Chapman-Kolmogorov test validates the Markov assumption by comparing the evolution of self-transition probabilities predicted by the MSM parametrized at lag time $\tau$ with direct estimates from the data at larger lag times $n\tau$.

Figure 14

(a) Implied time-scale test. Convergence of the largest relaxation time scale $t_2$ indicates a good Markov model fit; i.e., the slow eigenfunction of the associated dynamical operator is well approximated. (b) The Chapman-Kolmogorov test validates the Markov assumption by comparing the evolution of self-transition probabilities predicted by the MSM parametrized at lag time $\tau$ with direct estimates from the data at larger lag times $n\tau$. Values are obtained from an ensemble of short trajectories starting in the high-energy region utilizing the equilibrium distribution in the estimation of the maximum likelihood estimator (MLE) transition matrix; cf. Eq. (9).