The milestoning algorithm of Elber and co-workers creates a framework for computing the time scale of processes that are too long and too complex to be studied using simply brute force simulations. The fundamental objects involved in the milestoning algorithm are the first passage time distributions KAB(τ) between adjacent conformational milestones A and B. The method proposed herein aims to further enhance milestoning (or other interface based sampling methods) by employing an artificially applied force, akin to a wind that blows the trajectories from their initial to their final states, and by subsequently applying corrective weights to the trajectories to yield the true first passage time distributions KAB(τ) in a fraction of the computation time required for unassisted calculations. The re-weighting method is rooted in the formalism of stochastic path integrals. The theoretical basis for the technique and numerical examples are presented.

The task of calculating kinetic properties of complex systems from molecular dynamics simulations is challenging but also of considerable interest.1 In contrast to computational methods designed for thermodynamic equilibrium, in which the basic observables are averages over the conformational space with an invariant measure that does not need to obey the system’s exact dynamics, studies of kinetics require physically correct time-ordered trajectories to obtain time-correlation functions as the basic objects.2 Since each time-correlation function describes the relaxation under investigation as an average over all relevant trajectories, adequate sampling for accurate calculation of long-time kinetics can quickly become computationally intractable via direct simulation.3 This is because a direct and brute force method of this type would require simulation times that are long enough for the system to transition back and forth between all states of interest as many times as needed for the transition statistics to converge. In other words, the trajectory needs to be long enough to obtain a distribution of first passage times. Several computational methods have been developed to address the challenge of calculating kinetics, from the venerable transition state theory, Kramers rate theory, or Anderson’s reaction flux method,4–6 and continuing with more recent developments, such as transition path sampling (TPS),7 transition path theory (TPT),8 and a host of other recent methods based on path sampling reviewed, for example, in Ref. 9.

A common feature of several of these recent methods for the enhanced calculation of pathways is the assignment of interfaces that separate regions of conformational space. Trajectory segments are run between these interfaces, then the segments are combined all together to obtain the entire kinetics of interest. For example, in the weighted ensemble method10 (and its further developments11), one runs trajectories that are re-spawned at interfaces. Other prominent examples are the transition interface sampling (TIS) method of van Erp, Moroni, and Bolhuis12 and the milestoning method of Elber and co-workers.13,14 In both milestoning and TIS (as well as in the forward flux sampling method,15 which is related to TIS, or in the boxed molecular dynamics method16), one places interfaces and runs trajectories between the interfaces. The trajectory segments can be generated via integrating the natural dynamical equations of the system or can themselves be accelerated and reweighted to get the correct kinetics.17–21 Although for two-state systems transition state theory and related methods have been successfully used in the determination of the kinetics for many systems with well-defined reactant and product states, for which the dynamical bottleneck can be identified,22 there are many interesting problems in biophysics, and elsewhere, for which these assumptions do not hold. By contrast, transition path sampling approaches require no intuition for reaction mechanisms or advance knowledge of transition state, although the requirement of a dynamical bottleneck does persist.7,23 In the same category of methods is the milestoning algorithm of Elber and co-workers, which is a method for calculating kinetic properties, where the fundamental objects are the first passage time distributions KAB(τ) between adjacent conformation-space milestone states (configurations A and B in this case), where the milestone states do not necessarily need to be meta-stable states as in transition state theory. The key feature of the milestoning method (and also of the other interfacing methods) is that long trajectory pathways for large scale configuration changes can be broken up into shorter trajectories for which a linear network of transition probabilities between milestones can be devised. The aforementioned linear networks of transition probabilities can then be solved for such quantities as first passage time between any pair of milestones, including those at the extreme ends of the space and the flux through a given milestone, s, as a function of time, written as Ps(t) in Eq. (1). Some of the key gains from this treatment are that breaking up these long trajectory pathways into a network of shorter trajectories leads to increased sampling of the would-be under-sampled areas and that speed-up of computational efficiency is possible due to the capacity to run these short trajectories in parallel.13 In practice, previous milestoning calculations have been limited to calculating the constant flux values representative of the system at equilibrium, which can be thought of as the long time flux values limtPs(t). A method for calculating the time-dependent flux through a given milestone Ps(t) can be found in Paper II.24 The aim of the technique we present in Paper I is to increase the computational speed of the milestoning method via the addition of an artificial constant force (Fwind) along the vector pointing from the initial state to the final state for each pair of milestones in the simulation, causing the system to arrive at the destination configuration in far fewer time steps than if it were left to Brownian dynamics alone. The key idea which makes this possible is the use of a re-weighting function we have introduced previously,19,20,25,26 which generates a re-weighting coefficient for each trajectory, thus allowing the true distribution of first passage times to be recovered from the artificially accelerated trajectories. Preliminary calculations conducted on model systems, described in Sec. III, have demonstrated a computation time speedup by a factor of nearly 40 using this method.

The central quantity in milestoning is the flux through a given milestone;13 it is prescribed by the probabilities

Ps(t)=0tQs(t)10ttKs(τ)dτdt,
Qs(t)=2δ(t)Ps(0)+0tQs±1(t)Ks±1(tt)dt,
(1)

where Ps(t) is the probability of being at milestone s at time t (or, more specifically, arriving at any time t′ < t and not leaving before time t13), Qs(t) is the probability of a transition to milestone s at time t, and Ks(τ) is the probability of transitioning out of milestone s after an incubation time of τ. Thus 0ttKs(τ)dτ is the probability of an exit from milestone s anytime between 0 and tt′, which makes 10ttKs(τ)dτ the probability of there not being an exit from milestone s over that same time period. Since the probability of two independent simultaneous events happening concurrently is the product of the two events, the equation for Ps(t) is simply integrating the concurrent probabilities of arriving at milestone s and not leaving over the time frame from time 0 to t. In dissecting the meaning of Qs(t), the first term, 2δ(t)Ps(0), simply represents the probability that the system is already occupying milestone s at time t = 0, where the factor of 2 is present since the δ function is centered at zero, meaning only half of its area would be counted without this factor. Qs±1(t″) is the probability that the system transitioned into one of the two milestones adjacent to s at an earlier time t″. Ks±1(tt) is the probability of a transition from milestones s ± 1 into milestone s. Thus the second term of the second line of Eq. (1) is another concurrent probability: the probability of the system entering an adjacent milestone at an earlier time, and then transitioning into milestone s between time t and 0. It is important to note that all functions Ps(t) and Qs(t) are calculated using the respective values of Ks(τ) between adjacent milestones, and thus the set of Ks(τ) between all milestones of interest contains all the information needed to calculate kinetics using the milestoning method.

The essential connection to make in regard to combining the milestoning method with re-weighting of artificially accelerated trajectories is that a K function between two milestones located at x = A and x = B, KAB(τ), is nothing more than a probability distribution as a function of lifetime describing the conditional probability that a system found in state A at time t = 0 will be found, for the first time, in state B at time t = τ,

KAB(τ)=P(xB,τ|xA,0).
(2)

Given this relationship, we can now begin to make the connection between milestoning and re-weighting of artificially accelerated trajectories. Assuming Langevin dynamics with the addition of a wind force,

mẍ=γmẋV(x)+ξ(t)+Fwind,
(3)

where γ is the friction coefficient, V(x) is the potential, ξ(t) is the random force, and Fwind is a constant force applied in the direction of the goal milestone for each run; conditional probabilities reflecting first passage transitions from milestone A to B can be expressed as

P(xB,τ|xA,0)=DξW[ξ]δ(x(τ)xB).
(4)

In this equation, W[ξ(t)] is the probability distribution representing the joint probabilities of all possible series of random kicks, so multiplying by the delta function δ(x(τ) − xB) and integrating select for only the portion of the distribution which represents a series of random kicks, which results in a transition from state A to state B given an incubation time τ. It then follows suit that the integral in this equation is simply the expectation value for the probability of a transition from A to B for each incubation time point τ, which, again, is equivalent to KAB(τ). Because of fluctuation dissipation, ⟨ξ(t)ξ(t′)⟩ = 2kBTmγ, and assuming a normally distributed stochastic term, the random force in Langevin dynamics is a Gaussian distribution with variance w ≡ 2kBTmγ. Thus it can be shown that

W[ξ(t)]=exp(12w0tξ(t)2dt).
(5)

With W[ξ(t)], our weighting function for joint probabilities of random kick sequences in terms of our random force ξ(t) in hand, we can now write the noise history ξ(t) in terms of the trajectory x(t), it generates assuming overdamped Langevin dynamics

ξ(t)2=(γmẋ+V(x)Fwind)2.
(6)

We are ultimately interested in measuring conditional probability distributions in configuration space, x(t), not random force space, ξ(t). Additionally considering the fact that the functional Jacobian for the transformation is constant27 and can be thought of as built into the measure of Dx, we can define S[x(t)] thusly

S[x(t)]ξ(t)2dt=(γmẋ+V(x)Fwind)2dt.
(7)

We can then express our conditional probability using the Wiener formalism of path integrals28 as

P(xB,τ|xA,0)=(xA,0)(xB,τ)DxexpS[x(t)]2w,
(8)

where

W[x(t)]expS[x(t)]2w
(9)

is the weight of each trajectory x(t). We can now define the re-weighting factor for obtaining the true weight of an artificially accelerated trajectory as

W[x(t)]Wf[x(t)]=expS[x(t)]Sf[x(t)]2w,
(10)

where the f subscript indicates a function generated under the influence of the artificial Fwind force. In practice, once a trajectory x(t) is generated (in the presence of the wind force), the actions are calculated in discrete numerical form using

Sf[x(t)]imγΔxiΔt+ViFwind2Δt,
S[x(t)]imγΔxiΔt+Vi2Δt,
(11)

using the Ito discretization scheme29 for stochastic calculus. Although in the continuum limit, depending on the rules of stochastic calculus (i.e., Ito vs. Stratonovich), one may need to account for an extra term involving the second derivative of the potential in the integrand of the action employed in Eq. (7), the discrete Ito version utilized in this work does not require it and can be used as expressed in Eq. (11).30 

The re-weighting factor is calculated from Eq. (10) and stored in an array. When post-processing to compute the KAB(τ) distribution for a particular pair of milestones A and B by histogramming trajectories by lifetime τ, instead of adding 1 to a particular bin each time the lifetime of a particular trajectory falls within the bounds of that bin, the weight W[x(t)] corresponding to that trajectory is instead added. It is clear from Eq. (10) that as Sf[x(t)] for the artificially accelerated trajectory approaches S[x(t)], the weight of the trajectory approaches unity, and thus the method reduces to an unweighted histogram for Fwind=0 as it should.

In practice, given a well-chosen wind coordinate that significantly accelerates trajectories from the initial milestone to the final milestone, it is not surprising that shorter and lower weight trajectories are sampled with higher frequency. This can lead to re-weighted transition time distributions that have excellent convergence and accuracy in the regime of fast transitions and poorer convergence and accuracy in the tails (the integrand of the action integral is positive definite so the weights of trajectories increase exponentially with time), which can be seen in some of the following numerical demonstrations.

A simple two well potential (see the inset of Fig. 1) of the equation y = (x − 1)2(x + 1)2 was chosen to be the model system upon which the wind-assisted milestoning methodology could be developed. In running wind-assisted milestoning, the potential to which the particle is being subjected is first divided into any number of milestones, in this case, 7 milestones, thus 6 separate inter-plane spaces. Next, numerous Langevin trajectories are run both from left to right and right to left between each pair of adjacent milestones. The number of time steps required to go from the starting milestone to the destination milestone for each trial of each pair and the weight of each trajectory is then stored in an array as mentioned in Sec. II. As shown in Figs. 1 and 2, this method has demonstrated a roughly fivefold speedup with no discernible sacrifice to accuracy and a more than tenfold speedup if one is willing to sacrifice some accuracy in the tail of the distribution.

FIG. 1.

Calculation time as a function of the Fwind force for all K(τ) distributions in both directions for six subspaces ranging from −2 to 2 on the bistable harmonic potential (inset). The calculation took 507 s for Fwind=0pN and just 44 s for Fwind=12pN.

FIG. 1.

Calculation time as a function of the Fwind force for all K(τ) distributions in both directions for six subspaces ranging from −2 to 2 on the bistable harmonic potential (inset). The calculation took 507 s for Fwind=0pN and just 44 s for Fwind=12pN.

Close modal
FIG. 2.

The transition probability distribution K23(τ), i.e., the transition probability from milestone 2 to milestone 3 as a function of lifetime, calculated using Fwind forces ranging from 0 to 12 pN. The plots indicate that the rapid decrease in computation time due to the added Fwind force has little effect on accuracy.

FIG. 2.

The transition probability distribution K23(τ), i.e., the transition probability from milestone 2 to milestone 3 as a function of lifetime, calculated using Fwind forces ranging from 0 to 12 pN. The plots indicate that the rapid decrease in computation time due to the added Fwind force has little effect on accuracy.

Close modal

Thus far, we have approached the wind-assisted reweighted milestoning (WARM) method from the standpoint of speeding up the calculation by pushing Fwind until the K(τ) functions begin to distort. Here we will explore the possibility that even slightly distorted K(τ) functions can yield useful information, allowing for even greater computational speedup. The flux value for a given milestone s, Ps(t), should approach the probability predicted by the Boltzmann distribution generated from the configurational partition function as time approaches infinity. Given a discrete space in x, subject to our 1D potential y = (x − 1)2(x + 1)2, the Boltzmann distribution function can be obtained in the usual way, shown in Eq. (12) below,

limtPs(t)=eβU(xs)n=1NSeβU(xn),
(12)

where Ns is the total number of milestone configurations and xn signifies the spatial position of each milestone. This discrete space approximation for the equilibrium flux values is utilized below as a test for accuracy in Fig. 3 (dashed lines). Section III consists of dividing the space for the bistable 1D potential between x = −2 and x = 2 into 11 subspaces bounded by 12 milestones. First hitting trajectories were run between each pair of adjacent milestones, and then each pair of K(τ) functions describing a transition away from each milestone was normalized (e.g., for milestone 3, all trajectories must terminate at either milestone 2 or milestone 4, therefore 0K32(τ)dτ+0K34(τ)dτ=1). The normalized K(τ) functions are then integrated over all time, and these values are placed in a matrix, K, of equilibrium transition probabilities. The equilibrium flux values for the vector representing the set of milestones, P, are then found by numerically solving for the eigenvector: P·K = P, with reflecting boundary conditions implemented at the extremes.31 By using this method to determine equilibrium flux values, it is demonstrated in Figs. 3 and 4, that even when Fwind is set to a value strong enough to distort the K(τ) functions, accurate equilibrium flux values can still be calculated.

FIG. 3.

The top of this figure shows plots for one of the transition probability distributions K(τ) for the bistable 1D potential with different values of Fwind implemented. Note that although the distributions distort considerably for higher values of τ when the system is pushed with high magnitude Fwind, the equilibrium flux values in the plot below remain fairly constant. The dashed line indicates the equilibrium probabilities of the system occupying each of the milestone configurations, as predicted by the Boltzmann distribution. Color scheme legend applies to both plots.

FIG. 3.

The top of this figure shows plots for one of the transition probability distributions K(τ) for the bistable 1D potential with different values of Fwind implemented. Note that although the distributions distort considerably for higher values of τ when the system is pushed with high magnitude Fwind, the equilibrium flux values in the plot below remain fairly constant. The dashed line indicates the equilibrium probabilities of the system occupying each of the milestone configurations, as predicted by the Boltzmann distribution. Color scheme legend applies to both plots.

Close modal
FIG. 4.

Showcase of the effects of applying higher magnitude Fwind which are strong enough to significantly distort the K(τ) functions. This figure facilitates a direct comparison of gain in computational speed with the accuracy of the equilibrium flux values (measured as X2). Note that while there is no appreciable change in accuracy, calculation time drops from 1109 s to 26 s, a speedup by a factor of nearly 40.

FIG. 4.

Showcase of the effects of applying higher magnitude Fwind which are strong enough to significantly distort the K(τ) functions. This figure facilitates a direct comparison of gain in computational speed with the accuracy of the equilibrium flux values (measured as X2). Note that while there is no appreciable change in accuracy, calculation time drops from 1109 s to 26 s, a speedup by a factor of nearly 40.

Close modal

Two additional test systems for the WARM technique were implemented for further validating the method in two dimensions. Both systems have double well shapes; however for one well, the barrier to transition from one well to the other is primarily energetic, while the other is primarily entropic (see Fig. 5). The potential with the energetic barrier is a generalization of the 1D potential from Sec. III B to two dimensions, and the potential with the entropic barrier is the same potential implemented by Elber and Faradjian in their original paper on milestoning.13 As can be seen in the data below, the WARM method successfully re-weighted first passage time distributions (K(τ)) generated using artificially accelerated trajectories to yield the true first passage time distributions which would have resulted from trajectories in the absence of the wind force (Figs. 6 and 7). In both cases, the method achieved more than 60% faster computation times with very little sacrifice in terms of accuracy (Figs. 8 and 9).

FIG. 5.

Potentials used in the 2D WARM calculations. In the first case, the primary barrier to crossing from one well to the other is the height of the barrier relative to the strength of the “kicks” from the random force in the Langevin equation. In the second potential,13 the barrier to crossing between wells is entropic, in that a trajectory which results in a transition between wells must find its way through the gap at the center; i.e., the likelihood of a transition is not limited by any sort of uphill battle but instead by decreased degeneracy in the number of possible trajectories which result in a transition.

FIG. 5.

Potentials used in the 2D WARM calculations. In the first case, the primary barrier to crossing from one well to the other is the height of the barrier relative to the strength of the “kicks” from the random force in the Langevin equation. In the second potential,13 the barrier to crossing between wells is entropic, in that a trajectory which results in a transition between wells must find its way through the gap at the center; i.e., the likelihood of a transition is not limited by any sort of uphill battle but instead by decreased degeneracy in the number of possible trajectories which result in a transition.

Close modal
FIG. 6.

The transition probability distribution K12(τ), i.e., the transition probability from milestone 1 (the line x = −1) to milestone 2 (the line x = 0) on the 2D potential with the energetic barrier as a function of lifetime, calculated using Fwind forces ranging from 0 to 12 pN. The plots indicate that the rapid decrease in computation time due to the added wind force has almost no effect on accuracy.

FIG. 6.

The transition probability distribution K12(τ), i.e., the transition probability from milestone 1 (the line x = −1) to milestone 2 (the line x = 0) on the 2D potential with the energetic barrier as a function of lifetime, calculated using Fwind forces ranging from 0 to 12 pN. The plots indicate that the rapid decrease in computation time due to the added wind force has almost no effect on accuracy.

Close modal
FIG. 7.

The transition probability distribution K12(τ), i.e., the transition probability from milestone 1 (the line x = −0.5) to milestone 2 (the line x = 0) on the 2D potential with the entropic barrier as a function of lifetime, calculated using Fwind forces ranging from 0 to 1 pN. The plots indicate that the rapid decrease in computation time due to the added Fwind force has almost no effect on accuracy.

FIG. 7.

The transition probability distribution K12(τ), i.e., the transition probability from milestone 1 (the line x = −0.5) to milestone 2 (the line x = 0) on the 2D potential with the entropic barrier as a function of lifetime, calculated using Fwind forces ranging from 0 to 1 pN. The plots indicate that the rapid decrease in computation time due to the added Fwind force has almost no effect on accuracy.

Close modal
FIG. 8.

Calculation time as a function of the Fwind force for all K(τ) distributions in both directions for two subspaces ranging from −1 to 1 on the x axis of the 2D potential with the energetic barrier. All trajectories were run using β = 0.123. The highest value of Fwind yielded a faster computation time by a factor of 4.17 than the unassisted calculation with very little distortion to the K(τ) function.

FIG. 8.

Calculation time as a function of the Fwind force for all K(τ) distributions in both directions for two subspaces ranging from −1 to 1 on the x axis of the 2D potential with the energetic barrier. All trajectories were run using β = 0.123. The highest value of Fwind yielded a faster computation time by a factor of 4.17 than the unassisted calculation with very little distortion to the K(τ) function.

Close modal
FIG. 9.

Calculation time as a function of the Fwind force for all K(τ) distributions in both directions for two subspaces ranging from −0.5 to 0.5 on the x axis of the 2D potential with the entropic barrier. All trajectories were run using β = 3.0 so as to ensure that transitions over the barrier instead of through the small gap were highly unlikely. The highest value of Fwind yielded a faster computation time by a factor of 4.78 than the unassisted calculation with almost no distortion to the K(τ) function.

FIG. 9.

Calculation time as a function of the Fwind force for all K(τ) distributions in both directions for two subspaces ranging from −0.5 to 0.5 on the x axis of the 2D potential with the entropic barrier. All trajectories were run using β = 3.0 so as to ensure that transitions over the barrier instead of through the small gap were highly unlikely. The highest value of Fwind yielded a faster computation time by a factor of 4.78 than the unassisted calculation with almost no distortion to the K(τ) function.

Close modal

In order to demonstrate that the WARM method possesses no inherent limitations due to scaling, the method was applied to an 11 dimensional hyperspace. For this model, the 11D potential was defined as

V(x1,x2,,x11)=(x11)2(x1+1)212n=211xn2x12+n=211xn4,
(13)

where the first term is the same bistable potential in x1 used in the first one dimensional example, the second term couples motion in the 10 dimensions orthogonal to barrier height in x1, and the third term simply confines the system to a reasonably sized configurational space using a quartic potential. In order to develop some intuition for this potential, see Fig. 10, then just imagine that there are nine other dimensions which have the same effect as y on barrier height in x1.

FIG. 10.

2D representation of the 11D coupled potential. The y in the second term (red) has been left as a parameter in this plot. The surfaces shown are for values for the parametric y of 0, ±1, and ±1.5, where the deepest well corresponds to y = 1.5 and the shallowest corresponds to parametric y = 0. Just as the well becomes deeper, the further from the system wanders from the origin in the y direction in this 2D model, the 11D system also encounters deeper wells in the xn dimensions the further it wanders from the origin in each xn dimension.

FIG. 10.

2D representation of the 11D coupled potential. The y in the second term (red) has been left as a parameter in this plot. The surfaces shown are for values for the parametric y of 0, ±1, and ±1.5, where the deepest well corresponds to y = 1.5 and the shallowest corresponds to parametric y = 0. Just as the well becomes deeper, the further from the system wanders from the origin in the y direction in this 2D model, the 11D system also encounters deeper wells in the xn dimensions the further it wanders from the origin in each xn dimension.

Close modal

Accordingly, the milestones must be defined as hyperplanes, given the general definition of a hyperplane

a1x1+a2x2+a3x3++anxn=b.
(14)

To keep things as simple as possible, we set a2 through a11 equal to zero, and a1 = 1, allowing us to define two hyperplanes as x1 = −1 and x1 = 1. In this scenario, the features of interest are the transitions between the wells at x1 = −1 and x1 = 1, and thus the 11D Fwind is applied with zero components in all dimensions except for x1 where it is used to push the system over the barrier between wells. Trajectories were initialized at the point {−1, 0, 0, …, 0} and terminated upon reaching their hyperdimensional “finish line,” i.e., when the hyperplane a1x1 + a2x2 + a3x3 + ⋯ + anxn < 1 evaluates to false. The WARM method was successfully applied to this 11 dimensional potential, and a speedup by a factor of 4.5 was observed (Figs. 11 and 12).

FIG. 11.

CPU time as a function of the magnitude of the Fwind in 11D. The maximum speedup measured was a factor of 4.5.

FIG. 11.

CPU time as a function of the magnitude of the Fwind in 11D. The maximum speedup measured was a factor of 4.5.

Close modal
FIG. 12.

Plot of the K(τ) functions generated for each data point in the CPU time vs. Fwind plot for the 11D system.

FIG. 12.

Plot of the K(τ) functions generated for each data point in the CPU time vs. Fwind plot for the 11D system.

Close modal

In all of the preceding examples, Fwind was applied to the system as a constant force applied in a straight line, perpendicular to the parallel milestone hyperplanes, but Fwind can be defined any way we choose. This is particularly useful when the milestones are carved using Voronoi polyhedra or other more complicated tesselations. This section demonstrates a method whereby the directionality of Fwind is defined by a vector field which allows Fwind to blow in a curved path between two nearly orthogonal milestones (see Figs. 13 and 14); i.e., our wind has become a tornado. In order to define this vector field, the point of intersection between the two planes was determined, then a function was created which finds the straight line connecting the current position to this point of intersection, and then defines Fwind at that point to be a vector both orthogonal to that straight line and pointing in a clockwise direction. When first passage times were calculated going from the milestone shown in red toward the milestone shown in green (Fig. 13), the vector field is simply multiplied by −1 to cause our tornado to spin counterclockwise. Using a wind force defined in this manner, we obtain an efficient directionality for Fwind which biases the system toward both leaving its initial milestone in the right direction and approaching its destination milestone, regardless of the positioning of the milestones in configuration space. Another advantage of this scheme is that our curved vector field of Fwind can be defined without any knowledge of the system itself, we only need to know the positions of the milestones, which are always known in milestoning calculations. Figures 13–18 illustrate the application and results of this method using two different 2D potentials, the Muller-Brown potential,32 and a simpler Muller-inspired potential with two Gaussian wells we will refer to as our Gaussian potential. The Gaussian potential is defined as

V(x,y)=exp[(2(x0.8)2+y2)]1.3exp(x+1)2+(y1.5)2+0.2x2+0.2y2
(15)

and the Muller potential is defined as

V(x,y)=hk=14exp[ak(xxk0)2+bk(xxk0)(yyk0)+ck(yyk0)2],
(16)

where

A=(200,100,170,15),a=(1,1,6.5,0.7),b=(0,0,11,0.6),c=(10,10,6.5,0.7),x0=(1,0,0.5,1),y0=(0,0.5,1.5,1),h=0.005.

A speedup factor of 4 was observed in both the Muller potential and the Gaussian potential (Figs. 15 and 17), although the K(τ) functions in the Gaussian potential example displayed less distortion than those produced in the calculations performed using the Muller potential.

FIG. 13.

Representation of the vector field approach to applying Fwind to push milestoning trajectories between two nearly orthogonal planes, subject to our Gaussian potential. The green milestone is defined as the plane for which y44x=0.7, and the red milestone is defined as the plane for which y = 1.5. The vector wind is configured to show the Fwind scheme for accelerating trajectories going from red to green.

FIG. 13.

Representation of the vector field approach to applying Fwind to push milestoning trajectories between two nearly orthogonal planes, subject to our Gaussian potential. The green milestone is defined as the plane for which y44x=0.7, and the red milestone is defined as the plane for which y = 1.5. The vector wind is configured to show the Fwind scheme for accelerating trajectories going from red to green.

Close modal
FIG. 14.

Plot of the same milestone placement and Fwind scheme as the Gaussian potential example applied to the Muller potential and with a directionality for accelerating trajectories from the green milestone to the red one.

FIG. 14.

Plot of the same milestone placement and Fwind scheme as the Gaussian potential example applied to the Muller potential and with a directionality for accelerating trajectories from the green milestone to the red one.

Close modal
FIG. 15.

CPU time as a function of Fwind magnitude for the Gaussian potential.

FIG. 15.

CPU time as a function of Fwind magnitude for the Gaussian potential.

Close modal
FIG. 16.

Plot of the K(τ) functions corresponding to different magnitudes of Fwind as applied to the Gaussian potential.

FIG. 16.

Plot of the K(τ) functions corresponding to different magnitudes of Fwind as applied to the Gaussian potential.

Close modal
FIG. 17.

CPU time as a function of Fwind magnitude for the Muller potential.

FIG. 17.

CPU time as a function of Fwind magnitude for the Muller potential.

Close modal
FIG. 18.

Plot showing the K(τ) functions corresponding to different magnitudes of Fwind as applied to the Muller potential.

FIG. 18.

Plot showing the K(τ) functions corresponding to different magnitudes of Fwind as applied to the Muller potential.

Close modal

We have proposed and tested WARM, a wind-assisted reweighting method for accelerating milestoning calculations, whereby the true probability density functions for first passage transition times between milestones, KAB(τ), are recovered from artificially accelerated trajectories via the re-weighting method described in Sec. II. These KAB(τ) functions are central to milestoning calculations, and thus the WARM method presented herein shows potential for broad applications.

Our method has been shown to be effective in tests on one and two dimensional potentials with both energetic and entropic barriers, as well as on multi-dimensional models, implying that the method should have no scaling limitations. Thus the next step will be to test the method on more complex biomolecular transitions. The simplest extension would be to apply a force vector to a single atom which steers the system toward the configurational change of interest. In this case, the re-weighting factors S[x(t)] and Sf[x(t)] could be calculated by including the 3d force components both with and without the components of the applied force, respectively.

It should be noted that our application of the WARM method to both the high-dimensional model and our vector field-based Fwind implementation demonstrate that this technique can be applied to systems that are too complex to intuit the placement of the artificial forces. One such example for which the vector-field method can be applied concerns rotational local dynamics in DNA and RNA.33–35 Given an initial and a final milestone configuration, one could determine the vectors curving around the field from each atom’s initial position to its final position. Artificial forces, Fwind, could then be placed upon all atoms in the system pointing in the direction of these vectors and with a magnitude proportional to the length of the vectors. A zero cutoff could also be added so as not to waste computational resources on applying and accounting for Fwind force atoms which are beginning at a position fairly close to their destination. Another way to construct clever wind force biases could involve guessing the directions of a slow manifold, e.g., by averaging out the momenta over times that are longer than fast periods but smaller than the slow periods, as in momentum-enhanced self-guiding methods.21,36,37 The wind force could then be applied only in the direction of the slow manifold.

The main limitation of the WARM method, regardless of the number of dimensions present, is obtaining good re-weighting in the longer τ range. This is simply a matter of under-sampling, which leads to non-overlap between the probability of the bundles of trajectories with and without wind (probability non-overlap, by the way, is also the limiting factor for umbrella sampling and related techniques in the realm of thermodynamics). If Fwind is pushing the system to the next milestone so quickly that longer values of τ, relevant to the true KAB(τ) distribution, are not being sampled, then there just is not enough density present (or even none at all, in numerical implementations) to re-weight. This is why the accuracy in the low tau regime is often still quite good when too high of an Fwind has caused the latter portion of the distribution to turn to noise. Thus far, the limitations on the WARM method appear to be solely dependent on whether or not we have pushed the force so hard that trajectories in the longer τ region of the true KAB(τ) are even being sampled. For this reason, systems whose true KAB(τ) distribution functions possess thick tails place the greatest limitations on the degree of computational speedup achievable by WARM. This issue can be addressed in a couple of different ways. One approach is to simply define more milestones in the space, the other is to combine the WARM method with some sort of artificial heating methods, both modifications which will yield KAB(τ) functions which decay more rapidly after their peak.

We believe that, upon implementation into molecular dynamics packages on top of modules that perform interface sampling, the WARM method has the potential to be a useful tool for the determination of the long-time relaxation properties of macromolecules.

Funding from the National Institutes of Health, Grant No. R01GM089846, and from the National Science Foundation, Award No. CMMI-1404818, is kindly acknowledged.

1.
C.
Dellago
,
P. G.
Bolhuis
,
F. S.
Csajka
, and
D.
Chandler
, “
Transition path sampling and the calculation of rate constants
,”
J. Chem. Phys.
108
(
5
),
1964
1977
(
1998
).
2.
S.
Kirmizialtin
and
R.
Elber
, “
Revisiting and computing reaction coordinates with directional milestoning
,”
J. Phys. Chem. A
115
(
23
),
6137
6148
(
2011
).
3.
H.
Berendsen
, “
Protein folding—A glimpse of the holy grail?
,”
Science
282
(
5389
),
642
643
(
1998
).
4.
H.
Eyring
, “
The activated complex in chemical reactions
,”
J. Chem. Phys.
3
(
2
),
107
115
(
1935
).
5.
E.
Wigner
, “
The transition state method
,”
Trans. Faraday Soc.
34
,
29
41
(
1938
).
6.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
USA
,
2001
).
7.
P.
Bolhuis
,
D.
Chandler
,
C.
Dellago
, and
P.
Geissler
, “
Transition path sampling: Throwing ropes over rough mountain passes, in the dark
,”
Annu. Rev. Phys. Chem.
53
,
291
318
(
2002
).
8.
P.
Metzner
,
C.
Schütte
, and
E.
Vanden-Eijnden
, “
Transition path theory for Markov jump processes
,”
Multiscale Model. Simul.
7
(
3
),
1192
1219
(
2009
).
9.
L. T.
Chong
,
A. S.
Saglam
, and
D. M.
Zuckerman
, “
Path-sampling strategies for simulating rare events in biomolecular systems
,”
Curr. Opin. Struct. Biol.
43
,
88
94
(
2017
). Theory and simulation macromolecular assemblies.
10.
G. A.
Huber
and
S.
Kim
, “
Weighted-ensemble Brownian dynamics simulations for protein association reactions
,”
Biophys. J.
70
,
97
110
(
1996
).
11.
M. C.
Zwier
,
J. L.
Adelman
,
J. W.
Kaus
,
A. J.
Pratt
,
K. F.
Wong
,
N. B.
Rego
,
E.
Surez
,
S.
Lettieri
,
D. W.
Wang
,
M.
Grabe
,
D. M.
Zuckerman
, and
L. T.
Chong
, “
WESTPA: An interoperable, highly scalable software package for weighted ensemble simulation and analysis
,”
J. Chem. Theory Comput.
11
(
2
),
800
809
(
2015
).
12.
T. S.
van Erp
,
D.
Moroni
, and
P. G.
Bolhuis
, “
A novel path sampling method for the calculation of rate constants
,”
J. Chem. Phys.
118
(
17
),
7762
7774
(
2003
).
13.
R.
Elber
and
A.
Faradjian
, “
Computing time scales from reaction coordinates by milestoning
,”
J. Chem. Phys.
120
(
23
),
10880
10889
(
2004
).
14.
R.
Elber
, “
Long-timescale simulation methods
,”
Curr. Opin. Struct. Biol.
15
,
151
156
(
2005
).
15.
R. J.
Allen
,
P. B.
Warren
, and
P. R. ten
Wolde
, “
Sampling rare switching events in biochemical networks
,”
Phys. Rev. Lett.
94
,
018104
(
2005
).
16.
D. R.
Glowacki
,
E.
Paci
, and
D. V.
Shalashilin
, “
Boxed molecular dynamics: A simple and general technique for accelerating rare event kinetics and mapping free energy in large molecular systems
,”
J. Phys. Chem. B
113
(
52
),
16603
16611
(
2009
).
17.
T. B.
Woolf
, “
Path corrected functionals of stochastic trajectories: Towards relative free energy and reaction coordinate calculations
,”
Chem. Phys. Lett.
289
,
433
441
(
1998
).
18.
D. M.
Zuckerman
and
T. B.
Woolf
, “
Efficient dynamic importance sampling of rare events in one dimension
,”
Phys. Rev. E
63
,
016702
(
2001
).
19.
C.
Xing
and
I.
Andricioaei
, “
On the calculation of time correlation functions by potential scaling
,”
J. Chem. Phys.
124
(
3
),
034110
(
2006
).
20.
J.
Nummela
and
I.
Andricioaei
, “
Exact low-force kinetics from high-force single-molecule unfolding events
,”
Biophys. J.
93
(
10
),
3373
3381
(
2007
).
21.
J.
MacFadyen
,
J.
Wereszczynski
, and
I.
Andricioaei
, “
Directionally negative friction: A method for enhanced sampling of rare events
,”
J. Chem. Phys.
128
,
114112
(
2008
).
22.
D. G.
Truhlar
,
B. C.
Garrett
, and
S. J.
Klippenstein
, “
Current status of transition-state theory
,”
J. Phys. Chem.
100
(
31
),
12771
12800
(
1996
).
23.
E.
Vanden-Eijnden
, “
Transition-path theory and path-finding algorithms for the study of rare events
,”
Annu. Rev. Phys. Chem.
61
,
391
420
(
2010
).
24.
G.
Grazioli
and
I.
Andricioaei
, “
Advances in milestoning. II. Calculating time-correlation functions from milestoning using stochastic path integrals
,”
J. Chem. Phys.
149
,
084104
(
2018
).
25.
D.
Jeong
and
I.
Andricioaei
, “
Reconstructing equilibrium entropy and enthalpy profiles from non-equilibrium pulling
,”
J. Chem. Phys.
138
(
11
),
114110
(
2013
).
26.
J.
Nummela
,
F.
Yassin
, and
I.
Andricioaei
, “
Entropy-energy decomposition from nonequilibrium work trajectories
,”
J. Chem. Phys.
128
(
2
),
024104
(
2008
).
27.
A. E.
Cardenas
and
R.
Elber
, “
Modeling kinetics and equilibrium of membranes with fields: Milestoning analysis and implication to permeation
,”
J. Chem. Phys.
141
(
5
),
054101
(
2014
).
28.
H.
Kleinert
,
Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets
, 3rd ed. (
World Scientific
,
2004
).
29.
K.
Itô
, “
On a stochastic integral equation
,”
Proc. Jpn. Acad.
22
(
1-4
),
32
35
(
1946
).
30.
A. B.
Adib
, “
Stochastic actions for diffusive dynamics: Reweighting, sampling, and minimization
,”
J. Phys. Chem. B
112
(
19
),
5910
5916
(
2008
).
31.
R.
Elber
, “
Milestoning theory and a simple example
,” http://clsb.ices.utexas.edu/web/Milestoning\%20theory\%20and\%20a\%20simple\%20example.pdf (Online; accessed September 30, 2008).
32.
K.
Müller
and
L.D.
Brown
, “
Location of saddle points and minimum energy paths by a constrained simplex optimization procedure
,”
Theor. Chim. Acta
53
(
1
),
75
93
(
1979
).
33.
E. N.
Nikolova
,
E.
Kim
,
A. A.
Wise
,
P. J.
O’Brien
,
I.
Andricioaei
, and
H. M.
Al-Hashimi
, “
Transient Hoogsteen base pairs in canonical duplex DNA
,”
Nature
470
(
7335
),
498
U84
(
2011
).
34.
H.
Zhou
,
I. J.
Kimsey
,
E. N.
Nikolova
,
B.
Sathyamoorthy
,
G.
Grazioli
,
J.
McSally
,
T.
Bai
,
C. H.
Wunderlich
,
C.
Kreutz
,
I.
Andricioaei
 et al, “
m1A and m1G disrupt A-RNA structure through the intrinsic instability of hoogsteen base pairs
,”
Nat. Struct. Mol. Biol.
23
(
9
),
803
(
2016
).
35.
Y.
Xu
,
J.
McSally
,
I.
Andricioaei
, and
H. M.
Al-Hashimi
, “
Modulation of hoogsteen dynamics on DNA recognition
,”
Nat. Commun.
9
(
1
),
1473
(
2018
).
36.
I.
Andricioaei
,
A. R.
Dinner
, and
M.
Karplus
, “
Self-guided enhanced sampling methods for thermodynamic averages
,”
J. Chem. Phys.
118
,
1074
1084
(
2003
).
37.
J.
MacFadyen
and
I.
Andricioaei
, “
A skewed-momenta method to efficiently generate conformational-transition trajectories
,”
J. Chem. Phys.
123
,
074107
(
2005
).