Utilizing time-series measurements for entropy production estimation in partially observed systems

Estimating the dissipation, or the entropy production rate (EPR), can provide insights into the underlying mechanisms of nonequilibrium driven processes. Experimentally, however, only partial information can be accessed, and the ability to estimate the EPR varies depending on the available data. Here, we test different degrees of observed information stemming from coarse-grained time-series trajectory data, and apply several EPR estimators. Given increasing amount of information, we show a hierarchy of lower bounds on the total EPR. Further, we present a novel approach for utilizing waiting times in hidden states to provide a tighter lower bound on the total EPR.

Estimating the dissipation, or the entropy production rate (EPR), can provide insights into the underlying mechanisms of nonequilibrium driven processes.Experimentally, however, only partial information can be accessed, and the ability to estimate the EPR varies depending on the available data.Here, we test different degrees of observed information stemming from coarse-grained timeseries trajectory data, and apply several EPR estimators.Given increasing amount of information, we show a hierarchy of lower bounds on the total EPR.Further, we present a novel approach for utilizing waiting times in hidden states to provide a tighter lower bound on the total EPR.
The entropy production, or energy dissipation, is a fundamental physical quantity necessary to characterize the thermodynamics of nonequilibrium processes [1,2].In living systems, for example, the dissipation rate is closely related to the consumption rate of chemical fuel molecules, such as Adenosine triphosphate (ATP), by molecular motors [3].The entropy production calculated along a single trajectory is a stochastic quantity that follows a set of mathematical relations, collectively known as the fluctuation theorems [4][5][6][7][8][9][10], which have been validated experimentally [11][12][13].
Estimating the total EPR from experimental data for a driven system is not always trivial, and can be challenging owing to the limited resolution and the huge number of degrees of freedom [2].While calculating the EPR is straightforward given complete information about the nonequilibrium degrees of freedom, practically, only partial information is available.Such coarse-grained observation, where only some of the degrees of freedom are monitored or resolved, often cannot be treated with the Markovian approximation [14][15][16][17], and can only provide a lower bound on the total dissipation [18][19][20][21].
Partial information can stem from an observed subsystem, such that the rest of the system is hidden, or from coarse-graining some of the microstates into several meso-states [21][22][23][24][25].The observed information may also be the transitions between states rather than the states themselves [26][27][28].Therefore, different levels of coarse-graining can be considered based on the partial information available about the system.
There are several estimators for partial EPR that do not require any prior information about the system, such as the number of states or the underlying topology.The thermodynamic uncertainty relations (TUR), for example, provide a lower bound on the entropy production from the fluctuations of the transition fluxes or first passage times [29][30][31][32][33][34][35][36][37][38][39][40].
Another approach for inferring a lower bound on the total EPR is based on an optimization problem, search-ing over systems with known topology and the same observed statistics, preserving the first and second-order mass transfer rates [41], or the waiting time statistics [42], or both [43].A similar approach was also demonstrated for discrete-time Markov chains, searching over the possible underlying hidden Markov models given the number of hidden states [44].
Building on the deep connection between the dissipation and the breaking of time-reversal symmetry [45], many estimators rely on the direct link between the EPR and the difficulty of distinguishing between forward and reverse processes, quantified by the relative entropy, or the Kullback-Leibler Divergence (KLD), between them [21,[46][47][48][49][50][51].Calculating the KLD between probability distributions of forward and reverse trajectories for stationary data series can be done using different approaches [46][47][48][49]52].For example, the Plug-in method requires estimating the probabilities of sequences of data, discarding the information about transition times [48].
Applied to semi-Markov processes, the KLD breaks into two contributions, one of which captures irreversibility in the sequence of states, and the other captures irreversibility in waiting time distributions (WTD) [46].For second-order semi-Markov processes, this KLD estimator can detect and quantify entropy production even in the absence of observable currents [46,53,54].Machine-Learning (ML) tools have also been used for entropy production estimation by exploiting the irreversibility of data series [55][56][57][58].The core idea is to optimize an objective function whose extremum is the KLD between the forward and reverse trajectories of sequences of states.For example, the recurrent neural network estimator for entropy production (RNEEP) estimates the EPR from coarse-grained data of partially observed systems, using a recurrent neural network to solve the optimization problem [55].
In this work, we focus on a continuous-time Markov chain (CTMC) model over a discrete set of states, in which a subset of the microstates are coarse-grained, or  "lumped", into a single macrostate.We consider different levels of observed statistics from different coarsegraining (CG) approaches, and infer the EPR from the observed data using the KLD estimator [46], the Plug-in estimator [48,49], and the RNEEP estimator [55], when applicable.These estimators do not require prior knowledge of the systems, and only use the observed statistics to infer and quantify time-irreversibility. First, we use the sequence of observed microstates and coarse-grained macrostates, and the transitions between them.Then, we include information about transitions between the hidden microstates within the coarse-grained macrostates (intra-transitions).Finally, we reformulate the trajectory data of observed states and transitions, and intratransitions within macrostates, by labeling the coarsegrained macrostates according to the number of times they are visited before jumping into an observed state.
We apply the CG approaches to two model systems, namely, a 4-state (Fig. 1) system in which two of the states are coarse-grained into a single hidden state, and the discrete Flashing Ratchet (Fig. 2) with time-varying potential, whose values cannot be observed.We provide a unifying comparison between the estimators on the different CG schemes, and show how additional information is exploited for inferring tighter lower bounds on the total EPR.We begin by explaining the different coarse-graining approaches, taking the 4-state model system as an example (Fig. 1a).In the first CG approach, termed full-CG, we lump together a subset of the microstates into a single observed state, giving rise to a second-order semi-Markov process (Fig. 1b), since the waiting time in the hidden state depends on the state visited before [46].In this example, states 1 and 2 are observed, whereas states 3 and 4 can not be distinguished and are recorded as a single state, H. Here, the waiting time in H is the sum of the corresponding waiting times in the microstates 3 and 4 before jumping to one of the observed states.
In the second CG approach, termed semi-CG, we assume an observer can record intra-transitions within the hidden states (Fig. 1c).For example, a sequence of 1 → 4 → 3 → 2 is recorded as 1 → H → H → 2, with the corresponding waiting times, i.e., the time spent in the first visit to H and the time spent in the second visit to H are recorded separately.In this case of observed intra-transitions within a coarse-grained state, the initial and final microstates are not known, as both are lumped together to the same macrostate.Still, the added information can be utilized for improving the lower bound on the total EPR.The Plug-in and the RNEEP estimators rely on the sequence of states, and can be directly applied to the semi-CG trajectory.However, in order to apply the KLD estimator to the semi-CG data, we first need to reformulate the trajectory to a second-order semi-Markov process, and harness the information of the WTD.The transformation, depicted in Fig. 1d, consists of two steps.First, we look for all the consecutive sequences of the hidden state H, and record their length.Second, all sequences with the same length are considered a new state, i.e., a sequence of n appearances of H is labeled H n .The waiting time associated with H n is now the sum of the individual waiting times in the n consecutive appearances of H.This new representation of the semi-CG observed data gives rise to a second-order semi-Markov process, from which we infer a tighter lower bound on the total EPR.
The Plug-in estimator, σ plug , was proposed for approximating the KLD rate between the forward and reverse sequences of discrete stationary time series, by counting sequences of data and calculating their probabilities [48,49].The approximated mth-order KLD between sequences of length m is: where p(x 1→m ) and p(x m→1 ) are the probabilities of a forward sequence x 1→m = (x 1 , ..., x m ) and the backward one x m→1 = (x m , ..., x 1 ).These probabilities can be estimated from the number of appearances of each sequence in a long trajectory.Based on the approach in [48], the slope of D x m as a function of m, gives the entropy production per step in the limit of large m.However, for a non-Markov process that cannot be described by a semi-Markov process of any order, calculating dx m is challenging for large values of m.Therefore, the following ansatz [59] has been proposed: where dx ∞ , c, and γ, are the fit parameters for dx m as a function of m.Our Plug-in estimator for the entropy production rate per time, is thus: where τ is the mean waiting time in each step.Note that this estimator can be directly used for both semi-CG and full-CG partial information framework without any modifications to the trajectory data.The KLD estimator, σ KLD , derived by calculating the KLD between forward and reverse trajectories in semi-Markov processes, has two contributions [46]: where the affinity, σ aff , stems from observed currents, and the σ WTD stems from time-asymmetries in WTD.
In order to apply Eq. 5 to second-order semi-Markov processes, the observed states are reformulated as doublets, [ij], where the first index is the previous state, and the second index is the current state [46].The affinity contribution is: where p (ijk) is the probability to observe the sequence ) being the probability to jump to state k after jumping from i to j, and R [ij] being the fraction of visits to [i, j].The affinity, σ aff , is governed by the relation between the forward and reverse transition probabilities.
The WTD contribution stems from the Kullback-Leibler divergence between WTD associated with forward (i → j → k) and backward (k → j → i) transitions: ) is the WTD in state j given that the previous state was i and the following is k, τ is the average waiting time per state, and D[u(x)||v(x)] is the Kullback-Leibler Divergence between two probability distributions, u(x) and v(x), defined as D[u(x)||v(x)] = x u(x) ln (u(x)/v(x)).See Supplemental Material (SM) [60] for details regarding WTD estimation.Note that for a fully observed system following Markovian dynamics, σ WTD vanishes, and one cannot infer non-zero EPR without observable currents.
The RNEEP estimator, σ RNEEP , is formulated as an optimization problem [55], with a specific objective function to be minimized using stochastic gradient descent.The input of the problem is the set of all sequences of length m from a single long trajectory, and the solution is the coarse-grained entropy production rate per step along the input trajectory.Similar to the Plug-in estimator, the RNEEP uses the discrete sequence of states and does not exploit the WTD data, so estimating the full probability distributions of the waiting times is not required.Intuitively, this estimator should yield similar results to the Plug-in estimator, Eq. 4, and to σ aff , Eq. 6, as it uses the same information (see SM for further discussion [60]).The RNEEP can be directly applied to both full-CG and semi-CG frameworks, and it can be implemented by different machine learning models, such as recurrent or convolutional neural networks [55,56].Following the approach of [55], we use a recurrent neural network, whose input is a sequence of some length m, x m t = (x t , x t+1 , ..., x t+m−1 ), and its output is h θ (x m t ), where θ represents the learnable weights of the network.The output of the RNEEP is [55]: FIG. 3. Entropy production rates for the 4-state system.Total EPR, σtot (dark red), KLD estimator, σKLD, for the semi-CG (dark blue) and full-CG (light blue) data, Plug-in estimator, σ plug , for the semi-CG (dark orange) and full-CG (light orange), RNEEP estimator, σRNEEP,m, for the semi-CG (light to dark purple for increasing sequence length m) and full-CG (light to dark green for increasing sequence length m) data, and the affinity contribution, σ aff , for the full-CG data (red).Rates: where xm t is the time-reversed sequence of x m t .The RNEEP estimator is the solution of the optimization problem of minimizing the following objective function over ∆S θ (x m t ) for all possible sequences of length m: where E t in the expectation over t, and E (x m t ) is the expectation over the observed sequences x m t .See SM [60] for a detailed explanation of the implementation of the estimators and the different numerical considerations.
We have evaluated the performance of the three EPR estimators on two coarse-grained systems, the 4-state system with 1 coarse-grained state, and the discrete flashing ratchet with the unobserved external potential.For each system, the two CG approaches, namely, the full-CG and semi-CG, were applied on trajectories of approximately N = 10 7 states, simulated using the Gillespie algorithm [61].The code was implemented in PyTorch is available in [62].
The Plug-in estimator was fitted by gathering statistics of sequences of various lengths according to Eq. 3, as done in [48].The RNEEP was calculated for different input sizes, where we have used the implementation of D.-K.Kim, et.al. [55] with some adjustments for our hardware, see [60].The KLD estimator was applied without modifications on the full-CG data [46], whereas the trajectory reformulation was used only for calculating σ KLD for the semi-CG statistics.
The results for the 4-states system (Fig. 1(a)) under the two CG schemes, semi-CG and full-CG, and the three estimators, RNEEP, Plug-in, and KLD, are presented in Fig. 3.The system has two observed states, 1, and 2, where the states 3 and 4 are coarse-grained into a single state H.The rates between the two observed microstates are tuned according to ω 12 = ω12 e x and ω 21 = ω21 e −x to mimic an external forcing, where the range of x was chosen to include the stalling force in which there is no observable current over the 1 − 2 link [21].
As expected, The bounds on the total EPR obtained from estimators applied to the semi-CG statistics are better compared with the same estimators applied to the full-CG trajectories.In the full-CG case, σ RNEEP,m is similar to the affinity, σ aff , as both estimators use the same data.The Plug-in estimator, σ plug , which also uses the same data of the full-CG trajectory, provides similar results to σ aff away from the stalling force.However, close to the stalling force, where σ aff vanishes, σ plug provides non-zero values that stem from the inherent bias of the method which assigns positive values to all the probabilities [49] (See [60]).The KLD estimator, σ KLD , provides the tightest lower bound for the full-CG data, as it is the only one that utilizes information of the irreversibility in WTD.Total EPR, σtot (dark red) of the full trajectory, KLD estimator, σKLD (dark blue), Plug-in estimator, σ plug (dark orange), and RNEEP estimator, σRNEEP,m (light to dark purple for increasing sequence length m) for semi-CG data.The transition rates for ∀i = j are: where states i (i ) are with turned-on (off) potential.
In the semi-CG case, the RNEEP estimator, σ RNEEP,m , provides a tighter bound for increasing sequence lengths m, and it converges to σ plug as in the full-CG case [60].The KLD estimator σ KLD , provides the tightest lower bound on the total EPR compared to the other estimators tested, given the reformulation of the trajectory data according to the Transformed semi-CG scheme (Fig. 1(d)).The difference in σ KLD between the Transformed semi-CG and the full-CG schemes reflects the additional information regarding irreversibility encoded in the intra-transitions between microstates in the hidden macrostate.Moreover, the irreversibility encoded in these intra-transitions in the semi-CG data is also reflected in the values of the 3 estimators that do not vary significantly near the stall force, in contrast to the estimators applied to the full-CG trajectories that strongly depend on the deviation from stalling conditions.
The results for the discrete flashing ratchet system [48,49] under the semi-CG scheme and the three estimators, RNEEP, Plug-in, and KLD, are presented in Fig. 4. The system is of a Brownian particle moving along a periodic one-dimensional line, under the influence of a linear potential V that can be switched on and off at a constant rate.The particle is described by its position in the "on", i, or "off", i , states.Under CG, the information on the potential is not accessible, and both the on and off states are lumped into a single macrostate H i .Here, we apply the 3 estimators to the non-Markovian semi-CG data, which includes the information of the intra-transition within the macrostates.Note that in [55], the RNEEP was compared to a semianalytical calculation of the KLD between trajectory distributions.However, it was shown that the Plug-in estimator yielded similar results to the semianalytical values for the semi-CG observed statistics [48].Similar to the 4-state system results, the RNEEP estimator, σ RNEEP,m , provides a tighter bound for increasing sequence lengths m.Moreover, the application of the KLD estimator to the Transformed semi-CG date yields the tightest lower bound on the total EPR, compared to the Plug-in [48,49] and the RNEEP [55].
In summary, we have compared time-irreversibilitybased EPR estimators for different coarse-graining schemes, focusing on KLD-based estimators, with (KLD) and without (Plug-in and RNEEP) WTD statistics, using two coarse-graining approaches.We have confirmed that the semi-coarse-graining framework, which includes intra-transitions data, yields tighter EPR bounds compared to the full-coarse-graining framework, as it exploits more information on time-irreversibility.In addition, we have proposed a novel approach for reformulating semi-CG trajectories, previously used for the Plug-in [48,49] and RNEEP estimators [55], for applying the KLD estimator [46].Using the Transformed semi-CG approach and the KLD estimator, we could distill time-irreversibility encoded in the intra-transitions within hidden microstates to achieve the tightest lower bound on the total EPR among the estimators we tested.Moreover, comparing the EPR bounds obtained from the full-CG and the semi-CG statistics provides a direct quantification of the time-irreversibility in the intra-transitions captured by each estimator.The proposed transformation and EPR estimators can be applied to other discrete-state continuous-time systems, to provide a lower bound on the total EPR, when only partial information is available.

FIG. 1 .
FIG.1.Illustration of the partial information frameworks for arbitrary 4-states system.(a) A fully-observed 4-states system.The trajectory (blue line) is described by the sequence of microstates and the corresponding waiting times (WT).(b) Full coarse graining (full-CG).States 3 and 4 cannot be resolved, and are lumped together to a single macrostate H (orange line) (c) Semi-coarse graining (semi-CG).States 3 and 4 cannot be resolved, but intra-transitions between the hidden microstates can be recorded, where consecutive visits in the hidden microstates (orange line with black markers) are recorded as sequence of H, with the corresponding WT of each hidden microstate between intra-transition events.(d) Transformed Semi-coarse graining (Transformed semi-CG).Each n consecutive visits to the hidden microstates in H are recorded as Hn (light orange, orange, and brown represent different sequence lengths) and the WT in Hn are the sum of the WT in consecutive visits in the hidden microstates 3 and 4.

FIG. 2 .
FIG.2.Discrete flashing ratchet of 3 states with periodic boundaries, and a potential that can be switched on (i) and off (i ) but is not accessible to the observer.State of the same position i, regardless of the potential, are coarse-grained into a single macrostate H1, H2, and H3.In the semi-CG framework, intra-transitions i ↔ i are recorded.

2 FIG. 4 .
FIG.4.Entropy production rates for the flashing ratchet.Total EPR, σtot (dark red) of the full trajectory, KLD estimator, σKLD (dark blue), Plug-in estimator, σ plug (dark orange), and RNEEP estimator, σRNEEP,m (light to dark purple for increasing sequence length m) for semi-CG data.The transition rates for ∀i = j are:ω ii = ω i i = ω i j = 1, ωij = e (V j −V i )/2 ,where states i (i ) are with turned-on (off) potential.

G.
Bisker acknowledges the Zuckerman STEM Leadership Program, and the Tel Aviv University Center for AI and Data Science (TAD).This work was supported by the ERC NanoNonEq 101039127, the Air Force Office of Scientific Research (AFOSR) under award number FA9550-20-1-0426, and by the Army Research Office (ARO) under Grant Number W911NF-21-1-0101.The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Office or the U.S. Government.