††thanks: E-mail<EMAIL_ADDRESS>
# Improved estimation for energy dissipation in biochemical oscillations

Zhiyu Cao  Zhonghuai Hou  Department of Chemical Physics & Hefei National
Laboratory for Physical Sciences at Microscales, iChEM, University of Science
and Technology of China, Hefei, Anhui 230026, China

###### Abstract

Biochemical oscillations, regulating the timing of life processes, need to
consume energy to achieve good performance on crucial functions, such as high
accuracy of phase period and high sensitivity to external signals. However, it
is a great challenge to precisely estimate the energy dissipation in such
systems. Here, based on the stochastic normal form theory (SNFT), we calculate
the Pearson correlation coefficient between the oscillatory amplitude and
phase, and a trade-off relation between transport efficiency and phase
sensitivity can then be derived, which serves as a tighter form than the
estimator resulting from the conventional thermodynamic uncertainty relation
(TUR). Our findings demonstrate that a more precise energy dissipation
estimation can be obtained by enhancing the sensitivity of the biochemical
oscillations. Moreover, the internal noise and amplitude powers effects have
also been discovered.

##  I Introduction

To achieve good performance of certain functions, living systems are
inherently nonequilibrium and dissipative. Recently, the relationship between
biochemical functions and nonequilibrium thermodynamics has been an active
area in statistical physics community (Bialek and Setayeshgar, 2005; Hu et
al., 2010; Lan et al., 2012; Lan and Tu, 2013; Skoge et al., 2013; Lang et
al., 2014; Sartori and Tu, 2015; Cao et al., 2015; Fei et al., 2018; Mora and
Nemenman, 2019; Zhang et al., 2020; Barato and Seifert, 2016; Lee et al.,
2018; Hasegawa, 2018; Hasegawa and Van Vu, 2019a; Marsland III et al., 2019;
Del Junco and Vaikuntanathan, 2020a; Guan et al., 2020; Del Junco and
Vaikuntanathan, 2020b; Cao et al., 2021). For instance, Lan $et$ $al.$ have
revealed a powerful trade-off relation between energy dissipation rate,
adaption speed and the maximum adaption accuracy underlying many sensory
systems (Lan et al., 2012; Lan and Tu, 2013). Lang $et$ $al.$ have
investigated the fundamental thermodynamic constraints on statistical
inference and learning of biochemical signaling networks (Lang et al., 2014).
Particularly, for biochemical oscillations which are essential in regulating
the timing of life processes, such as the cell cycle, circadian clocks, and
glycolysis, both accuracy of the period and sensitivity to external cues can
be ensured by dissipative processes simultaneously (Ferrell Jr et al., 2011;
Buzsáki and Draguhn, 2004; Nakajima et al., 2005; Novák and Tyson, 2008;
Goldbeter, 1997; Martiel and Goldbeter, 1987). Therefore, it is important to
measure the free energy dissipation in biochemical oscillation systems that
maintains the cyclic dynamics. However, in actual experiments, how to infer
the energy dissipation is of great challenge (Gnesotto et al., 2018; Seifert,
2019).

Recent progress in this topic is the thermodynamic uncertainty relation (TUR)
(Barato and Seifert, 2015; Gingrich et al., 2016; Pietzonka et al., 2016,
2017; Dechant and Sasa, 2018a, b; Agarwalla and Segal, 2018; Hasegawa and Van
Vu, 2019b; Horowitz and Gingrich, 2020), quantifying the trade-off between
energy dissipation $\Delta W$, the average $\left\langle R\right\rangle$ and
variance $\text{Var}\left(R\right)=\left\langle(R-\left\langle
R\right\rangle)^{2}\right\rangle$ of a time-integrated current observable $R$
in nonequilibrium steady states (here
$\beta={\color[rgb]{1,0,0}{\color[rgb]{0,0,0}1/k_{B}T}}$, $T$ is the
temperature of the environment and $k_{B}$ is the Boltzmann constant):

$\eta(R)=\frac{2\left\langle R\right\rangle^{2}}{\text{$\beta$Var}(R)\Delta
W}\leq 1,$ (1)

where $\eta(R)$ is the transport efficiency to properly quantify the
performance of living systems working with high accuracy, but low energy
dissipation (Dechant and Sasa, 2018b). Directly, TUR yields that the magnitude
of current fluctuation provides a lower bound of energy dissipation as $\Delta
W\geq\Delta W_{\text{TUR}}\equiv 2k_{B}T\left\langle
R\right\rangle^{2}/\text{Var}(R)$ with $\eta\left(R\right)=\Delta
W_{\text{TUR}}/\Delta W$. If $\eta(R)$ is close to $1$, the TUR acts as a
powerful tool for energy dissipation inference (Li et al., 2019; Manikandan et
al., 2020; Van Vu et al., 2020; Otsubo et al., 2020; Skinner and Dunkel,
2021a, b; Dechant and Sasa, 2021a). For instance, recently Li $et$ $al.$ have
showed that the fluctuations in nonequilibrium currents can be utilized to
infer the dissipation rate for the bead-spring model (Li et al., 2019). Otsubo
$et$ $al.$ have developed a framework for dissipation estimation by using the
TUR along with machine learning techniques (Otsubo et al., 2020), to list just
a few.

However, since the TUR is an inequality, only a rough bound can be provided in
many cases. For instance, it has been revealed by Hwang and Hyeon that the TUR
is generally not tight for several types of molecular motors (Hwang and Hyeon,
2018). Jack $et$ $al.$ have found that the TUR only yields a weak bound for
molecular-scale energy conversion (Jack et al., 2020). Also, in our recent
work (Cao et al., 2020), we have established the TUR for general biochemical
oscillations by calculating the transport efficiency
$\eta(\theta)=2\left\langle\theta\right\rangle^{2}/\text{$\beta$Var}(\theta)\Delta
W<1$, where the observable oscillatory phase
$\theta(\tau)=\int_{0}^{\tau}\dot{\theta}(t)dt$ is the current observable.
Both the analytical and numerical results have shown that the TUR is far from
tight for models of chemical oscillators, providing typically lower estimation
for energy dissipation than the actual value. Therefore, how to obtain a more
qualified estimation than the conventional TUR for biochemical oscillation
systems is still an interesting question.

In the presented paper, we try to address this question by revealing a trade-
off relation between transport efficiency and phase sensitivity (Hasegawa and
Arita, 2014a, b; Fei et al., 2018). The basic idea is to improve the
conventional TUR by considering the Pearson correlations between the chosen
current and another state-dependent observable, based on a strategy proposed
by Dechant and Sasa very recently (Dechant and Sasa, 2021b). For practical
purpose in biochemical oscillation systems, we choose the time integral of the
oscillatory amplitude $r$ as the state-dependent observable, which reads
$Q_{2}(\tau)=\int_{0}^{\tau}r^{2}(t)dt$. By using the stochastic normal form
theory (SNFT) we established before (Hou et al., 2006; Xiao et al., 2007; Ma
et al., 2008; Xiao et al., 2009), explicit theoretical expressions of the
Pearson correlations between $Q_{2}$ and $\theta$ can be derived, which allows
us to obtain the efficiency-sensitivity trade-off relation as
$\eta(\theta)\leq 1-2\alpha\kappa^{2}$ with $\kappa$ the phase sensitivity
characterizing the ability for biochemical circuits to respond to external
signals and $\alpha>0$ the control parameter denoting the distance to the
bifurcation point. Remarkably, this trade-off relation provides a tighter
dissipation estimator for biochemical oscillations than the conventional TUR,
and the precision of this estimator can be further improved by enhancing the
sensitivity. Finally, we demonstrate our statements by detailed numerical
simulations in a circadian clock model.

##  II Improved estimation of the energy dissipation

###  II.1 Stochastic Normal Form Theory (SNFT)

We consider a general biochemical system of size $V$ including $N$ well-
stirred species and $M$ reactions as $(R_{1},\dots,R_{M})$. Generally, the
reaction $R_{\rho}$ can be written as:

$\mathbf{X}\to\mathbf{X}+{\bf v}_{\rho}$

where $\mathbf{X}=\left(X_{1},X_{2},\dots,X_{N}\right)$ with $X_{j}$ the
number of species $j$, and
$\mathbf{v}_{\rho}=\left(v_{\rho}^{1},v_{\rho}^{2},\dots,v_{\rho}^{N}\right)$
with $v_{\rho}^{j}$ the stoichiometric change of species $j$ in $R_{\rho}$. In
a mesoscopic system wherein intrinsic noise cannot be neglected, with the
assumption of existence of a “macro-infinitesimal” time scale (Gillespie,
2000; Xiao et al., 2007), the system’s dynamics can be described by the
chemical Langevin equations (CLEs) as

$\dot{x}_{j}=\sum_{\rho=1}^{M}v_{\rho}^{j}w_{\rho}(\bm{x})+\frac{1}{\sqrt{V}}\sum_{\rho=1}^{M}v_{\rho}^{j}\sqrt{w_{\rho}(\bm{x})}\xi_{\rho}(t),\>j=1,...,N.$
(2)

where $\bm{x}=(x_{1},\dots,x_{N})^{\text{T}}=\mathbf{X}/V$ denotes the
concentration vector, $w_{\rho}(\bm{x})$ is the reaction rate of $R_{\rho}$ as
a function of the concentrations $\bm{x}$, and
${\normalcolor\bm{\xi}(t)=(\xi_{1},\dots,\xi_{M})^{\text{T}}}$ is a vector of
independent Gaussian white noises with zero means and correlations
$\left\langle\xi_{\rho}(t)\xi_{\rho^{\prime}}(s)\right\rangle=\delta_{\rho\rho^{\prime}}\delta(t-s)$.

In the thermodynamic limit with $V\to\infty$, the noise term disappears and
the dynamics is described by the deterministic equation

$\dot{x}_{j}=F_{j}\left(\bm{x}\right)\equiv\sum_{\rho=1}^{M}v_{\rho}^{i}w_{\rho}\left(\bm{x}\right)$
(3)

Generally, to the occurrence of biochemical oscillation, we assume that the
system undergoes a supercritical Hopf bifurcation (HB) with the change of a
certain control parameter $\mu$. Eq.(3) has a unique stable point $\bm{x}_{s}$
with $\bm{F}(\bm{x}_{s})\equiv 0$, which loses stability at the HB point
$\mu=\mu_{c}$, in the way that the Jacobian matrix $\bm{J}$ with components
$J_{ij}=\left(\partial f_{i}/\partial x_{j}\right)|_{\bm{x}=\bm{x}_{s}}$ has a
pair of conjugate eigenvalues $\lambda_{\pm}=\alpha(\mu)\pm i\omega$ with
$\alpha(\mu_{c})=0$ (henceforth we uses $\alpha=\alpha(\mu)$ to represent the
control parameter). In the so-called supercritical region $\alpha>0$
($\mu>\mu_{c}$), the deterministic system shows a stable oscillation with
frequency given by $\omega$ and amplitude growing from zero. In the
subcritical region with $\alpha<0$ ($\mu<\mu_{c}$), no deterministic
oscillation can be observed. In the case where the system size is not large
such that the internal noise term in Eq.(2) can not be ignored, such as for
intracellular biochemical oscillation systems considered here, an interesting
phenomenon known as noise induced oscillations (NIOs) has been observed even
in the subcritical region where $\alpha<0$, demonstrating the constructive
role of internal noise in mesoscopic chemical oscillation systems (Ko et al.,
2010). In addition, an optimal system size exists where the NIO shows best
performance, knows as internal noise coherence resonance (INCR) (Hou and Xin,
2003; Zhou and Kurths, 2002; Hänggi, 2002).

In our previous works (Hou et al., 2006; Xiao et al., 2007; Ma et al., 2008;
Xiao et al., 2009), we have developed a stochastic normal form theory (SNFT)
to successfully elucidate the mechanism underlying NIO and INCR. When the
system locates near the HB, the motion of the oscillatory mode is much slower
than the other $N-2$ stable modes due to time-scale separation. Hence, the
system’s dynamics will be dominated by the oscillatory motion on a 2D center
manifold. According to SNFT, the stochastic dynamics governing the evolution
of the oscillation amplitude $r$ and and phase angle $\theta$ can be described
by (see Appendix A for details)

$\dot{r}=\alpha
r+C_{r}r^{3}+\frac{\varepsilon^{2}}{2Vr}+\frac{\varepsilon}{\sqrt{V}}\eta_{r}(t),$
(4)

$\dot{\theta}=\omega+C_{i}r^{2}+\frac{\varepsilon}{r\sqrt{V}}\eta_{\theta}(t)$
(5)

wherein $C_{r}<0$ and $C_{i}>0$ are system-dependent constants determined by
the nonlinear terms of $\bm{F}\left(\bm{x}\right)$ at the stable point,
$\eta_{r}$ and $\eta_{\theta}$ are independent Gaussian white noises with zero
mean and unit variance, $\varepsilon$ denotes an effective noise intensity
determined by the details of $\bm{F}\left(\bm{x}\right)$. According to Eqs.(4)
and (5), the steady-state (SS) distribution of $r$ reads

$p_{ss}\left(r\right)=N_{r}\exp\left[-\frac{V}{4\varepsilon^{2}}\left(2\alpha
r^{2}+C_{r}r^{4}\right)+\ln r\right]$ (6)

and $\theta$ is uniformly distributed with $\left[0,2\pi\right]$. Therefore,
the system exhibits a stochastic oscillation with most-probable amplitude
given by

$r_{m}=\left(-\frac{\sqrt{\alpha^{2}-2C_{r}\varepsilon^{2}/V}+\alpha}{2C_{r}}\right)^{1/2}$
(7)

satisfying $\partial p_{ss}\left(r\right)/\partial r|_{r_{m}}=0$.

Clearly, in the deterministic limit ($V\to\infty$),
$r_{m}=\sqrt{-\alpha/C_{r}}$ corresponding to a stable limit cycle and
frequency
$\omega_{s}=\omega+C_{i}r_{m}^{2}=\omega+\alpha\left|C_{i}/C_{r}\right|$,
which only exists for $\alpha>0$ in the supercritical region. If the system
size is finite, however, the internal term $2C_{r}\varepsilon^{2}/V$ in the
square-root will take effect and $r_{m}$ is not zero even for $\alpha<0$
(subcritical region), corresponding to the occurrence of NIO. In the case
$\left|\alpha\right|\gg 2C_{r}\varepsilon^{2}/V$, one has for NIO
$r_{m}\simeq\varepsilon/\sqrt{-2\alpha V}$ which scales as $V^{-1/2}$, and the
frequency is approximately
$\text{$\omega_{s}$}\simeq\omega+C_{i}\varepsilon^{2}/\left(2\left|\alpha\right|V\right)$.
Therefore,

$\omega_{s}=\begin{cases}\omega+\alpha\left|C_{i}/C_{r}\right|&\left(\alpha>0\right)\\\
\omega+C_{i}\varepsilon^{2}/\left(2\left|\alpha\right|V\right)&\left(\alpha<0\right)\end{cases}$
(8)

###  II.2 Transport Efficiency and Phase Sensitivity

The purpose of the present work is to figure out a way to improve the
estimation of energy dissipation (or entropy production) related to the
stochastic oscillations. As mentioned in the introduction, one usually uses
the thermodynamic uncertain relation (TUR) as an inference of the real energy
dissipation via $\Delta W\geq\Delta W_{\text{TUR}}\equiv 2k_{B}T\left\langle
R\right\rangle^{2}/\text{Var}(R)$ wherein $R$ is some well-defined current
variable, and $Var\left(R\right)=\left\langle R^{2}\right\rangle-\left\langle
R\right\rangle^{2}$ denotes the variance of $R$. Correspondingly, the
transport efficiency for $R$ reads $\eta\left(R\right)=\Delta
W_{\text{TUR}}/\Delta W=2k_{B}T\left\langle
R\right\rangle^{2}/\text{Var}(R)\leq 1$. For the oscillatory dynamics
considered here, it is convenient to choose $R$ as the change of phase angle
within a given time interval $\left(0,\tau\right)$ , i.e.,
$R\left(\tau\right)\to\theta\left(\tau\right)=\int_{0}^{\tau}\dot{\theta}\left(t\right)dt$.
By simply rewriting and setting $k_{B}T=1$ from now on, the transport
efficiency can be expressed as
$\eta(\theta)=2k_{B}T\left\langle\theta\right\rangle^{2}/\text{$[$Var}(\theta)\Delta
W]=v_{\theta}^{2}/D_{\theta}\dot{W}$ where
$v_{\theta}=\lim_{t\to\infty}\left\langle\theta\right\rangle/t$ is the phase
speed,
$D_{\theta}=\lim_{t\to\infty}(\left\langle\theta^{2}\right\rangle-\left\langle\theta\right\rangle^{2})/2t$
is the phase diffusion constant, and $\dot{W}=\lim_{t\to\infty}\Delta W/t$ is
the dissipation rate.

By using the SNFT, the mean and variance of the phase
$\theta(\tau)=\int_{0}^{\tau}\dot{\theta}dt$ can be calculated as
$\left\langle\theta(t)\right\rangle\simeq\omega_{s}t$ and
$\left\langle(\theta(t)-\left\langle\theta(t)\right\rangle)^{2}\right\rangle\approx\varepsilon^{2}t/Vr_{m}^{2}$.
Hence the velocity $v_{\theta}$ is simply $\omega_{s}$ and the phase diffusion
constant is given by $D_{\theta}\simeq\varepsilon^{2}/2Vr_{m}^{2}.$ It is also
possible to obtain the theoretical expression for $\dot{W}$ by using the SNFT,
which is after some manipulation given by
$\dot{W}\simeq\left(L_{12}-L_{21}\right)V\omega_{s}^{2}r_{m}^{2}$, where
$L_{12}$ and $L_{21}$ are model-dependent parameters determined by the linear
transformation of $F\left(\bm{x}\right)$ at the fixed point $\bm{x}_{s}$ (see
Appendix A for more details), and being independent of the control parameter
$\alpha$ and system size. Consequently, the transport efficiency reads

$\eta_{\theta}\simeq\frac{\omega_{s}}{\varepsilon^{2}\left(L_{12}-L_{21}\right)}=\frac{\omega+\alpha\left|C_{i}/C_{r}\right|}{\varepsilon^{2}\left(L_{12}-L_{21}\right)}$
(9)

and the TUR asserts that $\eta_{\theta}\leq 1$. Although the expression of
$\eta_{\theta}$, Eq.(9), gives no hint that the TUR holds, we indeed
demonstrate numerically in our previous work that for the well-known
Brusselator model Cao et al. (2020), $\eta_{\theta}\sim 0.4$ which is far
below the upper bound 1.0 in the vicinity of the Hopf bifurcation.

For oscillation systems, another important quantity is the phase sensitivity
quantifying the ability of the biochemical circuits to respond to external
signal (Hasegawa and Arita, 2014b, a). Instead of dealing with the entire
system, we employ the phase reduction method (Kuramoto, 2003; Goldobin et al.,
2010) to reduce the whole state space to a single phase variable $\phi$
characterizing the timing of oscillation, and the phase sensitivity $\kappa$
can be obtained by comparing the phase shift after perturbations. The phase
$\phi$ in Eq.(6) is defined on the limit cycle of the unperturbed
oscillations, and the definition can be expanded into the entire
$\bm{x}$-space by introducing the isochron (the two states are assigned the
same phase if trajectories originated from two states converge onto the limit
cycle at the same time). Following this definition, the deterministic phase
evolution equation can be expressed as
$\dot{\phi}=\Omega=\nabla_{\bm{x}}\phi\cdot\bm{F}(\bm{x})$. For a weak
external signal $\bm{\beta}(t)$, the deterministic term reads
$\bm{F}_{\alpha}(\bm{x})=\bm{F}(\bm{x})+k\bm{\beta}(t)$ with $k$ the control
parameter, and the phase shift incurred by a parametric perturbation $k\to
k+\delta k$ can be obtained as $\dot{\phi}=\Omega+\delta
k\left[\nabla_{\bm{x}}\phi\cdot\bm{\beta}(t)\right]$. Then, the global phase
sensitivity parameter $\kappa$ can be defined as the normalized value of
signal-independent factor $\nabla_{\bm{x}}\phi$ along the limit cycle with
$r=r_{m}$. For oscillations near the Hopf bifurcation with
$\sqrt{-2C_{r}\varepsilon^{2}/V}<\left|\alpha\right|\ll\left|C_{r}/C_{i}\right|$,
the phase sensitivity $\kappa$ can be approximately calculated as
$\kappa\approx\partial\omega_{s}/\partial\alpha$ (Fei et al., 2018; Cao et
al., 2020), i.e. (see Eq.8)

$\kappa=\begin{cases}\frac{\left|C_{i}\right|\varepsilon^{2}}{2\alpha^{2}V}&\alpha<0\\\
\left|\frac{C_{i}}{C_{r}}\right|&\alpha>0\end{cases}.$ (10)

###  II.3 Pearson Correlation Coefficient

Here, we investigate the formulation of a scheme for the characterization of
correlations between the oscillatory amplitude and phase based on a
statistical measure known as the Pearson correlation coefficient, which has
been commonly used in the context of quantum entanglement (Maccone et al.,
2015; Pozsgay et al., 2017; Jebarathinam et al., 2020) and filtering theorem
(Benesty et al., 2009, 2008). The Pearson correlation coefficient for any two
random variables $R$ and $Q$ is defined as
$\chi(R,Q)=\text{Cov}(R,Q)/\sqrt{\text{Var}(R)\text{Var}(Q)}$ with
$\text{Cov}(R,Q)=\left\langle RQ\right\rangle-\left\langle
R\right\rangle\left\langle Q\right\rangle$ the covariance. The values of
Pearson correlation coefficient lie between $-1$ and $1$.

Then, we start to calculate the Pearson correlation coefficient
$\chi_{n}^{2}=\chi^{2}(r^{n},\theta)$ between the two observables
$R(\tau)=\theta(\tau)=\int_{0}^{\tau}\dot{\theta}(t)dt$ and
$Q_{r,n}=\int_{0}^{\tau}r^{n}\left(t\right)dt$, where the exponent $n$
quantifies the order of correlation between the oscillatory amplitude and
phase. By using the SNFE, we find that the change rate of the covariance,
$\lim_{\tau\to\infty}\frac{1}{\tau}\text{Cov}(r^{n},\theta;\tau)=\left\langle\left\langle
r^{n}\dot{\theta}\right\rangle\right\rangle_{ss}-\left\langle\left\langle\dot{\theta}\right\rangle\right\rangle_{ss}\left\langle\left\langle
r^{n}\right\rangle\right\rangle_{ss}$, between oscillatory phase and amplitude
is related to the higher-order moment of the amplitude as

$\displaystyle\lim_{\tau\to\infty}\frac{1}{\tau}\text{Cov}(r^{n},\theta;\tau)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\int_{0}^{\infty}drr^{n}\dot{\theta}P_{ss}(r)$
$\displaystyle\hphantom{}\hphantom{}-(\omega+C_{i}\left\langle
r^{2}\right\rangle_{ss})\left\langle r^{n}\right\rangle_{ss}$
$\displaystyle\approx C_{i}\left(\left\langle
r^{n+2}\right\rangle_{ss}-\left\langle r^{n}\right\rangle_{ss}\left\langle
r^{2}\right\rangle_{ss}\right).$ (11)

Particularly, we choose $n=2$ to calculate the covariance between oscillatory
phase and amplitude. According to Eq.(6), the change rate of covariance is
(see Appendix B for detailed derivation)

$\lim_{t\to\infty}\frac{1}{\tau}\text{Cov}(r^{2},\theta;\tau)=\begin{cases}-\frac{2C_{i}C_{r}r_{m}^{4}\varepsilon^{2}}{\alpha^{2}V}&\alpha>0\\\
0&\alpha<0\end{cases}.$ (12)

Here, we need to emphasize that our theoretical expression for normal
oscillations ($\alpha>0$) holds in the region near the Hopf bifurcation where
$\sqrt{-2C_{r}\varepsilon^{2}/V}<\alpha\ll\left|C_{r}/C_{i}\right|$. It can be
found that the phase and amplitude are highly decoupled with the covariance
$\lim_{t\to\infty}\text{Cov}(r^{n},\theta;\tau)/\tau\approx 0$ in the
subcritical region ($\alpha<0$). The highly decoupling feature is also the
reason why the sensitivity for noise-induced oscillations ($\kappa\sim
V^{\delta},\delta=-1$) is typically smaller than the normal oscillations
($\kappa\sim V^{\delta},\delta=0$), i.e., the oscillatory amplitude’s
adaptation to phase shift incurred by perturbation is much slower in the
subcritical region.

Then, we start to calculate the Pearson correlation coefficient
$\chi_{2}^{2}=\chi^{2}(r^{2},\theta)$, which reads as

$\displaystyle\chi_{2}^{2}$
$\displaystyle=\frac{\left[\text{Cov}(r^{2},\theta)\right]^{2}}{\text{Var}(r^{2})\text{Var}(\theta)}$
$\displaystyle\approx\frac{C_{i}^{2}\left|\left\langle
r^{4}\right\rangle_{ss}-\left\langle
r^{2}\right\rangle_{ss}^{2}\right|}{D_{\theta}}.$ (13)

From Eqs. (12) and (13), the Pearson correlation coefficient $\chi_{2}^{2}$
for normal oscillations ($\alpha>0$) can be obtained as

$\chi_{2}^{2}=\begin{cases}2\alpha\left(\frac{C_{i}}{C_{r}}\right)^{2}&\alpha>0\\\
0&\alpha<0\end{cases}.$ (14)

which is independent of the system size $V$. The Pearson correlation
coefficient $\chi_{n}^{2}=\chi^{2}(r^{n},\theta)$ for $n\neq 2$ can be
calculated numerically.

###  II.4 Improved TUR

Recently, it was proposed by Dechant and Sasa that increasing the number of
observables will achieve tighter bounds than the conventional TUR (Dechant and
Sasa, 2021b; Dechant, 2018). To be precise, they defined a generalized
transport efficiency as $\eta(R,Q)=\eta(R)+\chi^{2}(R,Q)$, where
$Q=\int_{0}^{\tau}dtq(\bm{x},t)$ is the time-integral of a state-dependent
(non-current) observable $q(\bm{x},t)$, and $\chi(R,Q)$ is the Pearson
correlation coefficient between $Q$ and the current observable $R$.
Interestingly, they found that the generalized transport efficiency,
$\eta(R,Q)$, is also smaller than $1$ just like the conventional one,
$\eta(R)$. Therefore, the generalized transport efficiency
$\eta(R,Q)\geq\eta(R)$ provides an improved estimator for energy dissipation
than the conventional one,

$\displaystyle\Delta W_{\text{TUR}}$ $\displaystyle=\frac{2k_{B}T\left\langle
R\right\rangle^{2}}{\text{Var}(R)}$ $\displaystyle\leq\Delta
W_{\text{I}}=\frac{2k_{B}T\left\langle
R\right\rangle^{2}}{\text{Var}(R)[1-\chi^{2}(R,Q)]}$ $\displaystyle\leq\Delta
W.$ (15)

It can be found that how much the estimation can be improved is directly
related to the value of Pearson correlation coefficient between the chosen
observables, and the two observables we chose above, the oscillatory amplitude
and oscillatory phase, meet the conditions of use.

Based on Eq.(14) and (15), the explicit expression for the generalized
transport efficiency can be obtained as
$\eta(r^{2},\theta)=\eta(\theta)+2\alpha(C_{i}/C_{r})^{2}\leq 1$ when
$\alpha>0$. Therefore, we eventually get an efficiency-sensitivity trade-off
relation for normal oscillations,

$\eta(\theta)+2\alpha\kappa^{2}=\frac{v_{\theta}^{2}}{D_{\theta}\dot{W}}+2\alpha\kappa^{2}\leq
1,$ (16)

which is the main result of our paper, showing that both phase accuracy
$D_{\theta}^{-1}$ and phase sensitivity $\kappa$ can be improved
simultaneously only by increasing the energy dissipation rate $\dot{W}$
without sacrificing the phase speed $v_{\theta}$ (Hasegawa and Van Vu, 2019a).
More importantly, such trade-off relation provides an improved estimator for
the dissipation rate,

$\Delta
W_{\text{I,2}}=\frac{2k_{B}T\left\langle\theta\right\rangle^{2}}{\text{Var}(\theta)(1-\chi_{2}^{2})},$
(17)

than the conventional TUR, and the improvement of it is

$\frac{\Delta W_{\text{I,2}}}{\Delta
W_{\text{TUR}}}=\frac{\dot{W}_{\text{I,2}}}{\dot{W}_{\text{TUR}}}=\frac{1}{1-2\alpha\kappa^{2}}>1$
(18)

with the TUR estimator $\dot{W}_{\text{TUR}}=\lim_{\tau\to\infty}\Delta
W_{\text{TUR}}/\tau$ and the improved estimator
$\dot{W}_{\text{I,2}}=\lim_{\tau\to\infty}\Delta W_{\text{I,2}}/\tau$.

Several conclusions can be obtained as follows. Firstly, according to the
trade-off relation Eq.(16), it can be found that the precision of dissipation
inference will be further improved by enhancing the phase sensitivity of
biochemical oscillations. In actual experimental design, a feasible strategy
to achieve a higher phase sensitivity of the networks is to enhance the phase-
amplitude coupling strength $C_{i}$ by maximizing the net flux of the phase-
advancing pathway relative to that of the phase-retreating pathway (Fei et
al., 2018). Thus, we believe that our analyses provide realizable guidelines
for improving the precision of dissipation estimation for biochemical
oscillations.

Secondly, we find that the generalized transport efficiency
$\eta(r^{2},\theta)=\eta(\theta)+2\alpha(C_{i}/C_{r})^{2}$ is independent of
the system size $V$. Since the magnitude of the internal noise is proportional
to $V^{-1/2}$, it can be revealed that our improved scheme is not negatively
affected by the internal noise in the system.

Thirdly, since the phase and amplitude are highly decoupled in subcritical
region ($\alpha<0$), such scheme cannot be applied to improve the estimation
of the energy dissipation for noise-induced oscillations.

In the following, we further highlight our motivation. As stated above, it has
been proposed that the TUR provides a powerful tool to estimate energy
dissipation. Recently, this bound has been optimized to provide a more
accurate estimation Li et al. (2019); Manikandan et al. (2020); Otsubo et al.
(2020); Van Vu et al. (2020); Busiello and Pigolotti (2019); Manikandan and
Krishnamurthy (2018); Gingrich et al. (2017); Kim et al. (2020); Busiello and
Fiore (2022) and even realize equality Dechant and Sasa (2021a); Manikandan et
al. (2020); Otsubo et al. (2020); Van Vu et al. (2020), which is of great
significance. Particularly, Manikandan $et$ $al.$ have found that the TUR
estimates entropy production exactly in the very short time limit, if the
observed current is optimally chosen, which provides a powerful strategy for
the dissipation inference Manikandan et al. (2020). Some optimization
procedure needs to be utilized to obtain the optimal current, where the
similar manipulations have been used to get the hyperaccurate currents Li et
al. (2019); Busiello and Pigolotti (2019); Busiello and Fiore (2022). Other
techniques such as the gradient ascent in machine learning have also been
applied to construct the short-time limit TUR estimator by Otsubu $et$ $al.$
Otsubo et al. (2020) and Vu $et$ $al.$ Van Vu et al. (2020). However, the
related procedure may be difficult to follow than measuring the dissipation
itself Dechant and Sasa (2021b). To be specific, in Ref Manikandan et al.
(2020), Manikandan $et$ $al.$ tested their inference scheme by numerically
calculating the optimal current rely on linear combinations of the basis. In
Ref Otsubo et al. (2020), Otsubu $et$ $al.$ demonstrated that their learning
protocol performs well by numerical experiments in nonlinear Langevin
dynamics. In addition to the TUR-based approach, Frishman and Ronceray
proposed a principled method, stochastic force inference, to evaluate the
corresponding entropy production based on approximating force fields and
diffusion coefficients Frishman and Ronceray (2020). To sum up, all the
dissipation estimations mentioned above require some specific preprocessing,
thus increasing the statistical effort. As a comparison, only the oscillatory
phase and amplitude need to be tracked by using the scheme proposed by us,
which is readily accessible, showing its benefit for experimentally
application. Also, the results of our manuscript have demonstrated that the
estimation of energy dissipation is considerably improved by considering the
correlations between observables, no matter how far from equilibrium the
system is. In addition, since biochemical functional benefits from operating
at the edge of instability, studying the oscillatory behavior near Hopf
bifurcation points can bring general inspiration Munoz (2018). Therefore, we
believe that our formulations provide an efficient estimator in terms of
experimentally accessible quantities.

Figure 1: (a) Pearson correlations $\chi_{2}^{2}$ as a function of the control
parameter $\alpha$ for the circadian clock model. The value of $\chi_{2}^{2}$
changes sharply near the critical point $\alpha=0$, due to the bifurcation
phenomenon. Line: theory. Dots: simulation. (b) The Pearson correlations
$\chi_{2}^{2}$ as a function of the phase sensitivity
$\kappa=\left|C_{i}/C_{r}\right|$ (green and pink dots). The slopes $\nu$ for
$\chi_{2}^{2}\propto\kappa^{\nu}$ have been calculated from fitting the
numerical data. (c) The slope $\nu$ and goodness $R^{2}$ for the linear fit
between $\ln(\chi_{2}^{2})$ and $\ln\kappa$ as a function of the control
parameter $\alpha$. The gray circle and blue dotted line represent the range
in which the scaling behavior $\chi_{2}^{2}\propto\kappa^{2}$ holds. The
establishment of the scaling behavior reveals how far from the Hopf
bifurcation the trade-off relation between dissipation and phase sensitivity
satisfies. The system size $V=1.6\times 10^{5}$.

##  III Simulations

Figure 2: (a) Estimations of the dissipation rate as a function of the exact
dissipation rate $\dot{W}$ for $\alpha>0$. It can be observed that the
improved estimator $\dot{W}_{\text{I,2}}$ outperforms the TUR estimator
$\dot{W}_{\text{TUR}}$. The parameters of circadian clock model can be found
in Appendix C. (b) The improvement of the tighter bound
$\dot{W}_{\text{I,2}}/\dot{W}_{\text{TUR}}$ as a function of the phase
sensitivity $\kappa=\left|C_{i}/C_{r}\right|$. The numerical results (dots)
verify our analytical expressions (lines). The system size $V=1.6\times
10^{5}$.

In this section, we illustrate the formal analytical results of the above
section within numerical simulations of the circadian clock model (Hou and
Xin, 2003), describing how living systems keep an internal sense of time. The
circadian clock model considered here incorporates the transcription of the
gene (G) involved in the biochemical clock and transport of the mRNA (R) into
the cytosol where it is translated into clock proteins ($\text{P}_{C}$) and
degraded. The protein can be degraded or transported into the nucleus
($\text{P}_{N}$) where it exerts a negative regulation on the expression of
its gene. For the parameters we examine (see Appendix C for details), the Hopf
bifurcation point locates at $v_{s}\simeq 0.25725$ with $v_{s}$ the
transcription rate of mRNA. In addition, parameter values used in the
stochastic normal form theory can be calculated from simulations as
$C_{r}\simeq-0.3474$, $C_{i}\simeq 0.5722$ and $\varepsilon^{2}\simeq 0.3556$.
By adjusting the transition rates, the values of $C_{i}$ and $C_{r}$ will
change, and can also be obtained.

By using the Euler methods, we numerically calculate Eqs.(4) and (5) with a
time step of $0.002$. Generally, after a long time $t_{ss}=10^{5}$ to ensure
the system reaches the steady state, $2\times 10^{5}$ trajectories with the
length $t_{0}=1$ are used to get the Pearson correlation coefficient,
$\chi_{n}^{2}=\chi^{2}(r^{n},\theta)=\left[\text{Cov}(r^{n},\theta)\right]^{2}/\text{Var}(r^{n})\text{Var}(\theta)$
for $Q_{r,n}=\int_{0}^{\tau}r^{n}\left(t\right)dt$ with $n>0$ the power, which
yields the corresponding improved estimator as

$\Delta
W_{\text{I,n}}=\frac{2k_{B}T\left\langle\theta\right\rangle^{2}}{\text{Var}(\theta)(1-\chi_{n}^{2})}.$

Then, the TUR estimator $\dot{W}_{\text{TUR}}=\lim_{\tau\to\infty}\Delta
W_{\text{TUR}}/\tau$ and the improved estimator
$\dot{W}_{\text{I,n}}=\lim_{\tau\to\infty}\Delta W_{\text{I,n}}/\tau$ can be
obtained numerically. The corresponding improvement reads as

$\frac{\dot{W}_{\text{I,n}}}{\dot{W}_{\text{TUR}}}=\frac{1}{1-\chi^{2}(r^{n},\theta)}>1.$

On the other hand, the exact dissipation rate $\dot{W}$ is obtained from the
simulation data of Eq.(2) (see Appendix A for details).

In Fig.1(a), the dependence of the Pearson correlations $\chi_{n}^{2}$
($n=1,2,3$) on the control parameter $\alpha$ are depicted for the circadian
clock model. The value range of the control parameter $\alpha$ ensures the
establishment of the SNFT. For noise-induced oscillations in the subcritical
region($\alpha<0$), the Pearson correlation coefficients are almost zero, and
they increase significantly after the control parameter crossing the critical
point $\alpha=0$ to reach the supercritical region for normal oscillations
($\alpha>0$). Those results verify our prediction that the correlation between
oscillatory phase and amplitude is highly decoupled and not sufficient to
improve the estimation of energy dissipation for noise-induced oscillations.
In addition, we notice that numerical results (dots) of the Pearson
correlations are in good agreement with our theoretical predictions, Eq.(14)
(line). Further, in Fig.1(b), we plot the Pearson correlations $\chi_{2}^{2}$
as a function of the phase sensitivity $\kappa=\left|C_{i}/C_{r}\right|$. The
scaling behaviors are consistent with our analytical result
$\chi_{2}^{2}\propto\kappa^{2}$, further confirming our theory. It can be
found that the slope $\nu$ for $\chi_{2}^{2}\propto\kappa^{\nu}$ is closer to
the analytical prediction $\nu=2$ for smaller $\alpha$, showing that our
theory is more accurate for near Hopf bifurcation region. To further explore
the extent to which our formulation holds generally away from a Hopf
Bifurcation, we plot the slope $\nu$ and goodness $R^{2}$ for the linear fit
between $\ln(\chi_{2}^{2})$ and $\ln\kappa$ as a function of the control
parameter $\alpha$ in Fig.1(c). The gray circle and blue dotted line represent
the range in which the scaling behavior $\chi_{2}^{2}\propto\kappa^{2}$ holds
($R^{2}\geq 0.97$). The establishment of the scaling behavior reveals how far
from the Hopf bifurcation the trade-off relation between dissipation and phase
sensitivity satisfies. For the circadian clock model, our formulation holds
for $\alpha\leq 0.055$.

In Fig.2(a), we show both the TUR estimator $\dot{W}_{\text{TUR}}$ and the
improved estimator $\dot{W}_{\text{I}}$ for the circadian clock model to
demonstrate how much the estimation of energy dissipation can be improved. The
conventional TUR, while a commonly used dissipation estimator, only provides a
trivial bound with the 0.6 efficiency, and the improved estimator
$\dot{W}_{\text{I}}$ is much closer to the exact value $\dot{W}$. On the other
hand, it can be found that $\dot{W}\geq\dot{W}_{\text{I}}$, which verifies the
efficiency-sensitivity trade-off relation we proposed [Eq.(16)]. In Fig.2(b),
we show the relationship between the improvement of the tighter bound
$\dot{W}_{\text{I}}/\dot{W}_{\text{TUR}}$ and phase sensitivity
$\kappa=\left|C_{i}/C_{r}\right|$. The theoretical predictions
$\dot{W}_{\text{I}}/\dot{W}_{\text{TUR}}=\frac{1}{1-2\alpha(C_{i}/C_{r})^{2}}$
are in good agreement with the numerical results, which demonstrates that the
estimation of the dissipation can be improved by enhancing the phase
sensitivity of biochemical oscillations.

Figure 3: (a) The improvements of the dissipation estimation
$\dot{W}_{\text{I}}/\dot{W}_{\text{TUR}}$ ($n=1,2,3$) as a function of the
system size $V$ for the circadian clock model in normal oscillation region
($\alpha>0$). The control parameter $\alpha=0.05$. (b) The improvements of the
dissipation estimation $\dot{W}_{\text{I}}/\dot{W}_{\text{TUR}}$ as a function
of the power of the amplitude, $n$, for the circadian clock model. The values
are independent of the power. The error bars represent the standard deviation
obtained from independent trials. The system size $V=1.6\times 10^{5}$.

In Fig.3 (a), we have shown that the improvements of the dissipation
estimation $\dot{W}_{\text{I,n}}/\dot{W}_{\text{TUR}}$ change little with the
system size $V$ in normal oscillations ($\alpha>0$). Moreover, we numerically
test whether the power of amplitude observables $n$ affect the improvement
$\dot{W}_{\text{I,n}}/\dot{W}_{\text{TUR}}$ in details. In Fig.3 (b), it can
be observed that values of $\dot{W}_{\text{I,n}}/\dot{W}_{\text{TUR}}$ change
little for different choices of the power of the amplitude observable
$Q_{r,n}(\tau)=\int_{0}^{\tau}r^{n}(t)dt$.

To further demonstrate the broad application of the proposed improved
estimation, we have also applied them to another well-known biochemical
oscillation system, the Brusselator model. Other details of the model and
parameters can be found in Appendix C. As shown in Fig.4, we find that our
main results, such as the theoretical expression of Pearson correlations [Eq.
(14)] and the improved estimation [obtained from Eq. (15)], still hold in the
Brusselator model. In Fig.4(a), it can be observed that numerical results
(dots) of the Pearson correlations are in good agreement with our theoretical
predictions, Eq.(14) (line). In Fig.4(b), both the TUR estimator
$\dot{W}_{\text{TUR}}$ and the improved estimator $\dot{W}_{\text{I}}$ has
been depicted for the Brusselator model. As expected, the conventional TUR
estimator only yields a loose bound, and our improved estimator
$\dot{W}_{\text{I}}$ is much more accurate.

Figure 4: (a) Pearson correlations $\chi_{2}^{2}$ as a function of the control
parameter $\alpha$ for the Brusselator model. The value of $\chi_{2}^{2}$
changes sharply near the critical point $\alpha=0$, due to the bifurcation
phenomenon. Line: theory. Dots: simulation. (b) Estimations of the dissipation
rate as a function of the exact dissipation rate for $\alpha>0$. It can be
observed that the improved estimator $\dot{W}_{\text{I,2}}$ outperforms the
TUR estimator $\dot{W}_{\text{TUR}}$. The parameters of the Brusselator model
can be found in Appendix C. The system size $V=1.6\times 10^{5}$.

##  IV Discussion

In this paper, we proposed an improved estimation for the energy dissipation
of biochemical oscillations by using the Pearson correlations between
oscillatory phase and amplitude, which are easily accessible in experimental
observations. Both the analytical and numerical results demonstrate that such
scheme can be further improved by enhancing the phase sensitivity of systems.
In addition, it has been revealed by us that the validity of our scheme is
independent of the system size and the power of oscillatory amplitude.

In our previous work (Cao et al., 2020), we have found that the dissipation
rate $\dot{W}\sim V^{\gamma}$, with $\gamma=1$ for supercritical region
($\alpha>0$), $\gamma=1/2$ for the critical point ($\alpha=0$) and $\gamma=0$
for subcritical region ($\alpha<0$), showing that biochemical oscillations
have a much lower energy dissipation for noise-induced oscillation.
Intuitively, one might think that less dissipation will lead to a easier
estimation, however, the estimator $\Delta W_{\text{I,n}}$ introduced by us is
not applicable for improving the estimation of dissipation due to the highly
decoupling of the phase and amplitude, which is deserved for further study.

Biomolecules, especially proteins, can act as tiny and highly functional
machines, such as kinesin Verhey et al. (2011) and ribosome Aitken et al.
(2010). To probe the operation of these bimolecular machines, it is not enough
just to know their structure, one needs to understand how the structure
generates specific conformational dynamics. Meanwhile, how much energy the
biological machine dissipates to perform certain functions is also a major
issue Lan et al. (2012). Particularly, molecular dynamics (MD) simulation is a
primary technique for studying bimolecular machines, producing information
about the conformational dynamics with spatial and temporal resolutions. For
the molecular systems, the key ingredients of the slow kinetics can be
obtained by using SNFT based on stochastic averaging Cao et al. (2020) or
variational approach with MD simulations Nuske et al. (2014)(Perhaps some
short time dynamics with energy input). Then, the findings of this manuscript
that correlations between observables can yield improved estimation could be
applied in these systems based on the accurate capture of the dominant motion.
However, it is still important to note that fast processes may also have
finite correction to the thermodynamic quantities Seifert (2019). As
stochastic normal form equations can be extended to other oscillatory systems
related to other types of bifurcations, such as relaxation oscillations, we
believe that our scheme may have a wider range of applications.

###### Acknowledgements.

This work is supported by MOST(2018YFA0208702), NSFC (32090044, 21790350,
21521001).

## Author declarations

The authors have no conflicts to disclose.

## Data availability

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

##  Appendix A The derivation of normal form and calculation of transport
efficiency

In this section, we introduce the derivation of stochastic normal form
equation and the calculation of conventional transport efficiency for self-
consistency.

###  A.1 Stochastic normal form theory

Firstly, we assume that the deterministic form of the chemical Langevin
equation, Eq.(2), has a unique stable point $\bm{x}_{s}$ with
$\bm{F}(\bm{x}_{s})\equiv 0$, which loses stability at the supercritical HB
$\mu=\mu_{c}$, where $\mu$ is the control parameter. Based on the Hopf theorem
(Hassard et al., 1981), the Jacobian matrix $\bm{J}$, whose components
$J_{ij}=\left(\partial f_{i}/\partial x_{j}\right)|_{\bm{x}=\bm{x}_{s}}$, has
a pair of conjugate eigenvalues $\lambda_{\pm}=\alpha(\mu)\pm i\omega$ with
$\alpha(\mu_{c})=0$. The other $N-2$ eigenvalues of $\bm{J}$,
$-\lambda_{j(\geq 3)}$, all have negative real parts with absolute values
considerably larger than $0$. Performing the variable transformation
$\bm{u}=\bm{T}^{-1}(\bm{x}-\bm{x}_{s})$, the linear part of Eq.(2) can be
transformed to Jordan form as
$\dot{\bm{u}}=\bm{\Lambda\bm{u}}+O(\bm{u}^{2})+\frac{1}{\sqrt{V}}\bm{\eta}\left(t\right),$where
$\bm{\Lambda}=\left(\begin{array}[]{cc}\alpha&-\omega\\\
\omega&\alpha\end{array}\right)\oplus\text{diag$\left(-\lambda_{1},\dots,-\lambda_{N}\right)$}$
and $\bm{\eta}=\bm{T}^{-1}\zeta(\bm{x}_{s},t)$ with
$\zeta(\bm{x},t)=\sum_{\rho}v_{\rho}^{j}\sqrt{w_{\rho}(\bm{x})}\xi_{\rho}(t)$.
The variances of $\bm{\eta}$ are
$\left\langle\eta_{i}(t)\eta_{j}(s)\right\rangle=2D_{ij}\delta(t-s)$ with
$\bm{D}=\bm{T}^{-1}\bm{G}(\bm{T}^{-1})^{\text{T}}$. The transformation is done
as follows. Firstly, we calculate the eigenvector $\bm{u}_{+}$ whose
eigenvalue $\lambda_{+}=\alpha+i\omega$, and normalize it to ensure the first
non-vanishing component is $1$. Secondly, we construct a matrix
$\bm{T}=(Re\bm{u}_{+},-Im\bm{u}_{+},\bm{r}_{3},\ldots,\bm{r}_{n})$ with
$(\bm{r}_{3},\ldots,\bm{r}_{n})$ are any set of real vectors which span the
union of the eigenspaces for $(\lambda_{3},\ldots,\lambda_{n})$. Finally, it
is allowed to perform the change of variables
$\bm{x}=\bm{x}_{s}+\bm{T}\bm{u}$.

When the system locates near the HB ($\left|\alpha\right|\ll 1$), the
evolution of the oscillatory mode related to $(u_{1},u_{2})$ is much slower
than the other $N-2$ stable modes due to the time-scale separation. Hence, the
system’s dynamics will be dominated by the slow motion on a 2D center manifold
spanned by the eigenvectors of $\lambda_{\pm}$. The oscillatory mode are ruled
by a normal form equation involving the time evolution of a complex variable
$Z=u_{1}+iu_{2}$, or a pair of coupled equations for the oscillation amplitude
$r$ and phase $\theta$ via $Z=re^{i\theta}$.We follow the standard procedure
to get the normal form,

$\displaystyle\frac{dZ}{dt}$
$\displaystyle=\left(\alpha+i\omega\right)Z+\left(C_{r}+iC_{i}\right)\left|Z\right|^{2}Z$
$\displaystyle\hphantom{\hphantom{}}\hphantom{}+\frac{1}{\sqrt{V}}\sum_{\rho}\left(v_{1\rho}^{\prime}+iv_{2\rho}^{\prime}\right)\sqrt{w_{\rho}}\xi_{\rho},$
(19)

where $\bm{v}_{j\rho}^{\prime}=\left(T^{-1}\bm{v}\right)_{j\rho}$, i.e.,

$\frac{dr}{dt}=\left(\alpha
r+C_{r}r^{3}\right)+\frac{1}{\sqrt{V}}\sum_{\rho}\chi_{r\rho}\circ\xi_{\rho},$
(20)

$\frac{d\theta}{dt}=\left(\omega+C_{i}r^{2}\right)+\frac{1}{\sqrt{V}}\sum_{\rho}\chi_{\theta\rho}\circ\xi_{\rho}$
(21)

with

$\chi_{r\rho}=\left(v_{1\rho}^{\prime}\cos\theta+v_{2\rho}^{\prime}\sin\theta\right)\sqrt{w_{\rho}},$
(22)
$\chi_{\theta\rho}=\frac{1}{r}\left(-v_{1\rho}^{\prime}\sin\theta+v_{2\rho}^{\prime}\cos\theta\right)\sqrt{w_{\rho}}.$
(23)

By using the “stochastic averaging” method (Arnold et al., 1996), the
following equation can be obtained

$\frac{dr}{dt}=\alpha
r+C_{r}r^{3}+\frac{K(r)}{V}+\frac{\varepsilon_{r}}{\sqrt{V}}\xi_{r},$ (24)

and

$\frac{d\theta}{dt}=\omega+C_{i}r^{2}+\frac{K(\theta)}{V}+\frac{\varepsilon_{\theta}}{r\sqrt{V}}\xi_{\theta}.$
(25)

Here,

$K(r)=\frac{1}{2\pi}\sum_{\rho}\intop_{0}^{2\pi}d\theta\left(\chi_{r\rho}\partial_{r}\chi_{r\rho}+\chi_{\theta\rho}\partial_{\theta}\chi_{r\rho}\right),$
(26)

$K(\theta)=\frac{1}{2\pi}\sum_{\rho}\intop_{0}^{2\pi}d\theta\left(\chi_{r\rho}\partial_{r}\chi_{\theta\rho}+\chi_{\theta\rho}\partial_{\theta}\chi_{\theta\rho}\right),$
(27)

which is related to the coupling effects between amplitude and phase.
$\varepsilon_{r}^{2}=\frac{1}{2\pi}\sum_{\rho}\intop_{0}^{2\pi}d\theta\chi_{r\rho}^{2}$
and
$\varepsilon_{\theta}^{2}=\frac{1}{2\pi}\sum_{\rho}\intop_{0}^{2\pi}d\theta\chi_{\theta\rho}^{2}$
are the averaged noise intensities. The main purpose of this method is to
approximate the system’s dynamics as the Markovian stochastic process when the
system reaches the steady state. Further, by expanding the reaction rates,
$w_{\rho}=\sum_{i+j=0}^{n}w_{\rho}^{ij}(r\cos\theta)^{i}(r\sin\theta)^{j}$,
$K(\theta)$ is zero (Hou et al., 2006). Thus, the averaged noise intensities
read as

$\varepsilon_{r}^{2}=\varepsilon_{\theta}^{2}=\frac{1}{2}\sum_{\rho}\left[\left(v_{1\rho}^{\prime}\right)^{2}+\left(v_{2\rho}^{\prime}\right)^{2}\right]w_{\rho}^{00}$
(28)

near the Hopf bifurcation point, i.e., the stochastic normal form equation can
be obtained as

$\dot{r}=\alpha
r+C_{r}r^{3}+\frac{\varepsilon^{2}}{2Vr}+\frac{\varepsilon}{\sqrt{V}}\eta_{r}(t),$
(29)
$\dot{\theta}=\omega+C_{i}r^{2}+\frac{\varepsilon}{r\sqrt{V}}\eta_{\theta}(t),$
(30)

where the $i+j\geqslant 2$ terms are neglected.

###  A.2 Steady state dissipation rate and conventional transport efficiency

In order to obtain the transport efficiency, we start to calculate the steady
state dissipation rate $\dot{W}$. Based on the framework of stochastic
thermodynamics (Seifert, 2005, 2012; Sekimoto, 2010; Jarzynski, 2011; Gaspard,
2004), the entropy balance equation reads as
$\dot{s}_{tot}(\tau)=\dot{s}_{m}(\tau)+\dot{s}(\tau)$, where $s_{tot}(\tau)$
is the total entropy production, $s(\tau)$ is the Shannon entropy and
$s_{m}(\tau)$ is the entropy flux. As $s(\tau)=-\ln
p\left(\bm{x},\tau\right)$, the change rate of the Shannon entropy is

$\displaystyle\dot{s}(\tau)$
$\displaystyle=\left[-\partial_{\tau}p(\bm{x},\tau)+\frac{2V}{p(\bm{x},\tau)}\sum_{i,j}\Gamma_{ij}J_{j}|_{\bm{x}(\tau)}\dot{x}_{i}\right]$
$\displaystyle\hphantom{}\hphantom{}-V\sum_{i}H_{i}\dot{x}_{i},$ (31)

where $H_{j}=2\sum_{k}\Gamma_{jk}f_{k}^{\prime}$ ($\bm{\Gamma}=\bm{G}^{-1}$)
with $\widetilde{f}_{k}=f_{k}-1/(2V)\sum_{j}(\partial G_{kj})/(\partial
x_{j})$. Then, the entropy production rate and entropy flux rate can be
identified as
$\dot{s}_{tot}(\tau)=-\partial_{\tau}p(\bm{x},\tau)+\frac{2V}{p(\bm{x},\tau)}\sum_{i,j}\Gamma_{ij}J_{j}|_{\bm{x}(\tau)}\dot{x}_{i}$
and $\dot{s}_{m}(\tau)=V\sum_{i}H_{i}\dot{x}_{i}$. As
$\dot{s}(\tau)=\lim_{t\to\infty}\left\langle\Delta s\right\rangle/t$ vanishes
in the steady state, the averaged entropy production rate can be obtained as

$\dot{S}_{tot}=\lim_{t\to\infty}\left\langle\Delta
s_{m}\right\rangle/t=V\sum_{i}\left\langle\left\langle
H_{i}\dot{x}_{i}\right\rangle\right\rangle_{ss}$ (32)

with $\left\langle\left\langle\cdot\right\rangle\right\rangle_{ss}$ denotes
the average over time and steady state (Xiao et al., 2009).

By using the variable transform, the theoretical expression of the entropy
production rate can be calculated in terms of $\bm{u}$, which reads
$\dot{S}_{tot}=2V\left\langle\left\langle\widetilde{\bm{f}}^{\text{T}}\bm{\Gamma}^{\text{T}}\dot{\bm{x}}\right\rangle\right\rangle_{ss}$.
By approximating $\widetilde{\bm{f}}(\bm{x)}\approx\bm{J}\bm{T}\bm{u},$ the
entropy production reads

$\dot{S}_{tot}=2V\left\langle\left\langle\bm{u}^{\text{T}}\bm{L}\dot{\bm{u}}\right\rangle\right\rangle_{ss}=2V\sum_{i,j}L_{ij}h_{ij}$
(33)

with $h_{ij}=\left\langle\left\langle
u_{i}\dot{u}_{j}\right\rangle\right\rangle_{ss}$.
$\bm{L}=\bm{T}^{\text{T}}\bm{J}^{\text{T}}\bm{\Gamma}^{\text{T}}\bm{T}$ are
model-dependent parameters taken the value at the stable point $\bm{x}_{s}$.
Note that in the steady state, $\frac{d}{dt}\left\langle\left\langle
u_{i}\dot{u}_{j}\right\rangle\right\rangle_{ss}=0$, thus we have
$h_{ij}=-h_{ji}$. Then, we have that

$\displaystyle h_{12}$ $\displaystyle=-h_{21}=\left\langle\left\langle
r\cos\theta\frac{d}{dt}\left(r\sin\theta\right)\right\rangle\right\rangle_{ss}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\dot{\theta}\cos^{2}\theta
d\theta\cdot\int_{0}^{\infty}r^{2}p_{ss}\left(r\right)dr\approx\frac{1}{2}\omega_{s}\left\langle
r^{2}\right\rangle,$ (34)

where the time average is substituted by averaging over $\theta$ due to
dominant oscillatory mode. $\omega_{s}=\omega+C_{i}r_{m}^{2}$ is the effective
phase angular velocity. Meanwhile, for $j>2,$we have
$h_{1j}=\left\langle\left\langle
r\cos\theta\dot{u}_{j}\right\rangle\right\rangle_{ss}\approx 0$ and
$h_{2j}=\left\langle\left\langle
r\sin\theta\dot{u}_{j}\right\rangle\right\rangle_{ss}\approx 0$. For $i,j>2$,
one can obtain that $h_{ij}=\left\langle\left\langle
u_{i}\dot{u}_{j}\right\rangle\right\rangle_{ss}=(\lambda_{i}-\lambda_{j})D_{ij}/[(\lambda_{i}+\lambda_{j})V]$.
Therefore, the averaged entropy production rate is

$\dot{S}_{tot}=V\left(L_{12}-L_{21}\right)\omega_{s}\left\langle
r^{2}\right\rangle+2\sum_{i,j>2}L_{ij}D_{ij}\frac{\lambda_{i}-\lambda_{j}}{\lambda_{i}+\lambda_{j}}.$
(35)

Here, $r_{m}$ is the most probable value of the amplitude in the steady state
with $\partial_{r}P_{ss}\left(r\right)|_{r=r_{m}}=0$. By going through our
derivation, the contributions from the remaining other $N-2$ stable modes can
also be identified as
$\dot{S}_{fast}=2\sum_{i,j>2}L_{ij}D_{ij}(\lambda_{i}-\lambda_{j})/(\lambda_{i}+\lambda_{j})$,
which is absent in the expressions obtained by the conventional steady state
formula Cao et al. (2015); Fei et al. (2018). Further, the steady state
dissipation rate (here we set $k_{B}T=1$)

$\dot{W}=k_{B}T\dot{S}_{tot}\approx
V\left(L_{12}-L_{21}\right)\omega_{s}r_{m}^{2}$ (36)

Now, we start to calculate the transport efficiency $\eta_{\theta}$. The mean
and variance of the phase $\theta(\tau)=\int_{0}^{\tau}\dot{\theta}dt$ can be
calculated as $\left\langle\theta(t)\right\rangle\approx\omega_{s}t$ and
$\left\langle(\theta(t)-\left\langle\theta(t)\right\rangle)^{2}\right\rangle\approx\varepsilon^{2}t/Vr_{m}^{2}$,
and the phase diffusion constant is given by
$D_{\theta}=\lim_{t\to\infty}\left\langle(\theta(t)-\left\langle\theta(t)\right\rangle)^{2}\right\rangle/2t\approx\varepsilon^{2}/2Vr_{m}^{2}.$
The transport efficiency reads as

$\eta_{\theta}=\frac{v_{\theta}^{2}}{D_{\theta}\dot{W}}\approx\frac{2\omega_{s}}{\varepsilon^{2}\left(L_{12}-L_{21}\right)}.$
(37)

In the main text, we use the Eq.(32) allows us to numerically calculate the
exact dissipation rate $\dot{W}$ in Fig.2. Since $\dot{x}_{i}$ and $H_{i}$ can
be obtained from the dynamics generating from Eq.(2), $\dot{s}_{m}$ can be
calculated numerically. By averaging over trajectories in steady states,
$\dot{S}_{tot}$ and $\dot{W}$ can then be obtained.

##  Appendix B Pearson correlation coefficient

In this section, we calculate the Pearson correlation coefficient
$\chi^{2}(R,Q)$ between the phase
$R\left(\tau\right)=\theta\left(\tau\right)=\int_{0}^{\tau}\dot{\theta}\left(t\right)dt$
and the amplitude $Q\left(\tau\right)=\int_{0}^{\tau}r^{2}\left(t\right)dt$.
The change rate of the covariance,
$C\left(r^{2},\theta\right)=\lim_{t\to\infty}\frac{1}{\tau}\text{Cov}_{r^{2},\theta}\left(\tau\right)$,
of these two variables can be calculated as

$\displaystyle\lim_{\tau\to\infty}\frac{1}{\tau}\text{Cov}(r^{n},\theta;\tau)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\int_{0}^{\infty}drr^{n}\dot{\theta}P_{ss}(r)$
$\displaystyle\hphantom{}\hphantom{}-(\omega+C_{i}\left\langle
r^{2}\right\rangle_{ss})\left\langle r^{n}\right\rangle_{ss}$
$\displaystyle\approx C_{i}\left(\left\langle
r^{n+2}\right\rangle_{ss}-\left\langle r^{n}\right\rangle_{ss}\left\langle
r^{2}\right\rangle_{ss}\right).$ (38)

Note that the integrals (averages) we are going to calculate all take the form
$I_{n}=\int_{0}^{\infty}r^{2n}\exp[V\varepsilon^{-2}(\frac{\alpha}{2}r^{2}+\frac{C_{r}}{4}r^{4})]dr^{2}=\int_{0}^{\infty}x^{n}\exp[V\varepsilon^{-2}(\frac{\alpha}{2}x+\frac{C_{r}}{4}x^{2})]dx$.
By setting
$y=\rho(x/A-1)=\sqrt{-\frac{VC_{r}}{4\varepsilon^{2}}(x+\frac{\alpha}{C_{r}})}$
with $\rho=\alpha/2\sqrt{-C_{r}\varepsilon^{2}/V}$ and $A=-\alpha/C_{r}\approx
r_{s}^{2}$ (for $\alpha>0$), we have

$I_{n}=\exp\left(\rho^{2}\right)\left(r_{m}^{2}/\rho\right)^{n+1}\int_{-\rho}^{\infty}\left(y+\rho\right)^{n}e^{-y^{2}}dy.$
(39)

For $\alpha>0$, integrals
$\int_{-\rho}^{\infty}p_{n}\left(y\right)e^{-y^{2}}$($p_{n}\left(y\right)$ are
polynomials of degree $n$) can be obtained by simple Gaussian integrals
$\int_{-\infty}^{\infty}p_{n}\left(y\right)e^{-y^{2}}$, for $\alpha=0$
integrals read $\int_{0}^{\infty}p_{n}\left(y\right)e^{-y^{2}}$ and for
$\alpha<0$ integrals are approximately zero. For $\alpha>0$, the covariance
reads

$C\left(r^{2},\theta\right)\approx
C_{i}\left[r_{m}^{4}\left(1+\frac{1}{2\rho^{2}}\right)-\left(r_{m}^{2}\right)^{2}\right]=-\frac{2C_{i}C_{r}r_{m}^{4}\varepsilon^{2}}{\alpha^{2}V}.$
(40)

We need to emphasize that such equation holds in the region where
$\sqrt{-2C_{r}\varepsilon^{2}/V}<\alpha\ll\left|C_{r}/C_{i}\right|$ due the
above approximation. Thus, for normal oscillation region ($\alpha>0$) , the
Pearson correlation coefficient $\chi^{2}\left(r^{2},\theta\right)$ can be
calculated as

$\chi^{2}\left(r^{2},\theta\right)\approx 2\alpha(C_{i}/C_{r})^{2}<1,$ (41)

which means that such scheme works well for oscillators with high value of
$\left|C_{i}/C_{r}\right|$ (independent of the system size).

Generally, for $Q_{r,n}=\int_{0}^{\tau}r^{n}\left(t\right)dt$ with $n$ the
power of oscillatory amplitude, the change rate of the covariance between
oscillatory phase and amplitude is related to the higher-order moment of
amplitude as

$\displaystyle\lim_{t\to\infty}\frac{1}{\tau}\text{Cov}(r^{n},\theta;\tau)$
$\displaystyle=\left\langle\left\langle
r^{n}\dot{\theta}\right\rangle\right\rangle_{ss}-\left\langle\left\langle\dot{\theta}\right\rangle\right\rangle_{ss}\left\langle\left\langle
r^{n}\right\rangle\right\rangle_{ss}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\int_{0}^{\infty}drr^{n}\dot{\theta}P_{ss}(r)$
$\displaystyle~{}-(\omega+C_{i}\left\langle
r^{2}\right\rangle_{ss})\left\langle r^{n}\right\rangle_{ss}$
$\displaystyle\approx C_{i}(\left\langle
r^{n+2}\right\rangle_{ss}-\left\langle r^{n}\right\rangle_{ss}\left\langle
r^{2}\right\rangle_{ss}).$ (42)

and the Pearson correlation coefficient $\chi_{n}^{2}=\chi^{2}(r^{n},\theta)$
can be calculated as

$\displaystyle\chi_{n}^{2}$
$\displaystyle=\frac{\left[\text{Cov}(r^{n},\theta)\right]^{2}}{\text{Var}(r^{n})\text{Var}(\theta)}$
$\displaystyle\approx\frac{C_{i}^{2}(\left\langle
r^{n+2}\right\rangle_{ss}-\left\langle r^{n}\right\rangle_{ss}\left\langle
r^{2}\right\rangle_{ss})^{2}}{D_{\theta}(\left\langle
r^{2n}\right\rangle_{ss}-\left\langle r^{n}\right\rangle_{ss}^{2})}.$ (43)

##  Appendix C Details of the models

###  C.1 The circadian clock model

Figure 5: Typical trajectories for the circadian clock model. (a)
$\alpha=0.05$. (b) $\alpha=-0.01.$ The system size $V=1.6\times 10^{5}$.

Here, we describe the details of the circadian clock model studied in the main
text. The vector $\bm{x}=\left(x_{1},x_{2},x_{3}\right)$ stands for the
concentrations of $\left(\text{R},\text{P}_{C},\text{P}_{N}\right)$. The
transcription rate of mRNA is chosen as the control parameter, represented by
$v_{s}$. As stated in the main text, the Hopf bifurcation point locates at
$v_{s}\simeq 0.25725$. The the deterministic reaction equations for the
current model reads with $\bm{w}$ the transition rates

$\text{$\frac{dx_{1}}{dt}=w_{1}-w_{2}$},$

$\text{$\frac{dx_{2}}{dt}=w_{3}-w_{4}-w_{5}+w_{6}$},$

$\text{$\frac{dx_{2}}{dt}=w_{5}-w_{6}$}.$

The descriptions of the reaction channels and values of parameters are listed
in Table I. Typical trajectories for the concentrations of
$\left(\text{R},\text{P}_{C},\text{P}_{N}\right)$ in this model have been
shown in Fig.5. In Fig.6(a), we plot the the transport efficiencies as a
function of the control parameter $\alpha$ for the circadian clock model. It
can be observed that the estimator proposed by us yields a significant
improvement over the conventional TUR.

Table 1: Descriptions of the circadian clock model | Reaction | Transition rate | Biochemical function
---|---|---|---
1 | G $\to$ R + G | $w_{1}=\frac{v_{s}k_{I}^{n}}{k_{I}^{n}+x_{3}^{n}}$ | Transcription
2 | R $\to$ | $w_{2}=\frac{v_{m}x_{z}^{n}}{k_{m}+x_{1}}$ | R degradation
3 | R $\to$ R + $\text{P}_{C}$ | $w_{3}=k_{s}x_{1}$ | Translation
4 |  $\text{P}_{C}$ $\to$ | $w_{4}=\frac{v_{d}x_{2}}{k_{d}+x_{2}}$ | Degradation of $\text{P}_{C}$
5 |  $\text{P}_{C}$ $\to$ $\text{P}_{R}$ | $w_{5}=k_{1}x_{2}$ | Transport of $\text{P}_{C}$ into the nucleus
6 |  $\text{P}_{N}$ $\to$ $\text{P}_{C}$ | $w_{6}=k_{2}x_{3}$ | Transport of $\text{P}_{N}$ out of the nucleus
$k_{I}=2.0$ nM, $n=4$, $v_{m}=0.3$ nM $\text{h}^{-1}$, $k_{m}=0.2$ nM,
$k_{s}=2.0$ $\text{h}^{-1}$, $v_{d}=1.5$ nM $\text{h}^{-1}$, $k_{d}=0.1$ nM,
$k_{1}=k_{2}=0.2$ $\text{h}^{-1}$

###  C.2 The Brusselator model

Here, we introduce the Brusselator model, involving two distinct biochemical
species $X$, $Y$, whose time evolution is governed by the following
deterministic kinetic equations:

$\frac{dX}{dt}=A-(B+1)X+X^{2}Y,$

$\frac{dY}{dt}=BX-X^{2}Y.$

In the deterministic limit, the system has a stable point $X_{s}=A$,
$Y_{s}=B/A$, which loses stability when the control parameter $B$ exceeds the
Hopf bifurcation point $B_{c}=1+A^{2}$. The normal biochemical oscillation
happens for $B>B_{c}$. By choosing $A=0.3$, we calculate the parameters in
stochastic normal form theory as $C_{r}\simeq-2.9028$, $C_{i}\simeq 5.2506$,
$\varepsilon^{2}=4$, $\omega=1$, and $\alpha=(B-1-A^{2})/2$. In Fig.6(b), the
the transport efficiencies as a function of the control parameter $\alpha$ for
the Brusselator model have been depicted. As expected, the improved estimator
proposed by us is much more accurate than the TUR bound.

Figure 6: The transport efficiencies as a function of the control parameter
$\alpha$. The black squares correspond to the TUR,
$\eta(\theta)=\dot{W}_{TUR}/\dot{W}$, for the oscillatory phase $\theta$ only,
and the red circles show $\eta(r^{2},\theta)=\dot{W}_{I,2}/\dot{W}$ including
the correlations between the phase and amplitude $r$. (a) The circadian clock
model. (b) The Brusselator model. The system size $V=1.6\times 10^{5}$.

## References

  * Bialek and Setayeshgar (2005) W. Bialek and S. Setayeshgar, Proceedings of the National Academy of Sciences 102, 10040 (2005).
  * Hu et al. (2010) B. Hu, W. Chen, W.-J. Rappel, and H. Levine, Physical review letters 105, 048104 (2010).
  * Lan et al. (2012) G. Lan, P. Sartori, S. Neumann, V. Sourjik, and Y. Tu, Nature physics 8, 422 (2012).
  * Lan and Tu (2013) G. Lan and Y. Tu, Journal of The Royal Society Interface 10, 20130489 (2013).
  * Skoge et al. (2013) M. Skoge, S. Naqvi, Y. Meir, and N. S. Wingreen, Physical review letters 110, 248102 (2013).
  * Lang et al. (2014) A. H. Lang, C. K. Fisher, T. Mora, and P. Mehta, Physical review letters 113, 148103 (2014).
  * Sartori and Tu (2015) P. Sartori and Y. Tu, Physical review letters 115, 118102 (2015).
  * Cao et al. (2015) Y. Cao, H. Wang, Q. Ouyang, and Y. Tu, Nature physics 11, 772 (2015).
  * Fei et al. (2018) C. Fei, Y. Cao, Q. Ouyang, and Y. Tu, Nature communications 9, 1434 (2018).
  * Mora and Nemenman (2019) T. Mora and I. Nemenman, Physical review letters 123, 198101 (2019).
  * Zhang et al. (2020) D. Zhang, Y. Cao, Q. Ouyang, and Y. Tu, Nature physics 16, 95 (2020).
  * Barato and Seifert (2016) A. C. Barato and U. Seifert, Physical Review X 6, 041053 (2016).
  * Lee et al. (2018) S. Lee, C. Hyeon, and J. Jo, Physical Review E 98, 032119 (2018).
  * Hasegawa (2018) Y. Hasegawa, Physical Review E 98, 032405 (2018).
  * Hasegawa and Van Vu (2019a) Y. Hasegawa and T. Van Vu, Physical Review E 99, 062126 (2019a).
  * Marsland III et al. (2019) R. Marsland III, W. Cui, and J. M. Horowitz, Journal of the Royal Society Interface 16, 20190098 (2019).
  * Del Junco and Vaikuntanathan (2020a) C. Del Junco and S. Vaikuntanathan, Physical Review E 101, 012410 (2020a).
  * Guan et al. (2020) S. Guan, L. Xu, Q. Zhang, and H. Shi, Physical Review E 101, 012409 (2020).
  * Del Junco and Vaikuntanathan (2020b) C. Del Junco and S. Vaikuntanathan, The Journal of chemical physics 152, 055101 (2020b).
  * Cao et al. (2021) Z. Cao, H. Jiang, and Z. Hou, The Journal of Chemical Physics 155, 234901 (2021).
  * Ferrell Jr et al. (2011) J. E. Ferrell Jr, T. Y.-C. Tsai, and Q. Yang, Cell 144, 874 (2011).
  * Buzsáki and Draguhn (2004) G. Buzsáki and A. Draguhn, science 304, 1926 (2004).
  * Nakajima et al. (2005) M. Nakajima, K. Imai, H. Ito, T. Nishiwaki, Y. Murayama, H. Iwasaki, T. Oyama, and T. Kondo, science 308, 414 (2005).
  * Novák and Tyson (2008) B. Novák and J. J. Tyson, Nature reviews Molecular cell biology 9, 981 (2008).
  * Goldbeter (1997) A. Goldbeter, _Biochemical oscillations and cellular rhythms: the molecular bases of periodic and chaotic behaviour_ (Cambridge university press, 1997).
  * Martiel and Goldbeter (1987) J.-L. Martiel and A. Goldbeter, Biophysical journal 52, 807 (1987).
  * Gnesotto et al. (2018) F. Gnesotto, F. Mura, J. Gladrow, and C. P. Broedersz, Reports on Progress in Physics 81, 066601 (2018).
  * Seifert (2019) U. Seifert, Annual Review of Condensed Matter Physics 10, 171 (2019).
  * Barato and Seifert (2015) A. C. Barato and U. Seifert, Physical review letters 114, 158101 (2015).
  * Gingrich et al. (2016) T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Physical review letters 116, 120601 (2016).
  * Pietzonka et al. (2016) P. Pietzonka, A. C. Barato, and U. Seifert, Journal of Statistical Mechanics: Theory and Experiment 2016, 124004 (2016).
  * Pietzonka et al. (2017) P. Pietzonka, F. Ritort, and U. Seifert, Physical Review E 96, 012101 (2017).
  * Dechant and Sasa (2018a) A. Dechant and S.-i. Sasa, Physical Review E 97, 062101 (2018a).
  * Dechant and Sasa (2018b) A. Dechant and S.-i. Sasa, Journal of Statistical Mechanics: Theory and Experiment 2018, 063209 (2018b).
  * Agarwalla and Segal (2018) B. K. Agarwalla and D. Segal, Physical Review B 98, 155438 (2018).
  * Hasegawa and Van Vu (2019b) Y. Hasegawa and T. Van Vu, Physical review letters 123, 110602 (2019b).
  * Horowitz and Gingrich (2020) J. M. Horowitz and T. R. Gingrich, Nature Physics 16, 15 (2020).
  * Li et al. (2019) J. Li, J. M. Horowitz, T. R. Gingrich, and N. Fakhri, Nature communications 10, 1 (2019).
  * Manikandan et al. (2020) S. K. Manikandan, D. Gupta, and S. Krishnamurthy, Physical review letters 124, 120603 (2020).
  * Van Vu et al. (2020) T. Van Vu, Y. Hasegawa, et al., Physical Review E 101, 042138 (2020).
  * Otsubo et al. (2020) S. Otsubo, S. Ito, A. Dechant, and T. Sagawa, Physical Review E 101, 062106 (2020).
  * Skinner and Dunkel (2021a) D. J. Skinner and J. Dunkel, Proceedings of the National Academy of Sciences 118 (2021a).
  * Skinner and Dunkel (2021b) D. J. Skinner and J. Dunkel, Physical review letters 127, 198101 (2021b).
  * Dechant and Sasa (2021a) A. Dechant and S.-i. Sasa, Physical Review Research 3, L042012 (2021a).
  * Hwang and Hyeon (2018) W. Hwang and C. Hyeon, The journal of physical chemistry letters 9, 513 (2018).
  * Jack et al. (2020) M. Jack, N. López-Alamilla, and K. Challis, Physical Review E 101, 062123 (2020).
  * Cao et al. (2020) Z. Cao, H. Jiang, and Z. Hou, Physical Review Research 2, 043331 (2020).
  * Hasegawa and Arita (2014a) Y. Hasegawa and M. Arita, Journal of The Royal Society Interface 11, 20131018 (2014a).
  * Hasegawa and Arita (2014b) Y. Hasegawa and M. Arita, Physical review letters 113, 108101 (2014b).
  * Dechant and Sasa (2021b) A. Dechant and S.-i. Sasa, Physical Review X 11, 041061 (2021b).
  * Hou et al. (2006) Z. Hou, T. J. Xiao, and H. Xin, Chemphyschem: a European journal of chemical physics and physical chemistry 7, 1520 (2006).
  * Xiao et al. (2007) T. Xiao, J. Ma, Z. Hou, and H. Xin, New Journal of Physics 9, 403 (2007).
  * Ma et al. (2008) J. Ma, T. Xiao, Z. Hou, and H. Xin, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 043116 (2008).
  * Xiao et al. (2009) T. Xiao, Z. Hou, and H. Xin, The Journal of Physical Chemistry B 113, 9316 (2009).
  * Gillespie (2000) D. T. Gillespie, The Journal of Chemical Physics 113, 297 (2000).
  * Ko et al. (2010) C. H. Ko, Y. R. Yamada, D. K. Welsh, E. D. Buhr, A. C. Liu, E. E. Zhang, M. R. Ralph, S. A. Kay, D. B. Forger, and J. S. Takahashi, PLoS biology 8, e1000513 (2010).
  * Hou and Xin (2003) Z. Hou and H. Xin, The Journal of chemical physics 119, 11508 (2003).
  * Zhou and Kurths (2002) C. Zhou and J. Kurths, Physical review letters 88, 230602 (2002).
  * Hänggi (2002) P. Hänggi, ChemPhysChem 3, 285 (2002).
  * Kuramoto (2003) Y. Kuramoto, _Chemical oscillations, waves, and turbulence_ (Courier Corporation, 2003).
  * Goldobin et al. (2010) D. S. Goldobin, J.-n. Teramae, H. Nakao, and G. B. Ermentrout, Physical review letters 105, 154101 (2010).
  * Maccone et al. (2015) L. Maccone, D. Bruß, and C. Macchiavello, Physical review letters 114, 130401 (2015).
  * Pozsgay et al. (2017) V. Pozsgay, F. Hirsch, C. Branciard, and N. Brunner, Physical Review A 96, 062128 (2017).
  * Jebarathinam et al. (2020) C. Jebarathinam, D. Home, and U. Sinha, Physical Review A 101, 022112 (2020).
  * Benesty et al. (2009) J. Benesty, J. Chen, Y. Huang, and I. Cohen, in _Noise reduction in speech processing_ (Springer, 2009), pp. 1–4.
  * Benesty et al. (2008) J. Benesty, J. Chen, and Y. Huang, IEEE Transactions on Audio, Speech, and Language Processing 16, 757 (2008).
  * Dechant (2018) A. Dechant, Journal of Physics A: Mathematical and Theoretical 52, 035001 (2018).
  * Busiello and Pigolotti (2019) D. M. Busiello and S. Pigolotti, Physical Review E 100, 060102 (2019).
  * Manikandan and Krishnamurthy (2018) S. K. Manikandan and S. Krishnamurthy, Journal of Physics A: Mathematical and Theoretical 51, 11LT01 (2018).
  * Gingrich et al. (2017) T. R. Gingrich, G. M. Rotskoff, and J. M. Horowitz, Journal of Physics A: Mathematical and Theoretical 50, 184004 (2017).
  * Kim et al. (2020) D.-K. Kim, Y. Bae, S. Lee, and H. Jeong, Physical Review Letters 125, 140604 (2020).
  * Busiello and Fiore (2022) D. M. Busiello and C. Fiore, arXiv preprint arXiv:2205.00294 (2022).
  * Frishman and Ronceray (2020) A. Frishman and P. Ronceray, Physical Review X 10, 021009 (2020).
  * Munoz (2018) M. A. Munoz, Reviews of Modern Physics 90, 031001 (2018).
  * Verhey et al. (2011) K. J. Verhey, N. Kaul, and V. Soppina, Annual review of biophysics 40, 267 (2011).
  * Aitken et al. (2010) C. E. Aitken, A. Petrov, and J. D. Puglisi, Annual review of biophysics 39, 491 (2010).
  * Nuske et al. (2014) F. Nuske, B. G. Keller, G. Pérez-Hernández, A. S. Mey, and F. Noé, Journal of chemical theory and computation 10, 1739 (2014).
  * Hassard et al. (1981) B. D. Hassard, B. Hassard, N. D. Kazarinoff, Y.-H. Wan, and Y. W. Wan, _Theory and applications of Hopf bifurcation_ , vol. 41 (CUP Archive, 1981).
  * Arnold et al. (1996) L. Arnold, N. Sri Namachchivaya, and K. R. Schenk-Hoppé, International Journal of Bifurcation and Chaos 6, 1947 (1996).
  * Seifert (2005) U. Seifert, Physical review letters 95, 040602 (2005).
  * Seifert (2012) U. Seifert, Reports on progress in physics 75, 126001 (2012).
  * Sekimoto (2010) K. Sekimoto, _Stochastic energetics_ , vol. 799 (Springer, 2010).
  * Jarzynski (2011) C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2, 329 (2011).
  * Gaspard (2004) P. Gaspard, The Journal of chemical physics 120, 8898 (2004).