Inferring entropy production rate from partially observed Langevin dynamics under coarse-graining.

The entropy production rate (EPR) measures time-irreversibility in systems operating far from equilibrium. The challenge in estimating the EPR for a continuous variable system is the finite spatiotemporal resolution and the limited accessibility to all of the nonequilibrium degrees of freedom. Here, we estimate the irreversibility in partially observed systems following oscillatory dynamics governed by coupled overdamped Langevin equations. We coarse-grain an observed variable of a nonequilibrium driven system into a few discrete states and estimate a lower bound on the total EPR. As a model system, we use hair-cell bundle oscillations driven by molecular motors, such that the bundle tip position is observed, but the positions of the motors are hidden. In the observed variable space, the underlying driven process exhibits second-order semi-Markov statistics. The waiting time distributions (WTD), associated with transitions among the coarse-grained states, are non-exponential and convey the information on the broken time-reversal symmetry. By invoking the underlying time-irreversibility, we calculate a lower bound on the total EPR from the Kullback-Leibler divergence (KLD) between WTD. We show that the mean dwell-time asymmetry factor - the ratio between the mean dwell-times along the forward direction and the backward direction, can qualitatively measure the degree of broken time reversal symmetry and increases with finer spatial resolution. Finally, we apply our methodology to a continuous-time discrete Markov chain model, coarse-grained into a linear system exhibiting second-order semi-Markovian statistics, and demonstrate the estimation of a lower bound on the total EPR from irreversibility manifested only in the WTD.

Introduction

Irreversible processes in living systems lead to the production of entropy, which is a measure of energy dissipation and a signature of the arrow of time.[1-5] Quantifying the entropy production can shed light on the underlying nonequilibrium dynamics and provide insights into the thermodynamic burden of biological processes.[6-8] There are primarily two methods to infer that a system is out-of-equilibrium: (i) invasive methods,[9-12] and (ii) non-invasive[13-16] methods. In invasive methods, the system’s response to a perturbation is measured following an external manipulation, and the violation of the fluctuation–dissipation theorem (FDT)[9,17-22] confirms the nonequilibrium nature of the underlying process. On the other hand, non-invasive methods do not require a direct perturbation to a system, and can detect the nonequilibrium nature of the process from various system properties, such as broken time-reversal symmetry,[23,24] presence of net probability current of observables,[7,13,16,25-29] or asymmetric probability density function (PDF) of the timing of maximal observable values.[30]

One can estimate the EPR for discrete[31] and continuous systems[32-34] given that all out-of-equilibrium system variables are accessible; otherwise, the EPR estimation becomes challenging,[35-39] and the best estimate would be a lower bound on the total EPR value. Several studies focused on the fluctuations of the EPR calculated from partial information.[40-48] The mathematical relations that bound the EPR using the fluctuations of time asymmetric and generic variables are known as the thermodynamic uncertainty relation (TUR)[49-54] and kinetic uncertainty relation (KUR),[55] respectively. These relations have also been generalized for semi-Markov processes.[56,57] Recently, a unified relation considering both thermodynamic and kinetic quantities has been proposed.[58-60] For systems with partial information, estimators like the passive partial entropy production[47,61,62] and the informed partial entropy production[61-64] are helpful to get a dissipation bound; however, these fail to provide a tight bound on the total EPR for vanishing net current. These average partial entropy production estimators satisfy fluctuations theorems, and as such, they can be derived as a Kullback–Leibler divergence between the forward trajectory and the backward trajectory under auxiliary dynamics.[61]

The k-variable irreversibility measure is defined as, , where k is the Boltzmann constant, D[p∥q] denotes the Kullback–Leibler divergence (KLD)[65,66] between two probability distributions p and q, defined by D[p ∥ q] = ∫dxp(x)log(p(x)/q(x)) and calculated on the positive support.[66] It is a measure of distinguishability[67] between two probability distributions, being non-negative in general and zero for identical distributions. Γ denotes the forward path of k nonequilibrium variables for a time duration t, whereas denotes the corresponding backward path. Owing to the chain-rule of the relative entropy,[68] the more nonequilibrium variables (larger k) included in the path probability measure, the better the KLD bound is, i.e., 0 ≤ σ1 ≤ … ≤ σ ≤ σ ≤ … ≤ σtot where σtot is the total EPR calculated by the KLD between the forward and reverse trajectories with all the nonequilibrium degrees of freedom.[69,70] Obtaining a tight bound for a continuous variable system using the KLD estimator is challenging since some of the nonequilibrium variables may be inaccessible and sampling the distribution of paths becomes difficult.

In a recent study,[71] Roldan et al. transformed the forward and backward time series data of an observed variable of a continuous hair bundle system into two independent and identically distributed time series using a whitening approximation to estimate the KLD from two univariate distributions. They first calculated the EPR bound using only the observed degree of freedom, i.e., the tip position of the hair bundle. Moreover, they used the finite time thermodynamic uncertainty relation to obtain a lower bound on the total EPR using two observables, the tip position and the transduction current, and found a better lower bound on the EPR compared to the one calculated using only one variable, as expected. The EPR estimate calculated with only one observed degree of freedom was typically three orders of magnitude smaller than the total EPR. However, using two observables and the TUR, their measure was three orders of magnitude better than their single-variable result for the oscillatory regime and few fold smaller than the total EPR, but in the quiescent regime, the result was three orders of magnitude smaller than the total EPR.

An estimator based on the KLD between waiting time distributions of the time forward and the time backward transitions between discrete states was shown to provide a lower bound on the total EPR,[64] given that the time-reversal operator does not lead to kinetic hysteresis.[72-75] Applied to a second-order semi-Markov process, this KLD estimator of the EPR breaks into two contributions,[64] the affinity EPR, EPRaff, which accounts for the net flux and affinity or the thermodynamic force,[68,70] and the waiting-time-distribution (WTD) EPR, EPRWTD, which accounts for the broken time-reversal symmetry in the waiting time distributions.[64] For second-order semi-Markov processes, which naturally emerge when “lumping” adjacent states,[64,76] the EPRWTD can provide a lower bound on the total EPR, even when the system does not have any net current observed and EPRaff = 0. Describing processes by transitions instead of states,[77] the KLD estimator for the EPR was further applied to waiting times in between observed transitions.[75,78]

Skinner et al. presented new estimators to obtain the lower bound on the entropy production rates using optimization techniques.[79,80] They found an estimator given observables characterizing one-step transitions and two successive transitions, whereas in another publication the authors proposed an estimator given the observed waiting time distributions.[80]

There are several studies on the effect of coarse-graining (CG) on the EPR estimation[81-94] specifically discussing whether the CG procedure preserves the EPR fluctuations or not. In a recent study, using a Markovian model of a driven molecular motor, Hartich et al. compared different coarse-graining schemes, “milestoning” and “lumping”, and found that the “milestoning” method can restore Markovian dynamics in the case of time-scale separation and preserves local detailed balance.[76,86]

The quantitative effect of the coarse-graining on the EPR was estimated in an experimental system of steady-state trajectories of a microtubule length using an optimization procedure of a two-step estimator, where it was demonstrated that increasing the spatial and temporal resolution while coarse-graining leads to an improved EPR bound.[79] Moreover, a recent study by Tan et al.[95] has found that the time-irreversibility varies non-monotonically with the lag time, i.e., the time intervals between the position measurements, which determines the dissipation timescale.[95]

Here, we quantify the irreversibility using a non-invasive method to provide a lower bound on the total EPR in a partially observed model system with continuous variables following oscillatory dynamics, where one of its observables is coarsegrained into a few discrete states. We simulate an oscillating hair-bundle model in which the bundle’s tip position is experimentally observed, whereas the position of the molecular motor is hidden. The coarse-grained process follows second-order semi-Markov statistics in the reduced state space (tip position variable space). In this model, the affinity entropy production contribution vanishes; therefore, the irreversibility information can only be accessed from the asymmetries of the waiting time distributions of the forward and the reversed transitions. After the decimation, we exploit the underlying broken time-reversal symmetry stemming from the difference in the PDFs of the waiting times for the upward and the corresponding downward transitions among different coarsegrained states, to calculate the EPR bound, EPRWTD, by applying the KLD estimator. We show that the ratio of the means of the dwell time PDFs of the forward and reverse trajectories, termed the mean dwell-time asymmetry factor, can qualitatively detect the broken time reversal symmetry, and its variation with the number of coarse-grained states is studied. We further calculate the ratio between the EPRWTD and the total EPR as a function of the number of coarse-grained states to evaluate the tightness of the lower bound, and find that with finer resolution (larger number of coarse-grained states), the EPRWTD provides a better lower bound on the dissipation rate.

The paper is organized as follows. First, we introduce the model system and outline the calculation of the total EPR. Then, we describe our coarse-graining procedure, second-order semi-Markovian dynamics of the coarse-grained system, different contributions to the EPR, and mean dwell-time asymmetry factor in the next section. Subsequently, the effect of coarse-graining on the broken time-reversal symmetry, the EPR estimate, and the tightness of the lower bound as a function of the number of coarse-grained states are discussed. Finally, we summarize and provide a future outlook.

Model system

We estimate the entropy production rate in a partially observed system described by a Langevin equation. To do so, we consider a model which captures the experimental observation of spontaneous oscillations of mechanosensory hair bundles of auditory hair cells.[71,96-100] These oscillations help to amplify the sound stimuli in the ear of vertebrates, and provide sensitivity and frequency selectivity. Moreover, these oscillations are known as “active” oscillations, and they are distinct from “passive” oscillations that are obtained by blocking the corresponding transduction ion channels.[71] The activity originates from various molecular motors, which cannot be experimentally accessed. However, another degree of freedom coupled to the activity of the molecular motors – the tip position of the hair bundle (X1) is experimentally observed. Due to the presence of activity, the system is out-of-equilibrium, and its dynamics is governed both by a conservative force V(X1,X2), where X represents the position of the center of mass of the molecular motors, and a non-conservative driving force, Fact(X1,X2). The system can be described by the following coupled stochastic differential equations.[71,96-98] where λ1 and λ2 are the friction coefficients of the hair bundle tip and the molecular motor, respectively, T and Teff are the environment temperature and the effective temperature characterizing the motor fluctuations, respectively, with ratio Teff/T > 1. ξ1 and ξ2 are two independent white noise terms with zero-mean and correlation and k is the Boltzmann constant. The functional form of the conservative force, V(X1,X2), which is proportional to the difference between the positions of the coupled variables,[96-98] is: where kgs and ksp are the stiffness coefficients, ΔX = X1 — X2 is the separation between the position of the hair bundle and the molecular motors, D is the gating swing, and N is the number of transduction channels. A = exp[(ΔG + (kgsD2)/2N)/(kBT)], and ΔG is the energy difference between the open and closed states of the ion channel. The active non-conservative force exerted by the molecular motors is defined by Fact(X1,X2) = Fmax (1 — SP0(X1,X2)). The probability of the transduction channel being open is P0(X1,X2), and is defined by P0(X1,X2) = 1/[1 + A exp(—kgsDΔX/NkBT)]. The non-conservative force depends on the maximum motor force acting on the system (Fmax), and the calcium-mediated feedback strength (S). The main sources of the non-equilibrium drive come from the ratio Teff/T being greater than unity, and the maximal force (Fmax) exerted by the molecular motors. This model[71,96-98] was shown to agree well with experimental results.

First, we numerically solve the coupled differential equations (eqn (1) and (2)) for a fixed ratio between the effective temperature and the temperature of the environment (Teff/T = 1.5), and different values of S (0.5, 1, 1.5) and Fmax (70 pN, 80 pN, 90 pN) to obtain simulated trajectories of the hair bundle tip position and the motor position (see Fig. 1 for details on all the parameters used). Although there is clearly a directional current in the X1–X2 plane (Fig. 1a) manifesting the nonequilibrium nature of the process, its signature is not obviously present in the trajectories of X1 or X2 as a function of time, which oscillate around their respective mean values (as shown in Fig. 1b and c) for a particular set of the driving parameter values, and ESI,† Fig. S1 for additional realizations with different parameters).

Fig. 1

Simulated trajectories of the observed (X1, the tip position of the hair bundle) and the hidden (X2, the position of the molecular motors) variables and the coarse-grained trajectory of the observed variable after spatial coarse-graining (a) the trajectories in the X1–X2 plane for fixed values of driving parameters: (Fmax = 70 pN, S = 1, Teff/T = 1.5). The color of the curve represents time going from dark to bright (b) X2 = X2 — 〈X2〉, as a function of time for fixed values of driving (Fmax = 70 pN, S = 1, Teff/T = 1.5). (c) X1 = X1–〈X1〉, as a function of time for the same values of driving parameters, which does not show any sign of net flux, (d) the coarse-grained trajectory for 3 CG states at above mentioned parameter values. All the quantities plotted are calculated for the following additional parameter values: λ1 = 2.8 pN ms nm–1, λ2 = 10 pN ms nm–1, kgs = 0.75 pN nm–1, ksp = 0.6 pN nm–1, D = 61 nm, kBT = 4 pN nm, ΔG = 10kBT.

As the system is driven out-of-equilibrium by the non-conservative force and the effective temperature, there is a positive dissipation rate. The total entropy production rate can be calculated from the forces and their conjugated currents:[71,101] where 〈…〉 represents the steady state average. The steady state rate of the dissipated heat to the reservoir at temperature T is , with ° being the Stratonovich product, and is the rate of work done by the active force.

Coarse-graining, lower bound on the total entropy production rate, and the mean dwell-time asymmetry factor

We used two approaches for spatial coarse-graining to discretize the continuous variable space (the trajectories of the tip position of the hair bundle, X1) into discrete states: (i) dividing the continuous variable space equally into N (N = 3, 4, 5, 6, 7) coarse-grained states with the ratios 1:1:1, 1:1:1:1, 1:1:1:1:1, 1:1:1:1:1:1, and 1:1:1:1:1:1:1, respectively. This type of equal coarse-graining is only possible for a smooth trajectory for a particular choice of the driving parameter values (Fig. S1, ESI,† e.g. Fmax = 70 pN, S = 1, and Fmax = 80 pN, S = 1), (ii) diving the continuous variable space into unequal division, where N (N = 3, 4, 5, 6, 7) coarse-grained states correspond to dividing the trajectory with the ratios and , respectively, as shown schematically in Fig. S2 of the ESI.† This type of coarse-graining is better suited to track the irregular oscillations of the tip of the hair cell bundle for driving parameter values Fmax = 90 pN, S = 1, and Fmax = 80 pN, S = 1.5, etc. (see ESI,† Fig. S1).

We have two layers of coarse-graining: (I) one of the dynamical variables describing the system is decimated (in our example, the tip position of the hair bundle is observed, but the positions of the molecular motor are hidden) (II) we further coarse-grained the observed variable into a few discrete states.

Our system is coarse-grained such that the topology of the coarse-grained system is linear, without any cycles. The probability for a transition between the neighbouring states is non-zero, but the transition probability from one boundary state to the other boundary state is zero, and vice versa. For example, in a 3 coarse-grained system (N = 3, 1: 1: 1 spatial division), the probabilities of jumping from macro-state 2 to state 3 or 1 are both non-zero, whereas given the system is in state 1, the probability of finding it in state 3 in the next jump is zero, and vice versa. The waiting time distribution of the dwell time at state 2 depends, however, on the state visited before, whether it was state 3 or state 1, rendering the process a second order Markov process. Thus, we consider states composed of the current state, i, and previous state, j, i.e., [i,j] when applying the KLD estimator. Similarly, the approach can be generalized to higher order semi-Markov processes.

Estimating dissipation is non-trivial in the absence of the observable currents or flows, but as dissipated systems exhibit broken time-reversal symmetry, time irreversibility can be exploited to infer the out-of-equilibrium nature of the underlying process from the time series.[64] Martínez et al. developed an estimator based on the waiting time distributions containing information about irreversibility in hidden states even in the absence of visible transitions among the observed states. They applied the technique[64] for a partially hidden network where a subset of states are hidden, and a molecular motor system where the internal states are unresolved. In both cases, their estimator is able to predict a non-zero bound on the entropy production rate at the stalling driving force (the driving parameter value at which the current between the observed states vanishes).

To estimate the lower bound of the irreversibility, we use the KLD estimator,[101,102] which relies on the broken time-reversal symmetry of the underlying waiting-time distributions.[64] Due to the presence of coupled hidden degrees of freedom, the jump process in the observed variable space becomes a second-order[64] semi-Markov. The jump probability depends on the previous state, the time since the last jump, and the final state. The last two conditions make the system direction-time dependent,[91] which means that the joint distribution of times and transitions (ψ’(t)) cannot be written as a product of the probability distribution for a transition (Φ’) and the probability distribution of the time t the system waits at the initial state n (ψ. As proved earlier,[64] the KLD estimator of the EPR for a second-order semi-Markov process consists of two contributions: the affinity EPR (EPRaff) and the waiting-time-distribution EPR (EPRWTD). EPRaff accounts for the net current and the thermodynamic force of the system. It is sometimes called the “equivalent dissipation”.[91] A non-Markovian system and its memoryless counterpart – a system with the same network topology generating a Markovian sequence of states – have the same expression, but, the rate constants are replaced with the effective rate constants for the non-Markovian system. The affinity EPR is written as where p(ijK) = R → [jk]) is the probability to observe the sequence i → j → k. R denotes the normalized occupancy probability at the CG state j given the previous CG state was i. The numerator and the denominator of the argument of the logarithmic function are of the form p([ij] → [jk]), which denotes the probability that the system makes a transition from a CG state j to a CG state k, given that the previous CG state was i. τ is the mean step duration given by where τ[ is the mean time the system spends at a CG state j, given that the previous CG state was i. The sum is performed over all CG states (i,j, and k). For the active hair bundle system, there is no contribution to the EPR from the affinity EPR, since the coarse-grained system is a linear chain of states.

The other component of the KLD estimator comes from the broken time-reversal symmetry in the waiting-time distributions, and is obtained using the following equation: where Ψ(t|ijk) denotes the probability density function of the time t the system spends at a CG state j before jumping to another CG state k, given that the previous CG state was i, i.e., for i → j → k transition. The WTD estimator, EPRWTD, or the “memory dissipation”,[91] is the additional contribution that only exists for non-Markovian systems in contrast to their memoryless Markovian counterpart. It was shown that a semi-Markov process results in non-exponential waiting time distributions,[103] which is related to memory.[91]

Since there is no net current in the observed variable space, the position of the hair-bundle tip, X1, we use the KLD estimator[64] to calculate a lower bound on the total EPR. In order to apply this estimator, which was developed for discrete states, to a continuous variable system, we coarse-grain the observed variable into a few discrete states (a realization of 3 CG states is shown in Fig. 1d), from which the lower bound is estimated by EPRWTD, and study how the bound varies as a function of the number of coarse-grained states.

In order to demonstrate that a lower bound on the total EPR can be inferred from the WTD asymmetry in a system with second-order Markov process statistics with a linear topology having zero net current, we use a simple 6-state (i = 1, 2, 3 and i’ = 1′, 2′, 3′, where states i and i′ are indistinguishable) continuous-time Markov chain (CTMC) model coarse-grained into a 3-state linear continuous-time second-order semi-Markov system (observed states 1″, 2″, 3″) as shown in Fig. 2a. The net current in the 6-state model mimics the net current in the X1–X2 plane of the active hair bundle model Fig. 1a, whereas the coarse-grained 3-state system resembles the coarse-grained, observed hair-bundle position, X1. We simulated trajectories using the Gillespie algorithm[104] for 108 steps, where after the decimation, we were left with approximately 106 jumps. Fig. 2b shows the difference in the distribution of the times the system waits at state 2″ for an upward transition (1″ → 2″ → 3″) and the corresponding downward transition (3″ → 2″ → 1″). The non-exponential distribution originates from the non-Markovian statistics of the coarse-grained trajectory, whereas the difference between the distributions of the upward and downward waiting times originates from the nonequilibrium nature of the process.[64] Therefore, we can measure the irreversibility from the Kullback–Leibler divergence between the waiting time probability density functions EPRWTD, for the coarse-grained system with zero EPRaff to provide a lower bound on the total EPR.

Fig. 2

The X1 – X2 trajectory of the hair bundle system coarse-grained into a linear topology in X1 state space after decimation of the X2 states with zero net flux motivates to use KLD estimator of the waiting times (t [s]): (a) The circles with the lines represents a 6 state system, which after decimation is reduced to a linear 3 state system, (b) non-zero contribution from the Kullback–Leibler divergence of the waiting time distributions: the distribution of the waiting times (t [s]) the system waits at CG state 2″ for an (1″ → 2″ → 3″) upward transition (blue solid line) and (3″ → 2″ → 1″) the downward transition (red solid line) for the following parameter values: u1 = 10, u2 = 3, d1 = 2, d2 = 4, r1 = 3, r2 = 3, l1 = 1, l2 = 1.

For example, the waiting time distributions for the hair bundle system at equilibrium (Fmax = 0 pN, T = Teff) and at nonequilibrium conditions driven according to eqn (1) and (2) are shown in Fig. 3a and b, respectively. The distinguishability between the two WTD in the latter case (b), results in a positive KLD, which bounds the total EPR. The estimation of the EPRWTD improves with increasing the number of simulation steps (Fig. 3c) as evident from the decreasing error and the plateauing of the estimation value for the active hair bundle model governed by eqn (1) and (2).[64]

Fig. 3

Entropy production rate estimation from the Kullback–Leibler divergence between the waiting time distributions for 3 equally spaced coarse-grained states of the active hair bundle’s tip position: probability density functions of the time (t [s]) that the system stays at state 2 for an upward transition (blue solid line), and for a downward transition (red solid line) for two different parameter values: (a) Fmax = 0 pN, Teff = T, and S = 1.5, (b) Fma× = 70 pN, Teff = 1.5 T, and S = 1.5, (c) EPRWTD. as a function of length of the simulation for Fmax = 70 pN, Teff/T = 1.5, and S = 1.5. The error bar at each point describes the standard error of the mean.

The unimodal nature of the waiting time distributions also refers to the underlying network topology. If a network has internal cycles, the densities could exhibit multimodal behaviour.[75] For a second-order semi-Markov process, the waiting time distributions are direction-time dependent. Thus, the mean dwell-times that the system spends at a particular state for the forward and the reverse transitions are not necessarily identical, and a deviation of their ratio from one provides information regarding the irreversible nature of the process.[80] We calculate the mean dwell-time asymmetry factor (MDAF), i.e., the ratio between the means of the dwell time distributions (〈τ or 〈τkji〉) of times spent at a CG state j before transitioning to i, given that it arrived from k, k → j → i, to the mean time the system spends at a CG state j for a i → j → k transition, (〈τ〉 or 〈τ〉). The ratio between the mean times the system spends at a particular state before transitioning to another state and the mean times along the opposite direction (〈τ〉/〈τ〉) being not equal to unity indicates a broken time-reversal symmetry in the system. To obtain the total MDAF for a system with N coarse-grained states, we average the individual MDAF of different transitions among the N coarse-grained states. Therefore, the total MDAF equals N–1 ∑ (〈τ〉/〈τ〉). The ratios stemming from the transitions among different coarse-grained states are plotted in the ESI† (Fig. S3).

In the following, we calculate the contribution of the EPRWTD from eqn (6), and the effect of coarse-graining on the EPR and the MDAF, or the time-reversal symmetry breaking.

Effect of coarse-graining on the entropy production rate estimation and mean dwell-time asymmetry factor

We exploit the time-reversal symmetry breaking in the coarse-grained system to estimate the EPR. Since the affinity EPR vanishes, the signature of the irreversibility can only be tracked from the KLD between waiting time distributions, EPRWTD.

First, The EPR estimate (EPRWTD) values are calculated using eqn (6) by coarse-graining the X1 variable into N CG states (where N = 3, 4, 5, 6, 7) by equal partitioning of the state space, and plotted as a function of N (Fig. 4a), for Fmax = 70 pN, S = 1, and Teff/T = 1.5. The lower bound on the EPR estimate is improved with increasing resolution. The maximal value of EPRWTD/EPRtot = 0.0013 at 7 coarse-grained states. Moreover, the MDAF is plotted as a function of the number of the coarse-grained states (Fig. 4b).

Fig. 4

Entropy production rate estimation and the mean dwell-time asymmetry factor (MDAF) for equal spacing coarse-graining of X1 trajectory: (a) EPRWTD (s–1) (WTD estimate of the EPR) as a function of the number of CG states with equal spacing for parameter values Fmax = 70 pN, S = 1, and Teff/T = 1.5 (b) MDAF (mean dwell-time asymmetry factor) as a function of the number of CG states. The other parameter values are the same as mentioned in Fig. 1. The lines are drawn to as a guide to the eye. The total EPR for this set of parameter values is 7.3312 s–1.

Next, we calculate the EPRWTD for several driving parameter values (Fmax = 70 pN, Fmax= 80 pN, Fmax = 90 pN, and S = 0.5, 1, 1.5) and for unequal spatial spacing of the coarse-grained states (N = 3, 4, 5, 6, 7). Both the estimate of the EPR (Fig. 5a) and the mean dwell-time asymmetry factor (Fig. 5b) increase with increasing spatial resolution. Indeed, the EPR estimate is correlated with the MDAF (Fig. 5c), which is related to the non-Markovian nature of the process and the memory involved.[105]

Fig. 5

Effect of coarse-graining on the EPRWTD (s–1) and the mean dwell-time asymmetry factor (MDAF): (a) EPRWTD (s–1) as a function of number of CG states (3 CG: 1:1:1, 4 CG: 1:1/2:1/2:1, 5 CG: 1:1/3:1/3:1/3:1, 6CG:1:1/4:1/4:1/4:1/4:1, and 7CG: 1:1/5:1/5:1/5:1/5:1/5:1) for different parameter values, (left) Fmax = 70 pN for S = 0.5, 1, 1.5, (middle) Fmax = 80 pN for S = 0.5, 1, 1.5, and (right) Fmax = 90 pN, S = 0.5, 1, 1.5. (b) The MDAF as a function of number of CG states for different parameter values: (left) Fmax = 70 pN, S = 0.5,1,1.5; (middle) Fmax = 80 pN, S = 0.5, 1, 1.5; (right) Fmax = 90 pN, S = 0.5, 1, 1.5. In both panels: red circle symbols correspond to S = 0.5, blue square symbols correspond to S = 1, and magenta triangle symbols correspond to S = 1.5. (c) The values of EPRWTD as a function of the MDAF for all transitions and all parameter values as mentioned earlier. The other parameter values used in these figures are mentioned in Fig. 1.

As we mentioned, EPRWTD was calculated for equal (Fig. 4) and unequal (Fig. 5) partitioning of the observed trajectory. For certain driving parameter values at which the trajectories are not that smooth or regular. In that case, the equal partition of the trajectory space of the observed variable would lack enough statistics for the boundary states in the time series. Therefore, we consider unequal spatial partitioning of the trajectory.

To assess the tightness of the bound, we compare the ratio between EPRWTD estimates and the total EPR (EPRtot) calculated for different driving parameter values, Fmax = 70 pN, 80 pN, 90 pN, S = 0.5, 1, 1.5, and for different coarse-graining levels (Fig. 6), and find that the tightest bounds is obtained for 7 CG states (N = 7), where the EPRWTD values are between 1 to 2 orders of magnitude smaller than the total EPR (Fig. 6). The tightness of the bounds for unequal partitioning for 7 CG states are given in Table S1 of ESI.†

Fig. 6

Tightness of the EPR bound (EPRWTD) as a function of number of CG states: ratio between the EPR estimates from the waiting time distribution (EPRWTD (s–1)) and the total entropy production rate (EPRtot (s–1)) for different parameter values. The coarse-graining corresponds to unequal divisions of the X1 state space. The parameter values are Fmax = 70 pN, S = 0.5, 1, 1.5 (upper row), Fmax = 80 pN, S = 0.5, 1, 1.5 (middle row), Fmax = 90 pN, S = 0.5, 1, 1.5 (lower row). The other parameter values used in this figure are as mentioned in Fig. 1.

Discussion

Most of the previous studies on partially observed systems were performed on Markov chains where some nodes are observed, and the rest are either traced out or lumped together into a hidden state. These processes are performed with the constraints of preserving different quantities (depending on the applied coarse-graining method) like the transition flux among the observed states[93] or preserving the mean value and fluctuations of the entropy production rate at stationary state[87] before and after the coarse-graining. In this paper, we have discussed a different partially observed system where one of the coupled variables following the Langevin dynamics is observed experimentally, and the other is hidden. In addition, we have two layers of coarse-graining, where we preserve the equilibrium density of a particular state before and after the coarse-graining, but due to the linear topology, it cannot support current; therefore, it loses the net current of the original system. We have shown the benefit of using the waiting time distributions in estimating the dissipation rate using the hair bundle cell oscillations as an example. If the edge current vanishes in the observed states, the waiting time distributions may capture the broken time-reversal symmetry in the case of driven systems, depending on the network topology. We infer the irreversibility of the dynamics by coarse-graining the observed system variable into a few discrete states and applying the KLD estimator.[64] The coarse-grained linear system considered in our study is not Markovian, but rather a second-order semi-Markov system, and the breaking of time-reversal symmetry is manifested in the KLD between the non-exponential waiting time distributions of the forward and the reversed transitions among different coarse-grained states.[64] We show that instead of using the full probability distributions, the first cumulants of the dwell time distributions (easier to obtain in experimental scenarios), already provide predictions for the broken time-reversal symmetry and the dissipation rates. This quantity is much easier to quantify, both experimentally and theoretically, serving as a straightforward footprint for time-irreversibility. We further study the mean dwell time asymmetry factor variation with the number of the coarse-grained states.

Berezhkovskii et al.[105-108] discussed the case of low-resolution experimental observables in nonequilibrium systems, where the non-Markovian dynamics breaks time-reversal symmetry manifested in differences in the forward and backward waiting times. As suggested by several studies,[105-108] the time asymmetry in the active hair bundle system arises when the following two conditions hold: (i) the reduced variable system follows non-Markovian statistics, and (ii) the system is out-of-equilibrium. Using a 6-state CTMC model which is coarse-grained into a linear 3-state system (Fig. 2), (mimicking the hair cell bundle system with one degree of freedom is decimated), we demonstrate that the resulting waiting time distributions calculated by the Gillespie algorithm[104] show characteristics of second-order semi-Markov statistics, and break time-reversal symmetry under nonequilibrium driving, and thus KLD estimator would be the good choice for the estimation of the EPR. The 6-state network decimated into 3 states mimics the coarse-graining of the X1 trajectory into 3 coarse-grained states (Fig. 1d), in which a fundamental cycle is lost, and the contribution of the EPRaff vanishes. Indeed, we infer a lower bound on the total EPR, which can be calculated from the KLD between the distributions.

We calculate EPR estimates (EPRWTD) of the continuous-space model system, an oscillating hair cell bundle, after coarse-graining the observed X1 trajectory to equal (Fig. 4a) and unequal (Fig. 5a) spatial divisions. Comparing the results for a particular set of parameter values, Fmax = 70 pN, S = 1, and Teff/T = 1.5, for which the trajectory is rather smooth and regular (see ESI,† Fig. S1). For the equal and unequal coarse-graining, the lower bounds on the total EPR (i.e., EPRWTD/EPRtot) are 0.0013, and 0.0024, respectively at parameter values Fmax = 70 pN, S = 1, and Teff/T = 1.5.

The tightness of the lower bounds on the total EPR, i.e., EPRWTD/EPRtot, is found to be 0.0013 for equal spatial division (Fig. 4a) for N = 7 CG state at parameter value Fmax = 70 pN, S = 1, and Teff/T = 1.5. Whereas, for unequal spatial division (Fig. 6), EPRWTD/EPRtot equals to 0.1244 for N = 7 coarsegrained states at Fmax = 80 pN, S = 0.5, Teff/T = 1.5, respectively. The similar values of the EPRWTD/EPRtot ratio results from the smooth nature of the X1 trajectory at the chosen parameter set (as can be seen from Fig. 1c) in contrast to the other parameter values (ESI,† Fig. S1). Equal spatial division for N = 5, 6, 7 coarse-grained states becomes challenging for parameter values that lead to very rugged trajectories due to the lack of statistics for the boundary states.

The inferred time-irreversibility and the EPRWTD estimate increase with finer spatial resolution, i.e., larger number of CG states. Testing a wide range of parameter values, the EPRWTD is smaller by 1 to 2 orders of magnitude compared to the total ERP for the largest spatial resolution (N = 7) considered and unequal spacing of the observed X1 trajectory, where the tightest bound, EPRWTD/EPRtot ~ 0.1244, is obtained for Fmax = 80 pN, S = 0.5, and Teff/T = 1.5. All the ratios (EPRWTD/EPRtot) for 7 coarse-grained states are listed in Table S1 in the ESI.†

Conclusions

In summary, the hair bundle system was used as a model to study the effect of coarse-graining on the lower bound on the total entropy production rate, and the mean dwell-time asymmetry factor. The lower bound on the EPR was estimated using the underlying broken time reversal symmetry induced by the active force for a system with Langevin dynamics and zero net current along the reduced variable space. This approach can be applied to a system following Langevin dynamics with an arbitrary number of observed and hidden states carrying a net flux which vanishes on the observed state-space to quantify the deviation from thermal equilibrium manifested in the irreversibility of the observed degrees of freedom.