Inference of non-exponential kinetics through stochastic resetting

Transitions between metastable states are a key feature of many chemical and physical phenomena, such as chemical reactions and protein conformational changes. Molecular Dynamics (MD) simulations are an appealing computational tool for estimating the kinetic rates associated with such transitions. They track the dynamics of the system in microscopic detail, providing accurate predictions of both thermodynamic and kinetic properties. However, MD simulations require an integration time step on the order of $∼1fs$⁠, limiting the simulations to a timescale of $∼1μs$⁠. Processes that occur on longer timescales, such as protein folding and crystal nucleation, cannot be sampled efficiently in standard MD simulations.^1–5 

Different enhanced sampling methods were developed to overcome this well-known timescale problem. An incomplete representative list includes umbrella sampling,^6,7 milestoning,^8–10 replica-exchange MD,^11–13 Gaussian accelerated MD (GaMD),^14–16 and metadynamics.^1,17,18 These methods promote extensive sampling of the transitions between metastable states within the accessible timescales of MD simulations, allowing the calculation of both thermodynamic averages and dynamic properties. For kinetics inference, most schemes assume that the underlying process has an exponential first-passage time (FPT) distribution,^{1–4,13,16,19–22} relying on transition state theory (TST)^23,24 or Kramers’ theory.^4,25–27

Exponential FPT distributions arise when a high energy barrier separates two narrow metastable states. However, many processes in nature have non-exponential FPT distributions. For instance, protein dynamics are sometimes better described by a sum of exponents with different rates, or skewed exponents.^28–31 Power-law kinetics were also found experimentally for transitions between two stable conformers of an enzyme.³² Even when the FPT distribution has an exponential tail, it can behave differently at short to intermediate times, for instance, obeying a power-law.³³ Yet, kinetic inference schemes for processes with non-exponential FPT distributions are currently lacking. We propose such a scheme, designed for simulations enhanced by stochastic resetting (SR).

Resetting is the procedure of stopping random processes and restarting them subject to the sampling of independent and identically distributed initial conditions. It was shown to expedite different kinds of processes, including randomized computer algorithms,^34–36 queuing systems,^37,38 and various first-passage and search processes.^39–47 We recently demonstrated the power of resetting for enhanced sampling of MD simulations.^5,48 We showed that it can expedite rare events in molecular simulations when used as a standalone tool⁵ and in combination with metadynamics simulations, which also improved the kinetics inference.⁴⁸ 

The first two moments of the FPT distribution provide a sufficient condition for SR to be beneficial: if the ratio between the standard deviation and the mean, known as the coefficient of variation (COV), is greater than unity, introducing a small resetting rate is guaranteed to lower the mean FPT (MFPT).⁴⁹ The COV of the exponential distribution is exactly equal to one, so processes with a broader FPT distribution can be accelerated with SR. However, they cannot be properly treated by most kinetics inference schemes from accelerated simulations since these assume an exponential FPT distribution.^{2–4,20–22} Resetting, on the other hand, can potentially provide accurate predictions, as it minimally perturbs the natural dynamics of the system between restart events.

In this paper, we present a method to extract the unbiased MFPT of processes with non-exponential FPT distributions from simulations accelerated by SR. We first present the theory underlying our method and use analytical FPT distributions to discuss its advantages and disadvantages. Next, we employ it to study a three-state, two-dimensional model system and the dynamics of a nine-residue alanine peptide in explicit water. We show that our method can consistently predict the MFPT with high accuracy, with speedups of more than an order of magnitude over brute-force unbiased simulations.

Ideally, T* should be greater than t′ but not too large so that sampling by resetting would still lead to acceleration over standard MD simulations. We imagine the following use case: First, be conservative. Given a few trajectories, roughly estimate a T* that would cut the far tail of the FPT distribution, leading to acceleration while leaving enough to fit. Choose the functional form that best fits the data. In the second step, use Eq. (1) to predict whether an even shorter timer will lead to greater acceleration. In doing so, one must still leave enough of the tail to fit it. Similar to other kinetic inference methods,^3,4,22 there will be a trade-off between speedup and accuracy. In all the examples below, we find that it is possible to estimate the scaling of the distribution tails through this procedure while benefiting from significant speedups.

For this analytical example, we know the exact survival function and can use Eq. (1) to obtain the exact MFPT under resetting for any choice of timer T*. We can also estimate the MFPT with no resetting (=5.005) through Eqs. (3) and (4), using the exact values of ⟨τ|τ ≤ T*⟩ and S(T*), and evaluating $k=−d⁡log(S(t))dt|t=T*$ with the exact derivative. Figure 1(b) shows the resulting estimations as a function of T* (blue, right y-axis). It also presents the speedup (green, left y-axis), with the speedup defined as the ratio between the MFPT without and with resetting, respectively. The dashed black line marks the results with the specific timer T* shown in Fig. 1(a). We observe that the MFPT estimation is within 1% of the true value for speedups up to $∼40$⁠. This is expected since the slope of the survival at T* is within 1% of k₂ at T*. This means that, for a process with this FPT distribution, we could spend 40 times less computational time per first-passage event compared to tedious, unbiased simulations and obtain nearly the same accuracy.

However, in more realistic scenarios, the accuracy of the predicted MFPT would also be limited by the number of trajectories sampled and how well they capture the slope of the distribution tails. We investigate the sensitivity of our approach to sampling noise by numerical sampling from the distribution in Eq. (6) to estimate its survival function. To estimate the MFPT in this case, we sampled 1000 batches for every timer, each composed of 100 FPTs that were numerically sampled from the distribution. We constructed trajectories with resetting in the following way: We first sampled a time from the FPT distribution without resetting and compared it with the timer. If the timer was smaller, we tallied it up and sampled a new time from the FPT distribution without resetting. We repeated the process until the sample was smaller than the timer. In that case, we added the sample to the sum, and the sampling of the trajectory was completed. The total summed time was a sample of the FPT with resetting at that value of T*.

Figure 2 shows the MFPT estimation as a function of the selected timer T* (top panel). The boxes show the range between the first and third quartiles (interquartile range, IQR), and the whiskers show extreme values within 1.5 IQR below and above these quartiles. The circles and horizontal lines give the mean and median, respectively. The associated speedups are plotted with green squares in the bottom panel. We find that a timer of 0.115 gives a speedup of 40, as anticipated, but provides a poor estimate of the unbiased MFPT due to insufficient sampling in the proximity of t′. Longer timers provide better estimations, with accurate averages (⁠$<0.6%$ absolute deviation) for speedups of up to $∼6$⁠. Increasing the number of samples improves the estimations. Figure 1 of the supplementary material shows results equivalent to those of Fig. 2 but sampling 1000 first-passage events for each timer. We find that a timer of 0.2 already gives $<4%$ absolute deviation, with a speedup of $∼24$⁠.

So far, we assumed that all simulations are performed serially, each one initiated only when the former is done. This would be the case if only one processor was available. However, we can usually perform simulations in parallel on several processors. This becomes increasingly affordable with the improvement in computer power, as demonstrated, for instance, by the parallel replica dynamics (ParRep) method.^54–56 As explained in Ref. 55, the improvement in computer power is mainly reflected in ever-greater levels of parallelization and not in per-processor speed. While it does not directly solve the timescale problem, it greatly benefits parallelizable enhanced sampling methods.

Consider, for instance, running 100 trajectories on 100 processors simultaneously; for unbiased statistics, one has to wait for all trajectories to end. Thus, the walltime of the sampling is the FPT of the longest trajectory, which is always larger than the empirical MFPT. For the hyperexponential distribution above, using 100 trajectories, it is $∼45$⁠, i.e., almost an order of magnitude greater than the MFPT [see Eq. (S1)]. The larger the COV, i.e., the larger the fluctuations in the FPT distribution, the larger the walltime compared to the MFPT. Resetting lowers the COV,⁵⁷ reducing the ratio between the longest first-passage time and the MFPT. This, in turn, translates to greater speedups in parallelizable settings. We demonstrate this by plotting the walltime speedup with 100 processors for the hyperexponential distribution (orange triangles in Fig. 2). The walltime speedup is defined as the ratio between the sampling walltimes without and with resetting, respectively. We find that in this case, it is larger by up to 40% than the equivalent speedup with a single processor.

Figure 3 shows the estimations of the MFPT with 100–2000 FPT samples (top panel). The median of the MFPT estimations is close to the true MFPT even with 100 samples with resetting. The mean is overestimated due to a few outlier values not shown. With 500 and 2000 samples, the mean is only $∼12%$ and $∼2%$ off the true MFPT, respectively. The bottom panel shows the estimates of α from the same data used in the top panel, demonstrating that we can accurately obtain power-law tails from simulations with resetting.

We next apply our inference scheme for MD simulations of a particle diffusing over a two-dimensional three-state potential [Fig. 4(a) and Eq. (S2)]. This potential was introduced by Khan et al. to represent a simple kinetic model with two first-order reactions.²⁰ The particle is initially positioned at state A. One reaction path leads to state I, with a kinetic rate of k₁. It is a reversible reaction, with transitions from I to A governed by rate k₋₁ < k₁. The second path leads to the product state B with rate k₂ < k₋₁. As Khan et al. point out, the system should reach a quasi-equilibrium between states A and I at time $k1+k−1−1$⁠, after which transitions from states A and I to B follow an effective rate of $κ=k−1k1+k−1k2$⁠. We assessed the values of k₁, k₋₁, and k₂ using unbiased simulations in which only a single reaction path was available, with the other blocked by a strong repulsive potential. The survival function shown in Fig. 4(b) confirms that transitions to state B follow rate κ at times larger than some t′, while faster at t < t′.

Figure 4(c) presents estimations of the MFPT (top panel) and time t′ (center panel) as a function of T*, sampling 100 first-passage events per batch for each timer. The associated speedups are plotted in the bottom panel. The speedup is similar with 1 and 100 processors because the COV without resetting is $∼1.05$⁠, very close to unity. We find that even the shortest timer, providing speedups of $∼6$⁠, estimates the MFPT within an order of magnitude of the unbiased value. As expected, the accuracy improves with larger timers. In addition, unbiased simulations fail to provide any insight on the timescale t′: the IQR is very broad, and the mean is higher than the third quartile due to large outlier values. However, simulations with resetting correctly estimate its order of magnitude with almost any timer choice.

Many proteins are characterized by a rugged free-energy surface (FES), with multiple metastable states separated by low energy barriers.^31,58,59 An extreme case is downward folding, where there are no energy barriers along the folding path.^28–30 The assumption of Poisson statistics, which is usually acceptable in systems with high energy barriers, is often invalid in protein dynamics. Therefore, we expect protein dynamics to be a natural proving ground for our new method. To demonstrate it, we employ our method to study the dynamics of a nine-residue peptide of alanine [Fig. 5(a)], which is the shortest peptide known to form a stable α-helix.⁶⁰ It also forms a stable “loop,” similar to other alanine chains.³³ This system was chosen to benchmark our approach due to its multiple transitions between metastable states on a long timescale, but not too slow, allowing benchmarking against brute-force unbiased simulations.

Figure 5(b) shows the FES of the system along two collective variables (CVs), obtained from a continuous, 2 µs-long unbiased trajectory. One-dimensional FES along these CVs are provided in Fig. 2 of the supplementary material. The first CV is the end-to-end distance x, calculated using the center of mass of the two edge residues, which identifies the closed “loop” state (A, x < 0.52 nm). The second is the average of three distances between pairs of H-donor nitrogen and acceptor oxygen within the peptide. This average, denoted here as y, identifies the helix⁶⁰ (C, y < 0.37 nm). We also identify a broad basin of metastable open configurations (B, x > 1.25 nm and y > 0.8 nm) and a metastable intermediate state, where two of the three H-bonds are formed (D, 1.42 < x < 1.67 nm and 0.39 < y < 0.44 nm). Representative configurations of the states are presented in Fig. 5(a).

We first sampled 100 configurations of state A, obtained in time intervals of 0.5 ns from a trajectory restricted to state A using a strongly repelling potential (see Sec. V) along the end-to-end-based CV. We then ran independent simulations with resetting, uniformly sampling the initial conditions from those configurations. The first-passage was defined as settling into any other stable state: B, C, or D. For comparison, we also performed 1000 brute-force unbiased simulations to determine that this process has a MFPT of $∼102ns$ and a COV of $∼1.46$⁠. The tail of the FPT distribution is exponential; hence, we use Eq. (4) for the MFPT estimations.

Figure 5(c) shows the MFPT estimations from simulations with resetting with different timers (top). The associated speedups are plotted in the bottom panel. The shortest timer provides speedups of $∼11$ and $∼18$ using 1 and 100 processors, respectively. Using this timer, we estimate the MFPT as 7.9 ns, about an order of magnitude lower than the true value. Larger timers give more accurate results and still lead to accelerations. For instance, using a timer of 20 ns we estimate the MFPT as $∼92ns$ with speedups of $∼2.2$ and $∼3.3$ using 1 and 100 processors, respectively. To achieve even higher speedups, resetting could be combined with metadynamics⁴⁸ or other enhanced sampling methods. However, a kinetics inference scheme that does not assume an underlying exponential FPT distribution for metadynamics must be developed first and left for future work.

To conclude, we presented an inference scheme for non-exponential kinetics from MD simulations accelerated by resetting. Almost all kinetics inference methods for enhanced sampling simulations assume an underlying exponential FPT distribution, but in many cases of interest, this assumption does not hold. Resetting is an especially appealing tool for non-exponential systems since it is guaranteed to expedite processes whose FPT distribution is broader than exponential. Moreover, it does not affect the dynamics between restarts and samples the natural dynamics of the underlying process up to the resetting time. If the FPT distribution has a well-behaved asymptotic decay starting at times that are slightly shorter than the resetting time, we can estimate the tail of the distribution by its behavior at shorter timescales. This, in turn, allows estimating the unbiased MFPT of non-exponential processes using simulations accelerated by resetting.

Our approach becomes increasingly favorable as the number of available processors grows. With standard unbiased simulations, a few trajectories with long FPTs dominate the total walltime for all simulations to end, regardless of the number of available processors. By limiting the trajectories to short timescales using resetting, we use the available resources more efficiently. Compared to sets of unbiased trajectories, we obtain similar accuracy at much shorter simulation times. As acknowledged earlier by the developers of ParRep, parallelizability is an increasingly important quality: while the rapid improvement in computer power does not enable longer simulations per se, it makes trivially parallelizable methods much more appealing.^55,56

We applied our method to several analytical FPT distributions, a three-state model potential, and an alanine peptide in explicit water. We obtain speedups of more than an order of magnitude with small statistical errors in the predicted MFPT. In all systems, we rely on the fact that the FPT distribution has either an exponential or power-law tail. However, our approach is much more general and can be employed for other kinds of distributions.

Samples from analytical distributions were obtained using Python. The script is available on the associated GitHub repository. Simulations of the three-state model system were carried out in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS),⁶¹ following the motion of a single particle with a mass of 40 gr mol⁻¹. They were performed in the canonical (NVT) ensemble at a temperature of 300 K, using a Langevin thermostat⁶² with a friction coefficient of 0.01 fs⁻¹.

We used input files by Ayaz⁶³ for the alanine peptide. The simulations were carried out in GROMACS 2019.6,⁶⁴ using the Amber03 force field⁶⁵ for the peptide and the extended simple point-charge (SPC/E) model⁶⁶ for the solvent. They were performed in the NVT ensemble at 300 K using a stochastic velocity rescaling thermostat.⁶⁷ Additional simulation settings are as reported in Ref. 60.

The integration time step was 1 fs for both systems. The progress of the systems in the CV space was measured using PLUMED 2.7.1.^68–70 For the alanine peptide, we tested whether the system entered a new basin every 1 ps. A first-passage event was considered only if the system was observed in the new basin for at least 5 consecutive tests. We also used PLUMED to restrict the system to state A when obtaining initial configurations. We added a harmonic potential $0.5kx−0.522$ at x > 0.52 nm, with k = 100 k_BT/nm² and x being the end-to-end-based CV.

Finally, in Figs. 2–4, we report statistics over 1000 independent sample batches. In Figs. 5 and 6, we present statistics over 1000 bootstrapping sets. The total number of simulations used for bootstrapping and the total number of trajectories ending in first-passage for each timer are given in Tables 1 and 2 of the supplementary material.

The supplementary material contains the following: Results for the hyperexponential distribution with extensive sampling, discussion on the expected walltime with multiple processors, details of the three-state potential, one-dimensional FES of the alanine peptide, and metadata for Figs. 4 and 5.

B.H. acknowledges the support from the Israel Science Foundation (Grant Nos. 1037/22 and 1312/22) and the Pazy Foundation of the IAEC-UPBC (Grant No. 415-2023). This project has received the funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 947731).

The authors have no conflicts to disclose.

Ofir Blumer: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Visualization (lead); Writing – original draft (lead). Shlomi Reuveni: Conceptualization (equal); Funding acquisition (supporting); Supervision (supporting); Writing – review & editing (equal). Barak Hirshberg: Conceptualization (equal); Funding acquisition (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal).

Example input files, source data for all plots, and a Python class for the proposed inference method are available in the GitHub repository: github.com/OfirBlumer/nonExponentialKinetics.