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

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–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–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–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–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–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–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.

Zoom In Zoom Out Reset image size

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.

2.1. 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 ${\mathbf{\Gamma }}_{[0,t]}\equiv {\left\{({x}_{1}(s),{x}_{2}(s),\dots \enspace )\left.\right)\right\}}_{s=0}^{t}$ a given trajectory described by the system variables and its corresponding time-reversed trajectory as ${~{\mathbf{\Gamma }}}_{[0,t]}\equiv {\left\{({\theta }_{1}{x}_{1}(t-s),{\theta }_{2}{x}_{2}(t-s),\dots \enspace )\right\}}_{s=0}^{t}$, where θ_α = ±1 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

$$\begin{equation}{\sigma }_{\text{tot}}={k}_{\mathrm{B}}\enspace \underset{t\to \infty }{\mathrm{lim}}\enspace \frac{1}{t}D\left[\mathcal{P}\left({\mathbf{\Gamma }}_{[0,t]}\right)\right\vert \left\vert \mathcal{P}\left({~{\mathbf{\Gamma }}}_{[0,t]}\right.\right)\left.\right],\end{equation} \tag{ 1 }$$

where k_B is the Boltzmann constant and $\mathcal{P}$ denotes the steady-state path probability [31–34]. Here $D[\mathcal{Q}{\Vert}\mathcal{R}]{\geqslant}0$ is the KLD between the probability measures $\mathcal{Q}$ and $\mathcal{R}$, which quantifies the distinguishability between these two distributions. For measures of a single random variable x the KLD is given by $D[\mathcal{Q}(x){\Vert}\mathcal{R}(x)]\equiv \int \mathrm{d}x\enspace \mathcal{Q}(x)\mathrm{ln}[\mathcal{Q}(x)/\mathcal{R}(x)]$. 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₁, ..., X_k are known. We define the k-variable irreversibility measure in terms of path probabilities of k mesoscopic variables

$$\begin{equation}{\sigma }_{k}\equiv {k}_{\mathrm{B}}\enspace \underset{t\to \infty }{\mathrm{lim}}\enspace \frac{1}{t}D\left[\mathcal{P}\left({\mathbf{\Gamma }}_{[0,t]}^{(k)}\right)\right\vert \left\vert \mathcal{P}\left({~{\mathbf{\Gamma }}}_{[0,t]}^{(k)}\right)\right],\end{equation} \tag{ 2 }$$

where ${\mathbf{\Gamma }}_{[0,t]}^{(k)}\equiv {\left\{({x}_{1}(s),\dots ,{x}_{k}(s))\right\}}_{s=0}^{t}$ and ${~{\mathbf{\Gamma }}}_{[0,t]}^{(k)}\equiv {\left\{({\theta }_{1}{x}_{1}(t-s),\dots ,{\theta }_{k}{x}_{k}(t-s))\right\}}_{s=0}^{t}$ 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:

$$\begin{equation}0{\leqslant}{\sigma }_{1}{\leqslant}\dots {\leqslant}{\sigma }_{k}{\leqslant}{\sigma }_{k+1}{\leqslant}\dots {\leqslant}{\sigma }_{\text{tot}}.\end{equation} \tag{ 3 }$$

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–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 σ₁ 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 σ₁ 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 σ₁ is a herculean task due to the difficulties in sampling the whole phase space of paths [38–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₁, ..., x_n ), i = 1, ..., n, that transforms the original time series and its time reversal into two new time series ${\xi }_{i}^{\text{F}}={\xi }_{i}({x}_{1},\dots ,{x}_{n})$ and ${\xi }_{i}^{\text{R}}={\xi }_{i}({x}_{n},\dots ,{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 ${\xi }_{i}^{\text{F}}$ and ${\xi }_{i}^{\text{R}}$, 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₁, ..., x_n ) and the backward (x_n , ..., x₁) time series. In that case, the removal of the correlations in the backward series is enough to provide a lower bound for σ₁:

$$\begin{equation}{\sigma }_{1}{\geqslant}{k}_{\mathrm{B}}{f}_{\mathrm{s}}D[p(\xi ){\Vert}q(\xi )]\equiv {\hat{\sigma }}_{1},\end{equation} \tag{ 4 }$$

where f_s = (Δt)⁻¹ 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(\xi ){\Vert}q(\xi )]\simeq \gamma {\sum }_{i}\;{\hat{p}}_{i}\enspace \mathrm{ln}({\hat{p}}_{i}/{\hat{q}}_{i})$ where $\hat{p},\hat{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.

2.2. Estimates based on TURs

We now consider the case in which one can measure a set of k ⩾ 2 mesoscopic variables X₁, ..., 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₁, ..., 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–46]. For example, consider an isothermal overdamped Langevin system described by two mesoscopic variables,

$$\begin{equation}{\dot {X}}_{1}={\mu }_{1}{F}_{1}({X}_{1},{X}_{2})+\sqrt{2{D}_{1}}{\eta }_{1}\end{equation} \tag{ 5 }$$

$$\begin{equation}{\dot {X}}_{2}={\mu }_{2}{F}_{2}({X}_{1},{X}_{2})+\sqrt{2{D}_{2}}{\eta }_{2},\end{equation} \tag{ 6 }$$

where μ_k , D_k are the mobility and diffusion coefficient of the kth variable, F_k (X₁, X₂) is a generic state-dependent force, and ${\eta }_{1}$ and ${\eta }_{2}$ are two independent Gaussian white noises with zero mean $\langle {\eta }_{1}(t)\rangle =\langle {\eta }_{2}(t)\rangle =0$ and autocorrelation $\langle {\eta }_{1}(t){\eta }_{1}(t{{\prime}})\rangle =\langle {\eta }_{2}(t){\eta }_{2}(t{{\prime}})\rangle =\delta (t-t{{\prime}})$ with $\delta$ Dirac's delta function. For this example, the entropy production can be expressed in terms of forces F_k and currents (velocities) ${\dot {X}}_{k}$ as follows

$$\begin{equation}{\sigma }_{\text{tot}}={\sigma }_{2}=\frac{{k}_{\mathrm{B}}{\mu }_{1}}{{D}_{1}}\left\langle {F}_{1}{\circ}{\dot {X}}_{1}\right\rangle +\frac{{k}_{\mathrm{B}}{\mu }_{2}}{{D}_{2}}\left\langle {F}_{2}{\circ}{\dot {X}}_{2}\right\rangle .\end{equation} \tag{ 7 }$$

Here, ⟨⋅⟩ denote steady state averages, ◦ the Stratonovich product [47, 48] and we have used the shorthand notation F_k ≡ F_k (X₁, X₂). 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–51]

$$\begin{equation}{\sigma }_{\text{tot}}{\geqslant}\frac{2{k}_{\mathrm{B}}}{\tau }\frac{{\langle j(\tau )\rangle }^{2}}{\mathrm{Var}[j(\tau )]},\end{equation} \tag{ 8 }$$

where Var[j(τ)] = ⟨j²(τ)⟩ − ⟨j(τ)⟩² 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

$$\begin{equation}{j}^{G}(t)={\int }_{0}^{t}{G}_{1}(s){\circ}\dot {{X}_{1}}(s)+{\int }_{0}^{t}{G}_{2}(s){\circ}\dot {{X}_{2}}(s),\end{equation} \tag{ 9 }$$

where G₁(t) ≡ G₁(X₁(t), X₂(t)) and G₂(t) ≡ G₂(X₁(t), X₂(t)), with G₁(x, y) and G₂(x, y) two arbitrary functions. For this model, the TUR (8) implies that [52]

$$\begin{equation}{\sigma }_{2}{\geqslant}\frac{2{k}_{\mathrm{B}}}{\tau }\frac{{\langle {j}^{G}(\tau )\rangle }^{2}}{\mathrm{Var}[{j}^{G}(\tau )]}\equiv {\hat{\sigma }}_{2},\end{equation} \tag{ 10 }$$

which holds for any observation time τ and any choice of the functions G₁ and G₂ that enter in equation (9). Notably, the bound (10) also applies for currents that are obtained by coarse graining the phase space variables X₁ and X₂. To which extent the bound (8) can be tightened by finding optimal currents has become an active area of research (see e.g. [53–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].

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₁ 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₁, we can only estimate σ₁, 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 ${\xi }_{i}^{\text{F}}$ (${\xi }_{i}^{\text{R}}$) as the difference between the observed values of the forward (backward) time series and the forecast of that value based on an AR model of order m = 10. Parameters of the model are determined from fits of the AR-model to the time-reversed series of positions. The residual time series ${\xi }_{i}^{\text{F}}$ and ${\xi }_{i}^{\text{R}}$ (figures 1(B)–(D), center column) obtained from the whitening transformation are uncorrelated (figures 1(B)–(D), right column) and are therefore i.i.d. processes in good approximation.

We find that the irreversibility measure ${\hat{\sigma }}_{1}$ given by equation (4) distinguishes active hair-bundle fluctuations (${\hat{\sigma }}_{1}{ >}0$) from passive fluctuations of a bistable system (${\hat{\sigma }}_{1}\simeq 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 ${\hat{\sigma }}_{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].

Zoom In Zoom Out Reset image size

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

$$\begin{equation}{\hat{s}}_{1}(\xi )\equiv {k}_{\mathrm{B}}{f}_{\mathrm{s}}\left[p(\xi )\mathrm{ln}\enspace \frac{p(\xi )}{q(\xi )}+q(\xi )-p(\xi )\right],\end{equation} \tag{ 11 }$$

which obeys ${\hat{s}}_{1}(\xi ){\geqslant}0$ for all ξ [60], and ${\hat{\sigma }}_{1}=\int \mathrm{d}\xi \enspace {\hat{s}}_{1}(\xi )$. We find that for all the analyzed values of ξ, the local irreversibility of active oscillations is ∼10³ times larger than for passive hair-bundle fluctuations and experimental noise.

Zoom In Zoom Out Reset image size

We now relate the estimate ${\hat{\sigma }}_{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₁ and of the center of mass of a collection of molecular motors X₂ [2, 59, 61] (see appendix D):

$$\begin{equation}{\lambda }_{1}{\dot {X}}_{1}=-\frac{\partial V}{\partial {X}_{1}}+\sqrt{2{k}_{\mathrm{B}}T{\lambda }_{1}}{\xi }_{1}\end{equation} \tag{ 12 }$$

$$\begin{equation}{\lambda }_{2}{\dot {X}}_{2}=-\frac{\partial V}{\partial {X}_{2}}-{F}_{\text{act}}+\sqrt{2{k}_{\mathrm{B}}{T}_{\text{eff}}{\lambda }_{2}}{\xi }_{2}.\quad \end{equation} \tag{ 13 }$$

Here, λ₁ and λ₂ are friction coefficients and ξ₁ and ξ₂ 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

$$\begin{equation}V({X}_{1},{X}_{2})=\frac{{k}_{\text{gs}}({X}_{1}-{X}_{2}{)}^{2}+{k}_{\text{sp}}{X}_{1}^{2}}{2}-N{k}_{\mathrm{B}}T\enspace \mathrm{ln}\left[\mathrm{exp}\left(\frac{{k}_{\text{gs}}D({X}_{1}-{X}_{2})}{N{k}_{\mathrm{B}}T}\right)+A\right],\end{equation} \tag{ 14 }$$

where k_gs and k_sp are stiffness coefficients; D is the gating swing of a transduction channel; and A = exp[(ΔG + (k_gs D²)/2N)/(k_B T)], ΔG being the energy difference between open and closed states of the channels and N the number of transduction elements. The force

$$\begin{equation}{F}_{\text{act}}({X}_{1},{X}_{2})={F}_{\mathrm{max}}(1-S{P}_{\mathrm{o}}({X}_{1},{X}_{2})),\end{equation} \tag{ 15 }$$

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

$$\begin{equation}{P}_{\mathrm{o}}({X}_{1},{X}_{2})=\frac{1}{1+A\enspace \mathrm{exp}(-{k}_{\text{gs}}D({X}_{1}-{X}_{2})/N{k}_{\mathrm{B}}T)},\end{equation} \tag{ 16 }$$

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₁ − X₂. 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₁ 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₁(t) or X₂(t) do not reveal obvious signs of a net current, which here would correspond to a drift. However, trajectories in the (X₁, X₂) 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 ${\hat{\sigma }}_{1}$ to the total entropy production σ_tot.

Zoom In Zoom Out Reset image size

In the stochastic model of hair-bundle oscillations given by equations (12) and (13) we deal with only two variables, therefore σ_tot = σ₂. From the analytical expression of σ₂, we find that the steady-state entropy production rate can be written as (see references [46, 62, 63] and appendix E)

$$\begin{equation}{\sigma }_{\text{tot}}=\frac{1}{T}\left\langle {F}_{1}{\circ}\frac{\mathrm{d}{X}_{1}}{\mathrm{d}t}\right\rangle +\frac{1}{{T}_{\text{eff}}}\left\langle {F}_{2}{\circ}\frac{\mathrm{d}{X}_{2}}{\mathrm{d}t}\right\rangle =-\langle {\dot {Q}}_{1}\rangle \left(\frac{1}{T}-\frac{1}{{T}_{\text{eff}}}\right)+\frac{\langle {\dot {W}}_{\text{act}}\rangle }{{T}_{\text{eff}}},\end{equation} \tag{ 17 }$$

where

$$\begin{equation}{F}_{1}({X}_{1},{X}_{2})=-\frac{\partial V({X}_{1},{X}_{2})}{\partial {X}_{1}}\end{equation} \tag{ 18 }$$

$$\begin{equation}{F}_{2}({X}_{1},{X}_{2})=-\frac{\partial V({X}_{1},{X}_{2})}{\partial {X}_{2}}-{F}_{\text{act}}({X}_{1},{X}_{2}),\end{equation} \tag{ 19 }$$

are the conjugated forces to X₁ and X₂, respectively. We have also introduced the steady-state average heat $-\langle {\dot {Q}}_{1}\rangle =-\langle (\partial V/\partial {X}_{1}){\circ}{\dot {X}}_{1}\rangle$ dissipated to the thermal bath at temperature $T$ and the power $\langle {\dot {W}}_{\text{act}}\rangle =-\langle {F}_{\text{act}}{\circ}{\dot {X}}_{2}\rangle$ 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 ${\hat{\sigma }}_{1}$ for trajectories X₁(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, ${\hat{\sigma }}_{1}$ provides a lower bound to the actual dissipation rate. Actually, the rate of entropy production estimated from ${\hat{\sigma }}_{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–68].

Zoom In Zoom Out Reset image size

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₁ 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 ${\mathbf{\Gamma }}_{[0,\tau ]}\equiv {\left\{({X}_{1}(s),{P}_{\mathrm{o}}(s))\right\}}_{s=0}^{\tau }$ 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.

Zoom In Zoom Out Reset image size

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₁, 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)={\overline{X}}_{1}(t)+\mathrm{i}{\overline{P}}_{\mathrm{o}}(t)$ and measure θ(t) = ϕ(t) + 2πN_ϕ(t), where $\phi (t)={\mathrm{tan}}^{-1}({\overline{P}}_{\mathrm{o}}(t)/{\overline{X}}_{1}(t))\in [0,2\pi ]$ is the phase and N_ϕ(t) is the net number of counterclockwise turns—the winding number. Here, ${\overline{X}}_{1}(t)={X}_{1}(t)-\langle {X}_{1}\rangle$, ${\overline{P}}_{\mathrm{o}}(t)={P}_{\mathrm{o}}(t)-\langle {P}_{\mathrm{o}}\rangle$. 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 ${\hat{\sigma }}_{2}^{\phi }$ for the two case studies in (figure 6(G)). These estimates are obtained using the relative fluctuations of the current, i.e. ${\hat{\sigma }}_{2}^{\phi }$ is given by equation (10) applied to the current ϕ. For an example trajectory in the quiescent regime of the dynamics, ${\hat{\sigma }}_{2}^{\phi }\sim 1{\mathit{\text{k}}}_{\mathrm{B}}/$ s is of the same order of magnitude as ${\hat{\sigma }}_{1}$ (figure 6(G), blue squares). Remarkably, operating in the oscillatory regime instead yields an estimate ${\hat{\sigma }}_{2}^{\phi }\sim 1{0}^{3}{k}_{\mathrm{B}}/$ s (figure 6G, red circles), which is three orders of magnitude larger than ${\hat{\sigma }}_{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 ${\hat{\sigma }}_{2}^{\phi }$ 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, ${\hat{\sigma }}_{2}^{\phi }$ is not significantly different from ${\hat{\sigma }}_{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 ${\hat{\sigma }}_{2}^{\phi }$ 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₁ and I. In principle, one could develop an optimization procedure to find the current that provides the tightest bound to entropy production [54–57]. Although we did not attempt here to employ such a procedure, it is useful to compare ${\hat{\sigma }}_{2}^{\phi }$ to the estimate ${\hat{\sigma }}_{2}^{S}$, 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₁ = F₁(X₁, X₂)/T and G₂ = F₂(X₁, X₂)/T_eff. In figure 6(H) we show that ${\hat{\sigma }}_{2}^{\phi }$ and ${\hat{\sigma }}_{2}^{S}$ are similar in the oscillatory region (F_max ∼ 60 pN), but that ${\hat{\sigma }}_{2}^{S}$ 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.

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.

PM acknowledges support by the French National Agency for Research (ANR-11-BSV5 0011 and ANR-16-CE13-0015) and by the Labex Cell(n)Scale (ANR-11-LABX-0038, ANR-10-IDEX-0001-02). JMRP acknowledges financial support from the Spanish Government (Grant Contract, FIS-2017-83706-R). The electron micrograph of the hair bundles shown in figure 1(A) was obtained by Atitheb Chaiyasitdhi, a PhD student in P Martin's group. We acknowledge stimulating discussions with Roman Belousov, Izaak Neri, Andre C Barato, Simone Pigolotti, Johannes Baumgart, Jose Negrete Jr, Ken Sekimoto, Ignacio A Martínez, Patrick Pietzonka and A J Hudspeth.

The data that support the findings of this study are available upon reasonable request from the authors.

Here we prove a general lower bound for the KLD between two multivariate probability densities P_X (x₁, ..., x_n ) and Q_X (x₁, ..., x_n ) that fulfill the following: there exits a one-to-one map ξ_i = ξ_i (x₁, ..., x_n ) with i = 1, ..., n, such that

• (a)  
  The transformed variables ξ_i are identically distributed under both P and Q, that is, the distributions P_Ξ(ξ₁, ..., ξ_n ) and Q_Ξ(ξ₁, ..., ξ_n ) have, respectively, identical marginal distributions p(ξ) and q(ξ) for any ξ_i (i = 1, ..., n);
• (b)  
  The transformed variables ξ_i are i.i.d. under the distribution Q, that is, Q_Ξ(ξ₁, ..., ξ_n ) = Π_i q(ξ).

The first step in the derivation is a simple application of the invariance of the KL distance under a one-to-one map:

$$\begin{equation}D\left[{P}_{X}({x}_{1},\dots ,{x}_{n}){\Vert}{Q}_{X}({x}_{1},\dots ,{x}_{n})\right]=D\left[{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n}){\Vert}{Q}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n})\right].\end{equation} \tag{ A1 }$$

Second, we can rewrite the relative entropy as

$$\begin{align}D\left[{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n}){\Vert}{Q}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n})\right]& =\int \mathrm{d}{\xi }_{1}\dots \int \mathrm{d}{\xi }_{n}{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n})\mathrm{ln}\enspace \frac{{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n})}{{{\Pi}}_{i}\enspace q({\xi }_{i})}\\ & =\int \mathrm{d}{\xi }_{1}\dots \int \mathrm{d}{\xi }_{n}\left[{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n})\mathrm{ln}\enspace \frac{{{\Pi}}_{i}\enspace p({\xi }_{i})}{{{\Pi}}_{i}\enspace q({\xi }_{i})}\right.\\ & \quad \left.+{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n})\mathrm{ln}\enspace \frac{{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n})}{{{\Pi}}_{i}p({\xi }_{i})}\right]\\ & =nD[p(\xi ){\Vert}q(\xi )]+D\left[{P}_{{\Xi}}({\xi }_{1},\dots ,{\xi }_{n}){\Vert}{{\Pi}}_{i}p({\xi }_{i})\right].\end{align} \tag{ A2 }$$

Because the KLD between two distributions is always positive, equations (A1) and (A2) yield the bound

$$\begin{equation}D\left[{P}_{X}({x}_{1},\dots ,{x}_{n}){\Vert}{Q}_{X}({x}_{1},\dots ,{x}_{n})\right]{\geqslant}nD[p(\xi ){\Vert}q(\xi )],\end{equation} \tag{ A3 }$$

and the inequality saturates if the transformed variables ξ_i (i = 1, ..., n) are also i.i.d. under P_Ξ(ξ₁, ..., ξ_n ), i.e. when P_Ξ(ξ₁, ..., ξ_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 estimation of the KLD rate from single stationary trajectories of both discrete and continuous random variables have been previously discussed [40]. For continuous random variables, the most common strategy is to make a symbolization or discretization of the time series [38]. Then, the KLD is estimated from the statistics of substrings of increasing length [15, 16]. The main limitation of this method is that one easily reaches lack of statistics even for short substrings. If the observed time series is non-Markovian, this limitation could yield inaccurate bounds for the entropy production. For instance, the KLD between two data substrings can be zero in non-equilibrium stationary states without observable currents [15, 16, 40].

Here we introduce a new method to estimate the KLD rate

$$\begin{equation}\frac{{\sigma }_{1}}{{k}_{\mathrm{B}}}\equiv \underset{t\to \infty }{\mathrm{lim}}\enspace \frac{1}{t}D\left[\mathcal{P}\left({\left\{x(s)\right\}}_{s=0}^{t}\right){\Vert}\mathcal{P}\left({\left\{x(t-s)\right\}}_{s=0}^{t}\right)\right],\end{equation} \tag{ B1 }$$

that is valid for continuous and possibly non-Markovian stochastic processes X(t). First, in practice one has access to discrete-time observations of the process x_i ≡ X(iΔt), i = 1, ..., n, i.e. a time series containing n = t/Δt consecutive samples of the process with sampling rate f_s = 1/Δt. The time discretization implies a loss of information yielding a lower bound to the KLD rate:

$$\begin{equation}\frac{{\sigma }_{1}}{{k}_{\mathrm{B}}}{\geqslant}{f}_{\mathrm{s}}\enspace \underset{n\to \infty }{\mathrm{lim}}\enspace \frac{1}{n}D[{P}_{X}({x}_{1},\dots ,{x}_{n}){\Vert}{Q}_{X}({x}_{1},\dots ,{x}_{n})],\end{equation} \tag{ B2 }$$

where Q_X (x₁, ..., x_n ) = P_X (x_n , ..., x₁) is the probability to observe the reverse trajectory (x_n , ..., x₁).

We can now apply the inequality (A3) to the right-hand side in equation (B2). To do that, it is necessary to find a one-to-one map ξ_i = ξ_i (x₁, ..., x_n ) that transforms the reverse time series (x_n , ..., x₁) into a sequence of n i.i.d. random variables, that is, into a white noise. Such a transformation is usually termed whitening transformation.

An example of whitening transformation is the time series formed by the residuals of an autoregressive (AR) model, which is the transformation that we will use along this paper. A discrete-time stochastic process Y_i is called AR of order m, AR(m), when its value at a given time is given by a linear combination of its m previous values plus a noise term. Such process is univocally determined by m ⩾ 1 real coefficients, a₁, a₂, ..., a_m , a discrete-time white noise η_i and a set of initial values Y₁, Y₂, ..., Y_m . The values of Y_i for i > m are given by the linear recursion

$$\begin{equation}{Y}_{i}=\sum\limits _{j=1}^{m}{a}_{j}{Y}_{i-j}+{\eta }_{i}.\end{equation} \tag{ B3 }$$

Inspired by the AR(m) process, we introduce the following one-to-one map

$$\begin{equation}{\xi }_{i}=\begin{cases}_{i}\quad \text{if}\enspace \mathit{\text{i}}{\leqslant}\mathit{\text{m}}\quad \\ {x}_{i}-\sum\limits _{j=1}^{m}{a}_{j}{x}_{i-j}\quad \text{if}\enspace \mathit{\text{i}}{ >}\mathit{\text{m}}\quad \end{cases},\end{equation} \tag{ B4 }$$

which is a linear transformation defined by a unitriangular matrix with Jacobian equal to one. With an appropriate choice of the coefficients a_j , one can get a new process (ξ₁, ..., ξ_n ) which is approximately i.i.d.. A good choice is given by a maximum likelihood fit of the process to the AR(m) model. In that case, the elements ξ_i in this new time series for i > m are usually called residuals of the original time series (x₁, ..., x_n ) with respect to the AR(m) model. Notice also that, if (x₁, ..., x_n ) is indeed a realization of the stochastic process (B3), then the residuals are i.i.d. random variables and the process (ξ_m+1, ..., ξ_n ) has correlations ⟨ξ_i ξ_j ⟩ = δ_ij for all i, j > m.

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):

$$\begin{equation}\frac{{\sigma }_{1}}{{k}_{\mathrm{B}}}{\geqslant}{f}_{\mathrm{s}}D[p(\xi ){\Vert}q(\xi )].\end{equation} \tag{ B5 }$$

We can obtain empirical estimates of p(ξ) and q(ξ) from a single stationary time series (x₁, ..., x_n ) as follows. We apply the transformation (B4) to both the original time series (x₁, ..., x_n ) and to its time reversal (x_n , ..., x₁) obtaining, respectively, two new time series $({\xi }_{1}^{\mathrm{F}},\dots ,{\xi }_{n}^{\mathrm{F}})$ and $({\xi }_{1}^{\mathrm{R}},\dots ,{\xi }_{n}^{\mathrm{R}})$, 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₁, ..., x_n ) and its time reverse (x_n , ..., x₁), but the inequality (A3) only requires uncorrelated residuals in the reverse series. For this purpose, we calculate the coefficients a₁, ..., a_m by fitting the reverse time series (x_n , ..., x₁) 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 $({\xi }_{1}^{\mathrm{R}},\dots ,{\xi }_{n}^{\mathrm{R}})$. In that case, k_B f_s D[p(ξ)||q(ξ)] is an estimate of σ₁ 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 $({\xi }_{1}^{\mathrm{F}},\dots ,{\xi }_{n}^{\mathrm{F}})$ obtained from the forward time series.

To summarize, our theory provides an estimate ${\hat{\sigma }}_{1}$ for the KLD rate σ₁ which can be evaluated as follows:

• (a)  
  Estimate the coefficients, a₁, ..., a_m , by fitting the time-reversed series (x_n , ..., x₁) 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 $({\xi }_{1}^{\mathrm{R}},\dots ,{\xi }_{n}^{\mathrm{R}})$.
• (b)  
  Apply the whitening transformation (B4) to the original series (x₁, ..., x_n ) and to its time reversal (x_n , ..., x₁) to obtain, respectively, new time series $({\xi }_{1}^{\mathrm{F}},\dots ,{\xi }_{n}^{\mathrm{F}})$ and $({\xi }_{1}^{\mathrm{R}},\dots ,{\xi }_{n}^{\mathrm{R}})$. Note that the new processes are not each other's time reversal.
• (c)  
  Obtain the empirical distributions p(ξ) and q(ξ) from the time series $({\xi }_{1}^{\mathrm{F}},\dots ,{\xi }_{n}^{\mathrm{F}})$ and $({\xi }_{1}^{\mathrm{R}},\dots ,{\xi }_{n}^{\mathrm{R}})$, respectively.
• (d)  
  Calculate the KLD between p(ξ) and q(ξ)

  $$\begin{equation}D[p(\xi ){\Vert}q(\xi )]=\int \mathrm{d}\xi \enspace p(\xi )\mathrm{ln}\enspace \frac{p(\xi )}{q(\xi )},\end{equation} \tag{ B6 }$$

  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 estimate $\hat{\mathcal{D}}$, which is given by

  $$\begin{equation}\hat{\mathcal{D}}=\gamma \enspace \sum\limits _{i}{\hat{p}}_{i}\enspace \mathrm{ln}\enspace \frac{{\hat{p}}_{i}}{{\hat{q}}_{i}},\end{equation} \tag{ B7 }$$

  where ${\hat{p}}_{i}={n}_{i}^{\mathrm{F}}/({\sum }_{i}{n}_{i}^{\mathrm{F}})$ and ${\hat{q}}_{i}={n}_{i}^{\mathrm{R}}/({\sum }_{i}{n}_{i}^{\mathrm{R}})$ are the empirical probabilities, obtained from the number of times ${n}_{i}^{\mathrm{F}}$ and ${n}_{i}^{\mathrm{R}}$ that the sequences $({\xi }_{1}^{\mathrm{F}},\dots ,{\xi }_{n}^{\mathrm{F}})$ and $({\xi }_{1}^{\mathrm{R}},\dots ,{\xi }_{n}^{\mathrm{R}})$ lie in the ith bin, respectively. The sum in (B7) runs over all bins for which ${n}_{i}^{\mathrm{F}}{ >}0$ and ${n}_{i}^{\mathrm{R}}{ >}0$. For simplicity, we used 100 bins of equal spacing ranging from the minimum to the maximum values of the residual time series $({\xi }_{1}^{\mathrm{F}},\dots ,{\xi }_{n}^{\mathrm{F}})$.The value of the estimate $\hat{\mathcal{D}}$ 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₀: 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₀ (reversibility) and γ ≃ 1 there is a larger statistical evidence to reject H₀.
• (e)  
  Finally, our estimate of ${\hat{\sigma }}_{1}$ is thus given by the KLD estimate $\hat{\mathcal{D}}$ times the Boltzmann constant and the data sampling frequency:

  $$\begin{equation}\hat{{\sigma }_{1}}={k}_{\mathrm{B}}{f}_{\mathrm{s}}\hat{\mathcal{D}}.\end{equation} \tag{ B8 }$$

In this section, we analyze the dependency of our irreversibility measure on the sampling frequency f_s. For this purpose, we evaluate ${\hat{\sigma }}_{1}$ defined in equation (B8) for 30 s recordings of the 182 cells that showed spontaneous oscillations at different sampling frequencies, ranging from 125 Hz to 2500 Hz (the latter corresponding to the data shown in figure 2 (B) in the main text). Figure 7 shows that the distribution of the irreversibility measure depends strongly on the sampling frequency of the data. Notably, the distributions shift toward higher irreversibility when the sampling frequency is reduced, until there is too much filtering f_s < 250 Hz such that oscillations cannot be distinguished clearly.

Zoom In Zoom Out Reset image size

In figure 8 we report the distributions of the parameter γ for the 182 recordings of spontaneous oscillations at different sampling frequencies. For all the analyzed cases, the distributions are right-sweked toward values of γ close to 1. Interestingly, the number of cells that display large γ (i.e. KS p-value p_KS < 0.05) attains its maximum in an intermediate frequency range ∼200–600 Hz. To gain further insight on this result, we plot in figure 9 box plots of the distributions of γ and of the corresponding irreversibility estimate ${\hat{\sigma }}_{1}$ as a function of the data sampling frequency. Notably, the median of γ is above 0.95 for intermediate sampling frequencies ranging from 208 to 625 Hz (figure 9(A)). For values of f_s of this frequency band 208–625 Hz, the median of γ is above 0.95 indicating that more than half of the cells display irreversibility with 'significant' KS p-value p_KS < 0.05 at those frequencies. Within this band, the median of the irreversibility measure ${\hat{\sigma }}_{1}$ decreases monotonically with f_s from 3.5k_B/s (f_s = 208 Hz) to 2.6k_B/s (f_s = 625 Hz).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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²⁺ 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].

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

$$\begin{equation}{\sigma }_{\text{tot}}=\frac{1}{T}\left\langle {F}_{1}{\circ}\frac{\mathrm{d}{X}_{1}}{\mathrm{d}t}\right\rangle +\frac{1}{{T}_{\text{eff}}}\left\langle {F}_{2}{\circ}\frac{\mathrm{d}{X}_{2}}{\mathrm{d}t}\right\rangle ,\end{equation} \tag{ E1 }$$

where F₁ = F₁(X₁, X₂), F₂ = F₂(X₁, X₂) 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

$$\begin{equation}\left\langle F{\circ}\frac{\mathrm{d}X}{\mathrm{d}t}\right\rangle =\underset{t\to \infty }{\mathrm{lim}}\enspace \frac{1}{t}{\int }_{0}^{t}F({t}^{\prime }){\circ}\mathrm{d}X({t}^{\prime }),\end{equation} \tag{ E2 }$$

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:

$$\begin{equation}\left\langle F{\circ}\frac{\mathrm{d}X}{\mathrm{d}t}\right\rangle \simeq \frac{1}{\tau }\sum\limits _{i=1}^{n}\left(\frac{F({t}_{i})+F({t}_{i-1})}{2}\right)(X({t}_{i})-X({t}_{i-1})),\end{equation} \tag{ E3 }$$

where t_i = iΔt and n = τ/Δt.

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 $\tau$ = 300 s. We then partition the simulation into ten consecutive traces of duration T = $\tau$/10 = 30 s. For each of these traces {X_α (t)} (α = 1, ..., 10) we compute the power spectral density as ${C}_{\alpha }(f)=(1/T){\left\vert {\int }_{0}^{T}{X}_{\alpha }(s){e}^{2\pi ift}\enspace \mathrm{d}t\right\vert }^{2}$. We then calculate the average of the power spectral density over the ten different traces $~{C}(f)=(1/10){\sum }_{\alpha =1}^{10}{C}_{\alpha }(f)$ and fit the estimate $~{C}(f)$ as a function of f to the sum of two Lorentzian functions [4, 22, 59]

$$\begin{equation}~{C}(f)=\frac{A}{{({f}_{\mathrm{o}}/2Q)}^{2}+{(f-{f}_{\mathrm{o}})}^{2}}+\frac{A}{{({f}_{\mathrm{o}}/2Q)}^{2}+{(f+{f}_{\mathrm{o}})}^{2}},\end{equation} \tag{ F1 }$$

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).

Zoom In Zoom Out Reset image size