Design principles for enhancing phase sensitivity and suppressing phase fluctuations simultaneously in biochemical oscillatory systems

Biological systems need to function accurately in the presence of strong noise and at the same time respond sensitively to subtle external cues. Here we study design principles in biochemical oscillatory circuits to achieve these two seemingly incompatible goals. We show that energy dissipation can enhance phase sensitivity linearly by driving the phase-amplitude coupling and increase timing accuracy by suppressing phase diffusion. Two general design principles in the key underlying reaction loop formed by two antiparallel pathways are found to optimize oscillation performance with a given energy budget: balancing the forward-to-backward flux ratio between the two pathways to reduce phase diffusion and maximizing the net flux of the phase-advancing pathway relative to that of the phase-retreating pathway to enhance phase sensitivity. Experimental evidences consistent with these design principles are found in the circadian clock of cyanobacteria. Future experiments to test the predicted dependence of phase sensitivity on energy dissipation are proposed.

B iochemical processes are innevitably noisy due to the stochastic nature of reactions, the small number of molecules involved, and the thermal fluctuations from environment 1,2 . Various regulatory mechanisms have evolved to suppress effects of noise in order to process information accurately in vital life processes such as biomolecule synthesis 3 , cell cyle 4 , and development 5 . At the same time, some of these systems also need to have a high sensitivity to external stimuli. For example, for many biochemical oscillatory systems, such as glycolysis, cyclic AMP signaling, cell cycle, circadian rhythms, and neural activities 4,6-8 , besides being accurate in their rhythmic timing, they also need to respond sensitively to external cues. In fact, one of the most salient properties of circadian rhythms is their ability to be entrained by the daily cycle in the environment so that their endogenous 24 h cycle can quickly synchronize with environmental signals 9,10 .
However, these two requirements, high sensitivity and low fluctuation, are incompatible for equilibrium systems due to the Fluctuation Dissipation Theorem (FDT) 11 . Briefly, for a perturbation of intensity ϵ applied to the conjugate variable of an observable A at time t = 0, FDT establishes the fluctuationresponse relation (FRR) AðtÞ h i ϵ À A h i 0 = βϵ[C A (t, t) − C A (t, 0)], where C A ðt; sÞ = AðtÞAðsÞ h i 0 is the two-time autocorrelation and β = 1/k B T is the reverse thermal energy. We immediately see that the long time response ΔA ≡ Aðt ¼ 1Þ h i ϵ À A h i 0 is linearly proportional to the variance, i.e., ΔA ¼ Àβϵσ 2 A . This means that a higher sensitivity (ΔA/ϵ) would necessarily lead to a higher fluctuation σ 2 A À Á in any equilibrium system. Such a FRR was also found in certain biochemical systems in their linear response regime 12 .
To understand how living organisms solve the challenge of enhancing sensitivity (responsiveness) and reducing noise (fluctuation) at the same time, we studied the dynamics of a large class of biochemical oscillators in which limit cycles exist with the focus on non-equilibrium effects in the underlying biochemical reaction networks where FRR breaks down. Recently, the relationship between biological regulatory functions and their energy cost has attracted much attention in non-equilibrium statistical physics community [13][14][15][16][17][18] . A previous study found that the phase diffusion constant can be suppressed by a dissipative process that consumes free energy 19 . In this work, by studying different types of limit cycle oscillators analytically and numerically, we investigated whether dissipative processes can enhance sensitivity and reduce fluctuation at the same time. More importantly, our study uncovered the key design principles for biochemical circuits to achieve these two goals simultaneously.

Result
Reduced phase description for biochemical oscillations. The dynamics of a biochemical reaction system {X 1 , X 2 , …, X N }, with fixed volume V and constant temperature, is described by chemical Langevin equation (CLE) 20 . For a biological oscillator, the concentration variable of i th species x i (t) oscillates. Instead of dealing with the entire system, we employ the phase reduction method first developed by Kuramoto [21][22][23] , which reduce Ndimensional state space to a single phase variable ϕ characterizing the timing of oscillation. Specifically, ϕ x L À Á along the deterministic limit cycle L is chosen to progress with a constant speed Ω = 2π/τ for convenience, where τ is the period. This definition of the phase can be extended to the whole basin of attraction of L 24 . If trajectories originated from two states eventually converge onto the limit cycle at the same time, these two states are assigned the same phase. An isochron is a line formed by all points with the same phase (see Fig. 1a).
Clearly, geometrical structure of isochrons is crucial to the phase response property: larger ∇ x ϕ would produce larger phase shifts for the same deviation from limit cycle 25 . In biology literature, a phase response curve (PRC) is commonly used to characterize oscillators' responsiveness [26][27][28][29] . The PRC Δϕ(ϕ) is determined by delivering a perturbation at a given phase ϕ of the oscillation for a given duration of time and comparing the shift in peak times between the perturbed trajectories and the unperturbed ones to obtain Δϕ 30 . Indeed, as shown in Supplementary Note 1, ∇ x ϕ is the key signal-independent factor in determining the amplitude of PRC. At a given phase, we define a dimensionless phase gradient vector ∇ Á is a dimensionless state variable (normalized by the range of variation in x i ). We further use the maximal Euclidean norm of ∇ x Ã ϕ along the limit cycle to define a global phase sensitivity parameter χ: Throughout this work, we assume that change in other signaldependent factors (see Supplementary Note 1) does not overwhelm the effect of χ on phase shift.  Fig. 1 Illustration of the phase response in biological oscillators. a The circle is the assumed stable limit cycle. The gray dashed lines represent equally separated isochrons. An unperturbed system (purple) progresses on the circle, while a perturbed system (cyan) is driven away from the circle by an impulsive signal between time 1 and time 2, and then relaxes back to the limit cycle. At time 2 (end of perturbation), it is moved to an isochron different from the unperturbed one. The difference of their phases determines the phase shift. b, c Diagrams of the signal and phase evolution of the perturbed and unperturbed system. Phase shift is induced during the perturbation and sustains after the perturbation Intuitively, a more sensitive circuit can enhance entrainment as it can more readily change its oscillation phase to sync with external stimuli 31 . It is also straightforward to show analytically that a higher sensitivity widens the range of synchronizable frequencies known as the Arnold tongue 32 and therefore enhances entrainability of biochemical oscillators (see Supplementary Note 2 for details).
Due to the stochastic nature (Poisson process) of the underlying chemical transitions, biochemical oscillations are noisy and the phase fluctuates 19 . The variance of phase fluctuations σ 2 ϕ grows linearly with time, σ 2 ϕ ¼ D ϕ t, which can be used to define a phase diffusion constant D ϕ . It can be shown using the phase reduction method that the probability distribution of phase fluctuation (δϕ) indeed follows a diffusive dynamics ∂ t Pðδϕ; tÞ ¼ D ϕ ∂ 2 δϕ P (see Supplementary Note 3). It is usually hard to directly derive the phase equation. As shown in Supplementary Note 4, the finite correlation time τ c due to phase diffusion is inversely proportional to D ϕ . Thus, we use this relationship to infer the phase diffusion constant from the autocorrelation function, which follows a damped oscillation 19 We now introduce the specific biochemical oscillators we study here. It is known that nonlinear autocatalytic biochemical reactions in open systems can exhibit oscillatory behaviors. For clarity, we employ a simple model (Fig. 2a) derived from glycolysis oscillation 33 and the Brusselator 34,35 as follows: where A, B, and D have fixed concentrations. Here we neglect the inhomogeneous spatial distribution and focus on the dynamics and energetics of the "well-stirred" reaction system. To study thermodynamics of the system properly, we include the backward reactions (see Fig. 2a) in the Brusselator model. An equilibrium steady state can be achieved when k ðeqÞ À2 k À3 ¼ k À2 ¼ k 0 À2 ½D ðeqÞ are pseudo-first-order rate constants and equilibrium values are labeled by superscript (eq). However, if the concentrations of B and D are sustained at values different from their equilibrium values by active processes such as biochemical synthesis and/or active transport (pumping), the system is driven out of equilibrium with the chemical potential difference Δμ DB ¼ Àk B T ln k À2 k À3 =k 2 k 3 ð Þserving as the chemical driving force for the reaction cycle X → Y → X.
To characterize the nonequilibrium cycle dynamics, we introduce a (global) irreversibility parameter γ: which is related to the chemical driving force for the reaction cycle (see Fig. 2). For the special case of an equilibrium system with γ = 1, detailed balance is satisfied and the cycle is fully reversible without net flux. When γ ≠ 1, the system is driven out of equilibrium by external free-energy sources resulting in a nonzero net cycling flux.
In general, external free energy can be utilized to change the chemical reactions in two ways: enhancing forward reactions or suppressing backward reactions. To understand effects of these two distinct changes, we further introduce γ 1 and γ 2 to characterize the (local) irreversibility of the forward and backward reaction respectively, i.e., k 2 ¼ k In addition to the Brusselator model, we have studied another class of biochemical oscillators driven by the general activator-inhibitor (AI) mechanism (see Supplementary Fig. 1). In the AI model, oscillation is driven by the ATP hydrolysis energy and γ −1 can be expressed as ½ATP½ADP ðeqÞ ½Pi ðeqÞ ½ATP ðeqÞ ½ADP½Pi . The two "local" irreversibility parameters (γ 1 and γ 2 ) are also introduced in the AI model to characterize the non-equilibrium effects in different parts of the phosphorylation-dephosphorylation (PdP) cycle that dissipates energy to drive the oscillation. Dissipation outside of the PdP cycle is roughly independent of γ and does not have a direct role in controlling the oscillation (see Supplementary Note 5A and Supplementary Fig. 2 for details).
Effects of free-energy dissipation on phase dynamics. From the chemical reaction rates, we can compute the free-energy dissipation rate (in units of k B T) 36 : where J þ i and J À i are the forward and backward fluxes of the i th reaction. As the dissipation rate also oscillates, we use ΔW ¼ R τ 0 _ Wdt to measure free-energy cost per period per volume. Here we study how ΔW affects the performance of the oscillation as measured by its phase diffusivity and phase sensitivity.
The reversible Brusselator model. a The reaction with B and D can be mapped into a unimolecular reaction with rate constants: Together with the autocatalytic reaction with rates k 3 and k −3 , they form a reaction cycle with a reversibility parameter γ ≡ k −2 k −3 / k 2 k 3 . When γ ≠ 1, the system is nonequilibrium with free-energy dissipation driving the cyclic flux X → Y → X. b The steady-state probability density P(x, y) (color plot) and the state-space fluxes (J x , J y ) (vector field) of the model. The two small blue boxes highlight the regions around the two opposite phases ϕ = 0 and ϕ = π in the deterministic limit cycle. c Details of the chemical reactions in a local region (e.g., the small boxes in b) in the state space. Two microscopic states, (X − 1, Y + 1) and (X, Y), are linked by two distinct reversible reaction pathways, which form a microscopic reaction cycle (loop) First, we briefly summarize the effect of energy dissipation on phase diffusion. As shown in Fig. 3a, b, the minimum phase diffusion constant D min (ΔW) depends inversely on the freeenergy cost per period ΔW consistent with our previous work 19 and the general thermodynamic uncertainty relation for biomolecular processes 37 . Here we further dissect the separate dependence on γ 1 and γ 2 . As shown in Fig. 3a, decreasing γ 1 or γ 2 can both reduce phase diffusion. However, D eventually saturates to a nonzero value even in the absence of backward reactions (γ 2 = 0) due to the finite stochasticity in the forward reactions with a finitie γ 1 (see Supplementary Fig. 3b). On the other hand, D seems to decrease continuously with γ 1 (see Supplementary Fig. 3c). As the number of forward reactions per unit time N c increases with 1/γ 1 , the averaging effect (1= ffiffiffiffiffi ffi N c p effect) in reducing the phase fluctuation can persist for small γ 1 (or large N c ). In fact, this noise reduction strategy, i.e., taking average over multiple irreversible steps, is also commonly employed in other biological processes 38 .
We now turn to the main focus of this work, i.e., to understand how free-energy dissipation affects phase response. In particular, we ask whether low phase fluctuation and high phase sensitivity can coexist in a dissipative system. For different values of the forward and backward irreversible constants (γ 1 and γ 2 ), we calculated phase sensitivity χ of the biochemical oscillator to external perturbations (stimuli). Remarkably, we observed in Fig. 3c, d that the phase sensitivity χ is bounded by a maximum value χ max that increases linearly with ΔW in a wide range of ΔW for different combinations of γ 1 and γ 2 : where K W is a constant whose value is given in the legend of Fig. 3.
We have confirmed the generality of Eq. (6) for other implementations of γ 1 , γ 2 ( Supplementary Fig. 4) in the Brusselator model, as well as for the AI model ( Supplementary Fig. 5).
The relation between phase sensitivity and entrainment to external periodic driving. Many biochemical oscillators are exposed to an external periodic signal ϵp(Ω, t) (e.g., temperature, light, etc.) that entrains the internal oscillation. The external signal varies with the system's intrinsic frequency Ω and acts on an internal parameter μ. Incorporating this periodic driving force as time-dependent perturbation f(t) = (∂ μ F)ϵp(Ω, t) in the phase reduction description, we obtain the dynamical equation for the phase ϕ where Z μ (ϕ) is the infinitesimal PRC function for perturbing μ 25 . Applying a pulse perturbation with duration Δt = δ t τ and intensity Δμ = δ μ on μ, the PRC can be calculated as Δϕ(ϕ; δ t , τ, δ μ , μ) = Z μ (ϕ)ΔtΔμ. We now relate the normalized PRC C(ϕ) = Δϕ(ϕ)/(δ t δ μ ) = Z μ (ϕ)τ with the performance of entrainment to the external periodic driving signal.
To characterize the phase dynamics in the presence of an external periodic signal, we use the phase difference ψ = ϕ − Ωt between the internal oscillator and the external signal. By averaging over a period (as ψ is a slow variable when ϵ is small), we have Entrainment to the external signal occurs, because Eq. (8) has a stable fixed point ψ 0 , i.e., Γ(ψ 0 ) = 0, Γ´(ψ 0 ) < 0. Now, we consider a sudden phase shift Δψ of the external signal, which is equivalent to an initial perturbation Δψ in ψ away from its fixed point. When Δψ is small, the perturbation δψ = ψ − ψ 0 follows, to the leading order, _ δψ ¼ Àt À1 e δψ where t e ¼ ϵΓ′ ψ 0 À Á À1 τ is the entrainment time for the system to catch up the shifted signal. It is clear that  19 . The fitting parameters are: W c = 41.29, W 0 = 36.44, and C = 0.2305. The volume of Gillespie simulation is V = 100. d χ exhibits linear dependence on ΔW and the maximum sensitivity χ max = K W ΔW + const. with the proportional constant K W ≈ 6.5 × 10 −3 stronger signals (larger ϵ) shorten the entrainment process. At the same time, the internal phase response property also has an important role. In general, both the steady state phase difference ψ 0 and the entrainment time t e depend on the shape of PRC (C(ϕ)) and how the periodic signal is applied (i.e., the waveform p(Ω, t)); however, as long as the shape of PRC remains approximately the same, enhancing phase sensitivity χ (the amplitude of PRC) reduces t e : We have simulated the AI model to test the above relation between the phase sensitivity χ and the entrainment time t e for a periodic driving k ðpÞ 1 ðtÞ ¼ k 1 ð1 þ ϵ cos ΩtÞ. At different phases of the oscillation, a phase shift Δψ = 0.4π was applied to the external signal, triggering the entrainment process. The oscillator was considered "entrained" when the peak difference between a perturbed and unperturbed system is small δψ < 0.1Δψ. Minimum number of entrainment cycles min ϕ t e f g=τ was recorded for different γ 1 and γ 2 . The results are shown in Fig. 4a, b for ϵ = 0.1.
The external periodic signal also affects phase fluctuations, which is now bounded by a "potential" well established by the external signal (Fig. 4c). By approximating the potential well by its harmonic form, the stochastic phase evolution equation in the presence of noise (from stochastic chemical reactions) can be solved to determine the phase variance: where D is the previously defined dimensionless phase diffusion constant in the absence of the external signal. For time t ≪ t e , the phase variance follows diffusion δψ 2 h i ¼ Dt=τ. For t ≫ t e , the phase variance reaches a constant σ 2 ψ ¼ Dt e =τ. From Eq. (10), it is clear that phase fluctuation is suppressed by reducing D even in the presence of external driving signal.
By using the Gillespie algorithm 39 , we also computed the normalized auto-correlation function: C xx ðtÞ ¼ xðt þ sÞxðsÞ h i s / x 2 h i ¼ expðÀσðtÞÞcosðΩtÞ, where x is a state variable. As shown in Fig. 4c, the amplitude of a typical autocorrelation function C xx (t) first decreases exponentially before reaching a steady state constant A ¼ exp ÀDt e =τ ð Þ , which verifies Eq. (10). Based on the value of A from the simulations, we obtain the phase variance σ 2 ψ , which gives D ¼ σ 2 ψ τ=t e . As shown in Fig. 4d, the dimensionless diffusion constant D decreases with energy dissipation ΔW consistent with the case without external signal (Fig. 3b).
Energy-enhanced phase-amplitude coupling leads to higher phase sensitivity. To understand the dependence of phase sensitivity on free-energy consumption analytically, we study the normal form of limit cycle oscillation originated from a Hopf bifurcation. Applying variable transformation, we obtain from CLE the stochastic Stuart-Landau equation 40 where μ, ω, β 1 , β 2 are parameters associated with original reaction rate constants and ξ r,θ are unit variance white noise terms. We note that β 2 characterizes the phase-amplitude coupling strength that affects the phase dynamics. By using the stochastic averaging  41 , we can show the two noise strength q r and q θ to be the same q 2 r ¼ q 2 θ Q to the leading order in amplitude. See Supplementary Note 6 for detailed derivations.
For μ > 0, we have the mean amplitude r s ¼ ffiffiffiffiffiffiffiffiffi ffi μ=β 1 p and angular velocity Ω ¼ ω þ β 2 r 2 s . The isochron with the phase θ can be determined from the mean-field solution (see Supplementary Note 7 for details) Thus, ∇ x ϕ ¼ À =r e r þ ð1=rÞe θ and the phase sensitivity can be calculated The dissipation of Stuart-Landau oscillators can be determined by solving the Fokker-Plank equation for (11) 19 , which results in the steady state probability distribution where Δ = Q/V and A is the normalization coefficient. To have a well-defined phase variable requires that the amplitude fluctuation is much smaller than r 2 s , leading to a constraint ρ μ= ffiffiffiffiffiffiffiffiffiffi 2β 1 Δ p ) 1 19 . The entropy production rate computed from P ss (r, θ) yields the minimum free-energy cost per cycle 42 (k B T = 1, see also Supplementary Note 8 for details).
where κ ¼ ω=β 2 r 2 s is a dimensionless parameter. Finally, assuming that the phase sensitivity is dominated by the radial contribution, i.e., β 2 /β 1 ≫ 1, we have χ ≈ β 2 /β 1 and κ ≪ 1. Inserting this into Eq. (15), we obtain that to the leading order: where K W ≈ [4πρ 2 ] −1 is a constant independent of β 2 . The linear relation (16) agrees with our numerical results. More importantly, the analytical results reveal that the amplitude-phase coupling constant β 2 has an important role in relating energy dissipation (Eq. (15)) to phase sensitivity (Eq. (13)). We further develop a toy model of limit cycle oscillator (see Supplementary Fig. 6) where β 2 can be explicitly calculated in terms of the microscopic nonequilibrium parameter γ. See Supplementary Note 9 for details. We found, in this simple model, that free-energy cost of increasing forward rates and suppressing backward rates both strengthen the phaseamplitude coupling β 2 , therefore enhancing the phase sensitivity.
Design principles for enhancing oscillation functions. Now that we show it is in principle possible to increase phase sensitivity and to suppress phase diffusion simultaneously in a nonequilibrium system that consumes free energy, the next logical question is what are the design principles for a biochemical oscillator to optimize desirable oscillatory behaviors with a fixed energy budget. Here we search for possible design principles by studying the performance of an oscillator, characterized by χ and D, for different combinations of γ 1 and γ 2 that lead to the same free-energy dissipation per period. In particular, we look for rules for designing the kinetic rate parameters in reactions that form a microscopic loop like the one shown in Fig. 2c. Unlike in simplified models 43 where only one reversible reaction exists between two states, these microscopic loops, which are formed by two distinct reaction pathways between two nearby microscopic states, are the basic building blocks in realistic oscillatory biochemical networks.
Balance the forward-to-backward ratio in antiparallel pathways to suppress phase diffusion. Consider N discrete states with equally spaced phases ϕ(n) = 2πn/N along a limit cycle, we study the case in which there are two antiparallel pathways i = 1, 2 between two neighboring states. For each pathway, its forward and backward transition rates are given by w ± i , respectively. The overall effect of the i th transition is denoted by the net rate w i ðϕÞ ¼ w þ i À w À i . These two pathways are antiparallel in the sense that w þ 1 and w þ 2 are in the opposite direction and they together form a counter-clockwise loop (so do w À 1 and w À 2 in the clockwise direction; see Fig. 2c for the Brusselator model).
It is easy to show that the average change of phase and its variance are 44 : where the subscript f, b can be either i = 1 or i = 2 depending on whether a given transition pathway i contributes to forward or backward movement of the phase. As shown in Fig. 2c for the Brusselator model, a given reaction pathway, either the B þ X"D þ Y reaction (i = 1) or the 2X þ Y"3X reaction (i = 2), can lead to a forward or a backward phase change depending on which part of the limit cycle (or phase) the transition occurs. For example, f = 1, b = 2 at phase ϕ = 0, whereas f = 2, b = 1 at phase ϕ = π (see Fig. 2c). From the above Eqs. (17) and (18), we obtain a dimensionless phase diffusion constant Our goal is to minimize phase diffusion under the constraint of a fixed energy dissipation ln w þ 1 w þ 2 =w À 1 w À 2 À Á = ln γ À1 $ ΔW=N. For simplicity, we assume the ratio w þ i =w À i take the same value independent of the phase, i.e., w þ i ðϕÞ=w À i ðϕÞ g i . The fixed energy constraint is now given by g 1 g 2 = γ −1 . By using w f (ϕ) = w b (ϕ) + Ω N and the expressions for For an ideal symmetric clock, we can assume g b j ϕ ¼ g f ϕþπ , and the summation of the first term in Eq. (20) can be expressed as P , which leads to a new expression for D and consequently its lower bound: where the lower bound is reached when g 1 = g 2 = γ −1/2 . Here we have implicitly assumed that the average of normalized backward rate N À1 P ϕ w b ðϕÞ=Ω N has weak dependence on g 1 , g 2 . Our analysis clearly suggests that one design principle to minimize the phase diffusion of biochemical oscillators with a given energy budget is to equalize (balance) the ratio between forward and backward transition rates in the two antiparallel pathways that form the (non-equilibrium) dissipative loop.
This design principle obtained from an ideal clock is consistent with our previous findings in the AI model 19 . We have also directly tested it by numerical simulations of the Brusselator model, where two adjacent states (X − 1, Y + 1) and (X, Y) are linked by a non-equilibrium reaction loop (see Fig. 2c). We compute the ratio r of the forward and backward flux of two reactions, respectively: 3 . We can then define M ≡ max r 1 ðtÞ h i τ ; r 2 ðtÞ h i τ À Á /min r 1 ðtÞ h i τ ; r 2 ðtÞ h i τ À Á to measure the matching of these two pathways averaged over the period. As shown in Fig. 5a Minimize the net backward flux relative to the net forward flux to enhance phase sensitivity. To look for strategies in enhancing phase sensitivity, we calculate k ∇ x Ã ϕ k by projecting the transition fluxes in the state space onto the phase (angular) direction and the amplitude (radial) directions. In addition to the reactions that form part of the loop (such as w 1 and w 2 shown in Fig. 2c), we also include the reaction that is not part of the dissipative loop, e.g., the A"X reaction in the Brusselator model. The net flux for this non-dissipative reaction is denoted by w 0 , which is not regulated by the non-equilibrium parameters (γ 1,2 ).
Recall that phase sensitivity χ measures the largest k ∇ x Ã ϕ k along the limit cycle. It can be shown that the maximum phase gradient occurs at the point in phase space where the energydriven fluxes most closely align with the phase direction (see Supplementary Note 10 for a more detailed analysis). A corresponding approximation is that the transition rate in the radial direction near the most sensitive point is controlled by w 0 . Thus, we can estimate χ as: where q i ¼ Δwi=wi Δr=rs is the relative change of transition rates induced by a perturbation given by a relative change Δr/r s of the amplitude. q i characterizes the sensitivity of the i th transition determined by the nonlinearity in the underlying reaction rates. In biochemical oscillators, for example, reactions of (pseudo) first order have q i~1 ; the autocatalytic reactions in the Brusselator model give q i~3 .
The dominant effect of energy dissipation is to increase w f , which roughly takes the form c 1 /γ 1 − γ 2 c 2 , where c 1,2 are concentration-related coefficients for biochemical circuits. At a fixed energy dissipation or a fixed w f /w 0 , Eq. (23) suggests that reducing the net backward flux relative to the net forward flux, i.e., minimizing r w ≡ w b /w f can lead to higher phase sensitivity.
This design principle, based largely on heuristic arguments, has been tested directly in the Brusselator model where the two net fluxes are We calculated r w ¼ min w 1 ; w 2 ð Þ=max w 1 ; w 2 ð Þnumerically along the entire limit cycle and compared min r w f g for different combinations of γ 1 and γ 2 for a fixed dissipation ΔW. As shown in Fig. 5c, d, phase sensitivity reaches its peak when r w is small, which confirms the strategy (design principle) for maximizing phase sensitivity by a maximum separation of net transition fluxes between anti-parallel pathways in a dissipative loop. Similar results were confirmed in the AI model (see Supplementary Fig. 9 for details). Experimental evidence in the Kai system. To look for experimental evidence of the design principles described above, we investigate the circadian clock of cyanobacteria, i.e., the Kai system, by studying the experimentally measured kinetic rates in the reactions between different phosphorylation states (U, T, D, S) of KaiC along its PdP cycle [45][46][47] .
To test the first design principle, we calculate the ratio of forward to backward flux r X⇆Y = J X!Y =J Y!X on each of the four links in the KaiC phosphorylation cycle U → T → D → S → U based on a kinetic model of KaiC phosphorylation with experimentally measured rate parameters [45][46][47] (see Supplementary Note 11 and Supplementary Fig. 10a for details of the model and Supplementary Table S1 for the parameters). We then compare the periodaveraged ratio 〈r X⇆Y 〉 τ to check whether these ratios are matched along the cycle. The KaiC model with the experimentally determined parameters (see Supplementary Information for  details) yields 〈r U⇆T 〉 τ :〈r T⇆D 〉 τ :〈r D⇆S 〉 τ :〈r S⇆U 〉 τ = 1:1.16:0.97:1.08, clearly showing that the forward-to-backward ratio in different links (pathways) of the KaiC phosphorylation cycle are properly balanced, which is consistent with the first design principle for minimizing phase diffusion.
To test the second design principle, we calculate the net flux for each link J ðnetÞ X!Y ¼ J X!Y À J Y!X . The net phosphorylation and dephosphorylation flux are then approximated as J S!U , respectively. From direct simulation of the model, the backward-to-forward (in terms of phosphorylation rhythm) net flux ratio r w is smaller during the subjective day with an average r w h i day % 0:2 than during the night with average r w h i night % 0:5. According to the second design principle, this r w behavior leads to a higher phase sensitivity during the subjective day (phosphorylation phase) than during the night (dephosphorylation phase), which is consistent with the experimentally measured PRC reported in ref. 46 (see Supplementary Note 11 and Supplementary Fig. 10 for details). This result confirms that the Kai system controls phase sensitivity by modulating the relative strength of the phosphorylation and dephosphorylation fluxes.

Discussion
In this study, we investigated whether and how biochemical systems can achieve high sensitivity and low fluctuation at the same time in the context of biological oscillators. In nonequilibrium systems, the FDT is broken. There is no unique relationship between fluctuation and response-they could have a positive correlation or a negative one depending on which parameters are varied. For example, in the Stuart-Landau model, if we vary a single variable without changing the others, χ and D satisfy a positively correlated relation: χ 2 = D/T eff . However, even in this case, the effective temperature T eff could be lowered by dissipation ΔW in different ways depending on which variable (μ or β 2 ) is varied (see Supplementary Note 12 for details). For a realistic non-equilibrium biochemical reaction network (circuit) such as the Brusselator model, where many parameters are affected by the dissipation in a correlated manner, our study shows that an increase in free-energy dissipation can lead to both a suppression of phase fluctuation and an enhancement of the phase sensitivity if it is used properly (or the circuit is designed properly). This result should be generally applicable to other nonequilibrium systems.
Our study revealed two design principles to achieve optimal performance with a finite budget of free-energy consumption in oscillatory systems. The findings of our analysis can be demonstrated by the expression of dissipation rate In principle, high sensitivity and low phase fluctuations for a given dissipation rate can be achieved by favoring the net forward flux w f over the net backward flux w b and balancing the forward-to-backward ratio w þ i =w À i among antiparallel reaction pathways, respectively. As these two design principles act on different combinations of the reaction rates, they can be satisfied simultaneously leading to biochemical circuits that have both high sensitivity and low fluctuations. Strong evidence in support of these two design principles are found in the Kai system. Aside from helping us understand the structure and dynamics of naturally occurring biological pathways, these design principles may also serve as the best practice rules for constructing efficient synthetic biochemical circuits for oscillations.
Other evidence for our theory may be found in experiments measuring PRC at different energetic states of the system. For example, our theory predicts that the PRC amplitude should decrease when the background ATP/ADP ratio is decreased. As far as we know, no such experiment has been done. The closest are PRC measurements (to light and drug/chemical perturbations) at different temperatures. In our analysis, the dimensionless energy dissipation ΔW has a temperature-dependent component βΔG (0) , where β = (k B T) −1 is the inverse thermal energy and ΔG (0) is a concentration-independent free-energy difference 48 . According to Eq. (6), this means the phase sensitivity increases as the temperature decreases if we neglect other  Fig. 6 Experimental connection and prediction of the dependence of phase sensitivity on irreversibility γ. a By assuming k 1 is increased by 70% and k 4 by 30% during a temperature up perturbation (0.1 period), the PRC from the AI model (black line) can be fitted to that of a temperature up PRC (symbols) observed in 57 . Simulation trajectory of oscillation along with the circadian time scale are shown in the inset. The peak of activator rhythm is defined as CT 0. Experimental PRC data from ref. 57 is aligned with the model PRC to peak at the same CT. b PRCs for the same perturbation as in a is computed for different background levels of ATP/ADP ratio or equivalently different values of γ (given in the legend). The extra free-energy dissipation for each γ choice are: ΔW − W c ≈ 170 (blue); 250 (red); 350 (black, same as in a). Amplitude of PRC increases as sensitivity χ is enhanced by Àlog γ or ΔW as shown in the inset, where the hollow circles correspond to different choices of (γ 1 , γ 2 ). This prediction can be tested experimentally temperature dependence. Interestingly, in the experimental systems we found, i.e., Neurospora 49 , Gonyaulax polyedra 50 , and chick pineal cell 51 , except for the case of light response in chick pineal cell where no T-dependence is found, PRC amplitudes indeed increase with 1/T (see Supplementary Fig. 11). However, many internal variables in biological systems can depend on temperature, therefore these temperature dependence measurements in in-vivo systems may not serve as direct tests of our theory.
The best system to test our theory directly is the relatively simple cyanobacterial circadian clock, especially the in vitro Kai system. Experimental studies in cyanobacterial circadian clocks uncover that metabolism is the fundamental synchronizer 46,47,52 . However, there is yet no direct experiments investigating the relationship between phase sensitivity and the background metabolic state of the system. Away from the onset of oscillation, the relative amplitude change Δr/r caused by a modest phase resetting signal (e.g., temperature pulse) should has only a weak dependence on the metabolic state of system. Under this assumption, our theory shows that the phase response is dominated by χ, which increases with the background ATP/ADP ratio. In particular, we predict that the invitro KaiABC oscillators with a lower background ATP/ADP ratio would generate a smaller PRC amplitude in response to the same temperature perturbation. In Fig. 6, we show the PRC in response to a particular change of kinetic rate constants in a generic AI model, to which the Kai system belongs. By choosing the reaction rate changes during the temperature perturbation properly, we can reproduce the observed PRC to a temperature pulse as shown in Fig. 6a. We then use our model to compute the PRC to the same perturbation (rate changes) but with different values of the energy dissipation per period (ΔW) or different background ATP/ADP ratios. As shown in Fig. 6b, the amplitude of the PRC decreases as ΔW decreases. This result is robust as the enhancement of PRC amplitude by free-energy dissipation holds true generally regardless of the specific perturbation and the specific model used (see Supplementary Fig. 12). This prediction can be tested in future experiments, e.g., by measuring the PRC for changes in temperature in different background ATP/ADP ratios. These types of experiments would not only test our hypothesis that energy dissipation affects phase response, they may also shed light on its molecular mechanisms 53 and how cellular metabolic activities affects other crucial functions of biological clocks, such as temperature compensation [54][55][56] .
Data availability. The data that support the findings of this study are available from the corresponding author on request.