Thermodynamics of collective enhancement of precision
Yoshihiko Hasegawa
physics.bio-ph cond-mat.stat-mech







hasegawa@biom.t.u-tokyo.ac.jp




Department of Information and Communication Engineering, Graduate
School of Information Science and Technology, The University of Tokyo,
Tokyo 113-8656, Japan

The circadian oscillator exhibits remarkably high temporal precision,
despite its exposure to several fluctuations. The central mechanism
that protects the oscillator from fluctuations is a collective enhancement
of precision, where a population of coupled oscillators displays higher
temporal precision than that achieved without coupling. Since coupling
is essentially information exchange between oscillators, we herein
investigate the relation between the temporal precision and the information
flow between oscillators in the linearized Kuramoto model by using
stochastic thermodynamics. For general coupling, we find that the
temporal precision is bounded from below by the information flow.
We generalize the model to incorporate a time-delayed coupling and
demonstrate that the same relation also holds for the time-delayed
case. Furthermore, the temporal precision is demonstrated to be improved
in the presence of the time delay, and we  show that the increased
information flow is responsible for the time-delay-induced precision
improvement.



Thermodynamics of collective enhancement of precision
    Yoshihiko Hasegawa
    December 30, 2023
=====================================================





§ INTRODUCTION


The circadian oscillator is a biological clock that is prevalent in
biological organisms ranging from bacteria to humans <cit.>.
Since the circadian oscillator is induced by biochemical reactions
on a cellular level, it is subject to several fluctuations. Despite
such stochasticity, the circadian oscillator is known to exhibit incredibly
high temporal precision <cit.>, where the
temporal standard deviation is approximately 3 to 5 min in 24
h <cit.>. In a single oscillator level, a phase
response curve, which quantifies the dynamics of the oscillator, is
optimized so as to achieve high temporal precision <cit.>.
Still, this single oscillator level improvement does not seem to fully
explain the abovementioned temporal precision. Another significant
precision improvement arises at the population level <cit.>.
A previous study <cit.> experimentally showed that
when oscillators are coupled, each oscillator exhibits higher precision
than that realizable without coupling. This phenomenon is referred
to as collective enhancement of precision (CEP) <cit.>.
Since mammalian circadian oscillation occurs primarily in the suprachiasmatic
nucleus (SCN), which is a collection of 10^4 cells, CEP is considered
to be responsible for the high temporal precision in the circadian
oscillator <cit.>.

CEP is induced by coupling, which is essentially the exchange of
information. Information enables the violation of the second law of
thermodynamics through the Maxwell demon. High precision is achieved
at the cost of energy consumption, which is known as the thermodynamic
uncertainty relation <cit.>.
When we ignore the coupling between oscillators, CEP appears to achieve
thermodynamically impossible precision, motivating us to analyze CEP
from the viewpoint of information analogous to the Maxwell demon (Fig. <ref>).
Using stochastic thermodynamics <cit.>,
we find that the temporal variance of coupled oscillators is bounded
from below by the information flow conferred by coupling. Furthermore,
we generalize the obtained relation to incorporate time-delayed coupling
and show that the same relation also holds for the time-delayed case.
Intriguingly, the temporal precision in the time-delayed case is improved
as compared to the non-delayed case, and we can ascribe this time
delay induced improvement to the increase in the information flow.




§ MODEL


We consider N coupled identical oscillators, which are entrained
by an external periodic signal with the angular frequency Ω
and are subject to noise (Fig. <ref>(a)). Let
ϕ_i be the phase of the ith oscillator, and let ω
be the angular frequency of oscillators. Then, the phase dynamics
of each oscillator can be described by the following forced Kuramoto
model <cit.>:

    ϕ̇_i    =ω+Lsin(Ω t-ϕ_i)+f_i(t)+ξ_i(t),
    
    f_i(t)     =∑_j∈ V,j iK_ijsin(ϕ_j-ϕ_i),

where V={1,2,...,N}, and we aggregate all coupling effects, which
act on ϕ_i, into a single coupling variable f_i(t). The
sinusoidal term in Eq. (<ref>) represents
the periodic signal, and L>0 denotes its strength. A sinusoidal
function in Eq. (<ref>) is a coupling function
that represents the interaction between oscillators, where K_ij
is the coupling strength between the ith and jth oscillators.
Note that we allow both symmetric (K_ij=K_ji) and asymmetric
(K_ij K_ji) couplings. The phase reduction, with which the
phase equation of Eq. (<ref>) is derived,
assumes Langevin equations to be interpreted in the Stratonovich sense.
Moreover, ξ_i(t) is zero-mean white Gaussian noise with the
correlation ⟨ξ_i(t)ξ_i'(t')⟩ =2Dδ_ii'δ(t-t')
(D is the noise strength). The derivation of Eq. (<ref>)
can be found in Ref. <cit.>.
We hereinafter assume that all of the oscillators are synchronized
to the periodic signal. Although CEP is often studied for L=0,
we assume L>0, which enables an approximation of the synchronized
behavior as a nonequilibrium steady state (NESS). 







Since there are efferent projections from the SCN to other neurons,
the output of the SCN has been suggested to be generated by the average
of a subset of oscillators, and not by the average of all oscillators
<cit.>. Therefore, we observe the average phase
Φ(t)=(1/N_O)∑_i∈ V_Oϕ_i(t), where V_O is
a set of observable oscillators, and N_O is its element number
(Fig. <ref>(b)) <cit.>. We also
define V_L and N_L to be a set of latent variables (the
output of which is not observed) and its element number, respectively.
Every oscillator belongs to the whole set V=V_O∪ V_L, and
V_O and V_L are disjoint, i.e., V_O∩ V_L=∅.
See Fig. <ref>(b). Defining a relative phase x_i
observed from a rotating frame as x_i=ϕ_i-Ω t, the
synchronizing assumption is that x_i is well concentrated around
0. This condition can usually be satisfied when |ω-Ω|/L
is sufficiently small and D is sufficiently weak. Let μ_A=⟨ A⟩,
𝒱_A=⟨ A^2⟩ -⟨ A⟩ ^2,
and 𝒞_AB=⟨ AB⟩ -⟨ A⟩⟨ B⟩,
where A and B are arbitrary random variables, and let X(t)=Φ(t)-Ω t=(1/N_O)∑_i∈ V_Ox_i(t),
which is the average phase seen from a rotating frame. The temporal
variance of the average phase Φ(t) is given by the variance
𝒱_X. Since each phase x_i is well concentrated
around 0 (we assumed that all of the oscillators are synchronized
with the periodic signal), we apply a Taylor series expansion to Eqs. (<ref>)
and (<ref>) to obtain 

    Ẋ    =-L(X-c)+F(t)+Ξ_O(t),
    ẋ_i    =-L(x_i-c)+∑_j∈ V,j iK_ij(x_j-x_i)+ξ_i(t)(i∈ V),

where F(t) is defined by 

    F(t)=1/N_O∑_i∈ V_Of_i(t)=1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j-x_i),

and c=(ω-Ω)/L and Ξ_O(t)=(1/N_O)∑_i∈ V_Oξ_i(t)
[⟨Ξ_O(t)Ξ_O(t')⟩ =(2D/N_O)δ(t-t'),
⟨Ξ_O(t)ξ_i(t')⟩ =(2D/N_O)δ(t-t')
for i∈ V_O, and ⟨Ξ_O(t)ξ_i(t')⟩ =0
for i∈ V_L]. Here, F(t) includes all of the coupling effects
that act on x_i for i∈ V_O. Therefore, the information
flow between X and F is of interest. We calculate the information
flow between X and F based on Refs. <cit.>,
which is defined as

    ℐ̇_X(X;F)     =-∫ J_X(X,F)∂/∂ Xln P(F|X)dXdF
        =-1/dt⟨ln(P(F(t)|X(t+dt))/P(F(t)|X(t)))⟩ ,

where the sign is reversed compared to the original definition. This
quantity is equivalent to the learning rate <cit.>.
In order to calculate the information flow, we introduce a hypothetical
time delay h>0 in Eq. (<ref>), which yields F(t)=(1/N_O)∑_i∈ V_O∑_j∈ V,j iK_ij(x_j(t-h)-x_i(t-h)).
The introduction of time delay h is a mathematical trick, and,
when calculating several quantities, we set h→0 afterwards.
Using a Taylor expansion for a sufficiently small h, we obtain

    Ḟ=h^-1(-F+1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j-x_i)).

Then, we can treat Eqs. (<ref>), (<ref>),
and (<ref>) as coupled Langevin equations. Let P(X,x,F)
be a time-dependent probability density function of X, x=[x_1,...,x_N],
and F (for notational convenience, we dropped time t in the
argument in P(X,x,F)). Due to the linearity of Eqs. (<ref>),
(<ref>), and (<ref>), the steady-state
distribution is a Gaussian distribution with mean [μ_X,μ_x_1,..,μ_x_N,μ_F]=[c,c,...,c,0].
The corresponding Fokker–Planck equation (FPE) with respect to P(X,x,F)
can be represented as 

    ∂_tP(X,x,F)=𝕃(X,x,F)P(X,x,F),

where 𝕃(X,x,F) is an FPE operator, which
is defined by  

    𝕃(X,x,F)     =-∂/∂ X(-L(X-c)+F)+D/N_O∂^2/∂ X^2-∂/∂ F1/h{ -F+1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j-x_i)}
        -∑_i∈ V∂/∂ x_i{ -L(x_i-c)+∑_j∈ V,j iK_ij(x_j-x_i)} +D∑_i∈ V∂^2/∂ x_i^2+2D/N_O∑_i∈ V_O∂^2/∂ x_i∂ X.

We integrate out x in Eqs. (<ref>)
and (<ref>) to obtain 

    ∂_tP(X,F)     =-∂_XJ_X(X,F)-∂_FJ_F(X,F),

where J_X(X,F) and J_F(X,F) are probability currents, which
are defined as 

    J_X(X,F)     ={ -L(X-c)+F} P-D/N_O∂ P/∂ X,
    
    J_F(X,F)     =h^-1{ -F+ℋ(X,F)} P.

in which P=P(X,F) and ℋ(X,F)=∫(1/N_O)∑_i∈ V_O∑_j∈ V,j iK_ij(x_j-x_i)P(x|X,F)dx.
Although the FPE operator 𝕃(X,x,F) of Eq. (<ref>)
contains cross terms, such as ∂_X∂_x_i, due
to non-vanishing correlation ⟨Ξ_O(t)ξ_i(t)⟩0
for i∈ V_O, they disappear when integrating out x.
From Eq. (<ref>), we can obtain the generalized second
law following the conventional procedure in stochastic thermodynamics
<cit.>.
We first consider the Shannon entropy S(X,F)=-∫ P(X,F)ln P(X,F)dXdF.
The time derivative is Ṡ(X,F)=Ṡ_X(X,F)+Ṡ_F(X,F),
where Ṡ_X(X,F)=-∫ J_X(X,F)∂_Xln P(X,F)dXdF,
and Ṡ_F(X,F)=-∫ J_F(X,F)∂_Fln P(X,F)dXdF.
We can divide J_X(X,F) of Eq. (<ref>) into two contributions,
J_X^R(X,F) and J_X^I(X,F), which are defined as J_X^I(X,F)=-L(X-c)P-(D/N_O)∂_XP
and J_X^R=FP, respectively. In feedback cooling of Brownian
particles, J_X^I and J_X^R correspond to the irreversible
and reversible portions of the probability current <cit.>.
Then, Ṡ_X(X,F) is calculated as 

    Ṡ_X(X,F)=∫(N_O(J_X^I)^2/DP+N_OJ_X^IL(X-c)/D)dXdF,

where J_X^I=J_X^I(X,F), and P=P(X,F). In Eq. (<ref>),
the first term is non-negative, and the second term is calculated
as (1/D)∫ N_OJ_X^I(X,F)L(X-c)dXdF=-(N_OL^2/D)⟨ (X-c)^2⟩ +L.
According to Ref. <cit.>, we have ℐ̇_X(X;F)=-Ṡ(X)+Ṡ_X(X,F),
where S(X)=-∫ P(X)ln P(X)dX. Under the synchronization assumption,
X is in a steady state, and thus Ṡ(X)=0. Combining these
representations, the variance 𝒱_X=⟨ (X-c)^2⟩
satisfies 

    𝒱_X≥𝒱_X^LB=D/L^2N_O(L-ℐ̇_X(X;F)),

where 𝒱_X^LB is the lower bound of the variance
𝒱_X. Equation (<ref>) is our main result.
When there is no coupling, 𝒱_X=D/(LN_O). Therefore,
without information flow, the variance does not improve beyond D/(LN_O),
and Eq. (<ref>) shows that CEP is induced by the information
flow due to coupling between oscillators. Equation (<ref>)
indicates that CEP has the same mathematical property as the feedback
cooling by the Maxwell demon <cit.>,
where the feedback reduces the variance. The angular frequency Ω
affects the inequality of Eq. (<ref>) through c=(ω-Ω)/L,
showing that the inequality does not have Ω dependence when
ω=Ω.  Note that we cannot reach the desired inequality
of Eq. (<ref>) by using the positivity of ∫ J_X(X,F)^2/P(X,F)dXdF
in the conventional total entropy production rate <cit.>.
The transfer entropy rate <cit.> is a similar
quantity to the information flow. It has been shown that the transfer
entropy rate from X to F is greater than or equal to the information
flow <cit.>.
Therefore, the same inequality as Eq. (<ref>) holds for
the transfer entropy rate but the bound becomes weaker than the information
flow case. 




§ EXAMPLES





 §.§ Globally coupled model


We calculate 𝒱_X and 𝒱_X^LB
analytically for a global uniform coupling case K_ij=K. Defining
Y=(1/N_L)∑_i∈ V_Lx_i, we have 

    Ẏ    =-L(Y-c)-F/R+Ξ_L(t),
    Ḟ    =h^-1{ -F+KN_L(Y-X)} ,

where R=N_L/N_O and Ξ_L(t)=(1/N_L)∑_i∈ V_Lξ_i(t)
[⟨Ξ_L(t)Ξ_L(t')⟩ =(2D/N_L)δ(t-t')
and ⟨Ξ_L(t)Ξ_O(t')⟩ =0]. Then,
Eqs. (<ref>), (<ref>), and (<ref>)
constitute coupled Langevin equations of the global uniform coupling
model. Taking a limit h→0 for which the hypothetical
time delay in F vanishes, we can calculate 𝒱_X,
𝒱_F, and 𝒞_XF analytically due to the
linearity of the coupled equations. Then, the information flow ℐ̇_X(X;F)
can be calculated as follows (Appendix <ref>):

    ℐ̇_X(X;F)=LR.

We plot 𝒱_X^LB and 𝒱_X of
the global uniform coupling case for D=0.001 (Fig. <ref>(a))
and D=0.1 (Fig. <ref>(b)). The other parameters
are shown in the caption of Fig. <ref>. In Fig. <ref>,
the solid line denotes 𝒱_X^LB, and the long-dashed
and dashed lines show 𝒱_X for K=1 and K=∞,
respectively. We also carried out a Monte Carlo simulation (for details,
please see Appendix <ref>), where 𝒱_X
are denoted by circles (K=1) and triangles (K=10, which is intended
to simulate K=∞), and 𝒱_X^LB are
denoted by squares (K=1) and crosses (K=10). Note that the Monte
Carlo simulation is carried out for equations that do not use the
Taylor approximation. Apparently, 𝒱_X≥𝒱_X^LB
for R≥0, which verifies the inequality relation. For R>1,
𝒱_X^LB becomes negative, and the lower bound
does not have a practical meaning. The difference is 𝒱_X-𝒱_X^LB=L^-2∫ J_X^I(X,F)^2/P(X,F)dXdF
and the probability current |J_X^I| becomes larger for larger
R, which is responsible for the gap between 𝒱_X^LB
and 𝒱_X. The inequality saturates when F→0,
i.e., there is no coupling. Figure <ref> shows
that the Monte Carlo results agree with the analytical calculations,
even for a relatively large noise intensity of D=0.1 (Fig. <ref>(b)),
and this demonstrates the validity of the synchronization assumption
and the Taylor approximation. When R is sufficiently small, the
variance 𝒱_X is close to 𝒱_X^LB,
which shows that oscillators maximally exploit the information flow
to improve the variance. As shown in the next model, higher temporal
precision can be achieved if the information flow ℐ̇_X(X;F)
can be increased.










 §.§ Time-delayed model


Next, we generalize the model to include time delay in the coupling,
which exists in reality and has been reported to be able to cause
different dynamics in the Kuramoto model <cit.>.
In the context of feedback cooling, a time-delayed case was investigated
in Ref. <cit.>. Let τ≥0
be the time delay of the coupling. Following Ref. <cit.>,
we incorporate the time delay into Eq. (<ref>)
using the following f_i(t) instead of Eq. (<ref>):

    f_i(t)=∑_j∈ V,j iK_ijsin(ϕ_j^τ-ϕ_i),

where superscript τ hereinafter denotes a time-delayed variable,
i.e., ϕ_j^τ(t)=ϕ_j(t-τ). Similar to the non-delayed
case, we introduce a relative phase x_i=ϕ_i-Ω t. For
the non-delayed case, μ_X=0 when ω=Ω. However,
for the time-delayed case, μ_X deviates from 0, even when
ω=Ω. Therefore, in order for the Taylor approximation
to yield reliable results (i.e., x_i should be well concentrated
around 0), Ωτ should be sufficiently small. In such
cases, we can use the Taylor approximation to obtain F=(1/N_O)∑_i∈ V_Of_i=(1/N_O)∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ),
and the introduction of the hypothetical delay h leads us to Ḟ=(1/h){ -F+(1/N_O)∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)}.
Following Ref. <cit.>, an effective FPE, which
does not include time-delayed variables, can be derived for the time-delayed
case. Let x^τ=[x_1^τ,...,x_N^τ]
and P(X,x,F;x^τ) be the probability
density of X, x, and F at time t, and x^τ
at time t-τ. Calculating the marginal distribution P(X,F)=∫ P(X,x,F;x^τ)dxdx^τ,
the FPE for P(X,F) is similar to that given by Eqs. (<ref>)–(<ref>).
Therefore, the inequality as in Eq. (<ref>) holds for
the time-delayed case (Appendix <ref>).




 §.§ Globally coupled time-delayed model


Next, we calculate 𝒱_X and 𝒱_X^LB
for a globally coupled time-delayed case. For this case, we obtained
the following equations: 

    Ẋ    =-LX+F+Ξ_O(t),
    Ẏ    =K{(N_L-1)Y^τ-(N-1)Y+N_OX^τ}
        -LY+Ξ_L(t),
    
    F     =K{ (N_O-1)X^τ-(N-1)X+N_LY^τ} ,

where we have applied a parallel translation so that μ_X=μ_Y=μ_F=0
holds for the steady state. Since the system is the steady state under
the synchronization assumption, we can use the Fourier transform to
calculate the covariance matrix <cit.>.
Let ℱ[⋯] be the Fourier transform operator, and
let ν be the Fourier variable. We define 𝒮_F(ν)=ℱ[⟨ F(t)F(0)⟩],
𝒮_X(ν)=ℱ[⟨ X(t)X(0)⟩],
and 𝒮_XF(ν)=ℱ[⟨ X(t)F(0)⟩].
By virtue of the Wiener–Khinchin theorem, 𝒱_X, 𝒱_F,
and 𝒞_XF can be obtained by 𝒱_X=(2π)^-1∫_-∞^∞𝒮_X(ν)dν,
𝒱_F=(2π)^-1∫_-∞^∞𝒮_F(ν)dν,
and 𝒞_XF=(2π)^-1∫_-∞^∞𝒮_XF(ν)dν
(Appendix <ref>). Since the integration
is complicated, we evaluate these integrals numerically.

Figure <ref>(a) shows 𝒱_X and 𝒱_X^LB
as a function of R for time delay τ=0.02 and K=1. The
other parameters are the same as in Fig. <ref>(a).
In Fig. <ref>(a), the dotted and dot-dashed lines
denote 𝒱_X and 𝒱_X^LB for
τ=0.02, respectively. The solid and long-dashed lines represent
𝒱_X^LB and 𝒱_X for the non-delayed
case and are identical to those in Fig. <ref>(a)
(shown for comparison). We also carried out Monte Carlo simulations
for τ=0.02, and these data are shown by circles and triangles,
which correspond to 𝒱_X and 𝒱_X^LB,
respectively. We again confirm that the inequality of Eq. (<ref>)
also holds for the delayed case. Figure <ref>(a)
shows that 𝒱_X for τ=0.02 is lower than that
of the non-delayed case. Furthermore, 𝒱_X^LB
is also lower for the delayed case. In order to clarify the effect
of the time delay, we plot 𝒱_X and 𝒱_X^LB
as a function of τ for N=N_O=10 and K=1 in Fig. <ref>(b).
The solid and dashed lines denote 𝒱_X and 𝒱_X^LB,
respectively, obtained analytically, and the circles and triangles
denote 𝒱_X and 𝒱_X^LB, respectively,
obtained by Monte Carlo simulations. For τ<0.03, the analytical
and Monte Carlo results agree, whereas they do not agree for τ>0.03.
This is because ϕ_j^τ-ϕ_i in the coupling deviates
from 0 when τ increases, which degrades the reliability of
the Taylor approximation. Intriguingly, Fig. <ref>(b)
shows that when τ increases to τ∼0.04, both 𝒱_X
and 𝒱_X^LB decrease. Since the lower bound
𝒱_X^LB is smaller in the presence of a time
delay, the precision improvement can be ascribed to the increase in
the information flow ℐ̇_X(X;F). Although the time
delay in the Kuramoto model is often studied from the viewpoint of
coherence and incoherence <cit.>, we show
that the time delay is beneficial to temporal precision of oscillators.




§ CONCLUSION


We obtained the inequality relating the temporal variance in the Kuramoto
model with the information flow conferred by coupling between oscillators.
CEP is a universal phenomenon and is not limited to coupled oscillators.
For instance, the precision of cellular concentration inference is
substantially improved when the concentration is measured by a population
of cells <cit.>. A similar calculation can be
applied in such cases, which will be left for future studies. 





§ GLOBALLY COUPLED MODEL


We calculate 𝒱_X^LB and 𝒱_X
for a uniformly globally coupled model, where K_ij=K. The Langevin
equations for X and Y are given by Eqs. (<ref>)
and (<ref>), respectively. The coupling variable F
is 

    F=KN_L(Y-X).

Introducing hypothetical delay h in Eq. (<ref>),
we have 

    F(t)=KN_L(Y(t-h)-X(t-h)).

Assuming that h is sufficiently small, the coupling variable obeys
Eq. (<ref>). Equations (<ref>), (<ref>),
and (<ref>) constitute coupled Langevin equations.
The corresponding FPE is 

    ∂/∂ tP(X,Y,F)     =-∂/∂ X{ -L(X-c)+F} P(X,Y,F)+D/N_O∂^2/∂ X^2P(X,Y,F)
        -∂/∂ Y{ -L(Y-c)-F/R} P(X,Y,F)+D/N_L∂^2/∂ Y^2P(X,Y,F)
        -∂/∂ F1/h{ -F+KN_L(Y-X)} P(X,Y,F).

The steady-state distribution of Eq. (<ref>)
is a Gaussian distribution: 

    P(X,Y,F)=𝒩(X,Y,F|μ_XYF,Σ_XYF),

where 𝒩(·) denotes a multivariate Gaussian distribution,
and μ_XYF and Σ_XYF are
its mean vector and covariance matrix, respectively. We can obtain
the mean and covariance in the steady state. For h→0,
the hypothetical delay vanishes, and F reduces to Eq. (<ref>).
The mean vector is μ_XYF=[c,c,0], and the covariance
matrix Σ_XYF is 

    Σ_XYF=[[   (KN_O+L)D/N_OL(K(R+1)N_O+L)             DK/L(K(R+1)N_O+L)              -KRD/K(R+1)N_O+L;             DK/L(K(R+1)N_O+L) (KN_OR+L)D/RN_OL(K(R+1)N_O+L)                DK/K(R+1)N_O+L;              -KRD/K(R+1)N_O+L                DK/K(R+1)N_O+L     N_O(R+1)RK^2D/K(R+1)N_O+L ]].

With Eq. (<ref>), the information flow
is obtained as follows: 

    ℐ̇_X(X;F)=-𝒞_XF{𝒞_XF(D+N_O𝒞_XF)-N_O𝒱_F𝒱_X}/N_O𝒱_X(𝒞_XF^2-𝒱_F𝒱_X).

From Eqs. (<ref>) and (<ref>), we obtain
Eq. (<ref>). 




§ TIME-DELAYED MODEL


In the main text, we consider the time-delayed case. Assuming that
Ωτ is sufficiently small, we can apply the Taylor approximation
to obtain

    dX/dt    =ω-Ω-LX+F+Ξ_O(t),
    dx_i/dt    =ω-Ω-Lx_i+∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)
        +ξ_i(t)(i∈ V),
    dF/dt    =1/h{ -F+1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)} .

When deriving Eq. (<ref>), we again introduced
the hypothetical time-delay h as in Eq. (<ref>).
Suppose that the mean of x in the steady state is
[μ_x_1,μ_x_2,...,μ_x_N]=[c_x_1,c_x_2,...,c_x_N],
where [c_x_1,c_x_2,...,c_x_N] satisfies the following
relation:

    ω-Ω-Lc_x_i+∑_j∈ V,j iK_ij(c_x_j-c_x_i-Ωτ)=0(i∈ V).

 Equations (<ref>)–(<ref>)
can be written as

    dX/dt    =-L(X-c_X)+F-c_F+Ξ_O(t),
    dx_i/dt    =-L(x_i-c_x_i)+∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)-c_f_i
        +ξ_i(t)(i∈ V),
    dF/dt    =1/h{ -F+1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)} ,

where 

    c_f_i    =∑_j∈ V,j iK_ij(c_x_j-c_x_i-Ωτ),
    
    c_X    =1/N_O∑_i∈ V_Oc_x_i,
    
    c_F    =1/N_O∑_i∈ V_Oc_f_i.

Following Ref. <cit.>, we can obtain the time
evolution of P(X,x,F)=∫ P(X,x,F;x^τ)dx^τ
where P(X,x,F;x^τ) is the probability
density of X, x, and F at time t and of x^τ
at time t-τ. Let 𝒜(X,x,F) be an arbitrary
function of X, x, and F. The derivation in Ref. <cit.>
is constructed on the Ito interpretation. If Langevin equations of
interest obey the Stratonovich interpretation, they have to be converted
to equivalent Ito equations <cit.>. Since Eqs. (<ref>)–(<ref>)
include only additive noise terms, the calculations afterwards do
not depend on the stochastic integral. Taking the expectation for
d𝒜(X,x,F)/dt, we apply the Ito calculus
to obtain 

    ⟨d𝒜/dt⟩    =⟨∂𝒜/∂ X[-L(X-c_X)+F-c_F]+1/h∂𝒜/∂ F[-F+1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)].
        +∑_i∈ V∂𝒜/∂ x_i[-L(x_i-c_x_i)+∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)-c_f_i]
       .+D/N_O∂^2𝒜/∂ X^2+∑_i∈ VD∂^2𝒜/∂ x_i^2+∑_i∈ V_O2D/N_O∂^2𝒜/∂ X∂ x_i⟩ _P(X,x,F;x^τ)
        =⟨∂𝒜/∂ X[-L(X-c_X)+F-c_F]+1/h∂𝒜/∂ F[-F+⟨1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)⟩ _P(x^τ|X,x,F)].
        +∑_i∈ V∂𝒜/∂ x_i[-L(x_i-c_x_i)+⟨∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)⟩ _P(x^τ|X,x,F)-c_f_i]
       .+D/N_O∂^2𝒜/∂ X^2+∑_i∈ VD∂^2𝒜/∂ x_i^2+∑_i∈ V_O2D/N_O∂^2𝒜/∂ X∂ x_i⟩ _P(X,x,F).

Here, we explicitly write the probability density with
which the expectation is calculated. Letting P(B) be a probability
density function of arbitrary random variable(s) B (∫ P(B)dB=1),
we define ⟨⋯⟩ _P(B)=∫⋯ P(B)dB
(e.g., ⟨⋯⟩ _P(x^τ|X,F,x)=∫⋯ P(x^τ|X,x,F)dx^τ).
Since 𝒜(X,x,F) is an arbitrary function,
using the integration by parts, we obtain 

    ∂/∂ tP(X,x,F)=𝕃_τ(X,x,F)P(X,x,F),

where 𝕃_τ(X,x,F) is the following FPE
operator: 

    𝕃_τ(X,x,F)     =-∂/∂ X[-L(X-c_X)+F-c_F]
        -∑_i∈ V∂/∂ x_i[-L(x_i-c_x_i)+⟨∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)⟩ _P(x^τ|X,x,F)-c_f_i]
        -∂/∂ F1/h[-F+⟨1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)⟩ _P(x^τ|X,x,F)]
        +D/N_O∂^2/∂ X^2+∑_i∈ VD∂^2/∂ x_i^2+∑_i∈ V_O2D/N_O∂^2/∂ X∂ x_i.

For details of the derivations, please see Ref. <cit.>.
Integrating Eqs. (<ref>) and (<ref>)
with respect to x, we obtain

    ∂ P(X,F)/∂ t    =-∂/∂ X[-L(X-c_X)+F-c_F]P(X,F)
        +D/N_O∂^2/∂ X^2P(X,F)
        -∂/∂ F1/h{ -F+ℋ_τ(X,F)} P(X,F),

where 

    ℋ_τ(X,F)=⟨1/N_O∑_i∈ V_O∑_j∈ V,j iK_ij(x_j^τ-x_i-Ωτ)⟩ _P(x,x^τ|X,F).

Since Eq. (<ref>) is identical to Eqs. (<ref>)–(<ref>)
when replacing ℋ(X,F) with ℋ_τ(X,F)
and F→ F-c_F in a drift term, we obtain

    ⟨ (X-c_X)^2⟩≥D/N_OL^2(L-ℐ̇_X(X;F)),

which is the same inequality as the non-delayed model (Eq. (<ref>)).




§ GLOBALLY COUPLED TIME-DELAY MODEL


We consider a globally coupled case K_ij=K to calculate 𝒱_X^LB
and 𝒱_X. The calculations can be classified into two
cases: N_L≥1 and N_L=0. 

For N_L≥1, we obtain Eqs. (<ref>)–(<ref>)
for X, Y, and F. We need to calculate the covariance matrix
of X and F. These quantities can be obtained via the Fourier
transform <cit.>. The detailed procedures (known
as the Rice method) are shown in Ref. <cit.>.
For an arbitrary random variable A(t), we can define the Fourier
transform and its inverse as follows: 

    A(ν)     =ℱ[A(t)]=∫_-∞^∞e^-iν tA(t)dt,
    
    A(t)     =ℱ^-1[A(ν)]=1/2π∫_-∞^∞e^iν tA(ν)dν,

where i in a roman typeface denotes the imaginary unit.
Then, applying the Fourier transform, we obtain 

    [[ X(ν); Y(ν); F(ν) ]]=ℳ[[ Ξ_O(ν); Ξ_L(ν);      0 ]],

where ℳ is a regular matrix defined by

    ℳ=[[                               L+iν                                  0                                 -1;                        -e^-iντKN_O L+K{ -1-e^-iντ(N_L-1)+N_L+N_O} +iν                                  0;       K{ -1+N_L-e^-iντ(N_O-1)+N_O}                        -e^-iντKN_L                                  1 ]]^-1.

We can consider the Fourier transform of the correlation
function: 𝒮_X(ν)=ℱ[⟨ X(t)X(0)⟩],
𝒮_F(ν)=ℱ[⟨ F(t)F(0)⟩],
𝒮_XF(ν)=ℱ[⟨ X(t)F(0)⟩],
𝒮_O(ν)=ℱ[⟨Ξ_O(t)Ξ_O(0)⟩]=2D/N_O,
and 𝒮_L(ν)=ℱ[⟨Ξ_L(t)Ξ_L(0)⟩]=2D/N_L.
Let A and B be arbitrary random variables. According to the
Wiener–Khinchin theorem, the following relation holds <cit.>:

    ⟨A(ν)B^*(ν^')⟩    =2πδ(ν-ν^')𝒮_AB(ν),
    𝒮_AB(ν)     =ℱ[⟨ A(t)B^*(0)⟩],

where superscript * denotes the complex conjugate. From Eqs. (<ref>)–(<ref>),
we obtain 

    𝒮_X(ν)     =|ℳ_11|^2𝒮_O(ν)+|ℳ_12|^2𝒮_L(ν),
    𝒮_F(ν)     =|ℳ_31|^2𝒮_O(ν)+|ℳ_32|^2𝒮_L(ν),
    𝒮_XF(ν)     =ℳ_11ℳ_31^*𝒮_O(ν)+ℳ_12ℳ_32^*𝒮_L(ν),

where ℳ_ij is the i,j-th element of ℳ.
When we assume a Gaussian distribution for P(X,F), from Eq. (<ref>),
it is sufficient to calculate 𝒱_X, 𝒱_F,
and 𝒞_XF for information flow ℐ̇_X.
From Eq. (<ref>), the inverse Fourier transform yields variance
𝒱_X and 𝒱_F and covariance 𝒞_XF:

    𝒱_X    =⟨ X(t)^2⟩ =1/2π∫_-∞^∞𝒮_X(ν)dν,
    𝒱_F    =⟨ F(t)^2⟩ =1/2π∫_-∞^∞𝒮_F(ν)dν,
    𝒞_XF    =⟨ X(t)F(t)⟩ =1/2π∫_-∞^∞𝒮_XF(ν)dν.

Since it is unlikely that we can obtain the inverse Fourier transforms
of Eqs. (<ref>)–(<ref>) analytically, we calculated
the transforms numerically.

We next consider the case in which N_L=0. In a similar way, we
have the following coupled equations: 

    dX/dt    =-LX+F+Ξ_O(t),
    
    F     =K(N-1)(X^τ-X).

The Fourier transform yields 

    [[ X(ν); F(ν) ]]=ℒ[[ Ξ_O(ν);      0 ]],

where ℒ is a regular matrix defined as 

    ℒ=[[             L+iν               -1; K(N-1)(1-e^-iντ)                1 ]]^-1.

Then, the Fourier transforms of the correlation functions are 

    𝒮_X(ν)     =|ℒ_11|^2𝒮_O(ν),
    𝒮_F(ν)     =|ℒ_21|^2𝒮_O(ν),
    𝒮_XF(ν)     =ℒ_11ℒ_21^*𝒮_O(ν).

Again, the variance and covariance can be obtained by the inverse
Fourier transform of Eqs. (<ref>)–(<ref>).




§ MONTE CARLO SIMULATION


In order to solve Langevin equations numerically, we use the Euler
method shown in Ref. <cit.> with time step ϵ=0.001.
When calculating the information flow ℐ̇_X(X;F),
we use the following approximation from Eq. (<ref>):

    ℐ̇_X(X;F)≃-1/ϵ⟨ln(P(F(t)|X(t+ϵ))/P(F(t)|X(t)))⟩ .

Here, we use F(t) based on Eq. (<ref>) [non-delayed
case] or Eq. (<ref>) [time-delayed case].
In Eq. (<ref>), we need to calculate the probability
density P(X,F). Assuming a multivariate Gaussian distribution for
P(X,F), we empirically calculate the mean vector and covariance
matrix for P(X,F). Then, we can numerically calculate the information
flow ℐ̇_X(X;F) as the average of 2×10^6
samples. Regarding the number of oscillators, we use N_O=10 and
N_L=0–10 (0≤ R≤1). Note that when we employ an excessively
large N, the simulation suffers from an artifact due to the time
discretization, because discretized coupling terms f_i(t)ϵ
are on the order of O(Nϵ).



40
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Refinetti(2005)]Refinetti:2005:CircBook
authorR. Refinetti,
  titleCircadian physiology (publisherTaylor &
  Francis, year2005), edition2nd ed.

[Ukai and Ueda(2010)]Ukai:2010:MamCircadian
authorH. Ukai and
  authorH. R. Ueda,
  journalAnnu. Rev. Physiol. volume72,
  pages579 (year2010).

[Moortgat et al.(2000)Moortgat, Bullock,
  and Sejnowski]Moortgat:2000:Precision
authorK. T. Moortgat,
  authorT. H. Bullock,
  and authorT. J.
  Sejnowski, journalJ. Neurophysiol.
  volume83, pages971 (year2000).

[Enright(1980)]Enright:1980:CircPrec
authorJ. T. Enright,
  journalScience volume209,
  pages1542 (year1980).

[Hasegawa and
  Arita(2014a)]Hasegawa:2013:OptimalPRC
authorY. Hasegawa and
  authorM. Arita, journalJ.
  R. Soc. Interface volume11, pages20131018
  (year2014a).

[Hasegawa and
  Arita(2014b)]Hasegawa:2014:PRL
authorY. Hasegawa and
  authorM. Arita,
  journalPhys. Rev. Lett. volume113,
  pages108101 (year2014b).

[Winfree(2001)]Winfree:2001:GeoBiolTime
authorA. T. Winfree,
  titleThe Geometry of Biological Time
  (publisherSpringer, year2001), edition2nd
  ed.

[Clay and DeHaan(1979)]Clay:1979:CEP
authorJ. Clay and
  authorR. L. DeHaan,
  journalBiophys. J. volume28,
  pages377 (year1979).

[Needleman et al.(2001)Needleman,
  Tiesinga, and Sejnowski]Needleman:2001:CEP
authorD. J. Needleman,
  authorP. H. E. Tiesinga,
  and authorT. J.
  Sejnowski, journalPhysica D
  volume155, pages324 (year2001).

[Kori et al.(2012)Kori, Kawamura, and
  Masuda]Kori:2012:OscReg
authorH. Kori,
  authorY. Kawamura, and
  authorN. Masuda,
  journalJ. Theor. Biol. volume297,
  pages61 (year2012).

[Cao et al.(2015)Cao, Wang, Ouyang, and
  Tu]Cao:2015:ClockEnergy
authorY. Cao,
  authorH. Wang,
  authorQ. Ouyang, and
  authorY. Tu, journalNat.
  Phys. volume11, pages772
  (year2015).

[Barato and Seifert(2015)]Barato:2015:UncRel
authorA. C. Barato and
  authorU. Seifert,
  journalPhys. Rev. Lett. volume114,
  pages158101 (year2015).

[Barato and Seifert(2016)]Barato:2016:BrClo
authorA. C. Barato and
  authorU. Seifert,
  journalPhys. Rev. X volume6,
  pages041053 (year2016).

[Gingrich et al.(2016)Gingrich, Horowitz,
  Perunov, and England]Gingrich:2016:TUP
authorT. R. Gingrich,
  authorJ. M. Horowitz,
  authorN. Perunov, and
  authorJ. L. England,
  journalPhys. Rev. Lett. volume116,
  pages120601 (year2016).

[Horowitz and Gingrich(2017)]Horowitz:2017:TUR
authorJ. M. Horowitz and
  authorT. R. Gingrich,
  journalPhys. Rev. E volume96,
  pages020103 (year2017).

[Ritort(2008)]Ritort:NEArticle
authorF. Ritort, in
  booktitleAdvances in Chemical Physics, edited by
  editorS. A. Rice
  (publisherWiley publications, year2008), vol.
  volume137, pp. pages31–123.

[Seifert(2012)]Seifert:2012:FTReview
authorU. Seifert,
  journalRep. Prog. Phys. volume75,
  pages126001 (year2012).

[Van den Broeck and
  Esposito(2015)]VandenBroeck:2015:Review
authorC. Van den Broeck
  and authorM. Esposito,
  journalPhysica A volume418,
  pages6 (year2015).

[Sakaguchi(1988)]Sakaguchi:1988:KuramotoExt
authorH. Sakaguchi,
  journalProg. Theor. Phys. volume79,
  pages39 (year1988).

[Acebrón et al.(2005)Acebrón,
  Bonilla, Pérez Vicente, Ritort, and Spigler]Acebron:2005:KuramotoReview
authorJ. A. Acebrón,
  authorL. L. Bonilla,
  authorC. J. Pérez Vicente,
  authorF. Ritort, and
  authorR. Spigler,
  journalRev. Mod. Phys. volume77,
  pages137 (year2005).

[Pikovsky et al.(2001)Pikovsky,
  Rosenblem, and Kurths]Pikovsky:2001:SyncBook
authorA. Pikovsky,
  authorM. Rosenblem,
  and authorJ. Kurths,
  titleSynchronization: A universal concept in nonlinear
  sciences (publisherCambridge University Press,
  year2001).

[Kuramoto(2003)]Kuramoto:2003:OscBook
authorY. Kuramoto,
  titleChemical Oscillations, Waves, and Turbulence
  (publisherDover publications, addressMineola, New
  York, year2003).

[Horowitz and Esposito(2014)]Horowitz:2014:InfoFlow
authorJ. M. Horowitz and
  authorM. Esposito,
  journalPhys. Rev. X volume4,
  pages031015 (year2014).

[Horowitz(2015)]Horowitz:2015:Multipartite
authorJ. M. Horowitz,
  journalJ. Stat. Mech. volume2015,
  pagesP03006 (year2015).

[Barato et al.(2014)Barato, Hartich, and
  Seifert]Barato:2014:CellInfo
authorA. C. Barato,
  authorD. Hartich, and
  authorU. Seifert,
  journalNew J. Phys. volume16,
  pages103024 (year2014).

[Hartich et al.(2016)Hartich, Barato, and
  Seifert]Hartich:2016:SensCap
authorD. Hartich,
  authorA. C. Barato,
  and authorU. Seifert,
  journalPhys. Rev. E volume93,
  pages022116 (year2016).

[Brittain et al.(2017)Brittain, Jones,
  and Ouldridge]Brittain:2017:LR
authorR. A. Brittain,
  authorN. S. Jones, and
  authorT. E. Ouldridge,
  journalJ. Stat. Mech. volume2017,
  pages063502 (year2017).

[Munakata and Rosinberg(2013)]Munakata:2013:Feedback
authorT. Munakata and
  authorM. L. Rosinberg,
  journalJ. Stat. Mech. volume2013,
  pagesP06014 (year2013).

[Horowitz and
  Sandberg(2014)]Horowitz:2014:SecLawLike
authorJ. M. Horowitz and
  authorH. Sandberg,
  journalNew J. Phys. volume16,
  pages125007 (year2014).

[Ito and Sagawa(2015)]Ito:2015:Maxwell
authorS. Ito and
  authorT. Sagawa,
  journalNat. Commun. volume6,
  pages7498 (year2015).

[Tomé(2006)]Tome:2006:EPinFPE
authorT. Tomé,
  journalBraz. J. Phys. volume36,
  pages1285 (year2006).

[Schreiber(2000)]Schreiber:2000:TI
authorT. Schreiber,
  journalPhys. Rev. Lett. volume85,
  pages461 (year2000).

[Allahverdyan et al.(2009)Allahverdyan,
  Janzing, and Mahler]Allahverdyan:2009:TE
authorA. E. Allahverdyan,
  authorD. Janzing, and
  authorG. Mahler,
  journalJ. Stat. Mech: Theory Exp.
  volume2009, pagesP09011
  (year2009).

[Hartich et al.(2014)Hartich, Barato, and
  Seifert]Hartich:2014:Bipartite
authorD. Hartich,
  authorA. C. Barato,
  and authorU. Seifert,
  journalJ. Stat. Mech: Theory Exp.
  volume2014, pagesP02016
  (year2014).

[Yeung and Strogatz(1999)]Yeung:1999:DelayedKuramoto
authorM. K. S. Yeung and
  authorS. H. Strogatz,
  journalPhys. Rev. Lett. volume82,
  pages648 (year1999).

[Rosinberg et al.(2015)Rosinberg,
  Munakata, and Tarjus]Rosinberg:2015:Delay
authorM. L. Rosinberg,
  authorT. Munakata, and
  authorG. Tarjus,
  journalPhys. Rev. E volume91,
  pages042114 (year2015).

[Rosinberg et al.(2017)Rosinberg, Tarjus,
  and Munakata]Rosinberg:2017:Delay
authorM. L. Rosinberg,
  authorG. Tarjus, and
  authorT. Munakata,
  journalPhys. Rev. E volume95,
  pages022123 (year2017).

[Guillouzic et al.(1999)Guillouzic,
  L'Heureux, and Longtin]Guillouzic:1999:SDDE
authorS. Guillouzic,
  authorI. L'Heureux,
  and authorA. Longtin,
  journalPhys. Rev. E volume59,
  pages3970 (year1999).

[Risken(1989)]Risken:1989:FPEBook
authorH. Risken,
  titleThe Fokker–Planck Equation: Methods of Solution
  and Applications (publisherSpringer, year1989),
  edition2nd ed.

[Fancher and Mugler(2017)]Fancher:2017:Sensing
authorS. Fancher and
  authorA. Mugler,
  journalPhys. Rev. Lett. volume118,
  pages078101 (year2017).




