Improving thermodynamic bounds using correlations

Given a set of trajectory data, we may in principle optimize the function $z_{i}(\bm{x})$ in Eq. (23) explicitly by maximizing the Pearson coefficient between $J$ and $Z$. In the case of a pure jump process with $M$ states, this corresponds to optimizing $M-1$ parameters—since the Pearson coefficient is invariant under a global rescaling of $Z$, we may set $z_{1}=1$ without loss of generality. In the case of a diffusion process, we may choose a reasonable set of (say $K$) basis functions $f_{k}(\bm{x})$ and write $z(\bm{x})=\sum_{k=1}^{K}c_{k}f_{k}(\bm{x})$, which (setting $c_{1}=1$) again corresponds to optimizing $K-1$ parameters. Provided that $K$ and the data set are not too large, this explicit optimization is feasible and has the advantage of yielding the tightest possible bound within the accuracy of the optimization algorithm. However, as $K$ or the size of the data set increase, this procedure becomes increasingly unfeasible and less tight but easier to compute heuristic bounds may be desired.

Given that, in many physically relevant situations, the local mean current Eq. (30) is expected to give a reasonably tight bound, we may use an approximate expression $\tilde{\bm{\nu}}^{\text{st}}(\bm{x})$ for the local mean velocity to compute $Z$. Such an approximate expression may be obtained directly from the trajectory data. Even if the measurement resolution is not sufficient to observe every transition, we may still obtain a coarse-grained mean velocity profile. For a spatial resolution $\Delta x$ and a temporal resolution $\Delta t$, we may divide the observation volume into cells of size $\Delta x^{N}$ and obtain an estimate for the local mean velocity in a given cell by considering the displacement at time $t+\Delta t$, averaged over all particles starting out in the cell at time $t$. An approximate local mean velocity may also be obtained from theoretical considerations by considering a simplified version of the dynamics, in which the local mean velocity can be computed explicitly and is expected to resemble the one of the actual dynamics. A simplified model may be obtained, for example, by linearizing a given non-linear model or by constructing a lower-dimensional effective model. Even though both approaches are generally valid only in a certain limit, in many cases, the resulting local mean velocity still yields an improved bound via the CTUR even beyond the strict validity of the approximations. As a concrete example, we consider in the next section the flow through a two-dimensional channel, which we approximate by a one-dimensional flow in an effective potential, for which an explicit expression of the local mean velocity can be obtained.

Another application, for which we anticipate this approach to be useful, is for systems in the presence of spatial disorder. Since the disorder configuration is generally not known precisely, we also cannot obtain an explicit expression for the local mean velocity. In this case, we may compute the local mean velocity for the disorder-free system, and provided that the disorder is not too strong, we can still expect the CTUR with the corresponding observable $Z$ to provide a notable improvement over the TUR. For example, we can imagine a channel as in Fig. 5, through which a current is driven by an external bias. However, due to variances in the fabrication of the channel or the presence of contamination, the actual potential experienced by the particle may be more complicated. Nevertheless, if we can obtain an expression for the local mean velocity of the ideal situation, we can use it to obtain an improved estimate of the entropy production.

Finally, we provide a concrete formula for obtaining a candidate for the observable $Z$. This relies on the observation that, in many systems, the local mean velocity at some position $\bm{x}$ and the probability of finding the particle at position $\bm{x}$ are approximately inversely proportional. This relation is exact for Langevin dynamics in a spatially periodic one-dimensional system; in this case, we have $\nu^{\text{st}}(x)p^{\text{st}}(x)=v^{\text{d}}/L$, where $v^{\text{d}}$ is the drift velocity and $L$ the period of the system. The same relation holds for a Markov jump dynamics with on a ring, $V^{\text{st}}_{i\pm 1,i}p_{i}^{\text{st}}=\pm\mathcal{J}/N$ where $\mathcal{J}$ is the total current on the ring.

However, even in higher-dimensional settings, we often expect the probability density $p^{\text{st}}(\bm{x})$ to be small where the magnitude of $\bm{\nu}^{\text{st}}(\bm{x})$ is large and vice versa. The reasoning for this relation is that, for ergodic dynamics, $p^{\text{st}}(\bm{x})d\bm{x}$ is equal to the temporal occupation fraction, i. e. the fraction of time that the particle resides in the volume element $d\bm{x}$ around $\bm{x}$ in the limit of long observation times. On the other hand, if the local mean velocity at $\bm{x}$ is large, this means that the particle tends to move at a large velocity at $\bm{x}$ and thus we expect it to only spend a short time in the volume element $d\bm{x}$. Note that this qualitative argument may break down, e. g. near vortices in the local mean velocity, where particles may move fast yet still mostly reside in a small volume. However, since we only require an approximation to the local mean velocity to improve the TUR bound on the entropy production, using the time-integral of the inverse occupation fraction as the observable $Z$ is still expected to yield a notable improvement in most cases. Specifically, dividing the observation volume into $\mathcal{N}$ cells (for simplicity, we assume they have equal volume $v$) and the observation time into $\mathcal{M}$ steps, we define as the occupation number $\mathcal{K}_{i}$ the number of times that any trajectory takes a value in the $i$-th cell. Then, we may define the observable $Z$ as

where $i(m)$ is the cell the trajectory is located in at time step $m$. Note that, while $Z$ is evaluated for each individual trajectory, $\mathcal{K}_{i}$ is computed from the entire set of trajectory data. We do not need to normalize the occupation number into an occupation fraction, since the overall normalization factor does not impact the Pearson coefficient. While the observable Eq. (43) may not yield the tightest bound in Eq. (25), it has the advantage that it may be computed using only the trajectory data without any additional input.

Finally, we remark upon a possible issue with the definition Eq. (43). If the phase-space volume $\mathcal{V}$ that is accessible to the system is not finite, then it becomes necessary to introduce a lower cutoff on the occupation number $\mathcal{K}_{i}$. To see this, note that for an ergodic system in the limit of long time and/or a large number of trajectories, we have $\mathcal{K}_{i}\propto p^{\text{st}}(\bm{x}_{i})v$, where $\bm{x}_{i}$ is the location of the $i$-th cell and $v$ is the cell volume. Then,

So the average of $Z$ is proportional to the volume $\mathcal{V}$ and diverges if the latter becomes infinite. For example, if $p^{\text{st}}(\bm{x})\propto\exp(-\|\bm{x}\|^{2}/(2r_{0}))$, then $1/p^{\text{st}}(\bm{x})$ grows exponentially for large $\|\bm{x}\|$ and its average is not well-defined. Intuitively, this means that Eq. (43) is dominated by rare events, which are never well sampled no matter how many trajectories we record. To avoid this problem, we may introduce a lower cutoff on $\mathcal{K}_{i}$, that is, we only count cells that have been visited in at least a fraction $\epsilon$ of all data points. In terms of the probability density, this means that we discard all points with $p^{\text{st}}(\bm{x})<\frac{\epsilon}{v}$. In doing so, the average of $Z$ becomes finite and we may use the resulting observable as a qualitative estimate for the local mean velocity. Obviously, the resulting quantity is only an approximation to the actual occupation fraction, however, we may still expect it to yield a useful improvement over the TUR, provided that $\epsilon$ is sufficiently small and the statistics of the current are not dominated by rare events. This lower cutoff also avoids another possible issue: In Eq. (43), we use the same data set to evaluate the occupation fraction and the trajectory-dependent observable $Z$. While this is not problematic for a large number of trajectories, since the influence of any single trajectory on $\mathcal{K}_{i}$ is small, it may lead to unintended correlations if the number of trajectories is not large. In the latter case, it may be preferable to divide the data set into two parts, using one to determine $\mathcal{K}_{i}$ and the other to compute $Z$. In doing so, we may encounter the situation where a given cell is visited in the second part but not the first part of the data set. In this case, Eq. (43) is no-longer well defined since we divide by zero. Introducing a lower cutoff excludes this possibility. However, we stress that even if $\mathcal{K}_{i}$ and $Z$ are computed from the same data set, Eq. (3) remains valid, so the worst outcome of unintended correlations is to reduce the quality of the estimate on the entropy production.

To demonstrate how Eq. (25) may be used to obtain a tight bound on the entropy production, we consider the model for the $\text{F}_{1}$-ATPase molecular motor introduced in Ref. [21] and further studied in Ref. [22]. This model describes the motion of a probe bead coupled to a rotating molecular motor, which either consumes a molecule called ATP in order to generate rotational motion, or, conversely, generates ATP when driven by an external torque. The probe is considered to be trapped inside a potential $U_{i}(x)$, which is determined both the joint between the probe and the motor and the internal structure of the motor. As the motor rotates in steps of length $L$, the potential depends on the current state of the motor as $U_{i}(x)=U_{0}(x-iL)$. In the simplest form of the model, the trapping potential is harmonic, $U_{0}(x)=kx^{2}/2$, and the motor rotates in steps of $L=120^{\circ}$. The transitions between the states of the motor are described by the position-dependent transition rates

where $W_{i}^{+}$ ($W_{i+1}^{-}$) is the rate of transitions from $i$ to $i+1$ (from $i+1$ to $i$). Here, $W_{0}$ quantifies the overall activity of the motor, which is proportional to the concentration of ATP in the environment of the motor. $\Delta\mu$ is the chemical potential difference driving the rotation; it is the amount of energy gained by the motor when consuming a single molecule of ATP. The parameter $0\leq\alpha\leq 1$ characterizes the asymmetry of the position-dependence of the rates. The spatial motion of the probe is described by one-dimensional Brownian motion in the potential $U_{i}(x)$. In addition, an external torque may be applied to the probe, which enters the equation of motion in the form of a non-conservative bias force $F$. As a consequence, the drift coefficient in Eq. (15) is given by $a_{i}(x)=(-U^{\prime}_{i}(x)-F)/\gamma$, where $\gamma$ is the friction coefficient, $F$ is an external force acting on the probe, while the diffusion matrix is $G=\sqrt{2k_{\text{B}}T/\gamma}$. An experimentally accessible current is the total displacement of the probe,

In the steady state, we have $J=v^{\text{d}}\tau$, where $v^{\text{d}}$ is the drift velocity. Because the system is effectively one-dimensional, the local mean velocity in the steady state is given by $\nu^{\text{st}}(x)p^{\text{st}}(x)=v^{\text{d}}/L$ where $p^{\text{st}}(x)$ is the $L$-periodic steady state probability density. In this case, we can thus reconstruct the local mean velocity from the trajectory data of the probe by taking a histogram of the probe positions. However, it should be noted that this local mean velocity is a coarse-grained quantity, the true local mean velocity entering Eq. (16) also depends on the state $i$ of the motor. We define the observable

which is proportional to $\bar{J}$. Since the proportionality factor cancels in the Pearson coefficient, $Z$ and $\bar{J}$ are equivalent with respect to Eq. (25). To asses the tightness of the various inequalities, we introduce the transport efficiencies

Both of these quantities are smaller than unity and measure the magnitude of the average transport relative to its fluctuations and the dissipation.

The efficiencies Eq. (48) are shown for the molecular motor model in Fig. 1 as a function of various parameters. The top-left panel shows $\eta$ as a function of the base activity $W_{0}$, which corresponds to the concentration of ATP in the experiment. For small activity both $\eta_{J}$ and $\eta_{J,\bar{J}}$ are comparable and small; in this limit, the transitions between the different motor conformations are not translated efficiently into motion and the dissipation is not reflected in the motion of the probe [22]. For large activity, $\eta_{J}$ saturates at a value of around $0.4$. However, when we compute $\eta_{J,\bar{J}}$ in this regime, we find that it saturates at a value close to unity, i. e. the maximum possible value. The top-right panel shows $\eta$ as a function of the chemical potential difference $\Delta\mu$. While this value cannot be readily changed in experiment, it yields important insight into the nature of the bound Eq. (25). For small $\Delta\mu$ the system is almost in equilibrium, and both the TUR and CTUR are close to an equality. However, as we drive the system out of equilibrium, the TUR ratio quickly drops, while the CTUR remains close to unity. This suggests that, while the TUR is generically only saturated close to equilibrium [33, 4], the improved bound from the CTUR Eq. (25) can yield an accurate estimate of the entropy production even far from equilibrium. The bottom-left panel shows $\eta$ as a function of the external load force. Close to the stall condition $FL=\Delta\mu$ neither of the bounds is tight. This is reasonable, since when the motor stalls, also the probe stops moving, while the motor keeps changing its conformation and thus dissipating energy. Interestingly, both bounds are close to unity slightly above the stall condition, i. e. when the external load is just strong enough to turn the motor in the opposite direction. Away from the stall condition, we again find that the TUR is rather loose, while the CTUR remains tight. Note that in this panel we also included the results obtained by optimizing the observable $Z$. Specifically, we write

and then numerically optimize the parameters $a_{k}$, $b_{k}$ such that $\chi(J,\tilde{Z})^{2}$ is maximal for the given trajectory data. Here we use $K=10$; further increasing of the number of parameters provides no notable improvement. The numerical optimization is done using Mathematica’s NMaximize command. The result are the empty orange squares in the bottom-left panel of Fig. 1. As can be seen, the observable $\bar{Z}$ is not truly optimal, so that some improvement of the lower bound on entropy production is possible. However, the heuristic choice $\bar{Z}$ already provides a useful estimate without the need for any parameter optimization. Finally, the bottom-right panel shows $\eta$ as a function of the asymmetry parameter $\alpha$. As was shown in Ref. [22] the motor can operate without internal dissipation close to $\alpha=0$, whereas other values of $\alpha$ result in a finite amount of internal dissipation. Since this difference is most pronounced at low activity, we choose $W_{0}$ such that the velocity remains constant at $v=0.65v_{\text{max}}$ for all values of $\alpha$; this corresponds to $W_{0}\tau_{v}\approx 10^{-2}$ in the top-left panel. While in this regime of low activity, neither bound is saturated, Eq. (25) does yield a considerable improvement over the TUR. Interestingly, neither bound shows a pronounced dependence on $\alpha$, which indicates that their tightness, at least for this model, is not related to the amount of internal vs. external dissipation.

The molecular motor model is essentially one-dimensional: Since the energy scale of the coupling between the probe particle and motor are large compared to the temperature, the position of the probe and the state of the motor are tightly coupled, so measuring the displacement of the probe is almost equivalent to measuring the internal transitions of the motor. Indeed, considering the internal transitions of the motor as an additional current $J_{2}$, Eq. (21) yields no notable improvement over the TUR. Because of this, it is reasonable that the entropy production can be estimated accurately from a measurement of the probe position. However, this tight coupling breaks down for slow switching (see the top-left panel of Fig. 1): In this regime, while the asymmetry in the transitions between the motor states still gives rise to entropy production, these transitions increasingly fail to produce a directed motion of the probe and thus its displacement can no longer be used to obtain an accurate estimate of entropy production.

In summary, we find that the TUR is generally not tight for this model, which mirrors the behavior in other types of molecular motors [14]. Viewed on its own, this would suggest that $\text{F}_{1}$-ATPase is not efficient at saturating the bound set by the TUR. When we take into account the correlations between the velocity and its local mean value via the CTUR, the resulting inequality is almost saturated. Thus, in reality, the motor operates close to the limit permitted by the thermodynamic bound in the biologically relevant parameter regime. In contrast to the TUR, which represents a trade-off between dissipation and precision, the biological meaning of saturating the CTUR is not so clear. We note that, in the case of a particle moving in a flat potential, where the local mean velocity is independent of space, the TUR and CTUR are equivalent and both saturated. The presence of a spatially varying local mean velocity generally increases the amount of fluctuations in the current and also induces correlations between the former and the latter. As a consequence, the TUR becomes less tight and we explicitly need to take into account the correlations to counter this. However, in case of the $\text{F}_{1}$-motor, the spatial inhomogeniety of the dynamics is also an integral part of the mechanism creating the motion of the motor in the first place. In that sense, the CTUR may be interpreted as a trade-off between dissipation, precision and the spatial structure necessary for the motion of the motor.

An obvious question is how the CTUR fares in the presence of multiple independent currents. To investigate this issue, we introduce a simple Markov jump model consisting of $N=10$ sites forming two loops with a common edge, see Fig. 2. While such two-cycle models are also used in modeling the chemical state space of molecular motors [34], we use a simpler setup in order to better understand its behavior in relation to the CTUR.

We parameterize the transition rates between connected states as

where $\beta$ is the inverse temperature, $E_{i}$ the energy of state $i$ and $f$ is the bias, with the choice $+f$ corresponding to counter-clockwise and $-f$ to clockwise transitions. As indicated in Fig. 2, we keep the bias of the right cycle fixed at a value $f_{0}$ while varying the bias $f$ in the left cycle; the bias on the shared link $5\rightarrow 6$ is set to $f-f_{0}$. As the current $J_{\text{L}}$, we choose the number of counterclockwise transitions through the link $2\rightarrow 3$ (indicated in red in Fig. 2), i. e. $\omega_{32}=1=-\omega_{23}$ and $\omega_{ij}=0$ otherwise in Eq. (17).

As the state-dependent observable $Z_{\text{L}}$, we choose the inverse of the steady state occupation fraction in the left ring, $1/p^{\text{st}}_{i|L}$, which corresponds to measuring the ratio of the time spent in state $i$ and the time spent in the left ring. The motivation behind this choice is that, even though the system no longer consists of a single cycle, we still expect the occupation fraction to be approximately inversely proportional to the current through the site, as argued in Section V.3. This observable can be measured by observing only the states in the left ring, for example if the right ring cannot be directly observed. As in the previous example, we define the ratios

which are shown in Fig. 3. In the absence of bias, $f=0$, the average current $J_{\text{L}}$ in the left ring vanishes and thus neither the TUR nor the CTUR give us a non-trivial estimate on the entropy production. The contribution of $J_{L}$ to the total entropy production increases with increasing bias $f$ and thus we can get a non-trivial lower bound on the entropy production. Quantitatively, we observe that the estimate provided by the CTUR is considerably more accurate than the TUR, in particular for larger bias. This finding is independent of the precise potential landscape; indeed, the bound obtained from the CTUR for flat and random configurations of energies $E_{i}$ is qualitatively similar. However, as noted in Eq. (35), in case of a Markov jump process, we are actually estimating the pseudo entropy production Eq. (36) $R$. Since the latter is smaller than the entropy production $\Delta S^{\text{irr}}$ for large bias (i. e. asymmetry between forward and reverse transitions), the estimate becomes worse when we increase the bias further, even though most of the entropy production in the system now stems from the left ring. To quantify the difference between $R$ and $\Delta S^{\text{irr}}$, we also consider the modified efficiencies

As can be seen from Fig. 3, the CTUR saturates this bound in the limit of large bias; thus, we can accurately estimate the pseudo entropy production from a measurement of the current and the occupation probabilities in the left ring. By contrast, just as for the molecular motor model, the TUR estimate becomes less accurate in the limit of large bias.

In the case of small bias, it is expected that the (pseudo) entropy production cannot be estimated accurately, since it mostly originates from the right ring, and is not reflected in the current through the left ring. Thus, at most, we can estimate the contribution to the (pseudo) entropy production that originates from the left ring. To illustrate this, we compare the estimate of $R$ obtained from the TUR and CTUR with the pseudo entropy production $R_{\text{L}}$ of only the left ring, i. e. the system in Fig. 2 with $f_{0}=0$. There is no rigorous inequality between the two, since the bias in the right ring may affect the occupation probabilities and currents in the left ring. However, we do observe that the CTUR for the current in the left ring almost exactly reproduces $R_{L}$. This is an example of a more general statement: While the CTUR can potentially provide an accurate estimate of the (pseudo) entropy production, it can only capture the contributions to entropy production that correspond to the chosen current observable.

Both models in the previous example have a very simple geometry: transitions only occur between neighboring sites along one-dimensional structures. To illustrate the usefulness of the CTUR also in more complicated situations, we study the transport of a Brownian particle in a two-dimensional, soft-walled channel described by the potential

Along the transverse $y$-direction, the particle is confined in a parabolic trap, whose strength is periodically modulated in the longitudinal $x$-direction, see Fig. 5. The period $L$ of the modulation is described by the wave number $\lambda=2\pi/L$, its amplitude by the parameter $0\leq\alpha<1$, where the condition $\alpha<1$ ensures that the potential remains trapping in the transverse direction. The longitudinal motion in the potential Eq. (53) is qualitatively similar to the motion in a one-dimensional periodic potential. This analogy can be made explicit by assuming that the relaxation in the transverse direction is fast compared to the longitudinal motion, so that the transverse coordinate is distributed according to the equilibrium distribution,

where $p^{\text{eq}}(y|x)$ is the conditional probability density of finding the particle at transverse coordinate $y$ given the longitudinal position $x$ and $Z(x)$ is the normalizing partition function. This assumption is justified for a relatively narrow channel and in the absence of bias. If it holds, we can describe the motion in the longitudinal direction via an effective one-dimensional Langevin dynamics

with the potential

For the form Eq. (53) of the potential, the integral can be evaluated explicitly and we find

The effective potential governing the motion in the longitudinal direction appears only due to diffusion in the transverse direction—in the center of the channel, the potential is flat—and is thus proportional to the temperature, rather than the depth of the channel potential. If we assume Eq. (55) to describe the dynamics in the longitudinal direction, then we can use it to compute the local mean velocity in the presence of a constant bias $F_{0}$,

We remark that this expression is only an approximation, since, strictly speaking, the assumption of a transverse equilibrium distribution Eq. (54) breaks down in the presence of a bias. Thus, the true local mean velocity has to be computed using Eq. (53) and depends on both the longitudinal and transverse coordinate. However, the solution of the two-dimensional problem is considerably harder and does not permit a closed-form expression. Thus, we use Eq. (58) as an approximate expression and the corresponding observable

to compute the CTUR for the current given by the displacement in the longitudinal direction. Note that, in principle, we may obtain even simpler approximations for the local mean velocity. Since we expect the dynamics to be described by an effectively one-dimensional motion, the local mean velocity and probability density are related by $\nu^{\text{st}}(x)\propto 1/p^{\text{st}}(x)$, see the discussion in Section V.3. Provided that the bias is not too strong, we may further expect the probability density to be qualitatively similar to the equilibrium density. Thus, a crude approximation of the local mean velocity is given by

where we evaluated the potential at $y=\sqrt{k_{\text{B}}T/k}$, which is a measure of the average distance of the particle from the center of the channel. The resulting bounds Eq. (48) obtained from the TUR and CTUR are shown in Fig. 6. For small bias and thus close to equilibrium, both the TUR and CTUR are close to an equality [33]. As we drive the system further away from equilibrium by increasing the bias, the TUR ratio decreases sharply. In this regime, the CTUR with either approximate local mean velocity yields a tighter bound than the TUR and remains close to an equality. For even stronger bias, the equilibrium argument behind the approximate expressions Eq. (58) and Eq. (60) as well as the one-dimensional approximation become increasingly invalid, and thus the tightness of the CTUR likewise starts to decrease. However, we stress that for moderately strong bias, the estimate on the entropy production from the CTUR with Eq. (58) is around $80\%$ of the true value, which is considerably tighter than the TUR estimate of less than $20\%$.