Mapping the dynamics of complex multi-dimensional systems onto a discrete set of states conserving mean first passage times: a Projective Dynamics approach
MMapping the dynamics of complex multi-dimensional systemsonto a discrete set of states conserving mean first passage times:a Projective Dynamics approach
Katja Sch¨afer ∗ and M. A. Novotny Department of Physics and Astronomy,Mississippi State University, Mississippi State, MS 39762-5167
Abstract
We consider any dynamical system that starts from a given ensemble of configurations and evolvesin time until the system reaches a certain fixed stopping criterion, with the mean first-passage timethe quantity of interest. We present a general method, Projective Dynamics, which maps the multi-dimensional dynamics of the system onto an arbitrary discrete set of states { ζ k } , subject only to theconstraint that the dynamics is restricted to transitions not further than the neighboring states ζ k ± . We prove that with this imposed condition there exists a master equation with nearest-neighbor coupling with the same mean first-passage time as the original dynamical system. Weshow applications of the method for Brownian motion of particles in one and two dimensionalpotential energy landscapes and the folding process of small bio-polymers. We compare resultsfor the mean first passage time and the mean folding time obtained with the Projective Dynamicsmethod with those obtained by a direct measurement, and where possible with a semi-analyticalsolution. PACS numbers: 05.40.-a, 05.10.Gg, 05.10.-a, 02.50.Ga a r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r NTRODUCTION
The mean first-passage time (MFPT) is the time required for an ensemble of identicalsystems starting from a given initial distribution of configurations to reach a fixed criterion.Being able to obtain the MFPT is important in a wide array of physical applications. Theseinclude the MFPT when a star of a certain mass goes supernova, to neuron firing dynamics,spreading of diseases, chemical kinetics, folding processes of polymers, the death of livingcreatures, and the decay of metastable states via nucleation and growth processes [1, 2]. TheProjective Dynamics (PD) method to calculate MFPT for systems with a discrete state space[3–6] mapped the MFPT problem for metastable decay of magnetic states associated withIsing-type (discrete spin states) models onto a nearest-neighbor coupled master equation.However, the PD method has not previously been shown to be valid for other dynamicalsystems, for example for systems with continuous state spaces [7]. The calculation of theMFPT for systems with long time scales has previously required that the dynamic systemsatisfy certain conditions. For example, for the milestoning method [8] to work the systemmust be close to equilibrium or have separation of time scales or committer surfaces. Weprove in this letter that the PD method has no such constraints, and hence it is generallyapplicable to obtain the MFPT of any dynamical system.We prove that for any MFPT problem (with a finite MFPT) there exists a discrete masterequation with nearest-neighbor coupling that has the same MFPT as the original system.This master equation has time-independent rates, even if the underlying system is non-Markovian. Every point in phase space in the original system must be mapped onto a statein the discrete master equation. We show that any mapping that leads to only nearest-neighbor coupling in the master equation preserves the MFPT provided the growing ( g k )and shrinking ( s k ) rates between all master equation states k are correct. Therefore, whetherthe g k and s k are known theoretically or are measured, the correct MFPT is obtained. Formetastable decay in the ferromagnetic d =3 Ising case such knowledge enabled calculationsof MFPT of about 10 Monte Carlo steps per spin [5]. It is anticipated that the knowledgethat the PD method gives the correct MFPT will enable long-time-scale simulations forother physical models. An example of where the PD method may be useful for long-timesimulations is studies of the intercalation time of Li ions in molecular dynamics studiesmodeling charging of batteries [9]. 2
HEORY
Consider any ensemble of systems evolving in phase space Γ under some dynamic. Theystart with some given distribution, and we desire to calculate the MFPT until each systemstops after it meets some fixed stopping criterion. The dynamics of the ensemble can bemapped onto a set of discrete non-overlapping states { ζ k } , with 0 ≤ k ≤ S , which cover theentire domain. We impose the restriction that the states ζ k can make transitions only to thestates ζ k ± . This condition can be met for many systems by adjusting the ‘length’ ∆ ζ of thestates, so that no transitions within a single time step d t further than to the neighboringstates exist. The states that form the stopping citerion will all be mapped to ζ .The time-evolution of the system can be described by a master equation with nearestneighbor coupling, which describes the change of occupation P k ( t ) of state k by the proba-bility flows between adjacent states. Namely,dd t P k ( t ) = g k +1 P k +1 ( t ) − [ g k + s k ] P k ( t ) + s k − P k − ( t ), (1)where the growing g k (shrinking s k ) rates represent the directional rates at which the systemtransits from state k to state k − k + 1). The initial conditions { P k (0) } are specified bythe ensemble of initial conditions of the original system.Let N ( k → k ± t, t + d t ) denote the number of times the system makes a transitionfrom state k to its neighboring states ( k ±
1) in the time interval ( t, t + d t ] and N k ( t ) thenumber of times the system resides in state k at time t . Then the growing rate for the timeinterval ( t o , t e ] can be obtained by measuring the system at time intervals d t . Explicitly g k d t = (cid:82) t e t o d tN ( k → k − t, t + d t ) (cid:82) t e t o d tN k ( t ) . (2)A similar expression exists for the shrinking rate s k .Assuming that the ensemble of systems has been fully absorbed at t e = ∞ after startingat t o = 0, we write in shorthand the master equationdd t P k ( t ) = (cid:88) k (cid:48) W ( k, k (cid:48) ) P k (cid:48) ( t ) , (3)where the elements of W ( k, k (cid:48) ) are the growing and shrinking rates, which have been mea-sured from t o = 0 up to some point in time when all systems in the ensemble are fully3bsorbed. By using the fundamental relation [10] (cid:80) k τ ( k ) W ( k, k (cid:48) ) = − , the general ex-pression for the MFPT τ ( S o ) for a system starting fully in state S o can be obtained as τ ( S o ) = S o (cid:88) i =1 g i + S − (cid:88) (cid:96) =1 min( S o , | S − (cid:96) | ) (cid:88) i =1 (cid:32) (cid:81) i + (cid:96) − j = i s j (cid:81) i + (cid:96)j = i g j (cid:33) (4)where S o denotes the starting state and S the total number of states within the domain.We prove that Eq. (4) gives the same MFPT independent of the particular choice ofthe states. For this purpose, let m ≡ k (cid:83) k + 1 be the joined interval of state k and k + 1, then for a single time step it follows that the transitions from state m to k − k + 2) in the new system are equal to the transitions from state k to k − k + 1 to k + 2) in the old system namely ˜ N ( m → k − , t, t + d t ) = N ( k → k − , t, t + d t ) and˜ N ( m → k + 2 , t, t + d t ) = N ( k + 1 → k + 2 , t, t + d t ), whereas the probability of occupationof state m rescales as ˜ P m ( t ) = P k ( t ) + P k +1 ( t ) = [ N k ( t ) + N k +1 ( t )] / (cid:80) k N k ( t ). Thus thegrowing rate at time t in the old and new system are related as˜ g m = P k ( t ) P k ( t ) + P k +1 ( t ) g k (5)and similarly for the new shrinking rate ˜ s m .That this procedure leaves the MFPT unchanged can be understood by noting the effectof joining two states on the master equation. The joining of adjacent states corresponds tothe addition of terms of the master equation of the same adjacent states, namelydd t P k ( t ) = g k +1 P k +1 ( t ) − [ g k + s k ] P k ( t ) + s k − P k − ( t ). (6)and dd t P k +1 ( t ) = g k +2 P k +2 ( t ) − [ g k +1 + s k +1 ] P k +1 ( t ) + s k P k ( t ). (7)It is easily seen that the cross terms vanish and we are left with an expression of the formdd t ˜ P m ( t ) = g k +2 P k +2 ( t ) − g k P k ( t ) − s k +1 P k +1 ( t ) + s k − P k − ( t ). (8)Using (5) and similarly for ˜ s m in (8) for a single time step, leads todd t ˜ P m ( t ) = g k +2 P k +2 ( t ) − [˜ g m + ˜ s m ] ˜ P m ( t ) + s k − P k − ( t ). (9)Hence the general form [compare Eq. (1) and (6)] is restored in the reduced system ofequations. Performing this joining repeatedly, finally only one state is left. This gives thesame MFPT as the original ensemble of systems [4].4 IG. 1: (color online) Kramer’s and a rough potential. The states are chosen as intervals alongthe x -coordinate of some fixed length ∆ l . COMPUTATIONAL STUDY/SIMULATION
To illustrate the theoretical framework, we demonstrate diffusion in both dimensions d =1and d =2, as well as the folding process of a model for a small linear polymer chain. Wechoose these examples to specifically address the independence of the method from anyexistence or knowledge of a free energy landscape or most probable escape pathway.As a d =1 historical example, we studied the diffusion process of particles subjected to aKramers or to a rough [11] potential (Fig. 1). The trajectories were integrated along 5120individual paths starting at x = − . −∞ ,
0] atthe absorbing boundary x = 0. The MFPTs were calculated using Eq. (4). In order to showthe freedom in choosing the length of the intervals, we show examples of bins with different5engths ∆ l . Starting at the absorbing boundary k =0 for x ≥
0, the system resides in state kwhen ( k − l ≤ x < k ∆ l . In both cases (smooth and rough), we chose ∆ l = 0 . , . . l = 0 . l = 0 . l = 0 .
4) of the starting state.The results do not change with any of these starting positions, which further indicates theindependence of the PD method.As a d =2 example we applied the PD method to the diffusion process at the entropicbarrier U entr ( x, y ) of Ref. [8]. We chose the coordinate ζ in two different ways. The trajec-tories are started at ( x o , y o )=( − . ,
0) and terminated once they reach x>
0. We obtainedresults for defining the states both with the coordinate ζ = x and ζ = r [Fig. (2)]. The firstrepresents a binning along the x-coordinate, interpreted as the progress towards the absorb-ing boundary, while the circular binning has no such interpretation. For ζ = x the states were( k − x ≤ x 065 and 0 . 13. For ζ = r , the states are defined as( S − k )∆ r ≤ r< ( S − k + 1)∆ r for 2 ≤ k ≤ S , where r denotes the distance from ( x o , y o ). Wechose ∆ r =0 . 025 and 0 . 05. For the circular binning, the bin k =1 is comprised of all pointswith x ≤ HP ) H , where H standsfor the hydrophobic and P for the polar monomers. The structure of the individual chainis maintained by three interatomic potentials: a bond-stretching potential, a bond-anglebending potential, and a rejecting potential avoiding an overlap of HP or P P monomers.The chain undergoes a folding process from the fully elongated state to a compact near nativestate caused by attractive Lennard-Jones type of interactions between HH monomers. ThePD method is given by choosing the coordinate ζ as the potential energy of the nativecontacts (the HH interaction energy). Although this choice seems apparent, it does notcorrespond to an actual dynamical coordinate. Starting with the fully elongated chain, for6 B T Projective Dynamics Brownian Riemann∆ l = 0 . l = 0 . l = 0 . . a . ± . 2) 18 . ± . 2) 18 . ± . 2) 18 . ± . 2) 18 . . a . ± . 3) 24 . ± . 3) 24 . ± . 3) 24 . ± . 3) 24 . . a . ± . 4) 36 . ± . 4) 36 . ± . 4) 36 . ± . 5) 36 . . a ± 1) 62( ± 1) 62( ± 1) 62( ± 1) 61 . . b . ± . 6) 46 . ± . 6) 46 . ± . 6) 46 . ± . 6) 46.20 . b ± 1) 80( ± 1) 80( ± 1) 80( ± 1) 810 . b ± 1) 179( ± 1) 179( ± 1) 178( ± 2) 178∆ r = 0 . 025 ∆ r = 0 . 05 ∆ x = 0 . 065 ∆ x = 0 . . c . ± . 1) 15 . ± . 1) 15 . ± . 1) 15 . ± . 1) 15 . ± . . c . ± . 2) 17 . ± . 2) 17 . ± . 2) 17 . ± . 2) 17 . ± . . c . ± . 3) 20 . ± . 3) 20 . ± . 3) 20 . ± . 3) 20 . ± . . c . ± . 1) 24 . ± . 1) 24 . ± . 1) 24 . ± . 1) 24 . ± . . c . ± . 5) 30 . ± . 5) 30 . ± . 5) 30 . ± . 5) 30 . ± . U/C h ≈ U/C h ≈ . d . ± . 03) 9 . ± . 03) 9 . ± . . d . ± . 03) 9 . ± . 03) 9 . ± . . d . ± . 03) 10 . ± . 03) 10 . ± . . d . ± . 04) 12 . ± . 04) 12 . ± . . d . ± . 05) 15 . ± . 05) 15 . ± . TABLE I: MFPTs for (a) Kramer’s potential, (b) rough 1 dim. (c) 2-dim. entropic barrier and (d)the folding of a polymer chain. The results of the Projective Dynamics method are given for theircorresponding intervals ∆ ζ = ∆ l Kramer’s, ∆ r and ∆ x U/C h polymer. Also shownresults obtained by direct measurement of the Brownian motion process and the semianalyticalsolution for the d =1 absorption processes. different temperatures of the external bath, we simulated the folding process due to theBrownian motion. We compared the MFPTs of the chain of reaching a potential energydeep enough in the basin of attraction of the energy ground state (i.e. the native state) sothat no escape (i.e. unfolding) from there is likely to occur. The results of the PD method7 IG. 2: (color online) The entropic barrier, the states are defined along the coordinate ζ = x (representing a progress toward the absorbing boundary) and ζ = r (representing circular intervals)of some fixed interval length ∆ ζ = ∆ x and ∆ r . as well as a direct measurement are shown in Table I. For all simulated temperatures thePD method gives the same results (within statistics) as a direct measurement.8 ONCLUSION AND DISCUSSION The Projective Dynamics (PD) method, as a tool for projecting multi-dimensional sys-tems onto one set of states { ζ k } and mapping the time-evolution onto a discrete masterequation with nearest neighbor coupling, was shown to be independent of the choice in thestates ζ k and the ‘length’ of the intervals which define the states. Provided only transi-tions between adjacent states occur, the PD method correctly obtains MFPTs. This resultmakes the PD method generally applicable to the dynamics of any system (discrete or con-tinuous, Markovian or non-Markovian, stochastic or deterministic). Subject only to thenearest-neighbor coupling constraint, the states { ζ k } can be chosen in multiple ways, with-out changing the mean first-passage time (MFPT). In many cases of interest, the averageoccupation probability for every state ζ k is also the same as for the original dynamical sys-tem. The calculation of the growing and shrinking rates required in the PD method mayfor some models and dynamics be obtained by using other methods, such as calculations ofhistories [12] or forward flux sampling [13] or accelerated molecular dynamics [14, 15].Useful discussions with Jerzy Blawzdziewicz, Weinan E, Ron Elber, Miroslav Kolesik,Cory O’Hern, Per Arne Rikvold, and Eric Vanden-Eijnden are gratefully acknowledged,as well as useful interactions at a 2008 summer school at the Aspen Center for Physics.This research was supported in part by the National Science Foundation through TeraGridresources provided by the Texas Advanced Computing Center (TACC) under grant numberTG-DMR090092. ∗ Electronic address: [email protected][1] S. Redner, A Guide to First-Passage Processes (Cambridge Univ. Press, Cambridge, UK,2001); errata http://physics.bu.edu/ ∼ redner.[2] S. Condamin, O. B´enichou, V. Tejedor, R. Voituriez, and J. Klafter, Nature , 77 (2007).[3] M. Kolesik, M.A. Novotny, and P.A. Rikvold, Phys. Rev. Lett. , 3384 (1998).[4] M.A. Novotny, M. Kolesik, and P.A. Rikvold, Comput. Phys. Commun. , 330 (1999)[5] M. Kolesik, M.A. Novotny, and P.A. Rikvold, Inter. J. Mod. Phys. C , 121 (2003).[6] M.A. Novotny, ‘A Tutorial on Advanced Dynamic Monte Carlo Methods for Systems withDiscrete State Spaces’, in Annual Reviews of Computational Physics IX , ed. D. Stauffer (World cientific, Singapore 2001), pp. 153–210.[7] K. Sch¨afer and M.A. Novotny, Computer Simulation Studies in Condensed Matter PhysicsXXI, eds. D.P. Landau, S.P. Lewis, and H.B. Sch¨uttler (Springer Verlag, Heidelberg, Berlin,2008).[8] A.K. Faradjian and R. Elber, J. Chem. Phys. , 10880 (2004).[9] I. Abou Hamad, M.A. Novotny, D.O. Wipf, and P.A. Rikvold, Phys. Chem. Chem. Phys. ,2740 (2010).[10] R. Zwanzig, A. Szabo, and B. Bagchi, Proc. Natl. Acad. Sci. , 20 (1992).[11] R. Zwanzig, Proc. Nat. Acad. Science , 2029 (1988).[12] N. Gulbahce, F.J. Alexander, and G. Johnson, Phys. Rev. E 321 (2002).321 (2002).