Reciprocity Between Robustness of Period and Plasticity of Phase in Biological Clocks

Circadian clocks exhibit the robustness of period and plasticity of phase against environmental changes such as temperature and nutrient conditions. Thus far, however, it is unclear how both are simultaneously achieved. By investigating distinct models of circadian clocks, we demonstrate reci- procity between robustness and plasticity: higher robustness in the period implies higher plasticity in the phase, where changes in period and in phase follow a linear relationship with a negative coef- ficient. The robustness of period is achieved by the adaptation on the limit cycle via a concentration change of a buffer molecule, whose temporal change leads to a phase shift following a shift of the limit-cycle orbit in phase space. Generality of reciprocity in clocks with the adaptation mechanism is confirmed with theoretical analysis of simple models, while biological significance is discussed.

Such biological clocks often work as pacemakers, to adapt to periodic events. One of the most prominent examples of such oscillators is a circadian clock [1,2]. To respond to periodic events, the following two criteria are generally imposed on a biochemical oscillator. 1. Robustness of period: If the period of an oscillator strongly depends on external conditions, the oscillator would not accurately predict time. For example, if the period of a circadian clock is sensitive to temperature, the clock malfunctions depending on the temperature. To avoid such error, the period of pacemakers should not be affected by external conditions such as temperature and nutrient compensation [3,4]. 2. Plasticity of phase: The period of the circadian clock of most organisms is known not to correspond precisely with 24 hours [5], and biological clocks are entrained with the external 24-hr cycle [6], so that the phase difference between the two does not increase with time. This entrainment is also necessary to adapt an abrupt change in the environment that may cause temporal misalignment between the internal and external cycles. For such entrainment, plasticity of the phase of the internal clock against external stimuli, e.g., changes in temperature and/or brightness, is needed.
Indeed, biological clocks satisfy both robustness and plasticity to changes in factors such as temperature and nutrient conditions, which change in the daily cycle. For example, circadian clocks of in vivo Drosophila [7], Neurospora [8], and in vitro cyanobacteria [9,10] show temperature compensation of a period and are entrained by cyclic temperature changes. Robustness of period is also important to stable entrainment since it can reduce the difference between the period of inner clock and external cycle. In spite of some studies discussing the compatibility between the two properties [7,[11][12][13], however, little is known about the quantitative relationship between the two properties.
To answer how the robustness of the period and plasticity of phase are compatible with each other, we first analyze two major models of a circadian clock, i.e., post-translational oscillator (PTO) [9,14,15] and transcription-translation-based oscillator (TTO) [12,13,15], which consists only of protein-protein interactions and both transcription and translation processes, respectively. Without imposing any special mechanism, we demonstrate that biological clocks with robustness of period against changes in an environmental factor generally exhibit phase entrainment against the cyclic change of that factor -reciprocity between the robustness of period and plasticity of phase: the plasticity increases with robustness.
For PTO model, we adopt the KaiC allosteric model [16], for in vitro cyanobacterial circadian clock system [9]. Here, KaiC protein consists of six monomers, each of which has a phosphorylation site. The protein has active and inactive forms. Active (inactive) KaiC are phosphorylated (dephosphorylated) step by step, respectively. Phosphorylation reactions are facilitated with KaiA as an enzyme and dephosphorylation reactions spontaneously progress without an enzyme. k p and k dp denote the rate of phosphorylation and dephosphorylation of KaiC, respectively, which depend on temperature as k p ∝ exp(−βE p ) and k dp ∝ exp(−βE dp ), where E p (E dp ) is the activation energy for phosphorylation (dephosphorylatiion), respectively, and β is the inverse temperature by taking the Boltzmann constant as unity. The temporal evolution of the concentration of each phosphorylated active (inactive) KaiC is given by rate equations (see model equations and Fig.1A of [17]).
This model shows a limit-cycle attractor in which the total phosphorylation level, i.e., the ratio of phosphory- Difference between periods at two temperatures (β1 = 1.0 and β2 = 1.5) (∆T /T , red circle) and the amplitude of the phase response curve against a transient jump of temperature from β1 to β2 (∆φ, green square) plotted against different values of E dp while Ep is fixed at 1.0. ∆T is normalized by the period at β1, and ∆φ is normalized by the duration of stimulus and difference between β1 and β2. ∆T /T and ∆φ are negatively correlated across the entire range of E dp . (B) Entrainability is plotted against various E dp with fixed Ep = 1.0. (C) Phase response curve against transient increase in β. As a stimulus, the inverse temperature β is increased from β1 to β2 for the duration of one unit of time. Lines of different colors represent the PRCs for different values of activation energy for the dephosphorylation reaction E dp . lated monomers, oscillates in time. We demonstrated that the robustness of the period against various environmental changes is achieved by enzyme-limited competition [18,19]: With the increase in temperature, the abundance of the active form of the KaiC molecule increases, which in turn decreases the abundance of the free KaiA molecule, and thus the increase in the rate of phosphorylation k p is canceled out, when the total KaiA amount, A total , is sufficiently small. This robustness in the period is achieved when E p is sufficiently larger than E dp . We use the difference in periods between two different temperature conditions (∆T /T ) as an indicator of the robustness of period. Its dependence upon E dp − E p is given in Fig.1A.
This clock, on the other hand, entrains against external periodic change, so that the phase of phosphorylation oscillator coincides with that of external cycle. By imposing external periodic change in temperature, we computed how many number of cycles are needed for the clock to entrain with this external cycle, and defined entrainability as the inverse of the number (see [17]). Dependence of the entrainability and ∆T /T upon E dp with fixed E p is plotted in Fig.1A (red circle) and B. As E dp − E p is smaller, ∆T /T becomes smaller and the entrainability is higher. In other words, if the period of the clock is more robust against temperature change, it is entrained faster with the external temperature cycle, i.e., the phase has higher plasticity.
Although this demonstrated the correlation between period robustness and phase plasticity, the entrainability here is a complicated indicator for the latter, as it can depend on the form of external cycle. Hence, we introduce a more tractable indicator for the plasticity of phase, by using a phase response curve (PRC) [20]. PRC is a function of phase and represents a phase shift introduced by a transient stimulus. When a transient stimulus is added to an oscillatory system, the period of oscillation is temporally altered depending on the phase when the stimulus was added. The period finally returns to its original value. In this time, the phase of the oscillator progresses (or is delayed) from the original phase because of the temporal shortening (or lengthening) of the period. PRC represents such a phase shift ∆φ as a function of the phase φ when the stimulus is applied. We computed PRC by transiently changing the inverse temperature from β 1 to β 2 for one time unit (see Fig.1C), by defining the phase of oscillation by the time when the total phosphorylation level takes maximum at φ = 0, 2π, · · · . As an indicator of the plasticity of phase, we measured the difference between maximum and minimum values of the phase change ∆φ in PRC [21] normalized by the magnitude of a stimulus by fixing its duration as one time unit. The dependence of ∆φ and ∆T /T on E dp with fixed E p is plotted in Fig.1A. When E dp is low, i.e., when the temperature dependence of dephosphorylation reaction is weak, ∆T /T is small and ∆φ is large. This reciprocity was also obtained against changes in other parameters, β 1 and A total (see Fig.3 of [17]). This indicates that a biochemical oscillator with a homeostatic period against an environmental change can easily shift its phase under the same environmental change. We also confirmed such reciprocity against change in ATP, i.e., the case of nutrient compensation (see Fig.5

of [17]).
Now, we examine if such reciprocity holds in the other class of circadian clocks, the TTO. In the TTO model, a clock-related gene is first transcribed and translated, and later such a translated protein represses the expression of its own gene with a time-delay. When the transcription rate decreases, the amount of such protein also decreases, which weakens the suppression of the clockrelated gene expression. Consequently, such genes are transcribed again, leading to the oscillation of the gene expression level. As a typical example of the TTO model, we choose here a model of a circadian clock of a fruit fly [22] (see model equations and Fig.1B of [17]). By varying the activation energy for mRNA degradation, E a , and fixing activation energies for other reactions, we measured ∆φ and ∆T /T using the same procedure as in the Kai model. Then, ∆T /T is low and ∆φ is high for a low E a value, and ∆T /T (∆φ) increases (decreases) with the increase in E a ( Fig.2 and  the reciprocity holds also in the TTO.
To discuss the reciprocity analytically, we then study the Stuart-Landau model, a minimal model for simple sinusoidal oscillation [23]. The model consists of the amplitude R and argument Θ, where R andΘ reach a constant value at the limit-cycle attractor. Indeed, this model is derived as a normal form close to the Hopf bifurcation point. We introduce an external parameter β: where, f 1 (β) is a response function of the first order term in complex Ginzburg-Landau equation, and f 2 (β) is that of the third order term (for choice of each functions, see [17]). Considering the stability of limit cycle, the relaxation of R after perturbation is assumed to be much faster than that of Θ [24]. Here, the period is given as: Thus, after the change β → β + ∆β, the dependence of period on β is given as: Here, we neglected higher order terms of ∆β, assuming that it is sufficiently smaller than β.
) is satisfied, the dependence of the period on f 2 (β) will be counterbalanced by f 1 (β), and the period is compensated against a change in β.
The argument Θ is defined only on a limit-cycle orbit, and we introduce the phase φ to extend the definition to the phase space out of the limit-cycle attractor, in particular to its basin. It is postulated that φ agrees with Θ on the limit-cycle orbit, i.e., different orbits from the same φ converge to the same point on the limit cycle having the same Θ value. Now, we will derive an isochrone, which is a set of points with the same φ on the phase space. φ is expected to have rotational symmetry, hence the isochrone of the Stuart-Landau equation against the parameter β is derived as (see a supplemental text of [17].) Then, we consider an operation that increases β from β 0 to β 0 +∆β and instantaneously reverses it to β 0 . By assuming that R instantaneously relaxes to R * (β 0 + ∆β) = (f 1 (β 0 + ∆β)) 1/2 while Θ remains unchanged, the phase after the above operation is derived as (5) Hence, when ∆β ≪ β, the change in phase is derived as: Therefore, from Eqs. (3) and (16), changes in the period and phase are represented by an equality. where which depend only on f 2 (β) and not on f 1 (β). Thus, when we construct f 1 (β), which compensates for the dependence of f 2 (β) on β according to Eq. (3), the phase is altered as ∆φ = c. On the other hand, when f 1 (β) is independent of β, the phase is also independent due to Eq. (16), while the period is strongly dependent on β as ∆ ln T = c/a = −∆f 2 (β)/(ω + f 2 (β)). We also confirmed the reciprocity is valid in the modified van der Pol oscillator [25] with strong nonlinearlity, i.e., beyond the neighborhood of Hopf bifurcation (See Fig.9 of [17]).
The origin of reciprocity is also understood from the viewpoint of adaptation motif. The standard minimal feedforward motif for adaptation consists of two components, x and y [26]. In the feedforward network in Fig.3A, an input changes both components x and y, while y gives an input to x. Here, the direct path to x and the indirect path via y from the input have opposite signs. Then, the response of the output x via the direct path is later canceled by y, and the adaptation behavior against the input is shaped. The degree of adaptation depends on the strength of the indirect regulation; weak regulation induces a partial adaptation and strong regulation leading to the cancellation of the two paths, induces perfect adaptation [27,28].
Our Stuart-Landau model also has a feedforward motif consisting of amplitude and angular velocity. When an environmental condition β is changed, the angular velocity and amplitude are altered by the terms f 1 (β) and f 2 (β). After a direct change in angular velocity, the  change is relaxed by the change in amplitude. The period is determined as the inverse of the angular velocity.
If changes in the amplitude are large, period is perfectly compensated and phase is plastic. In contrast, if the change in amplitude is small, the angular velocity shows partial adaptation leading to partial compensation of the period while the phase is only slightly altered. Therefore, the reciprocity is understood as the adaptation dynamics on a limit cycle. Indeed, the above argument of the adaptation on the limit cycle generally holds, for PTO and TTO models, where we can generally consider the scheme of Fig.3A. Environmental change directly influences the angular velocity while it is also buffered in the amplitude and then influences the phase. In a biochemical clock, the period mainly depends on the rate-limit reactions, which are slower than others. Environmental change will alter the speed of such rate-limit reactions, which is later counterbalanced by the change in the concentration of buffer molecules [29]. In fact, in the PTO model, the amount of free enzyme working as a buffer molecule can counterbalance the speed of the rate-limit reaction. Hence, the period of the oscillator is homeostatic against environmental changes. Likewise, in TTO model, mRNA plays the role of such buffer molecule.
In this time, the limit-cycle orbit of oscillators with compensation shifts in the phase space of chemical concentrations to change that of a buffer molecule (see Fig.3B). When homeostatic response is achieved, the concentration of a buffer molecule x should be changed with ∆x by the change in the external environment. Then the limit-cycle orbit will be shifted to change the con-centration of a buffer molecule, and the magnitude of such shift and the change in isocline will be O(∆x) considering that continuous change in the isocline against ∆x which is small. Then, ∆φ ∝ ∆x is expected. On the other hand, when the change in the concentration of a buffer molecule is not sufficient to counterbalance the environmental stimulus, the concentrations of other molecules will change. Let us represent the concentration of x needed for perfect adaptation as ∆x * . Then, the change in the concentration of the other molecule of the lowest order is proportional to ∆x * −∆x. The period also changes accordingly, so that ∆T /T ∝ ∆x * − ∆x is expected. By combining the two proportionally relationships, we obtain a∆φ + b∆T /T = ∆x * with coefficients of proportionality a and b.
We have shown that reciprocity exists in both the PTO and TTO models. The currently known mechanisms of circadian oscillation can be classified into the above two cases [15], and the reciprocity is expected to be achieved universally in circadian clocks [30]. In a circadian clock system of a mold, Neurospora crassa, it was reported that a loss-of-temperature-compensation mutant, frq-7, shows smaller phase shift against transient temperature change than the wild type [8,31,32]. Although the quantitative relationship between temperature compensation and phase plasticity was not investigated therein, we expect that a quantitative experiment will confirm our reciprocity, not only in Neurospora crassa but also in other organisms in which lossof-temperature-compensation mutants are isolated, e.g., fruit fly [33] and cyanobacteria [34]. Here, we demonstrated the reciprocity against changes in the temperature and the nutrient concentration, but from theoretical consideration, it is expected to hold generally against a variety of stimuli, such as the change in strength of light and transcription rate [35,36], as long as the adaptation mechanism works. Moreover, it is also expected that the reciprocity is not limited to the circadian clock; it holds generally as long as the adaptation mechanism with buffering molecules works [37]. Our reciprocity will give a general quantitative law for such adaptation systems.
This We introduce the KaiC allosteric model [16] (Fig.  1A). The KaiC protein has six monomers, and each monomer has multiple phosphorylation sites. Here, we assumed each KaiC monomers have only two phosphorylation states, phosphorylated and unphosphorylated. KaiC hexamer takes an active or inactive form. By denoting C i andC i as active and inactive forms with the i phosphorylated monomers, respectively, their temporal changes are given as: where A denotes the free KaiA protein that works as an enzyme for phosphorylation. A total is the total KaiA amount, which is a constant because the total amounts of both KaiC and KaiA are conserved quantities, and [x] denotes the concentration of x. K i is the dissociation constant between C i and A. k p and k dp denote the rate of phosphorylation and dephosphorylation of KaiC, respectively, which depend on temperature as k p ∝ exp(−βE p ) and k dp ∝ exp(−βE dp ), where E p (E dp ) is the activation energy for phosphorylation (dephosphorylatiion) and β is the inverse temperature by taking the Boltzmann constant as unity. For case of the nutrient compensation (Fig. 5), we considered the model where only the phosphorylation reaction speed, k p , is proportional to ATPto-ADP ratio and others are independent of it.

Transcription-translation-based oscillator (TTO) model
We introduce a model of a circadian clock of a fruit fly [22] (Fig.1B). The governing differential equations are C 0 C 1 C 2 C 3 C 4 C 5 C 6 C 0 C 1 C 2 C 3 C 4 C 5 C 6 AC 0 AC 1 AC 2 AC 3 AC 4  described below: where M is the mRNA of the clock-related gene (per mRNA); R is the protein precursor of a clock-related protein (PER protein); Q and P are an extranuclear protein and a nucleic protein, respectively; and [x] denotes the concentration of x. k, a, s, b, c, d, u, and v are rate constants, and h, a ′ , s ′ , b ′ , c ′ , d ′ , u ′ , and v ′ are dissociation constants.
In [22], it was reported that a and k are especially important to determine the length of the period. Following this report, we set the rate constant of each reaction to follow the Arrhenius equation, i.e., a kinetic constant of mRNA degradation as a ∝ exp(−βE a ) and that of transcription as k ∝ exp(−βE k ) where E a and E k are activation energies of mRNA degradation and transcription, respectively. with E dp = 0.1 (cyan line), 0.4 (green line), and 0.6 (orange line), plotted against the time, normalized by the period of temperature cycles. As shown with the red line, the temperature (to be precise, its inverse β) is periodically changed between β1 = 1.0 and β2 = 1.1 with the period at β1, where each interval at β1 and β2 is set identical. The phosphorylation oscillation is entrained with the external temperature cycle, at the time shaded in the figure. (B) Plots of phosphorylation levels per period of the external cycle, i.e., at the time when the temperature is switched from β1 to β2. Initially, the phosphorylation level changes per period, and then, after entrainment, phosphorylation level takes an almost constant, when the clock is entrained with the external cycle. For smaller E dp , the phosphorylation clock is entrained faster.

CALCULATION OF THE ENTRAINABILITY
To calculate the entrainability in Fig.1B in the main text, initially, 36 oscillators are set at same intervals of the phase, and the number of cycles needed for the oscillators to synchronize is computed. As for the numerical criteria, and we regard that the clock is entrained if the circular variance of phase from different initial conditions is smaller than 10 −9 . Here, the temperature cycle is applied as a square wave between β 1 = 1.0 and β 2 = 1.1 with an equal interval, with the period at the condition of β = β 1 . When entrainability is zero, the oscillator is never entrained. See also Fig.2.

ANALYSIS OF STUART-LANDAU EQUATION
We introduce the Stuart-Landau equation with an external parameter β: This form is derived from the complex Ginzburg-Landau equation, where f 1 (β) represents the change in the bifurcation parameter, and f 2 (β) in that for phase-amplitude coupling.
where, the first order term (f 1 (β)) and the third order term (f 2 (β)) should have different dependency upon β for the adaptation mechanism to work. Indeed the reciprocity holds generally for other forms of the third order term to satisfy adaptation. Considering the stability of limit cycle, the relaxation of R after perturbation is assumed to be much faster than that of Θ. Hence, R * , the steady-state value of R (> 0), is obtained as Here, the period is given as: Here, the argument Θ is defined only on a limit-cycle orbit, and we introduce the phase φ to extend the definition to the phase space out of the limit-cycle attractor, in particular to its basin. It is postulated that φ agrees with Here, we set initial ATP ratio x as 0.5. ∆T /T and ∆φ are calculated similarly to how they are calculated in Fig.1 in the main text.
Θ on the limit-cycle orbit, i.e., different orbits from the same φ converge to the same point on the limit cycle having the same Θ value. Now, we will derive an isochrone, which is a set of points with the same φ on the phase space. φ is expected to have rotational symmetry and is given as Moreover, the time evolution of φ should coincide with that of Θ. Thus, φ is derived as: Hence, Thus, Accordingly, we obtain Next, we will derive C. Because the phase φ(R, Θ, β) should coincide with Θ at R = R * , g(R * ) becomes zero, and then from Eq.(4) in the main text and Eq. (12), Thus, the isochrone of the Stuart-Landau equation against the parameter β is derived. Then, we consider an operation that increases β from β 0 to β 0 + ∆β and instantaneously reverses it to β 0 . At this time, by assuming that R instantaneously relaxes to R * (β 0 + ∆β) = (f 1 (β 0 + ∆β)) 1/2 while Θ remains unchanged, the phase after the above operation is derived as In contrast, the phase before the operation is given as Hence, when ∆β ≪ β, the change in phase is derived as: Therefore, from Eq.(6) in the main text and Eq. (16), changes in the period and phase are represented by an equality.

VAN DER POL OSCILLATOR MODEL WITH PARAMETER β
The van der Pol oscillator is one of the simplest nonlinear-oscillator models, and it is given by: The above equation can be decomposed into two ordinary differential equations as follows: where ǫ is the strength of nonlinearity. As ǫ increases, the system deviates more from a harmonic oscillator.
Here, we modify the van der Pol oscillator to show the change in the amplitude against a change in an external parameter as in the Stuart-Landau equation. We alter van der Pol oscillator as dx dt = y, where β is an environmental parameter. If ǫ is sufficiently small, the amplitude A can be derived by using perturbation calculation as where |A| * is a fixed point value of |A|. On the other hand, we introduce the dependence of the velocity on the environmental parameter as When ǫ is small, the above modified van der Pol oscillator is expected to demonstrate same behavior as the Stuart-Landau model, as described in the main text. When ǫ is large, however, the nonlinearity becomes large and the oscillatory behavior is altered from sinusoidal to relaxation.
When the intensity of nonlinearity, ǫ, is small, the magnitude of change in the period and magnitude of change in the phase are fitted well by a linear relationship a∆T /T + b∆φ = c (a, b, and c are constants) (see Fig.9A). Here, as the intensity of nonlinearity increases, the orbit of a limit cycle is deformed from a circle, and the dynamics shifts from sinusoidal oscillation to relaxation oscillation. Still, the reciprocity is valid as long as the magnitude of stimuli is not exceedingly large (see Fig.9B). Thus, the reciprocity is a universal feature beyond the neighborhood of Hopf bifurcation.