Quantifying entropy production in active fluctuations of the hair-cell bundle from time irreversibility and uncertainty relations

We introduce lower bounds for the rate of entropy production of an active stochastic process by quantifying the irreversibility of stochastic traces obtained from mesoscopic degrees of freedom. Our measures of irreversibility reveal signatures of time’s arrow and provide bounds for entropy production even in the case of active fluctuations that have no drift. We apply these irreversibility measures to experimental recordings of spontaneous hair-bundle oscillations in mechanosensory hair cells from the ear of the bullfrog. By analyzing the fluctuations of only the tip position of hair bundles, we reveal irreversibility in active oscillations and estimate an associated rate of entropy production of at least ∼3k B/s, on average. Applying thermodynamic uncertainty relations, we predict that measuring both the tip position of the hair bundle and the mechano-electrical transduction current that enters the hair cell leads to tighter lower bounds for the rate of entropy production, up to ∼103 k B/s in the oscillatory regime.


Introduction
Active systems are maintained out of equilibrium by processes that consume resources of energy and produce entropy. This is the case of living cells, where energy is provided in the form of a biochemical fuel such as adenosine triphosphate that drives active mesoscopic cellular processes. As discussed below, spontaneous oscillations of mechanosensory hair bundles of auditory hair cells provide an important example of active cellular fluctuations [1,2]. These oscillations have been proposed to amplify sound stimuli in the ear of many vertebrates, providing exquisite sensitivity and sharp frequency selectivity [3].
Active mesoscopic processes do not obey the fluctuation-dissipation theorem: measuring both the linear response of the system to weak external stimuli and spontaneous fluctuations provides a means to quantify deviations from thermal equilibrium [4][5][6][7][8][9][10]. A related important question is how entropy production can be estimated in active mesoscopic systems. In cases where active systems generate movement with drift, such as molecular motors moving along filaments [11][12][13], the rate of entropy production can be estimated from measurements of drift velocities and viscous forces [11,14]. More generally, to assess the entropy production in a stationary process, it is necessary to measure all the currents in the system and their conjugated thermodynamic forces. However, in many situations of interest, especially in biology, currents consist of non-equilibrium reactions, molecular conformational changes, or other types of flows that cannot be observed. Far from equilibrium, thermodynamic forces are also hard to identify even if the corresponding current is observable. As a consequence, it is relevant to develop methods for estimating entropy production when the available information about a system is limited.
Two strategies have been recently introduced. The first one is based on a quantitative measure of the irreversibility of a time series measured in the steady state. This quantitative measure is the relative entropy or Kullback-Leibler divergence (KLD), and it has been proved to provide a lower bound to the entropy production [15][16][17]. Because non-Markovian stochastic processes can be irreversible in absence of any drift, this approach can detect that a system is out of equilibrium even in the absence of any current or flow in the observed variables. This is the case of many experiments in hair cells, like the ones reported here. However, the method has been implemented mostly for discrete systems [15][16][17] and an extension to continuous variables is needed. The second method to estimate entropy production is based on the recently introduced thermodynamic uncertainty relations (TURs). These relations establish a lower bound to entropy production proportional to the relative fluctuations of currents [18,19]. Hence, the method is useful when currents can be observed but there is no access to the conjugated forces. Here, we introduce a novel technique to calculate the KLD for continuous variables. We analyze the feasibility and accuracy of the bounds provided by this new technique as well as the TUR by applying the two approaches to experimental and numerical data of the motion of hair bundles in the sensory cells of the ear.
Hair cells are the cellular microphones of the inner ear [20]. They transduce sound-evoked mechanical vibrations of their hair bundle-a cohesive tuft of cylindrical stereocilia that protrudes from their apical surface (figure 1(A))-into electrical signals that then travel to the brain. Fluctuations and response of the hair bundle provide a paradigmatic case study of nonequilibrium physics in biology [21][22][23]. Hair bundles from the ear of the bullfrog show noisy spontaneous oscillations [1]. Under periodic external stimulation, oscillatory hair bundles can actively amplify their response, resulting in a hysteretic behavior corresponding to a net energy extraction from the bundle [24]. Furthermore, it was shown that the fluctuation-dissipation theorem does not hold for oscillatory hair bundles, revealing that their spontaneous fluctuations are also active [4]. This finding demonstrates that hair bundle fluctuations must be described by nonequilibrium stationary states. Such behavior can be captured by a minimal two-variable stochastic model with nonlinear and non-conservative forces [1,2,[25][26][27]. In this model, one variable describes the tip position of the hair bundle, whereas the other variable describes the dynamics of a collection of molecular motors that power the bundle oscillations. Although the tip position can be measured, motors' fluctuations are hidden and hence can only be estimated from stochastic simulations [2,21,25]. Experimental and theoretical evidence led to the proposal that spontaneous fluctuations of the hair bund are akin to noisy limit-cycle oscillations close to a Hopf bifurcation [27]. Hence, as any active system, hair-bundle spontaneous fluctuations are characterized by probability fluxes in suitable phase spaces and by entropy production. Whether, and to what extent, tools from the emerging field of stochastic thermodynamics [28,29] can be used to estimate entropy production from measurements of active hair-bundle fluctuations remains an open question.
In this paper, we quantify irreversibility of the hair-bundle dynamics using experimental measurements of only the hair-bundle position X, an observable for which there is no measurable drift. From this limited information, we obtain lower bounds for the entropy production of its spontaneous fluctuations. In principle, one could also observe the mechano-electrical transduction current I produced by oscillatory hair bundles, revealing currents in the (X, I) plane [30]. Here, to test the accuracy of the bounds based on TURs [18,19], we use a well-tested dynamical model for oscillatory hair bundles to simulate trajectories in the (X, I) plane. Using this procedure, we obtain tighter bounds of the rate of entropy production of the hair bundle.
The paper is organized as follows. In section 2, we discuss generic properties of irreversibility and dissipation of mesoscopic nonequilibrium stationary states, as well as methods to quantify irreversibility from the statistics of few stochastic variables. In section 3, we provide estimates of irreversibility from experimental measurements of only the tip position X of hair bundles that show noisy oscillations. In section 4, we perform numerical simulations of a stochastic model of active hair-bundle fluctuations to compare the one-variable irreversibility estimates with the actual entropy production by the active system. In section 5, we use TURs to predict how much entropy production could be estimated by having access to both the tip position X and the transduction current I experimentally. Finally, in section 6 we discuss our main findings and conclude the paper. Mathematical derivations, details on experimental data analysis and on biophysical modeling are provided in the appendices.

Estimates based on time irreversibility
We first discuss the relation between entropy production and irreversibility for generic nonequilibrium stationary processes. Consider a physical system described by a set of variables labeled as X α , with α = 1, 2, . . . . In a stationary nonequilibrium process of time duration t, the physical system traces a trajectory in the phase space described by the stochastic processes X α (t). We denote by Γ [0,t] ≡ {(x 1 (s), x 2 (s), . . . ) )} t s=0 a given trajectory described by the system variables and its corresponding time-reversed trajectory asΓ [0,t] is the time-reversal signature of the αth variable. Assume now that X α are the variables that may be out of equilibrium, i.e. we do not include in Γ [0,t] those variables corresponding to thermal reservoirs, chemostats, etc. In that case, the steady-state rate of entropy production σ tot is given by where k B is the Boltzmann constant and P denotes the steady-state path probability [31][32][33][34]. Here D[Q R] 0 is the KLD between the probability measures Q and R, which quantifies the distinguishability between these two distributions. For measures of a single random variable x the KLD is given by Note that for isothermal systems, σ tot T is equal to the rate of heat dissipated to the environment at temperature T.
Often in experiments, only one or several of the nonequilibrium variables can be tracked in time. Consider the case where only X 1 , . . . , X k are known. We define the k-variable irreversibility measure in terms of path probabilities of k mesoscopic variables where )} t s=0 denote paths described by k variables. The k-variable irreversibility measure increases with the number of tracked degrees of freedom, providing a set of lower bounds to entropy production: It can also be shown that the estimator σ k equals the physical entropy production σ tot if the missing variables, X with > k, are at thermal equilibrium [35][36][37]. When the missing variables are not at thermal equilibrium, which is often the case in active systems, the estimate σ k σ tot yields only a lower bound for the entropy production rate. We now introduce a method to estimate the irreversibility measure σ 1 for any nonequilibrium steady state from a single stationary time series x i = X(iΔt) (i = 1, . . . , n) of a single variable X that is even under time reversal. We describe the technique for a single variable, but it can be generalized to several variables X α (t). In discrete processes, the KLD in σ 1 can be accurately measured from the statistics of sequences of symbols [15,16] and non-Markovian waiting time distributions [17]. In continuous processes however, estimating σ 1 is a herculean task due to the difficulties in sampling the whole phase space of paths [38][39][40].
The key idea of our method is to exploit the invariance of the KLD under one-to-one transformations. Suppose that there exists a one-to-one map ξ i (x 1 , . . . , x n ), i = 1, . . . , n, that transforms the original time series and its time reversal into two new time series ξ F i = ξ i (x 1 , . . . , x n ) and ξ R i = ξ i (x n , . . . , x 1 ) that are independent and identically distributed (i.i.d.) processes. Such a procedure is often called a whitening filter [41,42]. Because the new series are i.i.d., the KLD is now simple to calculate: it is given by the KLD between two univariate distributions p(ξ) and q(ξ), corresponding to the stationary probability distribution of ξ F i and ξ R i , respectively [40]. In general, it is not possible to find a one-to-one map that fully eliminates the correlations of both the forward (x 1 , . . . , x n ) and the backward (x n , . . . , x 1 ) time series. In that case, the removal of the correlations in the backward series is enough to provide a lower bound for σ 1 : where f s = (Δt) −1 is the sampling frequency and D[p(ξ) q(ξ)] = dξp(ξ) ln[p(ξ)/q(ξ)] is the KLD between the univariate distributions p(ξ) and q(ξ). We estimate D[p(ξ) q(ξ)] γ ip i ln(p i /q i ) wherep,q are empirical densities, and the sum runs over the number of histogram bins. We introduce the prefactor γ = 1 − p KS 1, where p KS is the p-value of the Kolmogorov-Smirnov (KS) statistic between p(ξ) and q(ξ), to correct the statistical bias of our KLD estimate [43]. The proof of the bound (4) and further details of the estimate are found in appendices A and B.

Estimates based on TURs
We now consider the case in which one can measure a set of k 2 mesoscopic variables X 1 , . . . , X k that are all even under time reversal. In this case, a measurable signature of irreversible dynamics is the emergence of currents in the (X 1 , . . . , X k ) plane. For Markovian nonequilibrium stationary states, it has been shown that the rate of entropy production can be computed analytically and expressed in terms of forces and conjugated currents, see e.g. [44][45][46]. For example, consider an isothermal overdamped Langevin system described by two mesoscopic variables,Ẋ where μ k , D k are the mobility and diffusion coefficient of the kth variable, F k (X 1 , X 2 ) is a generic state-dependent force, and η 1 and η 2 are two independent Gaussian white noises with zero mean For this example, the entropy production can be expressed in terms of forces F k and currents (velocities)Ẋ k as follows Here, · denote steady state averages, • the Stratonovich product [47,48] and we have used the shorthand notation F k ≡ F k (X 1 , X 2 ). This result can be generalized to nonequilibrium overdamped Langevin dynamics with more than two degrees of freedom, and similarly to Markov-jump processes with an arbitrary number of states [44,46]. The TUR approach, introduced in references [18,19], provides a bound to entropy production which is useful when one has not access to the forces. For finite-time trajectories in stationary processes, the following TUR holds [49][50][51] where Var[j(τ )] = j 2 (τ ) − j(τ ) 2 is the finite-time variance of any current, and τ > 0 the observation time. Equation (8) allows us to bound from below the entropy production rate of any nonequilibrium stationary state with statistics of any current j which may contain partial information about the dynamics of the system. For the example of a 2D stochastic model described in equations (5) and (6), one can construct a family of currents as follows where G 1 (t) ≡ G 1 (X 1 (t), X 2 (t)) and G 2 (t) ≡ G 2 (X 1 (t), X 2 (t)), with G 1 (x, y) and G 2 (x, y) two arbitrary functions. For this model, the TUR (8) implies that [52] which holds for any observation time τ and any choice of the functions G 1 and G 2 that enter in equation (9). Notably, the bound (10) also applies for currents that are obtained by coarse graining the phase space variables X 1 and X 2 . To which extent the bound (8) can be tightened by finding optimal currents has become an active area of research (see e.g. [53][54][55][56][57]). In particular, it has been shown that currents containing 'footprints of irreversibility' (i.e. information about stochastic entropy production) provide tight, optimal bounds to the finite-time TURs [54,56].

One variable irreversibility in active hair-bundle fluctuations
We now discuss irreversibility and entropy production in active mechanosensory hair cells from the bullfrog's ear. In experimental recordings of spontaneous hair-bundle oscillations, only the tip position X 1 of the bundle is measured (figures 1(B) and (C)). Hair-bundle oscillations take the shape of relaxation oscillations corresponding to an alternation of fast jumps between two extreme positions interspaced by dwell times. Measuring X 1 , we can only estimate σ 1 , which provides a lower bound to the total steady-state entropy production rate σ tot . We later compare this estimate to that obtained for a passive bistable system in a thermal bath (figure 1(D)).
In the following, we make use of autoregressive (AR) models for the whitening transformation. More precisely, we obtain the transformed time series We find that the irreversibility measureσ 1 given by equation (4) distinguishes active hair-bundle fluctuations (σ 1 > 0) from passive fluctuations of a bistable system (σ 1 0). Note that the estimate saturates to a plateau when the time series is long enough, in practice here longer than 10 s. Using a population of 182 hair cells that showed spontaneous hair-bundle oscillations [59], we obtain an exponential distribution ofσ 1 with mean value 3k B /s ( figure 2(B)). Interestingly, this result depends on the sampling frequency f s (see appendix C): irreversibility is maximal in the range f s ∼ (200-600) Hz where its value goes up to 4.3k B /s. This frequency dependency may provide additional information about timescales of the underlying active process [48].
We further quantify differences in irreversibility in typical examples of: (i) active oscillatory hair bundles (figure 3(A), top); (ii) hair bundles that were brought to quiescence upon exposure to a drug (gentamicin) that blocks the transduction channels (figure 3(A), middle); (iii) noisy signals produced by the recording apparatus when there was no hair bundle under the objective of the microscope ( figure 3(A), bottom). To further characterize differences in irreversibility, we apply the local irreversibility measure defined as which obeysŝ 1 (ξ) 0 for all ξ [60], andσ 1 = dξŝ 1 (ξ). We find that for all the analyzed values of ξ, the local irreversibility of active oscillations is ∼10 3 times larger than for passive hair-bundle fluctuations and experimental noise.

Entropy production rate of active hair-bundle fluctuations
We now relate the estimateσ 1 of entropy production from experimental recordings (figure 2(B)) to the entropy production σ tot which we obtain from stochastic simulations of hair-bundle oscillations. Spontaneous hair-bundle oscillations are thought to result from an interplay between opening and closing of mechanosensitive ion channels, activity of molecular motors that pull on the channels, and fast calcium feedback. This interplay can be described by two coupled stochastic differential equations for the position of the bundle X 1 and of the center of mass of a collection of molecular motors X 2 [2,59,61] (see appendix D): Here, λ 1 and λ 2 are friction coefficients and ξ 1 and ξ 2 in (12)-(13) are two independent Gaussian white noises with zero mean ξ i (t) = 0 (i = 1, 2) and correlation ξ i (t)ξ j (t ) = δ ij δ(t − t ), with i, j = 1, 2 and δ ij the Kronecker's delta. T is the temperature of the environment, whereas the parameter T eff > T is an effective temperature that characterizes fluctuations of the motors. This model can be interpreted as a system affected by two nonequilibrium constraints: the active force F act and its fluctuations whose magnitude is described by an effective temperature T eff that differs from the actual temperature T. Recall, however, that the real system is at temperature T and that the effective temperature T eff has the same origin as the active force, which is the activity of the motors. The conservative forces derive from the potential associated with elastic elements and mechano-sensitive ion channels where  We remark that the hair-bundle tip position X 1 is measurable experimentally whereas the motors' position X 2 is a hidden variable that is not accessible experimentally. A rationale for the choice of parameter values can be found in references [2,25] and in appendix D. In (A-C) we substracted to X 1 and X 2 their respective mean values.
is an active nonconservative force exerted by the molecular motors with a maximum value F max . The parameter S quantifies calcium-mediated feedback on the motor force [25] and is the open probability of the transduction channels. Note that equation (16) is the open probability of a two-state equilibrium model of a channel with a free-energy difference between open and close states depending linearly on the distance X 1 − X 2 . The parameter values in the model are constrained by experiments as discussed in [2,59,61]. As shown earlier [2,25], this model can capture key features of noisy spontaneous oscillations of hair-bundle position X 1 that have been observed experimentally ( figure 4(A)). The oscillation of the motors' position ( figure 4(B)) is known in the model but hidden in experiments. Trajectories of only X 1 (t) or X 2 (t) do not reveal obvious signs of a net current, which here would correspond to a drift. However, trajectories in the (X 1 , X 2 ) plane show a net current which is a signature of entropy production ( figure 1(C)). In the following, we will use this stochastic model to compare the irreversibility measureσ 1 to the total entropy production σ tot . In the stochastic model of hair-bundle oscillations given by equations (12) and (13) we deal with only two variables, therefore σ tot = σ 2 . From the analytical expression of σ 2 , we find that the steady-state entropy production rate can be written as (see references [46,62,63] and appendix E) where are the conjugated forces to X 1 and X 2 , respectively. We have also introduced the steady-state average heat − Q 1 = − (∂V/∂X 1 ) •Ẋ 1 dissipated to the thermal bath at temperature T and the power Ẇ act = − F act •Ẋ 2 exerted by the active force on the motors. Equation (17) takes into account the nonequilibrium effects produced in the model by the activity of the motors: the difference of effective temperature and temperature, and the active force.
We performed numerical simulations of equations (12) and (13) for different values of the control parameters F max and S (figure 5) to explore entropy production throughout the state diagram of the system. The quality factor Q of the oscillation -given by the ratio between the oscillation frequency and the bandwidth at half the maximal height of the power spectrum (see appendix F)-and the average open probability P o at steady state are displayed in figures 5(A) and (B) in the state diagram. The irreversibility measureσ 1 for trajectories X 1 (t) of spontaneous oscillations is shown in figure 5(C). This measure can be compared to the quantification of total entropy production σ tot of the model, given by equation (17), which is shown in figure 5(D). Irreversibility of trajectories and total entropy production correlate strongly. As expected,σ 1 provides a lower bound to the actual dissipation rate. Actually, the rate of entropy production estimated fromσ 1 is here typically three orders of magnitude smaller than the total entropy production. Clearly, measuring other degrees of freedom additional to the hair-bundle position would be required to obtain tighter bounds to the rate of entropy production with our method or other estimation techniques [52,[64][65][66][67][68].

Thermodynamic uncertainty relation (TUR) in the ear of the bullfrog
Noisy limit-cycle oscillations in, for instance, a two-dimensional phase space can reveal irreversibility in the form of probability currents. As discussed in section 2.2, the so-called TURs provide lower bounds to the rate of entropy production in terms of the mean and the variance of empirical time-integrated currents (see e.g. [18,19]). Here, we apply the finite-time TUR (8) to predict how much entropy production one can assess by measuring two mesoscopic degrees of freedom: the tip position X 1 of the hair bundle and the transduction current, normalized to its maximum value, P o = I/I max (see equation (16)). Specifically, we analyze two-dimensional stochastic trajectories Γ [0,τ ] ≡ {(X 1 (s), P o (s))} τ s=0 obtained from simulations of equations (12) and (13) in the quiescent (figure 6(A)) and oscillatory region (figure 6(D)) of the state diagram shown in figure 5. These trajectories reveal a larger circulating probability current within the oscillatory region, as expected, but also a smaller relative uncertainty.
To quantify these effects, we make use of the estimate for entropy production (10) based on the finite-time TUR applied to two different currents in the (X 1 , P o ) plane. More precisely, we analyze the fluctuations of the counterclockwise current φ and of the environmental entropy change S env , as described below. First, we map the dynamics into the complex plane z(t) =X 1 (t) + iP o (t) and measure and N φ(t) is the net number  (10), symbols) and the total entropy production rate σ tot (equation (17), lines) as a function of the observation time τ in the quiescent (blue squares, blue dotted line) and oscillatory (red circles, red dashed line) regimes.
(H) Comparison between the total entropy production rate σ tot (black diamonds), the one-variable irreversibility measure σ 1 (green squares), and the two-variable irreversibility measure from the TUR (equation (10)), as a function of the maximum motor force F max . For the latter case, we employ the counterclockwise current φ ( σ φ 2 , orange circles) and the environmental entropy change ( σ S 2 , red circles). In (G) and (H) the lines are a guide to the eye. Simulations were run for a total duration of 300 s at the two operating points with maximum motor force F max = 31 pN (A)-(C) and F max = 62 pN (D)-(F), corresponding to quiescent and oscillatory regimes, respectively. The rest of the simulation parameters were set to the same values as in figure 4. of counterclockwise turns-the winding number. Here,X 1 (t) = X 1 (t) − X 1 ,P o (t) = P o (t) − P o . Using sample trajectories of duration τ = 2 s, we found that the counterclockwise current j φ (t) = φ(t)/τ displays both a larger absolute mean and a larger signal-to-noise ratio, corresponding to more accurate currents, when the system operates in the oscillatory (figure 6(D-E)) rather than in the quiescent regime of the dynamics ( figure 6(A-B)). We show estimatesσ φ 2 for the two case studies in (figure 6(G)). These estimates are obtained using the relative fluctuations of the current, i.e.σ φ 2 is given by equation (10) applied to the current φ. For an example trajectory in the quiescent regime of the dynamics,σ φ 2 ∼1k B /s is of the same order of magnitude asσ 1 (figure 6(G), blue squares). Remarkably, operating in the oscillatory regime instead yields an estimateσ φ 2 ∼ 10 3 k B /s ( figure 6G, red circles), which is three orders of magnitude larger thanσ 1 and only a few fold smaller than σ tot .
To get further insights on entropy production upon varying the operating point in the state diagram of the system, we plotσ φ 2 as a function of the maximal motor force F max which is a control parameter of the Hopf bifurcation at fixed S = 0.94 ( figure 6(H)). In the quiescent region, F max < F c 50 pN,σ φ 2 is not significantly different fromσ 1 , underestimating entropy production (∼1k B /s) by about one order of magnitude below σ tot ∼ 6k B /s. Increasing F max , the two-variable irreversibility measureσ φ 2 and the total entropy production σ tot both exhibit a jump when the system enters the oscillatory region of the dynamics, which is indicative of the underlying Hopf bifurcation, as also observed for other oscillatory systems e.g. in references [69,70]. We note here that there are multiple ways to define a current from measurements of X 1 and I. In principle, one could develop an optimization procedure to find the current that provides the tightest bound to entropy production [54][55][56][57]. Although we did not attempt here to employ such a procedure, it is useful to compareσ φ 2 to the estimateσ S 2 , which results from the environmental entropy current j S (t) = S env (t)/τ and has been being proposed to provide a near-optimal estimate [54,56]. Here, S env (t) is given by equation (9) for the choice G 1 = F 1 (X 1 , X 2 )/T and G 2 = F 2 (X 1 , X 2 )/T eff . In figure 6(H) we show thatσ φ 2 andσ S 2 are similar in the oscillatory region (F max ∼ 60 pN), but thatσ S 2 provides a tighter bound for small (F max < 40 pN) and large (F max > 120 pN) maximal motor force where the hair bundle does not oscillate.

Discussion
In this work, we have shown that fluctuations of active systems can reveal the arrow of time even in the absence of net drifts or currents. The hierarchy of measures of time irreversibility introduced here provides lower bounds for the entropy production of an active process. We have demonstrated the applicability of the approach by estimating entropy production associated with experimental noisy oscillations of a single degree of freedom in the case of mechanosensory hair bundles from the bullfrog's ear. We have shown that quantifications of the arrow of time can efficiently discriminate quiescent and oscillatory hair bundles, as well as reveal transitions between the two regimes in response to changes in a control parameter (e.g. the calcium concentration as in reference [2]). However, using a model of active hair bundle oscillations, we also showed that estimating the rate of entropy production with only one degree of freedom yields a lower bound that can be orders of magnitude smaller than the total entropy production rate in the system. In the case of hair-bundle oscillations, we predict that measuring a second degree of freedom, e.g. the transduction current, would add sufficient information to get a tight bound. With two degrees of freedoms, the current in the phase space and its fluctuations can be used to bound entropy production by means of thermodynamic uncertainty relations. Overall, our results show that irreversibility measures can quantify entropy production in active matter, including living systems, from fluctuations of only a few mesoscopic degrees of freedom.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.
The first step in the derivation is a simple application of the invariance of the KL distance under a one-to-one map: Second, we can rewrite the relative entropy as D P Ξ (ξ 1 , . . . , ξ n ) Q Ξ (ξ 1 , . . . , ξ n ) = dξ 1 . . . dξ n P Ξ (ξ 1 , . . . , ξ n ) ln Because the KLD between two distributions is always positive, equations (A1) and (A2) yield the bound and the inequality saturates if the transformed variables ξ i (i = 1, . . . , n) are also i.i.d. under P Ξ (ξ 1 , . . . , ξ n ), i.e. when P Ξ (ξ 1 , . . . , ξ n ) = Π i p(ξ i ). If one can find a one-to-one map that transforms the original random variables into i.i.d. variables under both distributions P and Q, then (A3) becomes an equality and the exact KLD between the two multivariate distributions P X and Q X can be reduced to the KLD between single variable distributions p(ξ) and q(ξ), which is much easier to evaluate from real data. This is the key idea of our method to estimate the irreversibility of experimental time series.
the linear recursion Inspired by the AR(m) process, we introduce the following one-to-one map We now apply the bound (A3) to the KLD in the right-hand side of equation (B2), using the transformation defined by equation (B4). Since the contribution of the first, possibly correlated, m values of the time series ξ i , vanishes in the limit n → ∞, we obtain the following lower bound to the KLD rate (equation (7) in the main text): We can obtain empirical estimates of p(ξ) and q(ξ) from a single stationary time series (x 1 , . . . , x n ) as follows. We apply the transformation (B4) to both the original time series (x 1 , . . . , x n ) and to its time reversal (x n , . . . , x 1 ) obtaining, respectively, two new time series (ξ F 1 , . . . , ξ F n ) and (ξ R 1 , . . . , ξ R n ), which are stationary at least for i > m. The empirical PDFs obtained from the data of each series are estimations of, respectively, p(ξ) and q(ξ). Note that the same transformation (B4) must be applied to both the original time series (x 1 , . . . , x n ) and its time reverse (x n , . . . , x 1 ), but the inequality (A3) only requires uncorrelated residuals in the reverse series. For this purpose, we calculate the coefficients a 1 , . . . , a m by fitting the reverse time series (x n , . . . , x 1 ) to the AR(m) model in equation (B3).
As indicated in the previous section, the inequality (B5) is tighter when the residuals are uncorrelated in the forward series as well. This is the case of the experimental series that we have analyzed (see, for instance, figure 1(B) in the main text) although, in principle, it is not guaranteed by this procedure. We remark that the inequality (B5) is a rigorous result if the transformation (B4) applied to the reverse time series yields an uncorrelated series (ξ R 1 , . . . , ξ R n ). In that case, k B f s D[p(ξ) q(ξ)] is an estimate of σ 1 with at least two sources of error: (i) the discrete sampling of the process X(t) and (ii) the remnant correlation time in the residuals (ξ F 1 , . . . , ξ F n ) obtained from the forward time series. To summarize, our theory provides an estimateσ 1 for the KLD rate σ 1 which can be evaluated as follows: (a) Estimate the coefficients, a 1 , . . . , a m , by fitting the time-reversed series (x n , . . . , x 1 ) to aAR(m) model of order m > 1. A reasonable choice is m = 10, but it should be tuned to minimize the correlation time in the residuals (ξ R 1 , . . . , ξ R n ). (b) Apply the whitening transformation (B4) to the original series (x 1 , . . . , x n ) and to its time reversal (x n , . . . , x 1 ) to obtain, respectively, new time series (ξ F 1 , . . . , ξ F n ) and (ξ R 1 , . . . , ξ R n ). Note that the new processes are not each other's time reversal. (c) Obtain the empirical distributions p(ξ) and q(ξ) from the time series (ξ F 1 , . . . , ξ F n ) and (ξ R 1 , . . . , ξ R n ), respectively. (d) Calculate the KLD between p(ξ) and q(ξ) which can be estimated from numerical integration of the right-hand side in (B6) using the empirical normalized histograms p(ξ) and q(ξ). We call this estimateD, which is given bŷ are the empirical probabilities, obtained from the number of times n F i and n R i that the sequences (ξ F 1 , . . . , ξ F n ) and (ξ R 1 , . . . , ξ R n ) lie in the ith bin, respectively. The sum in (B7) runs over all bins for which n F i > 0 and n R i > 0. For simplicity, we used 100 bins of equal spacing ranging from the minimum to the maximum values of the residual time series (ξ F 1 , . . . , ξ F n ). The value of the estimateD of the KLD (B6) is weighted by a prefactor γ 1 defined in terms of the probability to reject the null hypothesis p(ξ) = q(ξ). We use this procedure to correct the statistical bias in the estimation of the KLD that appears when two stochastic processes have similar statistics [16,43]. For this purpose, we use the KS statistical test under the null hypothesis H 0 : p(ξ) = q(ξ) which yields a p-value p KS for the two distributions to be equal. Here, small p KS means that there is stronger statistical evidence in favor of the alternative hypothesis p(ξ) = q(ξ), thus γ = 1 − p KS serves as a weight of irreversibility: γ 0 when it is hard to reject H 0 (reversibility) and γ 1 there is a larger statistical evidence to reject H 0 . (e) Finally, our estimate ofσ 1 is thus given by the KLD estimateD times the Boltzmann constant and the data sampling frequency:σ

Appendix D. Biophysics of mechanosensory hair bundles
Details of the experimental procedure have been published elsewhere [2]. In short, an excised preparation of the bullfrog's (Rana catesbeiana) sacculus was mounted on a two-compartment chamber to reproduce the ionic environment of the inner ear. This organ is devoted to sensitive detection of low-frequency vibrations (5-150 Hz) of the animal's head in a vertical plane; it contains about 3000 sensory hair cells that are arranged in a planar epithelium. The basal bodies of hair cells were bathed in a standard saline solution and the hair bundles projected in an artificial endolymph. The preparation was viewed through a ×60 water-immersion objective of an upright microscope. Under these conditions, spontaneous hair-bundle oscillations were routinely observed. The oscillations could be recorded by imaging, at a magnification of ×1000, the top of the longest stereociliary row onto a displacement monitor that included a dual photodiode. Calibration was performed by measuring the output voltages of this photometric system in response to a series of offset displacements. Here, we analyzed 182 spontaneously oscillating hair bundles from data previously published [59]. Spontaneous hair-bundle oscillations were described by a published model of active hair-bundle motility [2] that rests on a necessary condition of negative hair-bundle stiffness, on the presence of molecular motors that actively pull on the tip links, and on feedback by the calcium component of the transduction current. Hair-bundle deflections affect tension in tip links that interconnect neighboring stereocilia of the bundle. Changes in tip-link tension in turn modulate the open probability of mechano-sensitive ion channels connected to these links. Importantly, the relation between channel gating and tip-link tension is reciprocal: gating of the transduction channels affects tip-link tension. Consequently, channel gating effectively reduces the stiffness of a hair bundle, a phenomenon appropriately termed 'gating compliance', which can result in negative stiffness if channel-gating forces are strong enough. Active hair-bundle movements result from the activity of the adaptation motors. By controlling tip-link tension, adaptation motors regulate the open probability of the mechanosensitive channels. The force produced by the motors is in turn regulated by the Ca 2+ component of the transduction current which thus provides negative feedback on the motor force [2]. When the fixed point of this dynamical system corresponds to an unstable position of negative stiffness, the system oscillates spontaneously. The maximal force exerted by the motors F max and the calcium feedback strength S are control parameters of the system and fully determine its dynamics (oscillatory, quiescent, bi-stable) [25].

Appendix E. Quantification of entropy production in numerical simulations of hair bundle oscillations
In this section, we provide numerical results for the stochastic model of the ear hair bundle given by equations (12) and (13) in the main text. The steady-state entropy production rate of the model is given by where F 1 = F 1 (X 1 , X 2 ), F 2 = F 2 (X 1 , X 2 ) and • denotes the Stratonovich product. Using the definitions of the forces in equations (18) and (19) one obtains after some algebra equation (17) in the main text. In all our numerical simulations, we estimate the steady-state averages of the type for a generic force F(t) = F(X(t), Y(t)) from a single stationary trajectory of total duration τ = 300 s and sampling time Δt = 1 ms as follows: where t i = iΔt and n = τ /Δt.

Appendix F. Estimation of the quality factor of stochastic oscillations
We estimate the quality factor Q of spontaneous hair-bundle oscillations from numerical simulations of the hair-bundle stochastic model given by equations (12) and (13) in the main text. For this purpose, we generate a single numerical simulation of duration τ = 300 s. We then partition the simulation into ten consecutive traces of duration T = τ/10 = 30 s. For each of these traces {X α (t)} (α = 1, . . . , 10) we compute the power spectral density as C α (f ) = (1/T) T 0 X α (s)e 2πift dt 2 . We then calculate the average of the power spectral density over the ten different tracesC(f ) = (1/10) 10 α=1 C α (f ) and fit the estimateC(f ) as a function of f to the sum of two Lorentzian functions [4,22,59] where Q is the quality factor, f o is the oscillation frequency and A > 0 is an amplitude parameter. Figure 10 shows examples of numerical simulations for which we apply this procedure to determine the value of the quality factor by extracting the value Q from the fit of the data to equation (F1). Notably, equation (F1) reproduces power spectra of hair-bundle simulations for oscillations with values Q that are in a wide range (figure 10).