Thermodynamic precision of a chain of motors: the difference between phase and noise correlation

Biological processes continuously involve energy dissipation and, therefore, are inherently out of thermodynamic equilibrium [1]. A paramount case is that of molecular motors that typically convert chemical energy stored in adenosine triphosphate into mechanical energy to move a cargo or actuate the deformation of parts of a cell [2]. The cooperation of several molecular motors, for instance, during the excitation of traveling waves in flagella or cilia, is still an open problem that can obtain new insight from the application of theoretical results in non-equilibrium statistical physics [3]: in particular, mesoscopic fluctuations—under the lens of stochastic thermodynamics—can provide access to the underlying microscopic mechanisms producing such fluctuations [4]. This has been recently shown in an experiment with mammalian sperm, where the fluctuations of the tail-beating wave have been analyzed and found to be more irregular than expected [5]. The amount of noise in a driven system, such as the beating of a sperm tail, is bounded from below by the amount of dissipation, through the celebrated thermodynamic uncertainty relations (TURs), which we resume here [6–8].

The simplest TUR mainly concerns systems with a time-integrated current $\theta(t)$ that steadily grows and fluctuates, for instance—at not too small times—its average grows as $\langle \theta(t) \rangle = Jt$ and its variance grows as $\textrm{Var}(\theta) = 2 \sigma t$, where J and σ are the constant average current rate and diffusivity, respectively. This minimal picture is qualitatively meaningful for: (1) the position of Brownian motion with an external constant force; (2) a current in a Markov process with broken detailed balance, such as in several models for molecular motors. Of course, it encompasses all stochastic variables, which admit a large deviation principle. In those systems it is possible to define a precision

$$\begin{equation} p = \frac{J^2}{\sigma} \end{equation} \tag{ 1 }$$

which (1) in a steady state is a constant, (2) has the dimension of an inverse time, and (3) its inverse $t^* = 1/p$ represents the time that separates the small time regime dominated by fluctuations and the large time regime dominated by the average growth. A very succinct summary of the so-called TUR states that

$$\begin{equation} p \unicode{x2A7D} \frac{\dot S_p}{k_{\textrm{B}} } = \frac{\dot W}{k_{\textrm{B}} T} \end{equation} \tag{ 2 }$$

where $\dot S_p$ is the entropy production rate in the system, and the last equality holds for systems driven out-of-equilibrium by forces that inject power $\dot W$ and in contact with a thermal bath at temperature T [6–8]. It has also been demonstrated that there is a diffusion process with drift, whose large deviation rate constrains the large deviation rate of the observed current, with a diffusion coefficient that saturates the TUR [9, 10].

In a recent experiment, the maximum 'thermodynamic precision' of the sperm's tail has been measured to be $p \sim \frac{1}{N}\frac{\dot W}{k_{\textrm{B}} T}$, where $N \sim 10^5$ is a good approximation of the number of molecular motors (dynein) present along the axoneme that actuates the tail [5]. Additionally, it has also been seen that reducing $\dot W$ (by slowly depleting the sperm's nutrients) leads to a proportional reduction in p, leaving the $\frac{1}{N}$ factor substantially unaltered. This suggests that the sperm tail is behaving closely to a single dynein motor which, in fact, is semi-optimized [11], with a precision similar to that measured in the tail but an energy consumption N times smaller. Therefore, it seems that the current measured to define the sperm's tail precision is proportional to the current that can be measured to define the dynein's precision: these two 'Brownian clocks' are proportional, including their fluctuations. A comparison with other [12–17]—indirect—measurements of p with sperm in different setups and with C. reinhardtii (which swims using two flagella made of the same axonemal structure but with a smaller length) is coherent with the previous picture. All these observations could be explained by conjecturing a strong coupling of fluctuations in the dynamics of nearby molecular motors, which is reasonable by considering the high dense packing of motors in the axoneme, and by some direct evidence via electron microscopy [18–20].

Here, we address this problem on a more abstract and, possibly, general ground: in a mesoscopic system driven by many microscopic motors, can the direct coupling (spatial correlation) among fluctuations of the microscopic motors be detrimental to the precision of the mesoscopic system? We stress that models of multi-motor systems [21–23] already take into account an indirect coupling among the motors, also called cooperativity, which comes from the fact that each motor is coupled to the state of the axonemal backbone at the point where the motor is attached, so that adjacent motors feel a similar 'external field' and are, therefore, indirectly coupled to each other. Here, we neglect the feedback effect of the backbone onto the motors and, therefore, we ignore the indirect coupling between motors that would be mediated by the backbone. We only consider the direct coupling between adjacent parts of the backbone (through its elasticity k) and—in the second part of the paper—a direct coupling among the noises, which represents a direct (e.g. mechanical) interaction among the underlying motors, see figure 1.

Zoom In Zoom Out Reset image size

We briefly summarise some papers which address questions that are similar but not identical to the present one. In [24], the energetic cost of two oscillators coupled by a generic force is considered. It is shown that the potential part of the interaction always decreases dissipation, while the non-potential part may decrease or increase it. In the same spirit [25] considers coupled oscillators, but with a more complicated coupling, which happens to be purely non-potential, with the result of increasing the cost when the interaction (and therefore order) increases. In [26], the effect of a conservative (Kuramoto) coupling on the TUR is considered, but the TUR is replaced by a 'sub-system-TUR', which compares fluctuations of a sub-system with the energy dissipated by the same sub-system, a measure which, in many cases, gives simply $1/N$ of the total dissipated energy. Basically, the sub-system has a precision which is bounded by the same bound of the whole system. Finally, [27] considers the total thermodynamic cost in a mean-field Kuramoto model, somehow recalling and reproducing the general results from [24] for a single model, where a potential coupling leads to reducing dissipation (when interaction and correlation increase), while no discussion of precision and TUR are given.

A general formalization of this problem is the following: see figure 1. In all these models one has a set of phases $\theta_i(t)$ that represent a discretization (in space) of the backbone of the flagellum, i.e. each ith phase represents a local observable (e.g. the local angle with respect to a fixed axis or the local curvature). The local (ith) phase evolves under the effect of interaction and motors, $\dot \theta_i(t) = F_i(t)$, with $F_i(t) = F_{i,\mathrm{int}}(t)+F_{i,\mathrm{mot}}(t)$: $F_{i,\mathrm{int}}$ represents some non-linear couplings with adjacent parts of the backbone, e.g. $F_{i,\mathrm{int}}(\theta_i,\theta_{i+1},\theta_{i-1})$ (the coupling is typically non-confining, so that phase slips are allowed), while $F_{i,\mathrm{mot}} = \omega_i+\eta_i(t)$ represents the coupling with the ith molecular motor, which has an average driving force ω_i and associated fluctuations $\eta_i(t)$. In this paper, we consider $F_{i,\mathrm{int}}$ to be conservative, i.e. derived from a potential, and of the Kuramoto type.

A typical precision measure in these models is related to a single current: the simplest case given by θ_i with a given i, i.e. in the steady state $p = \lim_{t \to\infty}\frac{2}{t}\langle \theta_i \rangle^2/\langle (\delta \theta_i)^2 \rangle$ with $\delta \theta_i = \theta(t)-\langle \theta \rangle$. The comparison of such a single current, in the spirit of the TUR, is with some measure of energy cost or dissipation: the total dissipation rate in the system can be decomposed as $\dot W = \sum_i \dot W_i$, where $\dot W_i = \langle \dot \theta_i \circ F_i \rangle$ is the energy dissipation rate of the single variable θ_i , which is, in principle, influenced by couplings. Conservative couplings (which is the case for Kuramoto models) always have a non-increasing effect on the total dissipation rate, and usually leave it unchanged. Non-potential couplings can either increase or decrease the dissipation. We have verified that the total dissipation rate in our case is substantially unaltered by the value of k.

In section 2, we present the Kuramoto model with the two kinds of noise which we consider in our study. In section 3, the results with independent noises are discussed. In section 4, we report the results with spatial correlations in the noise. Our conclusions and perspectives are drawn in section 5.

We consider a chain of N coupled oscillators:

$$\begin{equation} \dot\theta_i = \omega_i -\frac{k}{2}\left[\text{sin}\left(\theta_i-\theta_{i-1}\right)+\text{sin}\left(\theta_i-\theta_{i+1}\right)\right]+\sqrt{2D}\eta_i\left(t\right) \end{equation} \tag{ 3 }$$

with $i = 1\ldots N$, $k \unicode{x2A7E} 0$ is the coupling constant, D is the bare diffusivity parameter, η_i is a white noise with $\langle \eta_i \rangle = 0$ and two possible recipes for spatial correlations. The first is the uncorrelated case, discussed in section 3, the second is the correlated one, discussed in section 4:

$$\begin{align} \langle \eta_i\left(t\right)\eta_j\left(t^{^{\prime}}\right)\rangle & = \delta_{ij}\delta\left(t-t^{^{\prime}}\right) \end{align} \tag{ 4a }$$

$$\begin{align} \langle \eta_i\left(t\right)\eta_j\left(t^{^{\prime}}\right)\rangle & = \sum_{n = 1}^{N} \mathrm{e}^{-|i-n|/\Delta} \mathrm{e}^{-|j-n|/\Delta}\delta\left(t-t^{^{\prime}}\right)\nonumber\\ & = \frac{1}{1-\mathrm{e}^{-2/\Delta}}\bigg\{\mathrm{e}^{-|d|/\Delta}\left[1+|d|+\left(1-|d|\right)\mathrm{e}^{-2/\Delta}\right]- \mathrm{e}^{-s/\Delta}\nonumber\\ &\quad -\mathrm{e}^{-\left(2N-s+2\right)/\Delta}\bigg\} \end{align} \tag{ 4b }$$

where $d\equiv j-i$ and $s\equiv j+i$. Formula (4b ) is chosen for its easiest implementation, see section 4. When N is large (already, with $N = 10^2$, this is true), and s is of the order N (e.g. if i is in the bulk of the chain) and $|d| \gg 1$, equation (4b ) reduces to $\langle \eta_i(t)\eta_j(t^{^{\prime}})\rangle \approx |i-j|\mathrm{e}^{-|i-j|/\Delta}$. In the present work, we mainly consider uniform driving, $\omega_i = \omega$, with a few exceptions, where the effect of a narrow distribution of ω_i is discussed.

When k = 0 we have an exactly overdamped Langevin equation for each variable θ_i with $J = \omega_i$; then the uncoupled value of the local precision is $p_{\textrm{unc},i} = \omega_i^2/D$ and, in the uniform case, we can define $p_{\textrm{unc}} = \omega^2/D$. It is immediately evident that in the uncoupled case the TUR, equation (2), is not saturated, since $\dot W = N\omega^2$ and $k_{\textrm{B}} T = D$.

We recall that the Kuramoto chain, in contrast with the original mean-field Kuramoto model, does not have a synchronization transition and a long-range order [28, 29], but short-range order can be revealed by a suitable order parameter, as discussed later [30, 31].

The aim of this section is to show the effect of chain short-range order upon the precision of phase fluctuations. We have simulated the Kuramoto chain for several choices of the parameters D, k, $\{\omega\}$ for long trajectories. From each long trajectory, we focus on the central oscillator (i = 50 in the middle of the chain of length N = 100) to reduce the effects of the boundaries. We have verified that—apart from the oscillators very close to the boundaries—the results do not strongly depend upon i. We have computed time averages to get values for the average phase drift velocity $\langle \dot \theta_i \rangle$ and the mean-squared displacement (MSD), which is shown for the site i = 50 in figure 2. The MSD is seen to reach the asymptotic normal regime $\text{MSD}_i \sim 2 D_i t$ after more or less long transients, depending on the parameters. In the case of uniform driving, the asymptotic MSD is bounded from above by the raw diffusion coefficient ${\sim}2 D t$, reached for strong $D/k$, and from below by the center-of-mass diffusion ${\sim}2(D/N) t$, reached for low $D/k$, as explained later. When the driving is not uniform, the situation is similar, but the bounds are not strict and one can observe strong differences between different choices of the values of ω_i within the same distribution.

Zoom In Zoom Out Reset image size

In all cases, however, it is possible to measure the diffusivity D_i and compute the precision

$$\begin{equation} p_i = \frac{\langle \dot \theta_i \rangle^2}{D_i} \end{equation} \tag{ 5 }$$

which is shown in figure 3. For uniform $\omega_i = \omega$ we get a neat collapse by plotting $k p_i/\omega^2$ vs $D/k$. The master curve shows two clear instances of limiting behavior, already anticipated by the two instances of limiting behavior of D_i

Zoom In Zoom Out Reset image size

$$\begin{align} p_i \approx 1/D \;\;\;\; \left(D/k \gg 1\right) \end{align} \tag{ 6a }$$

$$\begin{align} p_i \approx N/D \;\;\;\; \left(D/k \ll 1\right) .\end{align} \tag{ 6b }$$

The explanation for this limiting behavior is the following. When $D/k$ is large, the level of order is low (see later, where we characterize the degree of order in the chain), and therefore each oscillator is substantially independent from the others and follows a biased Brownian motion $\theta_i \sim \omega_i t+ W^D_t$, where $W^D_t$ is the Wiener process with diffusivity D. By contrast, when $D/k$ is small, there is a high level of order and all the oscillators fluctuate close to each other and, consequently, close to the average phase $\overline{\theta} = \frac{1}{N}\sum_{i = 1}^N \theta_i$. We note that the average motion is not influenced by the internal interactions, i.e.

$$\begin{equation} \overline{\dot\theta} = \overline{\omega}+ \frac{\sqrt{2D} }{N}\sum_{i = 1}^N \eta_i = \overline{\omega}+ \sqrt{\frac{2D}{N}} \eta \end{equation} \tag{ 7 }$$

where η is a white noise with 1 amplitude. The MSD of the single phase, asymptotically, cannot be slower or faster than the MSD of the average phase, i.e. $D_i \to D/N$, and this explains the behavior $p_i \approx N/D$ at small $D/k$.

In the case of uniform phase velocity $\omega_i = \omega = \overline\omega$, the rate of dissipated energy is trivially scaling with N, i.e. $w = \sum_i \langle \omega_i \dot \theta_i \rangle = N \omega^2$. Therefore, the maximum precision imposed by the TUR is $p_{\mathrm{max}} = w/D = N\omega^2/D$, and this maximum is exactly saturated in the strong coupling limit, $D/k \ll 1$.

Evidence of short-range order in the chain is shown in figure 4, where the short-range order parameters are analyzed, defined as

$$\begin{equation} Q = \frac{1}{N-1}\sum_{i = 1}^{N-1}|\text{sin}\left(\theta_{i+1}-\theta_i\right)|. \end{equation} \tag{ 8 }$$

Small values of Q indicate higher order. Here, Q smoothly increases with D, signaling a decreasing order, as expected, with the exception of the case with non-uniform driving, where the order parameter is somehow constant (or slightly decreasing with D) at small values of D. Somehow the distribution of ω_i acts as a kind of noise, disturbing the short-range order. It is also clear that non-uniform driving can produce quite different values of Q, from sample to sample, at least at small/intermediate values of D.

Zoom In Zoom Out Reset image size

To obtain noises with spatial correlation, as given in equation (4b ), we define

$$\begin{equation} \eta_i\left(t\right) = \sum_{j = 1}^N \mathrm{e}^{-\frac{|i-j|}{\Delta}}\xi_j\left(t\right) \end{equation} \tag{ 9 }$$

where $\xi_j(t)$ are independent white noises with amplitude 1. As discussed below, equation (4b ), Δ acts as an effective spatial correlation length for the noise. We note that the auto-correlation $\langle \eta_i(t)^2 \rangle$ of this kind of noise is not 1, but depends on both i (weakly) and Δ. Therefore, to compare the results of such a noise implementation with those in the previous section, we define $D_{\mathrm{eff}} = D \langle \eta_i(t)^2 \rangle$, where i = 50 is the site where we are measuring the precision. Using this definition, we plot, in figure 5, the precision as a function of $D_{\mathrm{eff}}/k$. The idea behind this re-scaling of D is that we have, on the x axis, a measure of the amount (amplitude) of noise on the elected site, and it is interesting to check if this amount is meaningful in some regimes. In fact, we see that when $D_{\mathrm{eff}}/k$ is large enough (on values comparable to the case where the same transition happens for independent noises, figure 3), any effect of coupling disappears and the precision is that of an uncoupled particle with diffusivity $D_{\mathrm{eff}}$, i.e. $p \sim 1/D_{\mathrm{eff}}$. By contrast, when $D_{\mathrm{eff}}/k$ is small, the precision saturates on a line which is smaller than $N/D_{\mathrm{eff}}$. In figure 5, we plot the lines corresponding to the precision of the center of mass, which is $p_{\mathrm{cm}} = \omega^2/D_{\mathrm{cm}}$, with $D_{\mathrm{cm}}$ discussed in the appendix. It is seen that one can define an effective number of degrees of freedom, $N_{\mathrm{eff}} = D_{\mathrm{eff}}/D_{\mathrm{cm}} \lt N$, which decreases when Δ increases (see figure 6). Summarizing the results in figure 5, the effect of spatial correlations on the noise is negligible for the precision at strong noise, but becomes relevant and reduces p at small noise. Interestingly, the effect of increasing Δ on the short-range order of the phase, when $D_{\mathrm{eff}}/k$ is small, is negligible, see figure 4. For larger values of $D_{\mathrm{eff}}/k$ some ordering effects are seen, but it is less clear because they depend upon the choice of the observable. To obtain some further insights into the question of spatial order in the chain and of its dependence with Δ, we plot, in figure 7, the power spectra of the field $\text{sin} \theta_i(t)$. These are defined as

$$\begin{equation} S\left(q\right) = \frac{2}{N} \left\langle \left| \sum_{i = 1}^N \mathrm{e}^{-2\pi\mathcal{I} q i/N} \text{sin}{\theta_i\left(t\right)} \right|^2 \right\rangle \end{equation} \tag{ 10 }$$

where $\mathcal{I}$ is the imaginary unit. The case $\Delta = 0$ is the most clear one; it goes from an almost flat spectrum at large D, with very small precision, to perfectly q⁻² behavior, marking high spatial order and very large precision. When $\Delta\gt0$ we observe two main differences: (1) the spectrum is non-flat (i.e. there is spatial order in the phase), even at very large $D_{\mathrm{eff}}$. This is remarkable if compared with the fact that the precision is the same as in the case of $\Delta = 0$, i.e. it is the uncoupled value ${\sim}1/D_{\mathrm{eff}}$; (2) an additional length scale, certainly related to both $D_{\mathrm{eff}}$ and Δ, appears at small values of $D_{\mathrm{eff}}$, representing a maximum length below which the q⁻² behavior is observed and above which it is lost.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In conclusion, we have presented evidence of the non-trivial effect—on the precision—of different sources of correlation in a chain of coupled phases. When the coupling is induced by a simple force aligning the phases of adjacent sites, the precision fairly reflects the degree of order in the phase field, i.e. it is larger when the spatial order is higher. By contrast, when there is a correlation in the noise, even if the degree of order in the phase field is not strongly influenced, the precision dramatically drops as if the number of effective degrees of freedom is reduced. We can summarize our findings as follows: the precision p, being influenced by the diffusion coefficient, reflects a degree of dynamical order, which is not entirely explained by the degree of spatial (static) order.

Our observations represent, in our opinion, a first step in order to understand the low precision observed in the fluctuations of the flagellar beating, such as for the tails of sperm and for the cilia of C. reinhardtii. Certainly, the connection with biophysics must go through more refined and realistic models of noisy flagellar beating, investigations of which are in progress.

The authors would like to thank Claudio Maggi and Alessandro Sarracino for useful discussions. The authors acknowledge the financial support of the MIUR PRIN 2017 project 201798CZLJ.

To obtain the expression in equation (4b ) we need to evaluate the various contributors for the different values of the quantities in the exponents. By first considering the case with j > i, we can split the sum into three parts depending on the j values, that is

$$\begin{align} \langle \eta_i\left(t\right)\eta_j\left(t\right)\rangle = \sum_{n = 1}^{i-1} \mathrm{e}^{-\left(i+j-2n\right)/\Delta}+ \sum_{n = i}^{j-1} \mathrm{e}^{-\left(j-i\right)/\Delta} + \sum_{n = j}^{N} \mathrm{e}^{-\left(2n-i-j\right)/\Delta}. \end{align} \tag{ A1 }$$

Using the expression for the sum of a geometric series, the terms in equation (A1) can be written as

$$\begin{align} \langle \eta_i\eta_j\rangle & = \frac{\mathrm{e}^{-\left(j+i\right)/\Delta}}{1-\mathrm{e}^{-2/\Delta}} \left( \mathrm{e}^{2\left(i-1\right)/\Delta} -1 \right) + \left(j-i\right)\mathrm{e}^{-\left(j-i\right)/\Delta}+ \frac{\mathrm{e}^{\left(j+i\right)/\Delta}}{1-\mathrm{e}^{-2/\Delta}} ( \mathrm{e}^{-2j/\Delta}- \mathrm{e}^{-2\left(N+1\right)/\Delta}) \nonumber\\ & = \frac{1}{1-\mathrm{e}^{-2/\Delta}}\Big\{\!\mathrm{e}^{-\left(j-i\right)/\Delta} \cdot \Big[j-i+1+ \left(1-j+i\right)\mathrm{e}^{-2/\Delta}\Big] -\mathrm{e}^{-\left(j+i\right)/\Delta} - \mathrm{e}^{-\left(2N+2-j-i\right)/\Delta}\!\Big\} .\end{align} \tag{ A2 }$$

The final expression in equation (4b ) is obtained by taking into account that the case j < i is given by exchanging the indices and using the definitions $d\equiv j-i$ and $s\equiv i+j$.

Putting j = i in equation (A2), we obtain the self-correlation of the noise η_i that is given by

$$\begin{equation} \langle \eta^2_i\rangle = \frac{1+\mathrm{e}^{-2/\Delta}-\mathrm{e}^{-2i/\Delta} - \mathrm{e}^{-2\left(N+1-i\right)/\Delta}}{1-\mathrm{e}^{-2/\Delta}}, \end{equation} \tag{ A3 }$$

that can be used to have the quantity $D_{\mathrm{eff}} = D\langle \eta^2_i\rangle$.

By considering equations (7) and (9), the noise of the average motion can be written as

$$\begin{equation} \eta_{\mathrm{cm}}\left(t\right) = \sum_{i = 1}^{N}\eta_i\left(t\right) = \sum_{i = 1}^{N}\sum_{j = 1}^{N} \mathrm{e}^{-|i-j|/\Delta} \xi\left(t\right) \end{equation} \tag{ A4 }$$

and considering the properties of noise $\xi(t)$ its amplitude is

$$\begin{equation} \langle \eta_{\mathrm{cm}}^2\rangle = \sum_{i = 1}^{N}\sum_{j = 1}^{N} \sum_{n = 1}^{N} \mathrm{e}^{-|i-j|/\Delta} \mathrm{e}^{-|n-j|/\Delta}. \end{equation} \tag{ A5 }$$

The final expression of $\langle \eta^2(t)\rangle$ can be obtained by evaluating the different contributors of the absolute values and, after long but simple calculations, we can write that

$$\begin{align} \langle \eta_{\mathrm{cm}}^2\rangle & = \frac{1}{\left(1+\mathrm{e}^{-1/\Delta}\right)\left(1-\mathrm{e}^{-1/\Delta}\right)^3} \Big\{N\left(1-\mathrm{e}^{-2/\Delta}\right)\cdot \Big[1+\mathrm{e}^{-2/\Delta}+2e^{-1/\Delta}\Big( 1+\mathrm{e}^{-N/\Delta} \Big)\Big] \nonumber\\ &\quad -4\mathrm{e}^{-1/\Delta}-6\mathrm{e}^{-2/\Delta}-4\mathrm{e}^{-3/\Delta}+4\mathrm{e}^{-\left(N+1\right)/\Delta} + 8\mathrm{e}^{-\left(N+2\right)/\Delta}+ 4\mathrm{e}^{-\left(N+3\right)/\Delta}-2\mathrm{e}^{-2\left(N+1\right)/\Delta} \Big\}. \end{align} \tag{ A6 }$$

The $D_{\mathrm{cm}}$ value in the main text is given by $D_{\mathrm{cm}} = D\langle \eta_{\mathrm{cm}}^2\rangle/N^2$.