Advancements in Milestoning I: Accelerated Milestoning via "Wind" Assisted Re-weighted Milestoning (WARM)
AAdvancements in Milestoning I: Accelerated Milestoning via “Wind” AssistedRe-weighted Milestoning (WARM)
Gianmarc Grazioli and Ioan Andricioaei ∗ Department of Chemistry, University of California, Irvine, CA 92697 (Dated: August 20, 2018)The Milestoning algorithm created by Ron Elber et al. is a method for determining the time scaleof processes too complex to be studied using brute force simulation methods. The fundamentalobjects implemented in the Milestoning algorithm are the first passage time distributions K AB ( τ )between adjacent protein configuration milestones A and B . The method proposed herein aims tofurther enhance Milestoning by employing an artificial applied force, akin to wind, which pushes thetrajectories from their initial states to their final states, and subsequently re-weights the trajectoriesto yield the true first passage time distributions K AB ( τ ) in a fraction of the computation timerequired for unassisted calculations. The re-weighting method, rooted in Ito’s stochastic calculus,was adapted from previous work by Andricioaei et al. The theoretical basis for this technique andnumerical examples are presented. PACS numbers: Valid PACS appear here
I. INTRODUCTION
The task of calculating kinetic properties from molecu-lar dynamics simulations is a complex problem of consid-erable interest [1] [2]. In contrast to computational meth-ods designed for equilibrium calculations, in which thebasic observables are thermodynamic averages over con-formational points (structures) generated over an invari-ant measure without the need to obey exact dynamicalequations, studies of kinetics require physically correcttime-ordered trajectories to obtain time-correlation func-tions as the basic objects [3]. Since each time-correlationfunction describes the relaxation under investigation asan average over all relevant trajectories, adequate sam-pling for accurate calculation of long-time kinetics canquickly become computationally intractable via directsimulation [4]. This is because a direct, brute forcemethod of this type would require sufficiently long simu-lation times such that the system would be able to transi-tion between states of interest enough times that a statis-tically significant distribution of first passage times couldbe generated. Several computational methods have beendeveloped to address the challenge of calculating chemi-cal kinetics, starting with the venerable transition statetheory (TST) [5] [6], and continuing with more recent de-velopments, such as, transition path sampling (TPS) [7],transition path theory (TPT) [8], and transition inter-face sampling (TiS) [9]. Although transition state theoryhas been successfully used in the determination of thekinetics for many systems with well-defined reactant andproduct states, for which the “dynamical bottleneck” canbe identified [10], there are many interesting problems inbiophysics, and elsewhere, for which these assumptionsdo not hold. In contrast, transition path sampling ap-proaches require no intuition for reaction mechanisms or ∗ Electronic address: [email protected] advance knowledge of transition state, although the re-quirement of a ”dynamical bottleneck” does persist [7][11]. In the same category of methods is the milestoningalgorithm created by Ron Elber et al., which is a methodfor calculating kinetic properties, where the fundamentalobjects are the first passage time distributions K AB ( τ )between adjacent protein configuration milestone states(configurations A and B in this case), where the mile-stone states do not necessarily need to be meta-stablestates as in transition state theory. The key feature of themilestoning method is that long trajectory pathways forlarge scale configuration changes can be broken up intoshorter trajectories for which a linear network of tran-sition probabilities between milestones can be devised.The aforementioned linear networks of transition proba-bilities can then be solved for such quantities as first pas-sage time between any pair of milestones, including thoseat the extreme ends of the space, and the flux througha given milestone, s , as a function of time, written as P s ( t ) (equation 1). Some of the key gains from this treat-ment are that breaking up these long trajectory pathwaysinto a network of shorter trajectories leads to increasedsampling of the would-be under-sampled areas, and thatgains in computational efficiency are possible due to thecapacity to run these short trajectories in parallel [12].In practice, previous milestoning calculations have beenlimited to calculating the constant flux values represen-tative of the system at equilibrium, which can be thoughtof as the long time flux values lim t →∞ P s ( t ). A methodfor calculating the time-dependent flux through a givenmilestone P s ( t ) can be found in the companion paperto this article, also in this publication [13]. The aim ofthe technique we present in this paper is to increase thecomputational speed of the milestoning method via theaddition of an artificial constant force ( F wind ) along thevector pointing from the initial state to the final statefor each pair of milestones in the simulation, causing thesystem to arrive at the destination configuration in farfewer time steps than if it were left to Brownian dynam- a r X i v : . [ c ond - m a t . s t a t - m ec h ] O c t ics alone. The key idea which makes this possible is theuse of a re-weighting function we have introduced previ-ously [14] [15] [16] [17], which generates a re-weightingcoefficient for each trajectory, thus allowing the true dis-tribution of first passage times to be recovered from theartificially accelerated trajectories. Preliminary calcu-lations conducted on model systems, described in theNumerical Demonstration section, have demonstrated acomputation time speedup by a factor of nearly 40 usingthis method. II. THEORY
The quantity of most fundamental importance in mile-stoning is the flux through a given milestone, for whichthe equation is [12]: P s ( t ) = (cid:90) t Q s ( t (cid:48) ) (cid:34) − (cid:90) t − t (cid:48) K s ( τ ) dτ (cid:35) dt (cid:48) ,Q s ( t ) = 2 δ ( t ) P s (0) + (cid:90) t Q s ± ( t (cid:48)(cid:48) ) K ∓ s ± ( t − t (cid:48)(cid:48) ) dt (cid:48)(cid:48) (1)where P s ( t ) is the probability of being at milestone s attime t , (or, more specifically, arriving at time t (cid:48) and notleaving before time t [12]), and Q s ( t ) is the probabilityof a transition to milestone s at time t. K s ( τ ) indicatesthe probability of transitioning out of milestone s givenan incubation time of τ , thus (cid:82) t − t (cid:48) K s ( τ ) dτ is the prob-ability of an exit from milestone s anytime between 0and t − t (cid:48) , which makes 1 − (cid:82) t − t (cid:48) K s ( τ ) dτ the proba-bility of there not being an exit from milestone s overthat same time period. Since the probability of two in-dependent simultaneous events happening concurrentlyis the product of the two events, the equation for P s ( t ) issimply integrating the concurrent probabilities of arriv-ing at milestone s and not leaving over the time framefrom time 0 to t . Dissecting the meaning of Q s ( t ) (14),the first term, 2 δ ( t ) P s (0), simply represents the proba-bility 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 itsarea would be counted without this factor. Q s ± ( t (cid:48)(cid:48) ) isthe probability that the system transitioned into one ofthe two milestones adjacent to s at an earlier time t (cid:48)(cid:48) . K ∓ s ± ( t − t (cid:48)(cid:48) ) is the probability of a transition from mile-stones s ± s . Thus the second term of thesecond line of equation 14 is another concurrent proba-bility: the probability of the system entering an adjacentmilestone at an earlier time, and then transitioning intomilestone s between time t and 0. It is important to notethat all functions P s ( t ) and Q s ( t ) are calculated using therespective values of K s ( τ ) between adjacent milestones,thus the set of K s ( τ ) between all milestones of interest contains all the information needed to calculate kineticsusing the milestoning method.The important connection to make in regard to com-bining the milestoning method with re-weighting of arti-ficially accelerated trajectories is that a K function be-tween two milestones x = A and x = B , K AB ( τ ), isnothing more than a probability distribution as a func-tion of lifetime describing the conditional probability thata system found in state A at time t = 0 will be found,for the first time, in state B at time t = τ : K AB ( τ ) = P ( x B , τ | x A ,
0) (2)Given this relationship, we can now begin to make theconnection between milestoning and re-weighting of arti-ficially accelerated trajectories. Assuming Langevin dy-namics with the addition of a wind force: m ¨ x = − γm ˙ x − ∇ V ( x ) + ξ ( t ) + F wind (3)where γ is the friction coefficient, V ( x ) is the potential, ξ ( t ) is the random force, and F wind is a constant forceapplied in the direction of the goal milestone for eachrun; conditional probabilities reflecting first passage tran-sitions from milestone A to B can be expressed as: P ( x B , τ | x A ,
0) = (cid:90)
DξW [ ξ ] δ ( x ( τ ) − x B ) (4)In this equation, W [ ξ ( t )] is the probability distributionrepresenting the joint probabilities of all possible seriesof random kicks, so multiplying by the delta function δ ( x ( τ ) − x B ) and integrating selects for only the portionof the distribution which represents a series of randomkicks which results in a transition from state A to state B given an incubation time τ . It then follows suit that theintegral in this equation is simply the expectation valuefor the probability of a transition from A to B for eachincubation time point τ , which, again, is equivalent to K AB ( τ ). Because of fluctuation dissipation, (cid:104) ξ ( t ) ξ ( t (cid:48) ) (cid:105) =2 k B T mγ , the random force in Langevin dynamics is aGaussian distribution with variance w ≡ k B T mγ . Thusit can be show that: W [ ξ ( t )] = exp( − w (cid:90) t ξ ( t (cid:48) ) dt (cid:48) ) (5)With W [ ξ ( t )], our weighting function for joint proba-bilities of random kick sequences in terms of our randomforce ξ ( t ) in hand, we can now write the noise history ξ ( t ) in terms of the trajectory x ( t ) it generates: ξ ( t ) = ( m ¨ x + γm ˙ x + ∇ V ( x ) − F wind ) (6)We are ultimately interested in measuring conditionalprobability distributions in configuration space, x ( t ), notrandom force space, ξ ( t ), but since the Jacobian is builtinto the measure, x , we can define S [ x ( t )] thusly: S [ x ( t )] ≡ ξ ( t ) = ( m ¨ x + γm ˙ x + ∇ V ( x ) − F wind ) (7)Then, using the Ito formalism for stochastic calculus, wecan express our conditional probability using the Wienerformalism of path integrals [18] as: P ( x B , τ | x A ,
0) = (cid:90) ( x B ,τ )( x A , Dx exp (cid:18) − S [ x ( t )]2 w (cid:19) (8)In this form, it is clear that the exponential functionrepresents the weighting function for the trajectory x ( t ): W [ x ( t )] ≡ exp (cid:18) − S [ x ( t )]2 w (cid:19) (9)With the weight of each trajectory, W [ x ( t )], now for-mally defined, we can define the re-weighting factor forobtaining the true weight of an artificially acceleratedtrajectory as: W [ x ( t )] W f [ x ( t )] = exp (cid:18) − S [ x ( t )] − S f [ x ( t )]2 w (cid:19) (10)where the f subscript indicates a function generated un-der the influence of the artificial F wind force. In practice,once a trajectory x ( t ) is generated (in the presence of thewind force), the actions are calculated in discrete numer-ical form using: S f [ x ( t )] ≈ (cid:88) i (cid:18) m ∆ ν i ∆ t + mγ ∆ x i ∆ t + ∇ V i − F wind (cid:19) ∆ tS [ x ( t )] ≈ (cid:88) i (cid:18) m ∆ ν i ∆ t + mγ ∆ x i ∆ t + ∇ V i (cid:19) ∆ t (11)The re-weighting factor is calculated from Eq. (10) andstored in an array. When post-processing to compute the K AB ( τ ) distribution for a particular pair of milestones A and B by histogramming trajectories by lifetime τ , in-stead of adding 1 to a particular bin each time the life-time of a particular trajectory falls within the bounds ofthat bin, the weight W [ x ( t )] corresponding to that tra-jectory is instead added. It is clear from equation 10 thatas S f [ x ( t )] for the artificially accelerated trajectory ap-proaches S [ x ( t )], the weight of the trajectory approachesunity, thus the method reduces to an unweighted his-togram for F wind = 0 as it should. III. NUMERICAL DEMONSTRATIONModel System in One Dimension
A simple two well potential (see inset of figure 7) ofequation y = ( x − ( x + 1) was chosen to be themodel system upon which the wind-assisted mileston-ing methodology could be developed. In running wind-assisted milestoning, the potential to which the particle isbeing subjected is first divided into any number of mile-stones, in this case, 7 milestones, thus 6 separate spaces.Next, numerous Langevin trajectories are run both fromleft to right and right to left between each pair of adja-cent milestones. The number of time steps required to gofrom the starting milestone to the destination milestonefor each trial of each pair and the weight of each trajec-tory is then stored in an array as mentioned in the theorysection. As shown in figures 1 and 2, this method hasbrought about a more than tenfold speedup in computa-tion time with very little sacrifice in terms of accuracy. (cid:72) pN (cid:76) (cid:72) s (cid:76) FIG. 1: Shown here is calculation time as a function of the F wind force for all K ( τ ) distributions in both directions forsix subspaces ranging from -2 to 2 on the bistable harmonicpotential. The calculation took 507 seconds for F wind = 0 pN and just 44 seconds for F wind = 12 pN . Model System in One Dimension with Distortion
Thus far, we have approached the WARM methodfrom the standpoint of speeding up the calculation bypushing F wind until the K ( τ ) functions begin to distort.Here we will explore the possibility that even slightlydistorted K ( τ ) functions can yield useful information,allowing for even greater computational speedup. Theflux value for a given milestone s , P s ( t ), should approachthe probability predicted by the Boltzmann distributiongenerated from configurational partition function as timeapproaches infinity. Given a discrete space in x , subjectto our 1D potential y = ( x − ( x + 1) , the Boltzmanndistribution function can be obtained in the usual way,shown in equation 12, below: (cid:72) time steps (cid:76) F r e qu e n c y (cid:72) t r a j ec t o r i e s (cid:76) F wind (cid:61) wind (cid:61) wind (cid:61)
12 pN
FIG. 2: Shown in this figure is the transition probability dis-tribution K ( τ ), i.e. the transition probability from mile-stone 2 to milestone 3 as a function of lifetime, calculatedusing F wind forces ranging from 0 to 12 pN. The plots indi-cate that the rapid decrease in computation time due to theadded F wind force has almost no effect on accuracy. lim t →∞ P s ( t ) = e − βU ( x s ) (cid:80) N S n =1 e − βU ( x n ) (12)where N s is the total number of milestone configura-tions, and x n signifies the spatial position of each mile-stone. This discrete space approximation for the equilib-rium flux values is utilized below as a test for accuracy infigure 3 (dashed lines). The numerical demonstration inthis section consists of dividing the space for the bistable1D potential between x = − x = 2 into 11 subspacesbounded by 12 milestones. First hitting trajectories wererun between each pair of adjacent milestones, and theneach pair of K ( τ ) functions describing a transition awayfrom each milestone were normalized (e.g. for milestone3, all trajectories must terminate at either milestone 2 ormilestone 4, therefore (cid:82) ∞ K ( τ ) dτ + (cid:82) ∞ K ( τ ) dτ = 1).The normalized K ( τ ) functions are then integrated overall time, and these values are placed in a matrix, K , ofequilibrium transition probabilities. The equilibrium fluxvalues for the vector representing the set of milestones, P ,is then found by numerically solving for the eigenvector: P · K = P [19]. By using this method to determine equi-librium flux values, it is demonstrated in figures 3 and4, that even when F wind is set to a value strong enoughto distort the K ( τ ) functions, accurate equilibrium fluxvalues can still be calculated. Model Systems in Two Dimensions
Two additional test systems for the WARM techniquewere implemented for further validating the method intwo dimensions. Both systems have double well shapes,however for one well, the barrier to transition from onewell to the other is primarily energetic, while the otheris primarily entropic (see figure 3 below). The poten-tial with the energetic barrier is a generalization of the K( (cid:7590) ) (cid:7590) (timesteps)xP(x) FIG. 3: The plot at the top of this figure shows plots for one ofthe transition probability distributions K ( τ ) for the bistable1D potential with different values of F wind implemented. Notethat although the distributions distort considerably for highervalues of τ when the system is pushed with high magnitude F wind , the equilibrium flux values in the plot below remainfairly constant. The color scheme legend applies to both plots.FIG. 4: Here we show effects of applying higher magnitude F wind which are strong enough to significantly distort the K ( τ ) functions. This figure facilitates a direct comparison ofgain in computational speed with the accuracy of the equi-librium flux values (measured as X ). Note that while thereis no appreciable change in accuracy, calculation time dropsfrom 1109 s to 26 s, a speedup by a factor of nearly 40.
1D potential from the previous section to two dimen-sions, and the potential with the entropic barrier is thesame potential implemented by Elber and Faradjian intheir original paper on milestoning [12]. As can be seenin the data below, the WARM method successfully re-weighted first passage time distributions ( K ( τ )) gener-ated using artificially accelerated trajectories to yield thetrue first passage time distributions which would have re-sulted from trajectories in the absence of the wind force.In both cases, the method achieved more than 60% fastercomputation times with very little sacrifice in terms ofaccuracy. FIG. 5: Show here are the potentials used in the 2D WARMcalculations. In the first case, the primary barrier to crossingfrom one well to the other is the height of the barrier relativeto the strength of the “kicks” from the random force in theLangevin equation. In the second potential [12], the barrier tocrossing between wells is entropic, in that a trajectory whichresults in a transition between wells must find its way throughthe gap at the center, i.e. the likelihood of a transition is notlimited by any sort of uphill battle, but instead by decreaseddegeneracy in the number of possible trajectories which resultin a transition. (cid:72) time steps (cid:76) F r e qu e n c y (cid:72) t r a j ec t o r i e s (cid:76) F wind (cid:61) wind (cid:61) wind (cid:61)
12 pNF wind (cid:61)
20 pN
FIG. 6: Shown in this figure is the transition probability dis-tribution K ( τ ), i.e. the transition probability from mile-stone 1 (the line x = −
1) to milestone 2 (the line x = 0) onthe the 2D potential with the energetic barrier as a functionof lifetime, calculated using F wind forces ranging from 0 to 12pN. The plots indicate that the rapid decrease in computa-tion time due to the added wind force has almost no effect onaccuracy. (cid:72) pN (cid:76) (cid:72) s (cid:76) FIG. 7: Shown here is calculation time as a function of the F wind force for all K ( τ ) distributions in both directions fortwo subspaces ranging from -1 to 1 on the x axis of the 2Dpotential with the energetic barrier. All trajectories were runusing β = . F wind yielded a fastercomputation time by a factor of 4.17 than the unassisted cal-culation with very little distortion to the K ( τ ) function. (cid:72) time steps (cid:76) F r e qu e n c y (cid:72) t r a j ec t o r i e s (cid:76) F wind (cid:61) wind (cid:61) .7 pNF wind (cid:61) FIG. 8: Shown in this figure is the transition probability dis-tribution K ( τ ), i.e. the transition probability from mile-stone 1 (the line x = − .
5) to milestone 2 (the line x = 0) onthe the 2D potential with the entropic barrier as a functionof lifetime, calculated using F wind forces ranging from 0 to 1pN. The plots indicate that the rapid decrease in computa-tion time due to the added F wind force has almost no effecton accuracy. Model System in Eleven Dimensions
In order to demonstrate that the WARM method pos-sesses no inherent limitations due to scaling, the methodwas applied to an 11 dimensional hyperspace. For thismodel, the 11 D potential was defined as: V ( x , x , ..., x ) = ( x − ( x +1) − (cid:88) n =2 x n x + (cid:88) n =2 x n (13)where the first term is the same bistable potential in x used in the first one dimensional example, the secondterm couples motion in the 10 dimensions orthogonal tobarrier height in x , and the third term simply confinesthe system to a reasonably sized configurational spaceusing a quartic potential. In order to develop some intu- (cid:72) pN (cid:76) (cid:72) s (cid:76) FIG. 9: Shown here is calculation time as a function of the F wind force for all K ( τ ) distributions in both directions fortwo subspaces ranging from -.5 to .5 on the x axis of the 2Dpotential with the entropic barrier. All trajectories were runusing β = 3 . F wind yielded a faster computation time by afactor of 4.78 than the unassisted calculation with almost nodistortion to the K ( τ ) function.. ition for this potential, see figure 10, then just imaginethat there are nine other dimensions which have the sameeffect as y on barrier height in x . V(x,y) = (x - 1) (x + 1) - .5 y x + y FIG. 10: Show here is a 2D representation of the 11D coupledpotential. The y in the second term (red) has been left as aparameter in this plot. The surfaces shown are for valuesfor the parametric y of 0 , ± , and ± .
5, where the deepestwell corresponds to y = 1 . y = 0. Just as the well becomes deeper, thefurther from the system wanders from the origin in the y direction in this 2D model, the 11D system also encountersdeeper wells in the x n dimensions the further it wanders fromthe origin in each x n dimension. Accordingly, the milestones must be defined as hyper-planes, given the general definition of a hyperplane: a x + a x + a x ...a n x n = b (14)To keep things simple, we set a through a equal tozero, and a = 1, allowing us to define two hyperplanesas x = − x = 1. In this scenario, the features ofinterest are the transitions between the wells at x = − x = 1, thus the 11D F wind is applied with zerocomponents in all dimensions except for x where it isused to push the system over the barrier between wells.The WARM method was successfully applied to this 11dimensional potential, and a speedup by a factor of 4.5was observed (figures 11 and 12). (cid:72) pN (cid:76) (cid:72) s (cid:76) FIG. 11: This plot shows CPU time as a function of the mag-nitude of the F wind in 11D. The maximum speedup measuredwas a factor of 4.5. (cid:72) time steps (cid:76) F r e qu e n c y (cid:72) t r a j ec t o r i e s (cid:76) F wind (cid:61) wind (cid:61) wind (cid:61) wind (cid:61) FIG. 12: Shown in this figure are the K ( τ ) functions gener-ated for each data point in the CPU time vs. F wind plot forthe 11D system. Wind Force as a Vector Field
In all of the preceding examples, F wind was appliedto the system as a constant force applied in a straightline, perpendicular to the parallel milestone hyperplanes,but F wind can be defined any way we choose. This sec-tion demonstrates a method whereby the directionalityof F wind is defined by a vector field which allows F wind to blow in a curved path between two nearly orthogo-nal milestones (see figures 13 and 14), i.e. our wind hasbecome a tornado! In order to define this vector field,the point of intersection between the two planes was de-termined, then a function was created which finds thestraight line connecting the current position to this pointof intersection, and then defines F wind at that point to bea vector both orthogonal to that straight line and point-ing in a clockwise direction. When first passage timeswere calculated going from the milestone shown in red to-ward the milestone shown in green (figure 13), the vectorfield is simply multiplied by − F wind which biases the system toward both leaving its initialmilestone in the right direction and approaching its des-tination milestone, regardless of the positioning of themilestones in configuration space. Another advantage ofthis scheme is that our curved vector field of F wind canbe defined without any knowledge of the system itself, weonly need to know the positions of the milestones, whichare always known in milestoning calculations. Figures 13through 18 illustrate the application and results of thismethod using two different 2D potentials, the Muller-Brown potential [20], and a simpler Muller-inspired po-tential with two Gaussian wells we’ll call our Gaussianpotential. The Gaussian potential is defined as: V ( x, y ) = − exp [ − (2( x − . + y )] (15) − . exp [ − (( x + 1) + ( y − . )]+ . x + . y and the Muller potential is defined as: V ( x, y ) = h (cid:88) k =1 exp [ a k ( x − x k ) (16)+ b k ( x − x k )( y − y k ) + c k ( y − y k ) ]where: A = ( − , − , − , , a = ( − , − , − . , . b = (0 , , , . , c = ( − , − , − . , . x = (1 , , − . , − , y = (0 , . , . , h = . K ( τ )functions in the Gaussian potential example displayedless distortion than those produced in the calculationsperformed using the Muller potential. IV. CONCLUDING DISCUSSION
We have presented and tested a method for accelerat-ing milestoning calculations, whereby the true probabil- (cid:45) (cid:45) (cid:45) (cid:45)
FIG. 13: Shown here is a representation of the vector fieldapproach to applying F wind to push milestoning trajectoriesbetween two nearly orthogonal planes, subject to our Gaus-sian potential. The green milestone is defined as the planefor which y − x = − . y = 1 .
5. The vector wind is configuredto show the F wind scheme for accelerating trajectories goingfrom red to green. (cid:45) (cid:45) (cid:45) (cid:45) FIG. 14: This plot shows the same milestone placement and F wind scheme as the Gaussian potential example applied tothe Muller potential and with a directionality for acceleratingtrajectories from the green milestone to the red one. ity density functions for first passage transition time be-tween milestones, K AB ( τ ), are recovered from artificiallyaccelerated trajectories via the re-weighting method de-scribed in the Theory section. These K AB ( τ ) functionsare central to milestoning calculations, thus the WARMmethod presented herein shows potential for broad ap- (cid:72) pN (cid:76) (cid:72) s (cid:76) FIG. 15: This plot shows CPU time as a function of F wind magnitude for the Gaussian potential. (cid:72) time steps (cid:76) F r e qu e n c y (cid:72) t r a j ec t o r i e s (cid:76) F wind (cid:61) wind (cid:61) .5 pNF wind (cid:61) FIG. 16: This plot shows the K ( τ ) functions correspondingto different magnitudes of F wind as applied to the Gaussianpotential. plication.Our method has been shown to be effective on one andtwo dimensional potentials with both energetic and en-tropic barriers, as well as an 11 dimensional hyperspace,implying that the method should have no scaling limita-tions, thus the next step will be to test the method onchemical systems. The simplest application would be toapply a single force vector to a single atom which pushesthe system toward a configurational change of interest.In this case, the re-weighting factors S [ x ( t )] and S f [ x ( t )]could be calculated by summing the force componentsin the x , y , and z directions both with and without thecomponents of the applied force, respectively.The main limitation of the WARM method, regardlessof the number dimensions are present, is obtaining goodre-weighting in the longer τ range. This is simply a mat-ter of under-sampling. If F wind is pushing the systemto the next milestone so quickly that longer values of τ ,relevant to the true K AB ( τ ) distribution, are not beingsampled, then there just isn’t enough density present (oreven none at all) to re-weight. This is why the accu-racy in the low tau regime is often still quite good whentoo high of an F wind has caused the latter portion of thedistribution to turn to noise. Thus far, the limitationson the WARM method appear to be solely dependenton whether or not we’ve pushed the force so hard that (cid:72) pN (cid:76) (cid:72) s (cid:76) FIG. 17: This plot shows CPU time as a function of F wind magnitude for the Muller potential. (cid:72) time steps (cid:76) F r e qu e n c y (cid:72) t r a j ec t o r i e s (cid:76) F wind (cid:61) wind (cid:61) wind (cid:61) wind (cid:61) FIG. 18: This plot shows the K ( τ ) functions correspondingto different magnitudes of F wind as applied to the Muller po-tential. trajectories in the longer τ region of the true K AB ( τ )are even being sampled. For this reason, systems whosetrue K AB ( τ ) distribution functions possess fat tails placethe greatest limitations on the degree of computationalspeedup achievable by WARM. This issue can be ad-dressed in a couple different ways. One approach is tosimply define more milestones in the space, the other isto combine the WARM method with some sort of artifi-cial heating method, both modifications which will yield K AB ( τ ) functions which decay more rapidly after theirpeak.It should be noted that our application of the WARMmethod to both the high dimensional model, and ourvector field-based F wind implementation of demonstratethat this technique can be applied to systems too com-plex to intuit the placement of the artificial forces. Givenan initial and a final milestone configuration, one coulddetermine the vectors pointing from each atom’s initialposition to it’s final position. Artificial forces, F wind ,could then be placed upon all atoms in the system point-ing in the direction of these vectors and with a magnitudeproportional to the length of the vectors. A zero cutoffcould also be added so as not to waste computationalresources on applying and accounting for F wind forcesatoms which are beginning at a position fairly close totheir destination. We believe that, upon implementationinto a molecular dynamics package such as MOIL [21],the WARM method has the potential to be a useful toolfor the determination of the kinetic properties of macro-molecules. V. ACKNOWLEDGMENTS
IA acknowledges funds from an NSF CAREER award(CHE-0548047). [1] Christoph Dellago, Peter G. Bolhuis, Flix S. Csajka, andDavid Chandler. Transition path sampling and the cal-culation of rate constants.
The Journal of ChemicalPhysics , 108(5):1964–1977, 1998.[2] Christoph Dellago, Peter G. Bolhuis, and David Chan-dler. On the calculation of reaction rate constants inthe transition path ensemble.
The Journal of ChemicalPhysics , 110(14):6617–6625, 1999.[3] Serdal Kirmizialtin and Ron Elber. Revisiting and com-puting reaction coordinates with directional mileston-ing.
The Journal of Physical Chemistry A , 115(23):6137–6148, 2011. PMID: 21500798.[4] Herman Berendsen. Protein Folding - A Glimpse of theHoly Grail?
Science , 282(5389):642–643, October 1998.[5] Henry Eyring. The activated complex in chemical reac-tions.
The Journal of Chemical Physics , 3(2):107–115,1935.[6] E. Wigner. The transition state method.
Trans. FaradaySoc. , 34:29–41, 1938.[7] P. Bolhuis, D. Chandler, C. Dellago, and P. Geissler.Transition path sampling: Throwing ropes over roughmountain passes, in the dark.
Annu. Rev. Phys. Chem. ,53:291–318, 2002.[8] P. Metzner, C. Schtte, and E. Vanden-Eijnden. Transi-tion path theory for markov jump processes.
MultiscaleModeling and Simulation , 7(3):1192–1219, 2009.[9] Titus S. van Erp, Daniele Moroni, and Peter G. Bol-huis. A novel path sampling method for the calcula-tion of rate constants.
The Journal of Chemical Physics ,118(17):7762–7774, 2003.[10] Donald G. Truhlar, Bruce C. Garrett, and Stephen J.Klippenstein. Current status of transition-state the-ory.
The Journal of Physical Chemistry , 100(31):12771–12800, 1996.[11] E. Vanden-Eijnden. Transition-path theory and path-finding algorithms for the study of rare events.
Annu.Rev. Phys. Chem. , 61:391–420, 2010.[12] Ron Elber and Anton Faradjian. Computing time scales from reaction coordinates by milestoning.
Journal ofChemical Physics , 120(23):10880–9, March 2004.[13] Gianmarc Grazioli and Ioan Andricioaei. Advancementsin milestoning II: Calculating autocorrelation from mile-stoning data using stochastic path integrals in discreteconfiguration spaces. arXiv , 2015.[14] Jeremiah Nummela and Ioan Andricioaei. Exact low-force kinetics from high-force single-molecule unfoldingevents.
Biophysical Journal , 93(10):3373 – 3381, 2007.[15] Chenyue Xing and Ioan Andricioaei. On the calculationof time correlation functions by potential scaling.
TheJournal of Chemical Physics , 124(3):–, 2006.[16] Daun Jeong and Ioan Andricioaei. Reconstructingequilibrium entropy and enthalpy profiles from non-equilibrium pulling.
The Journal of Chemical Physics ,138(11):–, 2013.[17] Jeremiah Nummela, Faten Yassin, and Ioan Andri-cioaei. Entropy-energy decomposition from nonequilib-rium work trajectories.
The Journal of Chemical Physics ,128(2):–, 2008.[18] H. Kleinert.
Path Integrals in Quantum Mechanics,Statistics, Polymer Physics, and Financial Markets, 3rdEd.
World Scientific, 2004.[19] Ron Elber. Milestoning Theory and a Simple Exam-ple. http://clsb.ices.utexas.edu/web/Milestoning%20theory%20and%20a%20simple%20example.pdf . [On-line; accessed 30-Sep-2008].[20] Klaus Mller and LeoD. Brown. Location of saddle pointsand minimum energy paths by a constrained simplex op-timization procedure.
Theoretica chimica acta , 53(1):75–93, 1979.[21] Ron Elber, Adrian Roitberg, Carlos Simmerling, RobertGoldstein, Haiying Li, Gennady Verkhivker, ChenKeasar, Jing Zhang, and Alex Ulitsky. Moil: A programfor simulations of macromolecules.