Arbitrarily accurate representation of atomistic dynamics via Markov Renewal Processes
Animesh Agarwal, Sandrasegaram Gnanakaran, Nicholas Hengartner, Arthur F. Voter, Danny Perez
AArbitrarily accurate representation of atomistic dynamicsvia Markov Renewal Processes
Animesh Agarwal
Theoretical Biology and Biophysics,Los Alamos National Laboratory,New Mexico 87545, United States
Sandrasegaram Gnanakaran
Theoretical Biology and Biophysics,Los Alamos National Laboratory,New Mexico 87545, United States
Nicholas Hengartner
Theoretical Biology and Biophysics,Los Alamos National Laboratory,New Mexico 87545, United States
Arthur F. Voter
Physics and Chemistry of Materials,Los Alamos National Laboratory,New Mexico 87545, United States
Danny Perez ∗ Physics and Chemistry of Materials,Los Alamos National Laboratory,New Mexico 87545, United States a r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug bstract Atomistic simulations with methods such as molecular dynamics are extremely powerfultools to understand nanoscale dynamical behavior. The resulting trajectories, by the virtueof being embedded in a high-dimensional configuration space, can however be difficult toanalyze and interpret. This makes low-dimensional representations, especially in terms ofdiscrete jump processes, extremely valuable. This simplicity however usually comes at thecost of accuracy, as tractable representations often entail simplifying assumptions that arenot guaranteed to be realized in practice. In this paper, we describe a discretization schemefor continuous trajectories that enables an arbitrarily accurate representation in terms of aMarkov Renewal Process over a discrete state space. The accuracy of the model convergesexponentially fast as a function of a continuous parameter that has the interpretation of alocal correlation time of the dynamics. ∗ danny [email protected] . INTRODUCTION Atomistic simulations using molecular dynamics (MD) are widely recognized asextremely powerful tools to investigate the behavior of complex dynamical systemssuch as materials and biomolecules. While the predictive power of MD is extremelyhigh, the complexity of atomistic trajectories in the full phase space can be such thatinterpreting and rationalizing the results poses a significant challenge. This has madethe development of compact and interpretable representations of complex dynamicsan area of intense interest over the past few decades [1–7].In this paper, we consider coarse models that are defined on a discrete state space.This type of representation is particularly interpretable because it describes the dy-namics in terms of a sequence of waiting times punctuated by discrete jumps betweenstates, i.e., in terms of a jump process. Developing a jump process that captures thedynamics of a continuous system is a two step process that involves i) developping adiscretization procedure that maps a continuous trajectory x ( t ) into a sequence of dis-crete states and transition times ( { X n , T n } ), and ii) constructing the conditional jumpprobability distribution P r ( T n +1 − T n < t, X n +1 = j | ( X , T ) , ( X , T ) , ..., ( X n = i, T n )) that statistically reproduces the dynamics of the discretized process. Forexample, a popular approach begins by introducing a number of sets in configura-tion/phase space (which can be taken to tesselate that space, for simplicity) andby assigning to each point x ( t ) the discrete index X of the set that contains x ( t ).One then assumes that the state-to-state (i.e., set-to-set) dynamics is Markovian,i.e., P r ( T n +1 − T n < t, X n +1 = j | ( X , T ) , ( X , T ) , ..., ( X n = i, T n )) = P r ( X n +1 = j | X n = i )(1 − exp( − l i t )) = p ij F i ( t ), where l i is the total escape rate out of set i . Notethat in the Markovian case, the distribution of escape times is independent of thefinal state j . This assumption gives rise to a representation in terms of a continuous-time Markov chain (CTMC). CTMC are extremely compact, can be easily analyzed3ormally, and allow for the efficient generation of new discrete trajectories using thecelebrated BKL [8] or Stochastic Simulation algorithms [9]. These favorable charac-teristics have made CTMC representations extremely popular in the computationalsciences [6, 10–12].For all of their powerful features, CTMC representations of atomistic dynamicsare however generally not exact because the mapping from a continuous to a discretestate-space is not guaranteed to preserve the Markov property. The mapping betweena continuous and a discrete Markovian dynamics has been recently revisited using thetheory of quasi-stationary distributions (QSD) [13–15]. This analysis has shown thatthe accuracy of the Markovian representation can be quantified based on the spectralproperties of the generators of the dynamics restricted to each set by absorbingboundary conditions. This representation can be shown to becomes exact in thelimit where the spectral gap of each set becomes infinitely large [13], so that therelaxation time within each set tends to zero, which is in general not the case forfinite sets [16]. That being said, in practice, if sufficiently metastable sets can bedefined such that relaxation within every single set is fast on timescales characteristicof escapes, the corresponding Markovian representation will be very accurate.The limitation of CTMCs stems from the assumption that the state-to-state dy-namics is strictly memory-less. One possible solution is to invoke a more generalform of the transition probabilities that includes a memory of the past [16, 17]. Inthis paper, we explore a simpler alternative where the dynamics is represented interms of Markov Renewal Processes (MRPs), where Markovian constraints on theform of the probability distributions are slightly relaxed to P r ( T n +1 − T n < t, X n +1 = j | X n = i ) = p ij F ij ( t ); i.e., the difference with CTMCs is that the cumulative distri-bution of residence times T n +1 − T n given an initial state X n = i is now a generalnon-decreasing function of t and of the final state X n +1 = j such that F ij (0) = 0 andlim t →∞ F ij ( t ) = 1 [18]. Building on an analysis in terms of quasi-stationary distri-4utions (QSDs), we show that a novel discretization scheme allows for an arbitrarilyaccurate representation via MRPs for any number and shape of sets. It can thereforebe expected to be a powerful representation to accurately represent the dynamics ofa wide range of systems. II. DERIVATION
In the following, we consider the problem of obtaining a simplified representa-tion of an overdamped Langevin dynamics in R N in terms of a jump process overa discrete state space. We are specifically interested in identifying a mapping fromcontinuous dynamics to discrete jump processes that can provide an arbitrarily im-provable approximation to the exact discretized dynamics for any number and shapeof sets. We first consider such a jump process in discrete time, and then generalizeto continuous time. Note that the results derived below also hold in the contextof the Langevin equation at finite friction. The mathematical analysis is howeverextremely involved [19], so this general case is not explicitly analyzed in the follow-ing. As we will show, it also holds for any process that converges to a unique QSDwhen confined by absorbing boundaries. The approach we propose here can there-fore be expected to be useful in a very wide range of conditions, well beyond that ofobtaining a reduced representation of atomistic dynamics. A. State Definition
The first task is to define the discrete state space in which the jump processevolves. To do so, consider a number of arbitrarily disjoint connected sets Ω i definedover the N -dimensional configuration space. In contrast to the common approachdiscussed above, we assign to each x ( t ) an integer X that indexes the set where the5 i W j FIG. 1. Illustration of the trajectory discretization process. Two sets Ω i and Ω j arerepresented by ellipses and the trajectory x ( t ) by a continuous line. Each point along thetrajectory is assigned a discrete index that corresponds to the last set where it remainedfor a time t c without escaping, as indicated by its color. process last remained for at least a time t c without escaping; t c can be chosen tobe different for each set, but it will in the following be taken as a global constantfor simplicity of notation. The index assigned to a given point x ( t ∗ ) therefore corre-sponds to the last set where the trajectory x ( t ≤ t ∗ ) remained for at least t c withoutescaping. See Fig. 1 for a schematic illustration of the mapping. The state definitionis therefore not purely geometric, but depends on the past history of the process x ( t ). As will be shown below, this representation becomes especially powerful when t c is chosen to be larger than a suitably-defined local ”correlation” or ”memory” timeof the trajectory. Intuitively, this definition acts to limit the coupling between theconsecutive jumps, which is key to making the representation local in time.Indeed, the main result of this paper is that the statistics of trajectories discretizedin this way can be arbitrarily approached by a MRP by increasing t c , no matter thenumber and shapes of the sets { Ω } used in the discretization. We first considerthe discrete-time case, which we then extend to continuous time. We finally use6umerical simulation of biological systems to show that this representation can bevery accurate even for relatively small values of t c .This definition is based on insights gained through a formal analysis [13] of theParallel Replica Dynamics [20, 21] and Parallel Trajectory Splicing [22, 23] acceler-ated molecular dynamics method [24, 25] that demonstrated that the statistics oftrajectories conditioned on having remained for a long time within a set becomeespecially simple. B. Discrete time jump process
The first important result is that the discrete-time jump process that describes thesequences of states visited by the discretized process approaches a simple discrete-time Markov chain as t c → ∞ . To show this, consider the evolution of an ensembleof trajectories initialized at some point x ∈ Ω conditioned on remaining inside Ω .This problem has been rigorously investigated in Ref. 13; only a broad outline isreproduced here for completeness and details can be found in the original reference.The (unnormalized) probability distribution in configuration space evolves accordingto a Fokker-Planck equation of the form: ∂p∂t = ( −∇ V · ∇ + β − ∆) p with p ( x,
0) = δ ( x ) p ( x, t ) = 0 on ∂ Ω , (1)where V ( x ) is the underlying potential energy landscape and β is the inverse temper-ature. This relation can be rewritten in the eigen-basis of the generator ( {− λ k , φ k } )7nd of its adjoint ( { ψ k } ) as: p ( x, t ) = (cid:88) k exp( − λ k t ) φ k ( x ) (cid:90) Ω p ( y, ψ k ( y ) dy = (cid:88) k exp( − λ k t ) φ k ( x ) ψ k ( x ) (2)where the eigenvalues are ordered such that 0 < λ < λ ≤ λ , ... , where we notethat, for this dynamics, λ is guaranteed to be non-degenerate [13].Taking the limit t → ∞ , the last equation becomes:lim t →∞ p ( x, t ) = exp( − λ t ) φ ( x ) ψ ( x ) + O (exp( − λ t )) , or, in normalized form,lim t →∞ ˜ p ( x, t ) = p ( x, t ) (cid:82) Ω p ( x, t ) dx = φ ( x ) + O (exp( − ( λ − λ ) t )) . (3)Conditional on not having escaped, the normalized distribution function inside Ω tends to φ ( x ), the so-called quasi-stationary distribution (QSD) of Ω for all x , up toterms that decay exponentially with t . What are then the first escape statistics out ofΩ after having spend a time t → ∞ within it? The previous elementary derivationdirectly implies an important result: the escape flux out of a given infinitesimalsurface element centered on point x s on ∂ Ω — which is proportional to ∇ φ ( x s ) · n ( x s ), where n ( x s ) is the local normal to ∂ Ω — is also independent of x and of theescape time, since the probability density tends to a unique distribution φ ( x ) for allinitial conditions x . See Proposition 3 in Ref. 13 for a rigorous derivation.With this result, return to the issue at hand: consider a long continuous trajectory x ( t ) that enters set Ω i at t = 0 and then remains within it for a time t c → ∞ withoutescaping, at which time the discrete trajectory made a transition to state X n = i .What is then the probability that the next set in which the trajectory spends t c isΩ j (and hence that the discrete trajectory next makes a transition to X n +1 = j )?8rom the previous derivation, the distribution of first escape points on ∂ Ω i , andhence of the future of the trajectory after it escapes, is independent of x (0) and ofthe time at which the escape ultimately occurs. It follows that this distribution isalso independent of x ( t <
0) and hence of { X k We now investigate the corresponding jump dynamics in continuous time. Asshown in Ref. 13, strong statements can be made of the distribution of residencetime in set Ω i given that x ( t ) remained within it for t c → ∞ without escaping. Mostnotably, the distribution of exit times is exponential and escape points and times are9ncorrelated random variables (c.f. Proposition 3 in Ref. 13). However, the samecannot be said of the corresponding distribution of residence time in state i except forthe fact that it also cannot depend on x ( t < 0) and hence on { X k The previous result is strictly exact only in the limit where t c tends to infinity.In this limit, the jump process is also unfortunately non-informative, as the resi-dence time within each state diverges, hence providing very little information on10hereabouts of the system between jumps. Key to the practical relevance of thisrepresentation is the fact that convergence is exponentially fast as a function of t c . Itwas indeed shown in Prop. 6 of Ref. 13 that the difference between the joint distribu-tions of first escape times and first escape points conditioned on remaining in set Ω i for an infinite and a finite amount of time t c , respectively, decays as exp( − ( λ − λ ) t c ),as suggested by Eq. 3. Note that this holds for any functions of the escape times andpoints out of Ω i . Hence, local deviations between the MRP representation and theactual statistics of the discretized Langevin process decay exponentially with t c , ona state-per-state basis. Convergence is especially fast when the sets can be chosen tobe sufficiently metastable, i.e., to have a sufficiently large spectral gap λ − λ . Then, t c can be chosen so as to guarantee an accurate representation ( t c > / ( λ − λ )),while still preserving a good time resolution ( t c < /λ ). Note that it was recentlyshown that exponential convergence to the QSD also occurs for Langevin dynamicsat finite friction [19] and for a number of other random processes [27]. However, ex-ponential convergence is not guaranteed for arbitrary dynamics, so slow convergenceof the MRP representation with increasing t c is possible in some cases.An important advantage of the mapping to MRPs is that the accuracy of a jumpprocess based on sub-optimal set definitions can nonetheless be made arbitrarilyhigh by simply increasing t c . The MRP representation can therefore be expected tobe useful in a wide range of conditions, even with a rather naive definition of thesets { Ω } , which we demonstrate using numerical examples below. It should howeverbe kept in mind that if sufficiently metastable sets cannot be defined, the trade-offbetween resolution and accuracy could be unfavorable, resulting in either an accurateyet uninformative representation or an inaccurate but informative one.11 V. NUMERICAL RESULTS In order to demonstrate the accuracy of the MRP representation, we concentrateon P ij ( t ), the probability that the system is found in state j at time t given that itentered state i at time 0, i.e.,: P ij ( t ) = P r { X ( t ) = j | X (0) = i, T = 0 } , i, j ∈ X ; t ≥ . (4)For MRPs, P ij ( t ) is given by the solution of so-called renewal equations: P ij ( t ) = (1 − (cid:88) j ∈{ X } p ij F ij ( t )) δ ij + (cid:88) k ∈{ X } ( p ik F ik ∗ P kj )( t ) (5)where the second term on the R.H.S. is a convolution:( A ∗ B )( t ) = (cid:90) t A ( t − s ) dB ( s ) (6)so that Eq. 5 can be rewritten as: P ij ( t ) = (1 − H i ( t )) δ ij + (cid:88) k ∈{ X } (cid:90) t P kj ( t − s ) p ik F ik ( ds ); t ≥ i ∈ { X } (7)In the following, renewal equations were numerically integrated using the algo-rithm proposed in Ref. 28.We calculate the state-to-state transition probabilities in two well studied biomolec-ular systems: alanine dipeptide and chignolin. We generated 14 1- µ s long trajectoriesof alanine dipeptide where snapshots were saved every 2 ps. The technical details ofthese simulations are reported in the supplementary material. The 1- µ s long Chigno-lin trajectory (where snapshots were saved every 0.2 ns) was provided by D.E. ShawResearch; technical details of the simulations are provided in Ref. [29]. For bothsystems, we employ Markov State Models (MSMs) [6] and Perron Cluster-Clusteranalysis (PCCA) [30] to define metastable sets. Details on the simulations and12onstruction of the sets are provided in the supplementary material. Note howeverthat, for the present purpose, these details are not essential as we will show that ahigh accuracy MRP model can be obtained even for naive set definitions.With the set boundaries defined, we discretized the continuous trajectory follow-ing the procedure discussed in Section II A, namely each point x ( t ) was assigned adiscrete state X according to the set where the trajectory last remained for least atime t c without escaping. In order to assess the accuracy of the MRP, we comparethe transition probabilities directly computed from the discretized trajectory (thegold standard in this context), to those obtained by solving the MRP using renewalequations. These equations were parametrized by fitting to the observed transitiontime distributions p ij F ij ( t ) between every pair of states again using the discretizedtrajectory. Fig. 2 shows the two sets of transition probabilities for alanine dipeptidefor three different correlation times, t c = 2 ps, 20 ps and 40 ps, respectively. It isworth noting that the reference results also depend on t c , and thus that both thetarget and the approximation change with that parameter. It can be seen that evenfor the shortest correlation time (2 ps), the transition probabilities obtained by solv-ing the renewal equations are in very good agreement with the directly measuredprobabilities, but small deviations can be observed. As expected, the results becomeindistinguishable from the direct reference at the larger values of t c . In this case, theresults show that the dynamics can be very well approximated by a MRP even forshort correlation times, pointing to the fast convergence of the model with respectto t c .This convergence can be further demonstrated by the following experiment: in-stead of carefully defining sets so that they are as metastable as possible (and hencefor which a Markovian approximation might be quite accurate), we instead arbitrar-ily construct 4 states by partitioning the φ − ψ space into four rectangular cells using φ = 0 ◦ , 180 ◦ and ψ = 0 ◦ , 180 ◦ as dividing lines. We then discretize the trajectory13ccording to this new set definition and repeat the procedure described above. Fig.3 shows the results for three correlation times (2 ps, 20 ps and 40 ps). It can be seenthat the two sets of probabilities now deviate significantly at the shortest correlationtimes ( t c = 2 and 20 ps), indicating that the new states are not as metastable as theold ones, and hence require significantly longer correlation times. As the correlationtime is further increased to 40 ps however, the two probabilities nicely converge intoan essentially perfect agreement. This again support the above discussion: whilemodels built upon strongly metastable sets are preferable because t c can then betaken to be short compared to a typical first escape time from the set ( ∼ /λ ),ultimate convergence is guaranteed for any set definition by simply increasing t c .Fig. 4 shows the same quantities for two PCCA states of Chignolin for three dif-ferent correlation times ( t c = 0.2 ns, 2 ns and 10 ns). For the shortest correlationtime of 0.2 ns, the MRP completely fails to capture the underlying dynamics. Asthe correlation time is increased, the transition probabilities computed from Eq. 5converge to the reference values, again showing that the approximation of the under-lying dynamics can be systematically improved by increasing t c . These encouragingresults show promise that this formalism can be extended to more complex biologi-cal systems. Importantly, it can be valuable when atomistic information needs to beupscaled into coarser models [31]. V. CONCLUSION We have shown that Langevin dynamics in a high-dimensional configuration spacecan be mapped to a jump process over a discrete state space through a combinationof a novel mapping procedure, whereby the discrete state of the system correspondsto the last set in which the continuous trajectory spent a time t c without escaping,and of a MRP representation of the jump probabilities. This representation is shown14o become exact as t c → ∞ , and to converge exponentially fast to that limit, nomatter the number and size of the sets used in the discretization, in contrast to con-ventional Markovian representations. We therefore expect the MRP representationto accurately reproduce the dynamics of the system in a wide range of conditions, assupported by numerical examples. While this work demonstrates the formal powerof this class of models, important questions related to efficient set definition andparametrization of MRPs from molecular dynamics simulations remain to be devel-oped; these will be addressed in future publications. VI. ACKNOWLEDGEMENTS We thank D.E. Shaw Research for providing the chignolin trajectory. A.A, S.G,and A.F.V acknowledge support from the Joint Design of Advanced ComputingSolutions for Cancer (JDACS4C) program established by the U.S. Department ofEnergy (DOE) and the National Cancer Institute (NCI) of the National Institutes ofHealth. D.P. ackowledges support from Los Alamos National Laboratory’s (LANL)LDRD program under grant 20190034ER. LANL is operated by Triad National Se-curity, LLC, for the National Nuclear Security Administration of U.S. Departmentof Energy (Contract No. 89233218CNA000001). Computing resources were madeavailable by LANL’s Institutional Computing program. [1] A. E. Torda and W. F. van Gunsteren, Journal of computational chemistry , 1331(1994).[2] A. Ma and A. R. Dinner, The Journal of Physical Chemistry B , 6769 (2005). 3] P. Das, M. Moll, H. Stamati, L. E. Kavraki, and C. Clementi, Proceedings of theNational Academy of Sciences , 9885 (2006).[4] G. A. Tribello, M. Ceriotti, and M. Parrinello, Proceedings of the National Academyof Sciences , 5196 (2012).[5] B. Nadler, S. Lafon, R. R. Coifman, and I. G. Kevrekidis, Applied and ComputationalHarmonic Analysis , 113 (2006), special Issue: Diffusion Maps and Wavelets.[6] G. R. Bowman, V. S. Pande, and F. No´e, An introduction to Markov state models andtheir application to long timescale molecular simulation , Vol. 797 (Springer Science &Business Media, 2013).[7] F. Sittel, A. Jain, and G. Stock, The Journal of chemical physics , 07B605 1(2014).[8] A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, Journal of Computational Physics ,10 (1975).[9] D. T. Gillespie, Journal of computational physics , 403 (1976).[10] A. F. Voter, in Radiation effects in solids (Springer, 2007) pp. 1–23.[11] F. No´e, I. Horenko, C. Sch¨utte, and J. C. Smith, The Journal of chemical physics , 04B617 (2007).[12] J. D. Chodera, N. Singhal, V. S. Pande, K. A. Dill, and W. C. Swope, The Journalof chemical physics , 04B616 (2007).[13] C. Le Bris, T. Lelievre, M. Luskin, and D. Perez, Monte Carlo Methods and Appli-cations , 119 (2012).[14] T. Leli`evre, Handbook of Materials Modeling: Methods: Theory and Modeling , 1(2018).[15] G. Di Ges`u, T. Leli`evre, D. Le Peutrec, and B. Nectoux, Journal de Math´ematiquesPures et Appliqu´ees (2019).