Good Rates From Bad Coordinates: The Exponential Average Time-dependent Rate Approach

Our ability to calculate rate constants of biochemical processes using molecular dynamics simulations is severely limited by the fact that the time scales for reactions, or changes in conformational state, scale exponentially with the relevant free-energy barrier heights. In this work, we improve upon a recently proposed rate estimator that allows us to predict transition times with molecular dynamics simulations biased to rapidly explore one or several collective variables (CVs). This approach relies on the idea that not all bias goes into promoting transitions, and along with the rate, it estimates a concomitant scale factor for the bias termed the “CV biasing efficiency” γ. First, we demonstrate mathematically that our new formulation allows us to derive the commonly used Infrequent Metadynamics (iMetaD) estimator when using a perfect CV, where γ = 1. After testing it on a model potential, we then study the unfolding behavior of a previously well characterized coarse-grained protein, which is sufficiently complex that we can choose many different CVs to bias, but which is sufficiently simple that we are able to compute the unbiased rate directly. For this system, we demonstrate that predictions from our new Exponential Average Time-Dependent Rate (EATR) estimator converge to the true rate constant more rapidly as a function of bias deposition time than does the previous iMetaD approach, even for bias deposition times that are short. We also show that the γ parameter can serve as a good metric for assessing the quality of the biasing coordinate. We demonstrate that these results hold when applying the methods to an atomistic protein folding example. Finally, we demonstrate that our approach works when combining multiple less-than-optimal bias coordinates, and adapt our method to the related “OPES flooding” approach. Overall, our time-dependent rate approach offers a powerful framework for predicting rate constants from biased simulations.


Introduction
A major challenge in biomolecular simulation is to be able to accurately assess the transition rates (inverse of the mean residence time in a state) of complex processes, including conformational transitions and the binding/unbinding of macromolecules and their ligands.Processes of interest often involve rare events, where the system spends a large amount of time in a metastable state and rarely transitions to another relevant one, so the transition path time is typically orders of magnitude shorter than the time spent in either state. 1 Because of this, extracting rates directly from unbiased simulation is out of reach for all but the simplest of systems.
Numerous methodologies have been developed to accelerate rare conformational transitions, with the primary purpose being to compute ensemble-averaged observables. 2,3A major subclass of such methods operate by adding an additional biasing potential to the sys-tem's Hamiltonian, usually in terms of a small set of collective variables (CVs) which are believed or determined to be good descriptors of states of interest, or the path between them. 2,3][6][7][8][9][10][11][12] All of these methods pay the price of distorting the system's dynamics to obtain a much more rapid estimate of the underlying freeenergy landscape as a function of the CVs.
Most methods that tackle the problem of computing rates of rare transitions seek to generate a set of unbiased trajectories, either through combining direct sampling of many short trajectories from different starting points 13,14 as in Markov State Modeling, through Monte Carlo in trajectory space as in Transition Path Sampling, 15 or by generating trajectories that progress in a particular coordinate as in Forward Flux Sampling, 16 Steered Transition Path Sam-pling, 17 Weighted Ensemble, 18 Transition Interface Sampling, 19 and Milestoning. 20][25] As such, we are interested in approaches that build on CV biasing methods, which have been used to probe conformational transitions with sufficient computational efficiency even for relatively complex biological assemblies.The challenge already mentioned is that these methods alter the dynamics, which prevents any obvious solutions to inferring the unbiased time scales of events.However, starting with the Hyperdynamics method of Voter, it was shown that the first passage time of rare events could be approximately predicted using biased simulations if bias is not applied during the actual crossing through the transition state, by formulating an ansatz for how the bias accelerates time. 26,27This approach was originally developed using a time-independent potential defined on the whole system of interest, but later in the Infrequent Metadynamics (iMetaD) approach the same ideas were extended to CV biasing.MetaD 8 works by updating an external bias with a Gaussian centered at the current position in CV space every ∆ time steps (see Sec. 5.2 for details).iMetaD solves the problem of not biasing the transition over the barrier by only rarely updating the bias potential, such that it is unlikely to add bias on a high barrier during a fast crossing.iMetaD also introduces an additional approximation since the system is experiencing a time-dependent bias rather than a static one.The difficulty of avoiding adding bias during barrier crossings can also be mitigated by MetaD variants that only add bias within a region or up to a certain energy level, 28,29 which is now particularly easy to implement in the OPES variant of MetaD. 30MetaD and similar approaches have now been used and benchmarked for many different problems, especially for protein-ligand unbinding problems, as detailed in Ref. 31.
To extract the transition rate, these methods assume that a 'good' CV is used, and validate the rate estimates using a Kolmogorov-Smirnov (KS) test between the empirical and theoretical survival distributions.Unfortunately, for large and complex transitions, this might not be the case because finding a good CV is challenging.Moreover, CV quality indicators, such as the committor, 32 are expensive or intractable to compute.The Kramers time-dependent rate (KTR) method 33 was re-cently developed to extract transition rates from biased simulations, such as those used for iMetaD, but with much less sensitivity to CV choice.It introduced a new parameter γ, called the CV biasing efficiency, to be extracted from the biased simulations which scales the bias.In that work, it was shown that γ had a lower value for a poor CV in a simple 2D double-well potential, and as such it was assumed to relate to the CV quality. 33However, this has not been systematically demonstrated, and KTR has not been benchmarked on a problem where many CVs could be tested.Moreover, a direct connection between the KTR and iMetaD estimators has not been established.
In this work, we introduce a more general framework for computing rates from time-dependent biasing protocols, which allows us to treat the iMetaD and KTR estimators on the same footing.We then use this framework to propose a revision to KTR termed the Exponential Average Time-dependent Rate (EATR) method that bridges the two approaches.The EATR approach is shown to give the correct Kramers' rate when γ = 1 for an idealized 1D potential.Then, we use a Gō-protein system as a model to show how the prediction of rates depends on the choice of bias coordinate, and compare EATR's results to the true intrinsic rate.Importantly, we find that γ correlates with the intuition of CV quality.We find that for the poor biasing coordinates, the original KTR and EATR results are comparable and they enable an accurate recovery of the unbiased rates.Surprisingly, this is often true even in the frequent-biasing regime.The paper is organized as follows.First, we present a general theory for rate calculations from time-dependent biased simulations.We relate it to iMetaD and KTR, and then formulate the EATR approach.Then, we show results for a 1D overdamped Langevin dynamics, and for the unfolding process of protein G.We also show that the method can be extended beyond one biasing coordinate, presenting accurate results on protein G unfolding when biasing two CVs simultaneously.We end with conclusions and future perspectives of the work.

Transition rate for rare events
The rare event problem constitutes the stochastic crossing of a single free-energy barrier, where typically the waiting time to cross the barrier is much longer than the transition time over it.For a high barrier, the survival function S (t), which is the probability of a transition not occurring before time t, is given by an exponential distribution characterized by a single transition rate constant k 0 , Note that this survival probability is related to the probability of a transition occurring at time t via and it is also related to the cumulative distribution function (CDF), the probability that a transition occurred by time t, For Brownian dynamics, Kramers' rate theory [34][35][36] can be used with several approximations to calculate the barrier crossing rate, k 0 , from the bottom of a well in a potential U(x) with a single high barrier, where D is the diffusion coefficient.However, for most systems of interest, the diffusion coefficient and underlying potential (or freeenergy landscape) are not known and, therefore, one cannot directly use Eq. 4 to estimate the rate.Instead, the transition rate can be calculated using the survival function and a set of simulations i = 1, . . ., N where M ≤ N have crossed the barrier and N − M have not.Let t i be the time the i-th simulation crossed the barrier, and t i = T i the total simulation time for simulations i = M + 1, . . ., N. For right censored transition times, the likelihood is given by which is the product of the probabilities of the transitions occurring at times {t i } M i=1 and the probabilities of not transitioning before the times {T i } N i=M+1 for N − M simulations. 37n estimate for the transition rate can be obtained by substituting Eq. 1 and Eq. 2 into Eq. 5 and maximizing the logarithm of the likelihood with respect to k 0 , Note that the summation in the denominator takes into account the simulations that did not transition.When all simulations have crossed the barrier, Eq. 6 reduces to the inverse of the average barrier-crossing time, k 0 = ⟨t i ⟩ −1 where ⟨•⟩ denotes the average over the simulations.
The transition rate can also be calculated by fitting the CDF (Eq.3).To do so for the same set of simulations i = 1, . . ., N, we construct an empirical CDF h(t i ) which is the number of simulations that have transitioned before t i over the total number of simulations.The theoretical CDF can then be fit to the empirical CDF with a least-squares method by optimizing k 0 .

Expression for the general timedependent rate
For time-dependent biased simulations, such as a MetaD simulation, the transition rate is no longer a constant, and so as in Ref. 33, we consider how this time-varying rate could be taken into account.As a general expression, we can write the time-dependent rate as the unbiased rate k 0 scaled by a function of time We consider f (t) to account for the effect of the timedependence of the bias on the rate.In the following, we will relate it to different rate-extraction theories, such as iMetaD.Assuming a quasi-adiabatic biasing process, the survival function for the general expression of the timedependent rate is given by which can be used in Eq. 3 for fitting to an empirical CDF.Substituting this survival probability for the timedependent rate into Eq. 5 results in a general likelihood given by which can be simplified by taking its logarithm (10) Similarly to the unbiased case, we can maximize this expression with respect to k 0 , to obtain k * 0 as an estimator for the true rate,

Relation to iMetaD
In the hyperdynamics method, 26,38 the rate from transition state theory is scaled by the acceleration factor α = e βV(x)

X
(note that here we use the subscript X to denote a configurational average and unlabeled brackets to represent an average over separate trajectories) where V(x) is a fixed bias function added to the system's Hamiltonian as a function of the system's full coordinates.This α arises by considering the average effect in many individual trajectories whose time is dilated by a factor e βV i (t) , where V i (t) is the bias experienced by the system during simulation i at time t.Hyperdynamics then corresponds to a rate scaling function of the form and replacing it in Eq. 8 results in the survival probability Using this expression and assuming all simulations transitioned (M = N), the likelihood maximization of the rate gives In iMetaD, the form of the bias is also changing in time along with the configuration of the system in a history-dependent manner.Therefore, in iMetaD an acceleration factor for each simulation is approximated as a time average over that simulation instead of calculating it as a configurational average, In the Supplementary Information, we show that if we use this expression in our time-dependent rate formulation, we recover the standard rescaled rate theory found for iMetaD, We note that this result is derived using the LM approach for the case where all simulations have transitioned.
In Ref. 39, it was shown that directly fitting the theoretical CDF (obtained from Eq. 13) is less sensitive to outliers in the tail of the distribution.The KS test can be used to assess whether the transition distribution is well-described by the theoretical CDF (see Sec. 5.5).Recently, it was also suggested that the short time information from the CDF can be fit to get a more robust estimate of the rate. 40We note that the results from likelihood maximization (LM) and CDF fitting need not coincide, as we will describe below.

Kramers time-dependent rate and the CV biasing efficiency
Most of the transition-rate methods for biased simulations, such as those described above, assume that bias is applied along a perfect CV, where all added bias accelerates the barrier crossing event.However, for large biomolecular systems, choosing a priori a perfect CV for accelerating transitions to another targeted state is almost impossible.In practice, the bias is applied along non-ideal CVs, which insert bias along useless directions that are not aligned with the transition path.
To overcome this issue, the KTR theory 33 introduces a parameter, γ ∈ [0, 1], to account for the efficiency of the biased CVs.In addition to the unbiased rate, γ will also be estimated from the simulation transition times, and it will inform about the quality of the CV with γ → 0 reflecting poor CVs and γ ∼ 1 good ones.In the KTR approach as previously implemented, the efficiency of CVs is accounted for by defining the scaling function as where ⟨max V i (t)⟩ is the maximum bias applied at any point up to time t in simulation i, averaged over all simulations (denoted V MB (t) in Ref. 33).This form of treating the biasing potential was inspired by rate-calculation methods developed for forcespectroscopy, 41,42 where the barrier is reduced due to the external force, and therefore, by using ⟨max V i (t)⟩ it was assumed that the bias only affects the barrier height.Inserting Eq. 16 into Eq.8 gives which can be used directly in a CDF fit.Substituting this expression into the log-likelihood from Eq. 10, and maximizing it with respect to k 0 results in a γdependent expression for the unbiased rate To obtain the maximum likelihood estimate for both γ and k 0 , Eq. 18 is substituted back into the loglikelihood function and it is maximized numerically with respect to γ.

Exponential average time-dependent rate (EATR)
While rates estimated by the KTR approach are accurate (as shown in Ref. 33 and in the following sections), we show below in Sec.3.1 that it has the undesirable property that it does not agree with the iMetaD estimator when γ = 1, whereas we expect the iMetaD estimator to be correct for an ideal coordinate with a very high barrier and slow deposition time.The reason for this discrepancy is the way in which the average effect of the bias is defined-for iMetaD e βV i (t) is averaged, whereas for KTR the maximum of the biasing potential is averaged.
To unify the two theories, we propose the following modification to the KTR method, which will have the desired property of producing the same rates as iMetaD in the case where γ = 1.To do so, we introduce the scaling function This gives the survival probability Substituting this expression in Eq. 10 results in a loglikelihood of the form In the case where all simulations transition (M = N), the optimal unbiased k 0 as a function of γ is given by, where we have benefited from the idempotence of averages to rewrite the average of an average as a single average.Observing that when γ = 1, the term within brackets in the denominator is equivalent to t i α i , this estimator is then identical to the standard iMetaD estimator in Eq 15.Similarly as with the KTR, we can replace this expression into the log-likelihood and numerical maximize it with respect to γ to obtain estimates for both the unbiased rate and the efficiency of CVs.
Importantly, we note that Eq. 20 also provides us the option to numerically fit the biased empirical CDF to find the best values of k 0 and γ.As initial guesses, we use the LM estimates for k 0 and γ, and optimize using the Levenberg-Marquardt algorithm implemented in the SciPy Python package 43,44 to fit the empirical CDF to the theoretical CDF obtained from Eq. 20.The same can be done for the KTR method using the theoretical CDF from Eq. 17.

Benchmarking on a 1D Potential
The rate methods were tested first on the onedimensional matched-harmonic potential illustrated in Fig. 1a given by where the subscript 0 corresponds to the well and the subscript 1 corresponds to the barrier, Θ(x) is the Heav-iside step function, and x i , which is needed to make the potential continuous.Full simulation details for this model are given in Sec.5.1.
To estimate the unbiased rate, many Langevin dynamics simulations were performed on the potential from Eq. 23, starting from the bottom of the well (Fig. 1a).The first barrier-crossing time for each simulation was recorded, and the empirical CDF was calculated as described in Sec. 2. The unbiased rate was extracted by fitting the distribution to the expected Poisson process.The 2-sample KS test was performed to assess whether this transition is accurately described by a Poisson process.The p-value for the KS statistic is 0.97, demonstrating that the transition times are likely Poisson-distributed data.This yielded a log-rate of −5.98 ± 0.04 where rates are in units of τ −1 , with τ as the time unit, and the error is the standard deviation obtained from bootstrap analysis 45 (see Sec. 5.4).The log-rate calculated with Kramers' rate theory using Eq. 4 is −6.02, which agrees with the empirical rate within error.
We then performed well-tempered metadynamics (WT-MetaD) simulations (see Sec. 5.2) for this system to predict the rates using iMetaD, KTR, and EATR using different bias-deposition times ∆ dt with dt the MD timestep (varying from 10 to 10000 τ, which corresponds to fractions of the mean-first passage time varying from 10 −6 to 10 −2 ).In Fig. 1, we compare the methods for both the LM and CDF-fit for the situation where γ = 1 is enforced, because (i) this allows us to only assess the quality of the different time-dependent rate metrics, and (ii) for a 1D potential, all bias should go into promoting the transition as there are no orthogonal degrees of freedom; results obtained from fitting both k 0 and γ are shown in Fig. S1.We find that the original KTR method does not give a rate consistent with the empirical unbiased rate when γ = 1.On the other hand, we find that both the iMetaD and EATR methods are consistent with the expected values for the rate.Moreover, in agreement with Ref. 39, we find that fitting the CDF provides more accurate rate estimates than LM at small ∆ for all three methods, with the discrepancy between the two fits negligible for large ∆.The rate estimates for each fitting procedure improve as ∆ increases, which is consistent with the principles of iMetaD.This is also consistent with the results of the KS test.The 2-sample KS test was performed on iMetaD and the 1-sample test was performed for KTR and EATR as explained in Sec.5.5.The KS tests failed for the CDF fits at ∆ dt = 10 1 τ and 10 2 τ and passed for the CDF fits at ∆ dt = 10 3 τ and 10 4 τ.These KS test failures are shown for EATR in Fig. 1c as open circles.Given that the rate estimates are more accurate for fitting the CDF (even when the KS test fails), we report results from the CDF-fitting procedure below.We now focus on a more complex system, the unfolding of the B1 domain of protein G using a Gō-like potential and MD simulations.This system has the advantage that it is possible to obtain an unbiased estimate of the unfolding rate, while having a rich unfolding landscape complexity, and many possible choices of CVs to characterize the transition. 46Below, we will evaluate the quality of several good and bad CVs for predicting rates.For this study, we first considered the fraction of native contacts Q and distance between the ends of the protein (R ee ) which were shown to be a good and bad coordinate, respectively, for characterizing the folding of this protein in Ref. 46.In addition to these two CVs, we will consider the radius of gyration (R g ), the root-mean-squared-deviation from the native state (RMSD), and a recently developed lineardiscriminant analysis coordinate maximally separating states as defined through a clustering analysis 47 (LD1, see Sec. 5.3.2,Fig. S2).

Protein G unfolding
We first performed a 120 µs-long unbiased simulation to study the system's behavior.For this unbiased trajectory, we computed the potential of mean force (PMF) along each CV by taking the negative log of the histogram of observed CV values (Eq.26).The PMF for Q is shown in Fig. 2, and along all CVs in the SI Fig. S3, revealing a range of potential profiles and apparent barrier heights.Although the PMFs of each CV exhibit two wells, we know that R ee is a poor CV for characterizing unfolding because the unfolded ensemble contains configurations with small values of R ee comparable to the folded state, resulting in an unusually shaped basin at small values of the CV.
In contrast, we expect Q to be a good CV for unfolding, 48 and so we used Q to define when our system has transitioned out of the folded state by using the average committor in the unbiased simulation. 32,49,50To do so, we computed for every frame whether the simulation next reached the metastable free-energy minimum on the left (unfolded) before reaching the global minimum on the right (folded), and computed the average by binning these values as a function of the corresponding value of Q in that frame.Based on this result, shown in Fig. 2, we defined the unfolded region to be where Q < 0.35 because the average committor to the left of this point is effectively 1.0.To estimate the unbiased rate, we ran 200 unbiased simulations of the protein with randomized initial velocities, and stopped these simulations when Q dropped below Q = 0.35 using the COMMITTOR function of PLUMED. 51The residence time for the folded state was recorded for each simulation.The unbiased rate was determined to be 1.4 ± 0.1µs −1 using the CDF-fitting procedure previously described.The KS test using a Poisson distribution passed with p = 0.65, demonstrating a good fit.
We then performed 100 biased simulations for each CV at various bias-deposition times to determine how sensitive each method was to biasing speed.The re-covered rates from the methods are shown in Fig. 3a for each of the biased CVs as a function of the biasdeposition time ∆ dt.Bias deposition times varied from 1 ps to 10 ns, corresponding to fractions of the mean first passage time ranging from ∼ 10 −6 to ∼ 10 −2 as in the case of the simple potential in the previous section.As expected, longer hill deposition times are observed to generally increase the accuracy of all rate calculations.However, for intermediate to fast-deposition times, KTR and EATR extract unbiased rates closer to the true rate than does iMetaD, especially for the three CVs shown on the left of Fig. 3a.We also performed a similar study using untempered MetaD and find that, similarly, all methods work well, with KTR and EATR slightly out-performing iMetaD in the fast biasing regime (SI Fig. S4).
The KTR and EATR methods also give a measure of the CV biasing efficiency γ, which is shown in Fig. 3b.We find that Q and RMSD typically give higher values of γ than the end-to-end distance (R ee ) and the radius of gyration (R g ).This coincides with the physical intuition of protein unfolding CVs, where the number of contacts and the similarity to the folded structure should be most relevant.This also agrees with Ref. 48, which found that Q is a good CV for this system.The CV obtained from linear discriminant analysis (LD1, Ref. 47) appears to have a large value of γ for slow biasing and a small value of γ for fast biasing.A similar trend appears for all CVs tested, but this is most prominent in LD1, and we are still investi- gating the reason γ for LD1 is so much more sensitive than the other CVs here, while still serving as a very good CV for distinguishing folded and unfolded states (as proposed in our previous study 47 ).We note that the discrepancy between iMetaD and the true rate is most pronounced when KTR or EATR predict lower values of γ.
Our intuitive expectation is that a bad CV would require significant amounts of extra bias to be deposited before the system can overcome the apparent barrier in the FES for that CV, necessitating a low value of γ to compensate in our rate calculation.To check this, we computed the histogram of the maximum bias across the different simulations at different deposition rates.CVs with high γ and slow deposition have maximum bias that do not exceed the apparent barrier, while fast biasing and poor CVs require substantial extra energy to be injected into the system to effect transitions.Here, we use max V i because it is a direct ingredient in the KTR method.We note that another interesting quantity to compute here would be the amount of non-equilibrium work performed by the MetaD bias, which has been recently exploited in another estimator of rates from time-dependent biased simulations. 52or all these CVs, the KTR and EATR methods performed comparably well and consistently performed better than the iMetaD method.Interestingly, the KS test seems not to be as sensitive for KTR and EATR as it was in iMetaD.Fitting the CDF for iMetaD results in failed KS tests even where the error in the rate is small, but KTR never failed the KS test for these CVs and EATR only failed for one condition.This may be an effect of introducing γ as an additional fitting parameter, so it is possible to get fits closer to the empirical CDF with worse rate estimates.In Fig. S5, we show that the introduction of γ allows us to make good fits to the CDF; indeed, many pairs of (k 0 , γ) can be used to fit these data; however, we note that the resulting predicted rates are still quite close to the most confident rate prediction, so this small amount of flexibility is not a problem here in practice.

2D MetaD
MD studies involving large biomolecules or molecular assemblies will typically have many slow degrees of freedom characterizing transitions between important states, and hence we expect the need to use multiple CVs to bias the system in order to promote transitions in a reasonably short amount of simulation time.We wanted to assess whether KTR and EATR still work for this case, despite the fact that the role of γ in characterizing CV quality is less direct.To do so, we performed MetaD simulations while simultaneously biasing the end-to-end distance R ee and the radius of gyration R g .For these two CVs, a combined PMF from the long unbiased trajectory is shown in Fig. 5a.The rates obtained from each rate method are shown in Fig. 5b, all methods recovered the rate equally well apart from EATR for the fastest bias-deposition time, where the KS test failed.The corresponding CDFs and EATR fits are shown in Fig. 5c, confirming that the method is able to predict the distribution of transition times.For intermediate ∆ dt ≥ 10 2 ps, the value for γ improved in this case over only biasing R ee and improved slightly over only biasing R g .This suggests that biasing multiple CVs simultaneously increases the efficiency of biasing without affecting the rate estimation.This might also be the reason for why iMetaD preforms well in this case, while it failed for the same bias-deposition times in the 1D biasing cases (see Fig. 3).

Conclusions
In this work, we have developed a general rate theory for time-dependent biased simulations that encompasses several of the existing methods by using the bias scaling factor f (t) (from Eq. 11) with different analytical formulations.In practice, LM and CDF-fitting are two different manners for estimating the unbiased rate from an ensemble of simulations launched from a single state.Although we demonstrate in this work that the previously proposed KTR approach works robustly for a realistic problem, it does not predict as accurate results as iMetaD for the case of an ideal 1D CV.Therefore, we proposed the EATR formulation, which exactly coincides with iMetaD in the ideal CV case but which gives good rate estimates for bad coordinates more robustly than does iMetaD, with the additional benefit of reporting on the efficiency of that parameter through the parameter γ.
We have validated the methods over the more complex landscape of G-protein unfolding where we could still have ground truth.We have found that both the KTR and EATR methods offer accurate rate measurements from biased simulations.The accuracy of the rates determined from these methods are surprisingly insensitive to biasing rate and CV quality, even for frequent-biasing regimes where the average maximum bias likely exceeds the true free energy barrier.We find that the γ computed from both KTR and EATR report CV efficiencies γ that correlate with our qualitative intuition of what is a good or bad biasing coordinates.Overall, we find that the CDF fit is better than LM to obtain the rate, however, many solution pairs of {γ, k 0 } could pass the KS test, so getting a good fit is not sufficient to guarantee that the optimal unbiased rate was obtained; however, in practice the predicted rates were all very close to the highest confidence prediction (as shown in SI Fig. S5).We find that for 2D biasing landscapes all methods perform reasonably well, which is an indication that biasing many directions might be helpful for barrier-crossing enhancement.In the future, we hope to test the methods on more complex systems using multiple biasing dimensions and, if necessary, extend the theory to multiple dimensions as was done for force-spectroscopy in Refs.53 and 54.Despite this success, there are several aspects of the EATR method that can still be improved through future work.For example, we could take into account the effect of the bias on the pre-exponential factor.Although the method works well in the case of our coarse-grained model of protein G for very fast biasing, we have not solved the general problem of how to compute rates in the over-biasing regime where γ times the bias could still be larger than the true barrier, as for example for LDA in Fig. 4 at fast deposition times, which could lead to the overestimated rate and small γ (Fig. 3).Addressing this issue will be crucial for systematically using the EATR method for large systems with many slow degrees of freedom.Going forward, we would like to determine whether it is possible to find a theoretical interpretation for γ in multiple dimensions, e.g.whether it can be derived considering projection operator approaches, and investigate whether it is rigorously connected to non-equilibrium estimators of the effect of a time-varying bias on the rates. 52Finally, we note that in this work we have concentrated on MetaD as a way of applying bias, but given our new generalized rate estimator, it would be interesting to compare whether there are other more efficient CV-based biasing protocols that give similarly accurate results.

Overdamped Langevin dynamics on a 1D potential
We ran overdamped Langevin dynamics over the potential given in Eq. 23 using x 0 = −3, x 1 = 3, and ∆U = 8 k B T. We used an integration time step of dt =0.01 τ where τ is the time unit.The friction parameter used was 0.02 τ, which corresponds to the fric-tion coefficient ζ = 50 τ −1 and the diffusion coefficient D = 0.02 τk B T. All simulations were started from x = −3.Because the diffusion coefficient and potential are known, we can derive the standard Kramer's expression in the Smoluchowski limit 36 from Eq. 24 and use that to determine the unbiased rate.For this specific system, the theoretical rate is 9.49 × 10 −7 τ −1 .These Langevin dynamics simulations were performed using the PESMD tool in PLUMED. 51Some of these simulations were biased using WT-MetaD with a starting hill height of 1 k B T, a σ of 0.5, and a biasfactor of 2.0.MetaD simulations were performed for four biasdeposition times (∆ dt): 10 1 τ, 10 2 τ, 10 3 τ, and 10 4 τ.

Well-tempered metadynamics
In WT-MetaD simulations, a history-dependent biasing potential V(ξ, t) is generated at a position ξ in CVspace.V(ξ, t) is formed as a sum of Gaussians with width σ (which can differ for each CV) and height h deposited every ∆ steps.For a one-dimensional CV, this can be written as, where t j = j∆ dt are the times where hills were deposited prior to time t, and N hills = ⌊t/(∆ dt)⌋.Here ∆T is a tempering factor which causes the heights to decrease proportionally to how much bias is already applied at that point, and is specified in PLUMED by setting a biasfactor of the form (T + ∆T )/T , where T is the simulation temperature.In the original untempered MetaD, the hills are of constant height, i.e. ∆T → ∞.
In iMetaD, the pace ∆ dt would be taken to be large, such that the frequency of deposition (∆ dt) −1 becomes small.For notation simplicity, we have omitted the explicit dependence on ξ from Eq. 25 in all equations in the Theory.

Gō-like model of protein G
A Gō-like coarse grained model of the B1 domain of protein G was prepared to assess the accuracy of the rate extraction methods, starting from PDB ID 1PGB.This system was selected because it was previously used as a paradigmatic example of a two state folder with known good and bad reaction coordinates. 46,55,56n a Gō-like model, each residue is modeled as a bead at the position of the α-carbon.The force field for this model treats pseudo-bonds and angles harmonically, and pseudo-dihedrals using a Fourier series.Non-covalent interactions, as in Ref. 57, depend on whether the residues are in contact in the native structure, which is determined by whether the side chains of two residues contain heavy atoms within 4.5 Å of each other.The force field parameters for the model in Refs.46,56 were provided by the authors.Our implementation of the potential in LAMMPS 58 and input files for all simulations are provided in the GitHub for this article (see Sec. 5.5).

Molecular dynamics simulations
The MD simulations of the Gō-like model were performed using the LAMMPS software package. 58The software was updated partway through the project and the version used for each set of simulations is shown in Table S1.All simulations used a time step of dt =10 fs, and the temperature was held constant at 312 K using the Nosé-Hoover chain thermostat 59 with a damping factor of 1 ps and a chain length of 3.All simulations started from the folded structure.
For the unbiased simulations, we ran 200 replicates and ended the simulations when the protein model unfolded, which was defined to be when the CV Q decreased past 0.35 as described above.The empirical CDF for the transition times to the unfolded state were fit to Eq. 1 as explained in the Theory section to obtain the observed unbiased rate for this system.

Collective Variables
A variety of CVs were analyzed for the Gō-like model, which were used for the biased simulations and during the rate analysis.The first of these is the fraction of native contacts (Q), which captures the degree to which the protein is folded.This CV is the fraction of the contacts present in the native structure which are still present, and was defined as in Ref. 56.The end-toend distance (R ee ) was also used, as it was previously determined to be a poor coordinate. 46The RMSD of the protein with respect to the native structure and the radius of gyration (R g ) were included to compare with the previously used CVs for this system.
To define the LD1 coordinate, first we performed a cluster scan on the unbiased trajectory using the shapeGMM clustering algorithm 60 with 50,000 frames for training, 3 training sets and 15 attempts each, for cluster sizes (K) = 2,. . .,6.The training curve with cross validation from the scan is shown in Supplementary Fig. S2.We used the positions of the beads as input features for shapeGMM.We did a 5 state shapeGMM fit on the entire trajectory (∼ 1.2M frames) with 15 at-tempts to identify the distinct clusters.We then performed an iterative global alignment of the trajectory to the global mean and covariance.Multi-state Linear Discriminant Analysis (LDA) was performed on the globally aligned trajectory with frames from all 5 clusters.Only the first coordinate (LD1) out of four resulting LD coordinates has been used in this study. 47n Supplementary Fig. S2b,c we show that this coordinate completely separates the folded and unfolded states, with the other states appearing as intermediates.

Biased simulations
The collective variables and biasing for protein G were handled using PLUMED. 51The version of PLUMED used for each set of simulations is shown in Table S1.As is the case for the 1D potential, WT-MetaD was used to bias the simulations.A set of untempered MetaD simulations were also performed, the results of which are provided in the Supplementary Information.The parameters used for the WT-and untempered MetaD simulations are given in Supplementary Tab.S2.The values of σ were chosen for WT-MetaD according to the standard deviation of the biased CV in the folded state, and for untempered MetaD σ was chosen to be less than that used in WT-MetaD.We performed simulations at eight different bias deposition times (∆ dt): 1 ps, 10 ps, 100 ps, 200 ps, 500 ps, 1 ns, 5 ns, and 10 ns. 100 simulations were performed for each ∆ dt.The simulations were halted when the protein was determined to have unfolded, or when either wall-clock time reached 48 hours or a total simulation time of 10 µs was reached.

Potential of mean force and committor analysis
A long simulation of protein G was performed and the potentials of mean force (PMFs) along various CVs were determined from the unbiased simulation data using where ξ is the CV along which the potential of mean force is computed and P(ξ) is the probability density of ξ obtained by computing a normalized histogram.Committor analysis 32,49,50 along Q was done on this long simulation by assigning either 0 or 1 to each frame of the trajectory depending on whether the system visits the folded or unfolded state next, then taking the average for all frames associated with each value of Q.In order to prevent incorrect assignments to either state, for the committor analysis the system was considered to be in the unfolded state when Q < 0.25 and to be in the unfolded state when Q > 0.85.From this, Q = 0.35 was decided to be the critical value for unfolding, as illustrated in Fig. 2.This was chosen to be less than the transition state to prevent counting cases where the system enters the transition region, but fails to unfold.

Bootstrap analysis
Errors were obtained from bootstrap analysis. 45For this analysis, a new set of transition times was constructed by choosing random simulations from the original set with replacement.Once the new set had the same size as the original set, the rate calculation was performed on the new set.This was repeated 100 times and the standard deviation of the log of the rate and γ across these new sets is reported.

Kolmogorov-Smirnov test
The KS test was performed to assess whether the transition distribution is accurately described by the expected theoretical distribution.For the case of unbiased or iMetaD, it is a Poissonian distribution.The 2-sample KS test was used for the unbiased and iMetaD analyses.This version of the test determines the maximum deviation of the observed CDF from two samples and gives a p-value which, when sufficiently low, allows us to conclude that the samples most likely did not come from the same underlying distribution.We consider the empirical and theoretical distributions to coincide if p > 0.05.The 1-sample KS test was used for the KTR and EATR analyses, as generating large random samples from their distributions took a significant amount of time.This version of the test determines the maximum deviation of the observed CDF for one sample from a theoretical CDF, and gave the same results as the 2-sample test in all the cases that were checked.us and offering ideas for future directions.We are grateful to R.B. Best for providing us with his noncovalent force field parameters for the protein G Gōlike model.NM, SS, and GMH were supported by the National Institutes of Health under award number R35GM138312.SS was also partially supported by a graduate fellowship from the Simons Center for Computational Physical Chemistry (SCCPC) at NYU (SF Grant No. 839534).This work was supported in part through the NYU IT High Performance Computing resources, services, and staff expertise, and simulations were partially executed on resources supported by the SCCPC at NYU. PC was supported by the Flatiron Institute, a division of the Simons Foundation.NM thanks the Center for Computational Mathematics and GMH the Center for Computational Biophysics at the Flatiron Institute for their hospitality while a portion of this research was carried out.S5 Rate estimates using untempered MetaD

Figure 1 :
Figure 1: (a) The potential energy profile of the 1D matched-harmonic potential from Eq. 23.(b) The unbiased rate from Kramers' rate theory and from unbiased simulations are shown as the dashed black line and the dotted blue line, respectively.We compare these with the predicted rates from iMetaD, and from the KTR and EATR methods by asserting γ = 1 as a function of the bias-deposition time (∆ dt).Maximum likelihood estimates are represented by open symbols while CDF-fit estimates are filled symbols.Error bars are from a bootstrap analysis as described in Sec.5.4.(c) The empirical CDFs of observed transition times over the barrier are shown with their EATR-CDF fits; different bias deposition times are indicated in the caption with curves showing fastest to slowest appearing from left to right.Fits that fail the KS test are represented with open symbols.The unbiased empirical CDF (black points) is shown together with the Poisson-process distribution fit (black line) and the predicted distribution using Kramers' analytical rate (Eq.24, cyan dashed line).

Figure 2 :
Figure 2: The potential of mean force along the fraction of native contacts Q for the G ō-like model of the B1 domain of protein G colored according to the average committor function.The value of Q where the average committor is 0.5 is marked with the dashed line, while the critical value for unfolding Q = 0.35 is marked with the solid line.

Figure 3 :
Figure 3: (a) The rates obtained from fitting the CDF for iMetaD (black), KTR (blue) and EATR (orange) for each CV at various deposition times (∆ dt).The horizontal dashed line represents the empirical rate obtained from unbiased simulations.Open shapes indicate where the KS test failed.(b) γ values obtained for the KTR and EATR methods.Horizontal dashed lines represent the bounds placed on γ.In both panels, the error bars are computed from a bootstrap analysis as described in Sec.5.4.

Figure 4 :
Figure 4: Histograms of the maximum bias deposited for the WT-MetaD simulations for the five CVs and different biasdepositions times.As the biasing efficiency γ increases, there is a decrease in the amount of bias needed for transition.

Figure 5 :
Figure 5: (a) The potential of mean force projected onto the 2D R g /R ee space.(b) The rates obtained from each method at various deposition times.The horizontal dashed line represents the true unbiased rate for this system.(c) The empirical and EATR-fit CDFs for each deposition time.The deposition times in the legend are given in ps.(d) The γ values obtained for KTR and EATR at various deposition times.In panels b, c, and d, open shapes indicate where the KS test failed.We performed some simulations with even slower biasing with the result that not every simulation transitioned; while the resulting rate estimate was unaffected, the estimated value of γ dropped significantly.The effects on the rate and γ are shown in SI Fig. S6.

Figure
Figure S2: (a) Log-likelihoods of data coming from a shapeGMM model trained on K clusters as described in Sec.5.3.2.The agreement between training and cross-validation prediction suggests that the data is not over-fit for any of these numbers of clusters.(b) 2D scatter plot of sampled conformations colored according to their cluster assignments from ShapeGMM using a K = 5 model.Based on the position of each cluster in (Q,RMSD) space, we assign the 0 to be the folded state and 1 the unfolded state (colors defined in c).(c) 1D PMF profile (solid black) along the LD1 coordinate.Histograms of LD1 corresponding to separate clusters, normalized individually and colored according to cluster id are shown in transparent colors.

Figure
Figure S4: (a) The rates obtained from the four CVs biased using untempered MetaD using the iMetaD, KTR, and EATR methods.These methods provide accurate results in this biasing scheme.(b) The values of γ obtained from the KTR and EATR methods for each of the CVs.S6 Multiple k 0 and γ pairs pass the KS test

Figure S5 :
Figure S5: Same data as Fig. 3. (a) Shaded regions represent the ranges of rates for each of the five CVs from the iMetaD, KTR, and EATR methods which pass the KS test for some value of γ.(b) Shaded regions represent the ranges of γ values from the KTR and EATR methods which pass the KS test for some value of the rate.(c) Several EATR CDFs which pass the KS test for the LD1 CV at ∆ dt = 10 ps.(d) Several EATR CDFs which pass the KS test for the RMSD CV at ∆ dt = 1 ns.