Estimating dissipation from single stationary trajectories
EEstimating dissipation from single stationary trajectories. ´Edgar Rold´an and Juan M.R. Parrondo
Dep. F´ısica At´omica, Molecular y Nuclear and GISC. Universidad Complutense de Madrid. 28040-Madrid, Spain (Dated: November 8, 2018)In this Letter we show that the time reversal asymmetry of a stationary time series providesinformation about the entropy production of the physical mechanism generating the series, even ifone ignores any detail of that mechanism. We develop estimators for the entropy production whichcan detect non-equilibrium processes even when there are no measurable flows in the time series.
PACS numbers: 05.70.Ln, 05.20.-y, 05.40.-a
The relationship between irreversibility and entropyproduction forms the core of thermodynamics and sta-tistical mechanics. However, it had not been formu-lated quantitatively until the recent introduction of theKullback-Leibler distance or relative entropy in the con-text of fluctuation and work theorems [1]. The relativeentropy between two probability distributions, p ( x ) and q ( x ) is defined as D ( p || q ) ≡ (cid:88) x p ( x ) log p ( x ) q ( x ) , (1)and is a measure of their distinguishability [2]. The aver-age entropy production associated with a process drivenby an external agent turns to be equal to the relativeentropy between the two probability distributions de-scribing the process running forward and backward intime [1, 3–6]. This relative entropy can be thought of asthe distinguishability between the process and its timereverse, i.e., as the irreversibility exhibited by the pro-cess. The relationship between entropy production andrelative entropy has been derived in different scenarios:Hamiltonian dynamics [1, 3] and Langevin dynamics [5],and has also been tested in experimental situations [5].When applied to non-equilibrium stationary states(NESS), the entropy production per unit time reads (cid:104) ˙ S (cid:105) k = lim t →∞ t D (cid:104) p (cid:16) { x ( τ ) } tτ =0 (cid:17)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) p (cid:16) { x ( t − τ ) } tτ =0 (cid:17)(cid:105) (2)where k is the Boltzmann constant and p (cid:16) { x ( τ ) } tτ =0 (cid:17) isthe probability of observing a given trajectory { x ( τ ) } tτ =0 in phase space. Since we focus on stationary trajecto-ries —where the external forcing, if any, is constant—,there is no need of reversing the driving in the backwardprocess. Moreover, a sufficiently long single trajectorycan provide all the necessary statistics to compute therelative entropy in Eq. (2) and consequently the entropyproduction rate.Fortunately, the full information of the trajectory inthe phase space is not always necessary. Eq. (2) followsimmediately from the Gallavotti-Cohen theorem [7], byreplacing the relative entropy between trajectories with D ( p S ( s ) || p S ( − s )), where p S ( s ) is the probability to ob-serve an entropy production s in a time interval [0 , t ]. In general, the relative entropy calculated using partial in-formation, { ˜ x ( τ ) } tτ =0 where ˜ x ( τ ) is a non-invertible func-tion of x ( τ ), only provides a lower bound on the averageentropy production [1, 6, 8]. For stationary trajectories,instead of Eq. (2) one obtains a lower bound, which is metif ˜ x ( τ ) univocally determines the entropy production s .For discrete stationary trajectories x , . . . , x n , we candefine the relative entropy of n -strings as D n ( p F || p B ) ≡ (cid:88) x ,...,x n p ( x , . . . , x n ) log p ( x , . . . , x n ) p ( x n , . . . , x )(3)Following the above arguments, we arrive at: (cid:104) ˙ S (cid:105) k ≥ d ( p F || p B ) ≡ lim n →∞ n D n ( p F || p B ) . (4)This equation reveals a striking connection betweenphysics and the statistics of a time series. The l.h.s.is a purely physical quantity (it is proportional to theaverage dissipated energy per step), whereas the r.h.sis a statistical magnitude depending solely on the data x , x , . . . , but not on the physical mechanism generat-ing those data. Such a connection is a generalizationof the Landauer’s principle relating entropy productionand logical irreversibility [1, 9, 10]. Eq. (4) extends thisprinciple and suggests that we can determine the entropyproduction of an arbitrary NESS by computing the rel-ative entropy of forward and backward trajectories. Wecould, for instance, determine whether a biological pro-cess is active or passive or even estimate, or bound, theamount of consumed ATP by measuring the relative en-tropy of data generated in the process.In this Letter we explore the feasibility of such a tech-nique by analyzing the validity of Eq. (4) and developingestimators of the relative entropy. Our approach is gen-eral, but we use a discrete flashing ratchet as a case study,wherein direct comparison between analytical and empir-ical values of the relative entropy and the entropy pro-duction is possible. There have been previous attemptsto distinguish between equilibrium and NESS. Martin etal checked the fluctuation dissipation relationship in ex-perimental data from hair bundles of hair cells [11], butthis approach needs two types of data: spontaneous andforced fluctuations. Amman et al analyzed the possibility a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p to discriminate between equilibrium and non-equilibriumin a three state chemical system [12]. Finally, Kennel in-troduced in [13] criteria based on compression algorithmsto distinguish between symmetric and asymmetric timeseries in the context of chaotic signals, without any con-nection to dissipation. As we show in this Letter, relativeentropy provides a more general and simpler frameworkfor the problem of distinguishing between equilibriumand NESS and, moreover, yields estimations and lowerbounds on the entropy production.Two strategies have been considered to estimate therelative entropy between stochastic processes: the firstis based on brute-force counting of n -strings, obtainingempirical estimates of p ( x , . . . , x n ), and computing D n using Eq. (3); the second is based on string parsing, thebasic procedure of the Lempel-Ziv compression algorithm[14].The first strategy is simpler and more effective forMarkov chains. Our results indicate that this is still thecase for some non-Markov process [15]. Consequently,we will restrict ourselves in this Letter to estimations ofrelative entropy from empirical probability distributions.If the process and its reverse are Markovian, p ( x , x , . . . , x n ) = p ( x ) p ( x | x ) ...p ( x n − | x n ), the rel-ative entropy rate d defined in Eq. (4) can be expressedin terms of the relative entropy between distributions ofsubstrings of size 2: d ( p F || p B ) = (cid:88) x ,x p ( x , x ) log p ( x | x ) p ( x | x ) = D − D . (5)In the specific case of a trajectory and its reverse, the one-time statistics are identical and D ( p F || p B ) = 0. Thenfor Markovian dynamics d ( p F || p B ) = D , which can becalculated by frequency counting if the number of statesand possible transitions is not large. In general, if onedefines d k ≡ D k − D k − (6)then d k → d for k → ∞ . The limit is reached for finite k for the so-called k -th order Markov chains, i.e. whenblocks of size k , X k ≡ ( x n , . . . , x n + k − ), are Markovian[16]. In this case d ( p F || p B ) = d k +1 = d k +2 = . . . . Formore general processes, we will use the following ansatz,proposed by in Ref. [17] for Shannon entropy estimation: d k = d ∞ − c log kk γ , (7)where c and γ are parameters that, together with d ∞ ,can be obtained by fitting the empirical values of d k vs. k . We have tested the accuracy of these estimators and ofthe bound (4) in a specific example: a discrete flashingratchet [18], consisting of a particle moving in a one di-mensional lattice. The particle is at temperature T and Vr2 2 0 12V0 12’ 0’ 1’ 2’ 0’ 1’
FIG. 1: Discrete ratchet scheme. Particles can jump betweenthe states i → j , i (cid:48) → j (cid:48) , and i → i (cid:48) in a flashing asymmet-ric potential of height 2 V with periodic boundary conditions.The switching rate of the potential is r . V/kT − . . . . FIG. 2: Average dissipation per step (in units of kT ) in theflashing ratchet ( r = 1) and different estimations of relativeentropy using a trajectory with n = 10 steps and full in-formation, as a function of V /kT : analytical calculation ofthe average dissipation (black line), d (blue circles), d (redsquares). moves in a periodic and asymmetric potential of height2 V , which is switched on and off at a rate r (see Fig. 1).Trajectories are described by two variables: the positionof the particle, x = { , , } , and the state of the poten-tial (on or off), y = { , } .To define the dynamics of the particle, we start witha continuous time description based on rates of spatialjumps and switching. We assume that the motion in eachpotential obeys detailed balance: k i → j = e − β ( Vj − Vi )2 , and k i (cid:48) → j (cid:48) = 1 for i, j = 0 , , i (cid:54) = j . The system isdriven out of equilibrium by imposing constant switchingrates k i → i (cid:48) = k i (cid:48) → i = r , i = 0 , ,
2, which do not obeydetailed balance.We will focus on the dissipation per step : from thecontinuous trajectory ( x ( t ) , y ( t )) we generate a series( x n , y n ) comprising the states visited by the system.That is, we drop the information of the times when jumpsor switches occur. ( x n , y n ) is a Markov chain with transi-tion probabilities given by p α → γ = k α → γ / (cid:80) γ k α → γ , with α, γ = 0 , , , (cid:48) , (cid:48) , (cid:48) . Introducing these probabilities inEq. (5), d ( p F || p B ) = β (cid:80) (cid:104) V α − V γ (cid:105) , where the sum runsover transitions mediated by the thermal bath, i → j , i (cid:48) → j (cid:48) . The relative entropy turns out to be the aver-age dissipation per step in units of kT and we recoverthe main result, Eq. (2) [22]. It is also interesting to ex-plore the relationship between d and the stationary flows . . V/kT . . . . .
25 10 − − − − − − FIG. 3: Average dissipation per step (in units of kT ) in theflashing ratchet ( r = 1) and different estimations of relativeentropy using a trajectory with n = 10 steps and partial in-formation (position) as a function of V /kT : analytical calcu-lation of the average dissipation (black line), d (blue circles), d (green diamonds), d ∞ in Eq. (7) (orange circles with errorbars), and Monte-Carlo semi-analytical calculation of d (pur-ple crosses). Inset . Estimators for weak potentials in a log-logplot. We have added in the inset the analytical calculation of d (blue solid line). J αγ = p αγ − p γα between states α, γ = 0 , , , (cid:48) , (cid:48) , (cid:48) . If J αγ (cid:28) p αγ , we have: d (cid:39) (cid:88) αγ ( J αγ ) p αγ = (cid:88) α<γ ( J αγ ) p αγ . (8)which is a well known expression of the entropy produc-tion in continuous Markov systems [19], where d = d .Fig. 2 shows the dissipation, calculated analyticallyby solving the six-state Markov chain in the station-ary regime, and the estimations discussed above. Dueto Markovianity, relative entropies, d k , immediately con-verge d = d = d = . . . and d is equal to the entropyproduction per step. As long as one has a good esti-mation of p ( x , . . . , x k ), our approach provides accuratevalues of the entropy production, which is the case forweak potentials V (cid:39) kT . If V (cid:29) kT , then uphill jumps,0 →
1, 0 →
2, and 1 →
2, are so unlikely that they do notoccur in a finite trajectory. The higher order the statis-tics, the earlier this problem arises, as shown in Fig. 2.The reason is that d involves probability distributions ofthree-step trajectories, the sampling space is bigger and itis easier that some transitions i → j → k do not appearwhile their reverse do. Although these jumps are veryunlikely, they contribute significantly to d , as shown inFig. 2, where d and d have been calculated by restrict-ing the sum in D k to strings satisfying p ( x . . . x k ) (cid:54) = 0and p ( x k . . . x ) (cid:54) = 0.In real applications, it is more likely that one has onlypartial information of the trajectories. To study the ac-curacy of our estimators and of the inequality (4) in thiscase, we remove the information of the state of the po- .
00 0 .
01 0 .
02 0 .
03 0 . FL/kT − − − − − − − − − − FIG. 4: Average dissipation per step (in units of kT ) in theflashing ratchet ( r = 2 , V = 2 kT ) with external force F and different estimations of relative entropy using a trajec-tory with n = 10 steps with partial information (position):analytical calculation of the average dissipation (black line), d (blue circles, analytical values in blue dashed line), d (redsquares), d (green diamonds), and semi-analytical calcula-tion of d (purple crosses). The minimum in d corresponds tothe stall force. tential and consider trajectories described only by theposition { x k } nk =1 , which in general are not Markovian.As a consequence, the estimation of the relative entropy d ( p F || p B ) is more difficult, but even a good estimation of d only provides a lower bound on the relative entropy. Itis known that the Gallavotti-Cohen symmetry does nothold in the continuous flashing ratchet if the state of thepotential is not considered [20]. In fact, the bound (4)can be quite loose. For instance, if r → ∞ , switching isvery fast and the particle moves in an effective potential(the average of on and off) which is periodic. The po-sition x k becomes Markovian and the current vanishes.Using Eq. (8) one arrives at d = d = 0, whereas thedissipation per step is non-zero.In most cases however the bound given by Eq. (4) pro-vides significant information. In Fig. 3 we show the esti-mation of d using the empirical values of d k for k = 2 , d ∞ resulting from the fit of theansatz in Eq. (7). The error bars in Fig. 3 correspondto the error in the fit with a confidence interval of 90%.Our estimations clearly distinguish between the equilib-rium case ( V = 0) and the NESS. The empirical d k with k > d , or the bound provided by Eq. (4) is nottight. To clarify this question we need an analytical cal-culation of the relative entropy between two non-Markovprocesses. In our case, the relative entropy D n reads: D n = (cid:42) log (cid:80) y ,...,y n p ( x , y ; . . . ; x n , y n ) (cid:80) y ,...,y n p ( x n , y n ; . . . ; x , y ) (cid:43) (9)where the average is taken over all possible trajecto-ries. The probability distribution p ( x , y ; . . . ; x n , y n ) = p ( x , y ) × p ( x , y | x , y ) × . . . × p ( x n , y n | x n − , y n − ) isknown, but Eq. (9) cannot be calculated exactly. For-tunately, the log in Eq. (9) is a self-averaging quantityfor large n [21] and we can compute the average using asingle long typical trajectory [15]. We show in Fig. 3 thevalue of d obtained by this Monte-Carlo semi-analyticalcalculation (purple crosses), which is very close to theestimation d ∞ based on the ansatz Eq. (7).Although the relative entropy d underestimates the ac-tual dissipation, it does reproduce its asymptotic behav-ior. Entropy production decreases as V for small V , sodo d ∞ and d (see inset of Fig. 3). On the other hand, d ∝ V , since the current is J ∝ V (see Eq. (8)).We have found in several instances a similar qualitativeimprovement on the estimation of relative entropy whenusing blocks of size bigger than two. In particular, d and above outperform d , which, as indicated by Eq. (8),is equivalent to the standard calculation of entropy pro-duction using the currents observable from the availabledata; in our case, the spatial current. For a striking il-lustration of this effect we add an external force F to theflashing ratchet and study dissipation and relative en-tropy close to the stalling force F stall , for which the spa-tial current and d both vanish. Jumping rates are nowbiased in the direction of the force, giving the followingdetailed balance condition k i → j /k j → i = e − β ( V j − V i − F L ij ) , L ij = 1 being the distance between i and j .We have plotted in Fig. 4 the real dissipation, the an-alytical value of d and d and the empirical values of d , d , and d , close to the stalling force F stall . Recall that,for F = F stall , the position of the particle does not ex-hibit any flow and its average position remains constant.Consequently, d or any other estimation of entropy pro-duction based on flows will fail. However, the relativeentropy calculated using blocks of size 3 captures thenon-equilibrium nature of the time series.In conclusion, we have shown that the statistical prop-erties of a time series impose a lower bound on theentropy produced in generating the series. This lowerbound is valid even if we do not have any access or infor-mation of the physical mechanism generating the data.Finally, we have shown that the bound can be non-trivial,predicting dissipation even when the data do not exhibitany measurable flow. Our techniques could be applied todata from different sources. In the case of biological sys-tems, they could help to distinguish between passive andactive processes, and even to estimate ATP consumption.On the other side, as in the case of Landauer’s principle,relative entropy can be used to ascertain the minimalentropy production associated with a specific behavior,such as spatiotemporal patterns, excitable systems, etc.This in turn may influence the design of optimal deviceswith functionalities given by these behaviors. We acknowledge financial support from Grant MO-SAICO (Spanish Government), Becas de la Caixa paraestudios de M´aster en Espa˜na , and from the Max-Planck-Institut f¨ur Physik komplexer Systeme (Dresden, Ger-many). The main idea of the paper, namely, estimat-ing dissipation by relative entropy in stationary trajecto-ries, was suggested to us by Frank J¨ulicher and BenjaminLindner. We also acknowledge fruitful discussions withLuis Din´ıs and Abigail Klopper. [1] R. Kawai, J. M. R. Parrondo, and C. V. den Broeck,Phys. Rev. Lett. , 080602 (2007).[2] T. M. Cover and J. A. Thomas, Elements of informationtheory (Wiley, Hoboken, New Jersey, 2006), 2nd ed.[3] J. M. R. Parrondo, C. V. den Broeck, and R. Kawai, NewJ. Phys. , 073008 (2009).[4] U. Seifert, Phys. Rev. Lett. , 040602 (2005).[5] D. Andrieux, P. Gaspard, S. Ciliberto, N. Garnier,S. Joubaud, and A. Petrosyan, Phys. Rev. Lett. ,150601 (2007).[6] J. Horowitz and C. Jarzynski, Phys. Rev. E , 021106(2009).[7] E. Cohen and G. Gallavotti, J. Stat. Phys. , 1343(1999).[8] A. Gomez-Marin, J. M. R. Parrondo, and C. V. denBroeck, EPL , 50002 (2008).[9] R. Landauer, IBM J. Res. Dev. , 261 (2000).[10] D. Andrieux and P. Gaspard, P. Natl. Acad. Sci. USA , 9516 (2008).[11] P. Martin, A. J. Hudspeth, and F. J¨ulicher, P. Natl.Acad. Sci. USA , 14380 (2001).[12] C. P. Amann, T. Schmiedl, and U. Seifert, J. Chem. Phys. , 041102 (2010).[13] M. B. Kennel, Phys. Rev. E , 056208 (2004).[14] J. Ziv and N. Merhav, IEEE T. Inform. Theory (1993).[15] E. Rold´an and J. M. R. Parrondo, in preparation.[16] Z. Rached, F. Alajaji, and L. L. Campbell, IEEE T. In-form. Theory , 917 (2004).[17] T. Schurmann and P. Grassberger, Chaos , 414 (1996).[18] A. Ajdari and J. Prost, C.R. Acad. Sci. Paris II ,1635 (1992).[19] J. M. R. Parrondo and B. J. de Cisneros, App. Phys. A , 179 (2002).[20] D. Lacoste and K. Mallick, Phys. Rev. E , 021923(2009).[21] P. Jacquet, G. Seroussi, and W. Szpankowski, Theor.Comput. Sci. , 203 (2008).[22] Each transition α → γ obeying detailed balance con-tributes to d as β (cid:104) V α − V γ (cid:105) which is the average entropyincrease in the bath due to the transition. This still ap-plies for systems in contact with several baths at differ-ent temperatures. In our case, the constant flashing rate r can be interpreted as a transition mediated by a bathat infinite temperature ββ