Waveform distortion for temperature compensation and synchronization in circadian rhythms: An approach based on the renormalization group method
Shingo Gibo Teiji Kunihiro Tetsuo Hatsuda Gen Kurosawa

Shedding Light on the Future: Exploring Quantum Neural Networks through Optics
    Dagomir Kaszlikowski
    September 9, 2024
==============================================================================

§ ABSTRACT 

Numerous biological processes accelerate as temperatures increase, but the period of circadian rhythms remains constant, known as temperature compensation, while synchronizing with the 24h light-dark cycle. 
We

theoretically explores 
the possible relevance of waveform distortions in circadian gene-protein dynamics to

the temperature compensation and synchronization.

Our analysis of the Goodwin model provides a

coherent explanation 
of most of temperature compensation hypotheses. Using the renormalization group 

method, we analytically
demonstrate that the decreasing phase of circadian protein oscillations should lengthen with increasing temperature, 

leading to waveform distortions to maintain a stable period. This waveform-period correlation also occurs in other oscillators like Lotka-Volterra and van der Pol models. A reanalysis of known data nicely confirms
our findings on waveform distortion and its impact on synchronization range. 

Thus we conclude that circadian rhythm waveforms are fundamental to both temperature compensation and synchronization. 



Keywords: waveform distortion, renormalization group method, circadian rhythms, temperature compensation, synchronization






§ AUTHOR SUMMARY

Our daily rhythms are underlain by gene regulatory and biochemical networks, called circadian clocks. Although most biochemical reactions accelerate as temperature increases, the period of circadian rhythms is almost constant even with increasing temperature. This phenomenon is called temperature compensation, and the mechanism is still unclear. By applying a method of theoretical physics, the renormalization group method to a biological problem for the first time, we revealed that the waveform of gene dynamics should be more distorted from sinusoidal wave at higher temperature when the circadian period is stable to changes in temperature. This prediction as for the importance of waveform in temperature compensation is verified by analyzing published experimental data. Notably, the correlation between period and waveform distortion holds for other oscillator models, indicating the waveform distortion is important for determining the period in various types of oscillatory systems. Another unsolved problem of circadian clocks is to synchronize with environmental light-dark cycles. By theoretically analyzing a circadian clock model, we found that the frequency range for synchronization becomes narrower when the waveform is distorted.




§ INTRODUCTION

Humans exhibit sleep-wake cycles with an approximate 24h period, and these cycles persist 
under constant environmental conditions, a phenomenon termed the circadian rhythm. This temporal regulation 
exists in both humans and various organisms such as molds, plants, and insects <cit.>. Recent advances in genetic research 
through  insects, molds, mammals, and plants have unveiled 
that genes and proteins 



are involved as integral components in the primary mechanism governing autonomous circadian rhythms <cit.>. 

Understanding circadian rhythms holds promise for deciphering a multitude of sleep patterns, including sleep disorders such as advanced sleep phase syndrome (characterized by early awakening around 4:00 am), delayed sleep phase syndrome (marked by late awakening), non-24h sleep-wake disorder, and narcolepsy <cit.>. Notably, advanced and delayed sleep phase syndromes are believed to be linked to the circadian rhythm period <cit.>. Ongoing studies explore possible correlations between genetic characteristics revealed by large-scale genetic analysis and various sleep patterns <cit.>. However, the nature of the system is so intricate that it remains a challenge to link sleep patterns to specific genes. In such a situation, it would be meaningful to have recourse to mathematical models and obtain possible hints for the linkage and hopefully suggestions for studies of genetic dynamics.

One unresolved fundamental issue in circadian rhythm research is temperature
compensation <cit.> 
, in which the period 

keeps constant despite temperature-induced changes in reaction rates. 
Despite the extensive experimental and theoretical research on 
temperature compensation, the mechanism has remained elusive. 
Hypotheses have been proposed to explain temperature compensation, 
including the balance hypothesis, critical-reaction hypothesis, 
temperature-amplitude coupling hypothesis, and waveform hypothesis. 
The balance hypothesis proposes that the stability of the circadian period 
with temperature arises from a balance between period-lengthening and 
period-shortening reactions <cit.>. 
The critical-reaction hypothesis assumes that there should be critical 
reactions that determine the circadian period. If these reaction 
rates are stable against temperature variations, 
then the circadian period will similarly remain stable
<cit.>. The temperature-amplitude 
coupling hypothesis suggests that temperature-sensitive amplitudes 
in gene activity rhythms should generate a stable period by generating 
larger amplitudes at higher temperatures <cit.>. 
Lastly, the waveform hypothesis proposes that temperature-sensitive 
waveforms in gene activity rhythms should be correlated with a stable 
period in a manner that their higher harmonic components become 
larger and the distortion of the waveform increases at higher
temperatures <cit.>.  

Another unresolved issue in circadian rhythm research is synchronization 
with 24h environmental light-dark cycles.
Previous theoretical and experimental studies on synchronization revealed 
that if the internal period of the oscillation closely matches the 
external period, then it is more likely to synchronize with the forced 
period <cit.>. 
Additionally, experimental studies on several species uncovered genes and 
proteins in circadian systems affected by a light pulse <cit.>. 
In reality, the circadian rhythm must adjust to the 24h light-dark cycle while 
maintaining a temperature-compensated circadian period. 
Therefore, multiple questions arise. (i) Given the significant temperature 
variations between seasons, how do organisms synchronize their circadian 
rhythms with the 24h light-dark cycle across various temperatures 
<cit.>? (ii) if the gene activity 
rhythm of the circadian rhythms becomes more distorted as temperatures increase 
to achieve temperature compensation, how does the ease of synchronization 
change with temperature variations?
Theoretical analyses incorporating the findings of light pulse experiments
might provide further insights into these questions.
 
In the present paper, we investigate possible roles of the waveform distortion 
in temperature compensation based on analytical and numerical analyses 
of the Goodwin model for circadian rhythms and clarify how the 
waveform in gene activity rhythms tends to be more distorted at higher 
temperatures (e.g., steeper rise, longer tail) for temperature compensation.

To this end, we employ the renormalization group (RG) method, a powerful 
tool for analyzing various non-linear systems described by ordinary and partial 
differential equations, to derive global solutions that are valid in a global 
time domain <cit.>. 

Combining an index for waveform distortion, namely non-sinusoidal 
power (NS) introduced by two of the present authors (KG) <cit.>, with the result of 
the RG method,
we can obtain both a unified picture of 

the above mentioned theoretical hypotheses 
(balance hypothesis, critical-reaction 
hypothesis, temperature-amplitude coupling hypothesis, and waveform hypothesis) 
and quantify previous experimental data on Drosophila <cit.>. 
 Our analyses demonstrate that the fundamental role of the waveform distortions 
 in temperature compensation from both theoretical and experimental 
 perspectives in accordance with 
 
 the previous findings <cit.>.
Moreover, we reveal for the first time the mechanism by which 
the synchronization of circadian rhythms changes with temperature if 
the waveform in gene activity rhythms is more distorted at higher temperatures. 
We theoretically prove that the frequency range of the external force 
that synchronizes circadian rhythms becomes narrower if the waveform 
of gene activity rhythms is more distorted. This indicates that it is more 
difficult to synchronize with light-dark cycles at higher temperatures. The 
present result of synchronization is consistent with the previous experimental 
and numerical studies demonstrating that the magnitude of the phase shift 

caused by light pulses was smaller at higher temperature <cit.>.



§ RESULTS



§ WAVEFORM DISTORTION IN CIRCADIAN RHYTHMS



 §.§ Index for waveform distortion


Let the time dependence of a certain variable in the circadian rhythm system be 
expressed in a Fourier series as 
x(t)=∑_j=-∞^∞a_jexp(i(2π/τ)jt)), 
with a_j being the Fourier coefficients of

the oscillatory time series. 
Then, we introduce an index for describing   
the distortion of x(t) from a sinusoidal shape as

    NS=[ ∑_j=1^∞|a_j|^2j^m/∑_j=1^∞|a_j|^2j^q]^1/2 (m>q≥ 0),

where m and q are integers. 
Termed the "non-sinusoidal power (NS)", 
this index is designed to emphasize higher harmonics (m>q) 
as discussed in a previous paper <cit.>. 
For instance, we have NS=1 when the time series has a sinusoidal waveform 
in which only the coefficients for the fundamental component are non-zero 
(a_± 1≠0). Conversely, for non-sinusoidal time series, 
the coefficients for higher harmonics are non-vanishing, resulting in NS>1. 
The previous theoretical work demonstrated that a more distorted waveform 
(larger NS) at higher temperature is necessary for temperature compensation 
in the four-variable negative-feedback model <cit.>. 
However, it is unclear whether the relevance of waveforms 
to temperature compensation 
found in previous research has general validity not 
restricted to some specific model.

To explore the possible relevance of the waveform characteristics to temperature 
compensation and synchronization, we consider

the simplest model for circadian rhythms, known as the Goodwin model (Fig. <ref>A)
<cit.>. This model incorporates negative-feedback regulation
of gene expression, a mechanism established as essential for 
transcriptional-translational oscillations. 
The three-component Goodwin model reads:

    dx_1/dt=f(x_3)-k_1x_1,
    
       dx_2/dt=p_1x_1-k_2x_2,
    
       dx_3/dt=p_2x_2-k_3x_3,

where x_1(t) represents mRNA abundance, and x_2(t) 
and x_3(t) denote protein abundance. 
The function f(x_3) in the model signifies transcriptional regulation, 
and the parameters p_1 and p_2 denote protein synthesis and 
phosphorylation rates, respectively,
and k_i (i=1, 2, 3) represent degradation rates (Fig. <ref>A). 
By applying signal processing methods, Forger derived the period of this model 
as follows <cit.>: 

    τ=2π/√(k_1k_2+k_2k_3+k_3k_1)[ ∑_j=1^∞|a_j|^2j^4/∑_j=1^∞|a_j|^2j^2]^1/2.

Subsequently, two of the present authors indicated that this formula implies 
that temperature compensation of the period in this model occurs 
only when the waveform (NS) is distorted as temperature 
increases <cit.>.
Suppose that all reactions become faster as temperature increases in the model. 
Then, one can numerically demonstrate that the waveform tends to be more 
distorted (larger NS) at higher temperatures for temperature compensation 
(Fig. <ref>B, magenta line). 
We call this mechanism the waveform hypothesis for temperature compensation. 





 §.§ Theories of temperature compensation

The waveform hypothesis and the other three hypotheses for temperature 
compensation (balance hypothesis I, critical-reaction hypothesis II, and 
amplitude hypothesis III) are not mutually exclusive, which can be understood in a 
unified way through Eq. (<ref>):

 
  * I.
The balance hypothesis, previously explored theoretically by Ruoff <cit.>, 
suggests that the temperature compensation of the period is caused by a balance 
between the effects of reactions that shorten the period and those that lengthen 
the period. Equation (<ref>) illustrates that the balance between the effect 
of shortening the period and that of lengthening the period can be caused by 
changing the distortion of the waveform. 

  * II.
Equation (<ref>) demonstrates that even if some of the governing 
reactions in circadian rhythms
are temperature-insensitive <cit.>
 as discussed in the critical-reaction hypothesis. 
 Still, the temperature compensation  requires waveform distortion at high 
 temperatures 
 if reactions other than some of the governing reactions accelerate at higher 
 temperatures.
 
  * III.
Our numerical analysis of the circadian model (Fig. <ref>B) 

show that when temperature compensation occurs, the waveform is more 
distorted at higher temperatures, 
and the amplitude of the oscillation is larger. This tendency for the amplitude to 
increase at high temperatures 

is consistent with the amplitude hypothesis. According to Eq. 
(<ref>), a greater distortion of the 
waveform at higher temperatures is necessary, but not sufficient, for 
temperature compensation. 




 §.§ Synchronization in circadian rhythms

If the waveform is more distorted at higher temperatures for temperature 
compensation, then it would be intriguing to explore whether the temperature-
dependent waveform also affects the synchronization of circadian rhythms with 
environmental light-dark cycles at various temperatures. 
Theoretical and experimental studies of synchronization 
and circadian rhythms illustrated that the oscillation is more likely to 
synchronize with forcing cycles if the internal period is sufficiently 
close to the period of the forcing cycles <cit.>. In mammals and Neurospora, a light pulse is 
known to increase Per1 and frq mRNA expression <cit.>. To incorporate gene activation during the light phase, 
we employ

a model in which
Eq. (<ref>) for the change in mRNA expression is modified to

    dx_1/dt=f(x_3)-k_1x_1+Icos(Ω t),

where I represents the light intensity and Ω is 
the angular frequency of the light-dark cycles. 



 §.§ Main results

In Sections 3 and 4, we provide detailed discussions on how the waveform 
in gene activity rhythms should be distorted at higher temperatures using the 
RG method, as well as the synchronization in circadian rhythms.
The main results are pictorially summarized in Fig. <ref>B. 
The magenta line indicates that the waveform of the gene activity rhythms 
should be more distorted 
at higher temperatures for temperature compensation, whereas the cyan line 
indicates that synchronization 
with the light-dark cycles should become more difficult at higher temperatures 
because of the larger waveform distortion.



§ WAVEFORM DISTORTION AND TEMPERATURE COMPENSATION



 §.§ Numerical simulation of the waveform-period correlation

Equation (<ref>) and numerical simulations indicate that NS tends 
to be larger when the period is relatively stable even with 

increased parameter values (see Fig. <ref>B). 
To quantitatively reveal the correlation between waveform and period, 
we conduct

numerical simulations using a circadian clock model. 
In the analysis of the circadian clock model, the transcription function 
f(x_3)=r/x_3^n was considered for simplicity. 
We first search

for parameter sets in which oscillations occur.

We define 
those parameter sets as the reference parameter sets. Because 
many biochemical parameters have not yet been measured, we prepared

100 random reference parameter sets 
for the oscillations. k_1, k_2, k_3, p_1, p_2, and r 


were assigned uniformly distributed random values ranging from 0 to 10, 
and n 


was assigned a uniformly distributed random 
integer ranging from 9 to 15. The period obtained with each reference 
parameter set 


was denoted as τ_1.
Next, reaction rates often follow the Arrhenius equation, which states that a 10°C 

rise in temperature increases 
the reaction rate by a factor of 2-3. To incorporate the effect of high 
temperature, instead of 
using the Arrhenius equation, each parameter in the model's reference parameter 
set 


was randomly multiplied by a factor of 1.1-1.9, and the period and waveform 


were examined when the oscillation behavior persists.

The period obtained by increasing the parameters from each reference parameter 
set 


was denoted as τ_2, and the ratio, τ_2/τ_1, 

was called the relative period. 

To quantitatively analyze the correlation between period and waveform 
when temperature compensation occurs, 
we consider

the case of the relative period ≥ 0.85 because it has been experimentally 
confirmed that the circadian rhythm frequency at high temperature divided by 

that at low temperature of the wild-type ranges 
between 0.85 and 1.15 when the temperature
is increased by 10°C and temperature compensation occurs <cit.>. In the 
present numerical analysis, the period 

is relatively stable (relative period ≥ 0.85) 
in 34 of the 4900 parameter sets,
in qualitative agreement with previous theoretical analyses 
that the period often shortens with increasing 
reaction rates <cit.>. Because the range of parameter 
variation 

is 1.1-1.9 and the average value 

is 1.5, the reaction rate 

is accelerated by a factor of 1.5 on average, and the average relative period 

is approximately 1/1.5 ≈ 0.67. Figure <ref> indicates 
that when temperature compensation occurs, there is a clear correlation between the period and waveform.






 §.§ RG analysis of waveform-period correlations

In the previous section, we demonstrated that the temperature compensation 
of the period 
in the Goodwin model is always accompanied by an increase of the index
NS and waveform distortion correspondingly occurs 
as temperatures increase. 
This raises the following question: 
Is there a universal law governing the waveform distortion occurring when 
temperature increases? To 

answer this question, we derive

an approximate solution for the time evolution of the Goodwin model 
for the circadian rhythm using 

a powerful reduction method, 

called "the renormalization-group(RG) method " <cit.>. 
The solution obtained using the RG method can be interpreted 
as the envelope of the set of solutions given in the perturbation theory, 
which has been applied to various models, including (but not limited to) ODE, PDE, 
discrete systems, and stochastic equations 
<cit.>. 
To apply the RG method, we again set the transcriptional regulation 
function f(x_3) to be r/x_3^n. 
In this function, n is the cooperativity of the transcriptional 
regulation, which is a Hopf bifurcation parameter.
The approximate solution of the phosphorylated protein of the circadian clock 
reads (see 

Supplementary Information A.2)

    x_3(t)=( p_1p_2r/s_3)^s_3/s_1s_2+ε A_0sin(ω t)
    +ε^2 A_1A_0^2sin(2ω t+α)+o(ε^2)

where we have

    ε =n-s_4/s_3,
    
       s_1=k_1+k_2+k_3, 
     s_2=k_1k_2+k_2k_3+k_3k_1, 
    
       s_3=k_1k_2k_3, 
      s_4=(k_1+k_2)(k_2+k_3)(k_3+k_1),

with the angular velocity and the phase parameter of the second-order term

    ω = √(s_2)-εs_1s_3s_4/6(2s_1s_2^2-
    (s_1^2+6s_2)s_3)√(s_2)+o(ε^2),
    
       α = arctan( s_1/2√(s_2)),

as well as the amplitudes

    A_0=√(4((s_1^2+s_2)^2-2ε s_1s_3)
    (s_1^2+4s_2)s_3^3/ε (2s_1s_2^2-(s_1^2+6s_2)s_3)
    (s_1^2+s_2)^2s_4s_1s_2)( p_1p_2r/s_3) ^s_3/s_1s_2,
    
       A_1=s_4s_1/12s_3√(s_1^2+4s_2)( 
    s_3/p_1p_2r) ^s_3/s_1s_2.


The RG method provides an approximate but globally valid solution, 

and thus enable us to make a 
detailed investigation of 
the waveform distortion when temperature compensation occurs 
in the Goodwin model. 
The numerical analysis using the same parameter sets as used in Fig. <ref>A, 
in which the relative period in the model remains stable
and within the interval (0.85, 1.0) against the temperature variations, 

shows that the phase parameter α in the 2nd-order frequency

tends to decrease with increasing reaction rates
(see Fig. <ref>A). 
When the increase in reaction rates is small, the change 
in the phase of the second-order frequency scatters around zero 
and is negligible.
However, with a significant increase in the reaction rates, 
the phase of the second order always tends 
to decrease as the reaction rates increase.



The significance of the phase parameter α given in the second-order term 
on the waveform of the time series can be understood intuitively as follows: 
when the phase α is large, the increasing duration tends to become 
longer because of
the less overlap of the time profiles given by sin(ω t) and 
sin(2ω t+α) (Fig. <ref>BC, blue line). 
Conversely, a smaller α tends to result 
in a shorter increasing duration because of an additive effect of the two terms, 
which leads to a steeper slope on total, as presented 
in Fig. <ref>BC, red line. 
Therefore, the numerical results in Fig. <ref>A suggest 
that the decreasing duration 
of the time series elongates with as the reaction rate 
increases when the period is relatively stable despite 
the increasing reaction rate.


For a theoretical confirmation of  the numerical result that the 
phase parameter α given in the second-order term
tends to become smaller as the temperature increases 
when the period is temperature-compensated, 
we analyze

the sensitivity of the angular frequency ω and the phase α 
to the reaction rates by utilizing

the results of the RG method.

With use of Eqs. (<ref>) and (<ref>), we calculate

dω =(∂ω/∂ k_1)dk_1+(∂ω/∂ 
k_2)dk_2+(∂ω/∂ k_3)dk_3 and dα =
(∂α/∂ k_1)dk_1+(∂α/∂ k_2)dk_2+
(∂α/∂ k_3)dk_3. 
In Fig. <ref>A, we present the parameter regions given 
by the constraints dω =0 (red surface) and 
dα =0 (yellow surface) for 
the cooperativity n=12 of the transcription regulation. 
We can see that the region for dα<0 (outside yellow surface) 
includes that given 
by the constraint dω =0 for all k_i (i=1, 2, 3).
This implies that if the period is robust against a change 

in the parameters 
k_i (i=1, 2, 3), then the phase α in the second-order term 
always becomes smaller with increasing parameters.



Next, let us examine how the parameter regions given by the constraints dω =0 and dα =0 
change with variations of the cooperativity n. 
The numerical calculation 

shows that the region corresponding to dα <0 
includes that given 
by dω =0 for n=13 and 14. 
However, in the case of an exceedingly high cooperativity of 
transcription regulation n, which is not biologically realistic, 
dα can be positive when dω =0. For instance, for n=20, 
although dα is negative for most of the parameter space, 
there is a region in which dα >0 when dω =0 
(see Fig. <ref>B).
These results indicate that if circadian rhythms are stable 

under temperature variations, 
the slope in the increasing phase of phosphoproteins should become sharper as 
temperature increases.
Thus, we conclude that 
(i) the waveform of the gene activity rhythm should be more distorted 
at higher temperatures, and (ii) the rate of the increase 
in phosphoprotein levels 
should be greater at higher temperatures if temperature compensation is achieved. 
In principle, these features can be tested experimentally.



 §.§ Verification of the theoretical analysis of temperature 
compensation using published experimental data

The period formula Eq. (<ref>) of the Goodwin model
indicates that the non-sinusoidal index NS of the waveform 
of the circadian rhythms becomes larger, 
implying greater distortion of the waveform when all reactions are faster at 
higher temperatures during temperature compensation. 
To test this theoretical prediction of circadian gene activity in actual 
organisms, we analyze

the waveform of the activity rhythms of the timeless gene in Drosophila at 18 and 
29°C using published experimental data <cit.>.

First, we extract the time series of the average curve 
from Fig. 3C in a prior study <cit.> 
using WebPlotDigitizer per hour. 
Second, we add uniformly distributed noise between 
-0.4 and 0.4 
to the extracted 
data to consider data errors. Then, we interpolate
the time series every 0.1 h 
using spline interpolation (Fig. <ref>A).
The interpolated data were detrended by multiplying 
an exponential 
function so that the position of the local minima of the oscillations 
are approximately reproduced. 

Then, the detrended time series is fitted with 
a sum of trigonometric functions up to the third harmonics using 
the generalized harmonic analysis (GHA) method <cit.>. 
The width of the window for the analysis is set to one period.

Using the Fourier coefficients of the fitting time series, 
we 

evaluate the distribution and average value of NS. 
The 

resultant NS of the activity rhythms 
of the timeless gene, as defined by Eq. (<ref>), 
at a higher temperature (29 ^∘C) tend to be larger than 
that at a lower temperature (18 ^∘C), whereas

the NS values are somewhat varied (Fig. <ref>B), 
which 

is consistent with 

the prediction.
Experimental studies have demonstrated that temperature compensation 
can be impaired by genetic mutations. In the Drosophila mutant perL, 
the period increases with increasing temperature <cit.>. 
Equation (<ref>) implies that if the period increases with temperature, 
then the waveform of the circadian rhythm should be non-sinusoidal and more 
distorted at higher temperatures. 

Thus, it is predicted that the waveform 
of the circadian gene activity in perL should become more non-sinusoidal 
with higher temperatures. Again, we can quantify the waveform of perL using 
experimental data <cit.> (Supplementary Fig. <ref>). The waveform 
of circadian gene activity tends to be more non-sinusoidal 
at higher temperatures in perL, as observed in the wild-type, 
whereas the NS values varied, in line with the prediction.

	



§ THEORETICAL ANALYSIS OF SYNCHRONIZATION IN THE CIRCADIAN RHYTHM MODEL

The numerical 

reslut in Fig. <ref>B 

shows that the range of synchronization into the light-dark cycles tends 
to decrease as the waveform becomes more distorted 
in the simple circadian clock model.
To 

clarify the condition for synchronization, we again 
consider

the Goodwin model but with an external 
force incorporated as follows:


    dx_1/dt=f(x_3)-k_1x_1+Icos(Ω t)
    
       dx_2/dt=p_1x_1-k_2x_2
       dx_3/dt=p_2x_2-k_3x_3

where Icos(Ω t) is a periodic environmental change, 
such as a light-dark cycle. By eliminating x_1 and x_2, 
Eqs. (<ref>)-(<ref>) is converted to the following single equation:

    d^3x_3/dt^3+s_1d^2x_3/dt^2
    +s_2dx_3/dt+s_3x_3
    =p_1p_2f(x)+p_1p_2Icos(Ω t).

If the model is to admit a synchronization 
to the external cycle Icos(Ω t) at all, 
then x_3(t) should be written as the Fourier series 
x_3(t)=∑_j=-∞^∞a_jexp(iΩ jt).
Multiplying Eq. (<ref>) by dx/dt and integrating 
that for the interval t to t+2π/Ω, we have the following equation:

    Ω^3-ω^2Ω = 
    1/2p_1p_2I Rsinβ

where 

    ω=2π/τ=√(s_2∑_j=1^∞|a_j|^2j^2/∑_j=1^∞|a_j|^2j^4)

is the natural angular frequency without the external force and 

    R=
    |a_1|/∑_j=1^∞|a_j|^2j^4

 with β being the argument of a_1 such that a_1=|a_1|exp(iβ). 
Because -1≤sinβ≤ 1, when x(t) synchronizes with the external 
cycles, the angular frequency of the external cycles (Ω) should satisfy
the inequality


    |Ω^3-ω^2Ω |
    ≤1/2p_1p_2I R.

We note that
R defined by Eq. (<ref>) becomes smaller when the components 
of higher harmonics become larger and the waveform exhibits greater
distortion. Therefore, if the waveform is more distorted 
by, say, an increasing temperature, the bounds of Eq. (<ref>) become smaller,
and accordingly, the allowed region of the middle term is narrower. 
The left hand side of Eq. (<ref>) is a cubic function of Ω, 
which is monotonically increasing near Ω=ω. 
If the waveform is more distorted, 
 Eq.(<ref>)
should be smaller, making the allowed region of Eq. (<ref>) narrower. 
Then, the range of Ω that causes synchronization becomes narrower, 
as presented in Fig. <ref>. 
This indicates that the range of synchronization into light-dark cycles always
decreases as the waveform becomes more distorted in the simple circadian model, 
which is consistent with the numerical simulation in Fig. <ref>B.





§ DISCUSSION



We theoretically explored the conditions for clarifying the temperature 
compensation of the biological clock and its synchronization to light-dark cycles 
with a particular focus on waveform distortion. The theoretical analysis 
of the Goodwin model, one of the most studied models of biological 
clocks, revealed that waveform distortion of gene activity rhythms with increasing 
temperature is necessary for temperature compensation. Furthermore, we derived an 
approximate but globally valid solution to the waveform of the time profiles using 
the RG method as a powerful tool for global analysis. This allowed us to 
investigate, based on the analytical solution, whether there is a universal law 
for the mechanism by which the waveform changes with temperature 
variation.

The results indicated that temperature compensation is more likely to 
occur if the waveform is distorted if the decreasing duration of 
circadian protein oscillation elongates as temperature increases. 
Although theoretical predictions based on a model might not 
always be realized in real organisms, we quantified the gene activity rhythms of 
published experimental data using Drosophila. This quantification confirmed that 
the waveform is distorted at high temperatures, in accordance with our 
theoretical predictions.

It is notable that the systematic wave distortion governed by Eq. (<ref>), 
which we have found to hold in the Goodwin model, also applies to a wide 
class of non-linear oscillators used for biological phenomena different from 
biological rhythms, including the Lotka-Volterra model <cit.>, 
which is commonly used in ecology, and the van der Pol model, as 
presented in the Supplementary Infomration A.3 and A.4 (see also <cit.>). This suggests that 
exploring the possible significance of waveform distortion in other mathematical 
models, such as the Fitz-Hugh-Nagumo model in neuroscience, would be intriguing 
<cit.>.

To the best of our knowledge, this is the first study to apply the RG method, 
a powerful resummation method of the perturbation series first developed in physics, 
to circadian rhythm problems. In the RG method, secular terms appearing in the 
naïve perturbation series are renormalized into the 'integral constants' 

, which thus acquire the nature of the slow modes, 
making it a powerful tool for global and asymptotic analysis. 
Unlike the naive perturbation theory, the solutions given by the RG method provide 
a time evolution close to numerical simulations in 
the relevant global domain of time, 
offering an approximate solution for the period and waveform. 
The analytical results predict that longer tails of gene activity rhythms 
at higher temperatures occur for temperature compensation.

This study

also investigated the synchronization 
with environmental light-dark cycles at various temperatures <cit.>. The numerical simulations and theoretical analysis predict 
that as the distortion 
of the gene activity rhythms for achieving a temperature-compensated period 
increases, it becomes more difficult to synchronize with 
the light-dark cycle. 
This prediction aligns with the reported temperature-dependent variation 
in response 
to light pulses in Drosophila and Neurospora, 
displaying smaller phase shifts 
at higher temperatures <cit.> <cit.>.

As mentioned in the Introduction, previous theoretical and experimental 
studies, such as those using Drosophila <cit.>, suggested 
that waveforms in gene activity rhythms do not change under temperature variations, 
although they are blurred by error bands. These findings are 
apparently inconsistent with our current conclusion. We believe that this 
discrepancy stems from two main factors: differing assumptions about the 
temperature sensitivity of degradation rates and differing interpretations of 
experimental results regarding gene activity rhythms at distinct temperatures. 
First, the previous study assumed that all degradation rates are 
temperature-insensitive. Thus, the waveform of gene activity rhythms does not need 
to change with temperature. By contrast, we assume that some degradation 
rates at least should accelerate with temperature, leading to a more distorted 
waveform of gene activity rhythms at higher temperatures. Second, the previous 
study <cit.> interpreted their experimental results as indicating that gene 
activity rhythms at different temperatures can be collapsed onto each other by 
rescaling, supporting their prediction that temperature compensation occurs because 
of rescaling and the stability of the temperature insensitivity of degradation 
rates. Conversely, our quantification of their experimental data 
indicated that the waveform tends to be more distorted at higher 
temperatures, whereas variation in NS values was noted. Thus, the 

present result is consistent 

with our theoretical 
prediction. We believe that further systematic quantification of the waveforms of 
gene activity and/or protein activity rhythms in various circadian organisms will 
be essential for clarifying the importance of the waveform in circadian 
rhythms in the future.



§ METHODS



 §.§ Computation of ODEs

The ordinary differential equations in this work were calculated by a fourth-order Runge-Kutta method with MATLAB (The MathWorks, Natick, MA). The time step size Δ t for the Goodwin model (Eqs. (<ref>)-(<ref>)) was 0.001, and that for the Lotka-Volterra model (Eqs. (<ref>)-(<ref>)) and the van der Pol oscillator (Eq. (<ref>)) was 0.0001.



 §.§ Numerical analysis of synchronization with light-dark cycles

The angular frequency of the light-dark cycle Ω in Eq. (<ref>) was changed by 0.0001. We considered that the model was synchronized with the light-dark cycle when the change in the amplitude was smaller than the rounding error as in Fig. 1B.



 §.§ Estimation of waveform distortion from simulated time-series

We estimated waveform distortion, i.e., non-sinusoidal power (NS, see Eq. (<ref>)) from the simulated time-series by using the generalized harmonic analysis (GHA) method as in Figs. 1B and 2A. The GHA method estimates the amplitudes A_j and B_j, and the frequencies f_j, which minimize the squared residual ∫_0^L[x(t)-∑_j=1^j_max{ A_jcos(2π f_jt)+B_jsin(2π f_jt) } ]^2dt. The detailed procedure of GHA can be found in <cit.>.



 §.§ Estimation of phases of 

Fourier components from simulated time-series

By using GHA method, we 

computed the Fourier series x(t)=∑_j=1^j_max{ A_jcos(2π f_jt)+B_jsin(2π f_jt)}. This form was converted into x(t)=∑_j=1^j_maxr_jsin(2π f_jt+α_j), where r_j=√(A_j^2+B_j^2) and α_j =arctan (B/A).



 §.§ Computation of parameter space for temperature compensation using the RG method

The frequency change dω =(∂ω/∂ k_1)dk_1+(∂ω/∂ k_2)dk_2+(∂ω/∂ k_3)dk_3 and the phase change dα =(∂α/∂ k_1)dk_1+(∂α/∂ k_2)dk_2+(∂α/∂ k_3)dk_3 were calculated from Eqs. (<ref>) and (<ref>) by using Maple. The parameter space for dω =0 and dα =0 in Fig. <ref>AB were numerically generated by Maple.



 §.§ Estimation of waveform distortion from published experimental data

The average curves of experimental time-series of timeless luciferase (tim-luc) reporter (shown to recapitulate the dynamics of timeless gene expression) at 18 and 29 ^∘C of Figs. 3C and 4B in the published literature <cit.> were extracted using WebPlotDigitizer at 1 h intervals. Uniformly distributed noise between -0.4 and 0.4 was added to the extracted data, generating 100 time series datasets so that the experimental noise 
is roughly reproduced as in Fig. 5A and Supplementary Fig. 3A. Spline interpolation was applied to set the sampling interval to 0.1 h to smooth time series data. The time series data with noise were detrended by multiplying an exponential function so that the positions of local minima of the oscillations are approximately

reproduced. The Fourier coefficients of the detrended time series were quantified using GHA. The NS values were estimated from the coefficients up to the third harmonics as in Fig. 5B and Supplementary Fig. 3B.



§ ACKNOWLEDGEMENTS

We thank H. Nakao, H. Chiba, Y. Kawahara, A. Mochizuki for useful comments on this study. This work was supported by grants from Japan Science and Technology Agency (JPMJCR1913 to G.K.), and from the Japanese Society for the Promotion of Science, and the Ministry of Education, Culture, Sports, Science, and Technology in Japan (JP21K06105 to G.K., 19K03872 to T.K.).









§ SUPPORTING INFORMATION FOR




§ WAVEFORM DISTORTION FOR TEMPERATURE COMPENSATION AND SYNCHRONIZATION IN CIRCADIAN RHYTHMS: AN APPROACH BASED ON THE RENORMALIZATION GROUP METHOD


Shingo Gibo^1*, Teiji Kunihiro^2, Tetsuo Hatsuda^1, and Gen Kurosawa^1*


^1Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS), RIKEN 

^2Yukawa Institute for Theoretical Physics (YITP), Kyoto University

^*Correspondence: shingo.gibo@riken.jp (S.G), g.kurosawa@riken.jp (G.K)

This PDF includes: 

Supplementary Text 

Supplementary Figures 1 to 3 

Supplementary Tables 1 to 2





§ SUPPLEMENTARY TEXT



 §.§ A.1 Brief introduction of the renormalization group (RG) method 
using a simple model with Hopf bifurcation

In this section, we introduce the RG method <cit.> in a geometrical manner as 
formulated in <cit.> with a simpler prescription 
without the redundant 'time-splitting' procedure.
In this aim, we 

use a generic model with a Hopf bifurcation. 
A more detailed account of
the method is given in <cit.>.

Let us consider the model equation 

    Lx(t)=F(x(t); ε)

where x(t) is a state variable of our dynamics, 
L=∑_n=1^Na_n(d/dt)^n
 is a linear differential operator, F(x(t); ε) is a nonlinear function of x(t), 
 and ε is an internal parameter of F, which acts as a bifurcation parameter of the system. 
We assume that the model has a fixed point x_0 satisfying the equation
F(x_0; 0)=0, which is destabilized for ε >0 through 
the Hopf bifurcation. We are interested in the derivation of the reduced equation 
and an approximate but valid solution in a global domain of time around the critical point 
of the Hopf bifurcation.
Thus,

we apply

the perturbation theory and express the solution around an arbitrary 
time t=t_0 belonging to a global domain in the asymptotic regime (see below) in a power series of ε as follows:

    x(t;t_0)=x_0+ε u_1(t;t_0)+ε ^2u_2(t;t_0)+ε^3u_3(t;t_0)+o(ε^3).

Substituting Eq. (<ref>) into Eq. (<ref>) and equating the terms with the same powers of ε, we obtained

    O(ε^1): L'u_1=0,
    
       O(ε^2): L'u_2=f_1(u_1),
    
       O(ε^3): L'u_3=f_2(u_1,u_2),

where 

    L'=L-(∂ F/∂ x)|_x=x_0,ε =0,

and f_1(u_1) and f_2(u_1,u_2) are nonlinear functions that depend on F(x_0). 
Because Hopf bifurcation occurs at ε =0, two of the eigenvalues of the linear differential operator L' are written as 
± iω_0, with ω_0 being a real number, and the others have negative real parts as
Re(λ_k)<0 (k=1,...,N-2). 

Then, the first-order solution can be expressed as


    u_1(t;t_0)=A(t_0)cos(ω_0 t+θ (t_0))+∑_k=1^N-2c_k(t_0)e^λ_kt,

where A(t_0), θ (t_0), and c_k(t_0) are the 
integral constants which are assumed to

depend on the initial time t_0. 

Next, we consider the asymptotic regime as t → ∞ so that 
the second term describing the transient behavior has virtually become negligible. 

Then,
the first-order solution in this asymptotic regime can be expressed only by the first term as

    u_1(t →∞ ;t_0)=A(t_0)cos(ω_0 t+θ (t_0)).



Next, we proceed to the second-order equation. 
Substituting Eq. (<ref>) into Eq. (<ref>), we have








    L'u_2
    =b_1Acos(ω_0t+θ)
    +b_2A^2cos(2(ω_0t+θ))
    +b_2A^2,

where b_k (k=1, 2) are constants depending on f_1(u_1). 
It is to be noted that the inhomogeneous part (r.h.s.) contains a term proportional to cos(ω_0 t+θ), which is a zero mode of the linear operator that gives rise to secular terms in the particular solutions of the inhomogeneous equation.
The general solution to Eq. (<ref>) is given as a sum of a particular solution to the inhomogeneous equation 
and the general solution to the homogeneous equation.
Now, it is possible and convenient to choose 
the coefficients of the latter so that all of the secular terms vanish
 at t=t_0 <cit.>, which leads to the second-order solution as


    u_2(t;t_0)=   (t-t_0)d_1Acos(ω_0 t+θ)
    +(t-t_0)d_2Asin(ω_0 t+θ) 
       +d_3A^2cos(2(ω_0t+θ))
    +d_4A^2sin(2(ω_0t+θ))
    +d_5A^2,

where d_k (k=1⋯ 5) are constants depending on the right-hand side of Eq. (<ref>). 

Similarly, the third-order solution takes the form of

    u_3(t;t_0)=   (t-t_0)(f_1aA^3+f_1bA)cos(ω_0t+θ)
    +(t-t_0)(f_2aA^3+f_2bA)sin(ω_0t+θ) 
       +(t-t_0)^2f_3Acos(ω_0t+θ)
    +(t-t_0)^2f_4Asin(ω_0t+θ)
    +f_5A^2cos(2(ω_0t+θ)) 
       +f_6A^2sin(2(ω_0t+θ))
    +(t-t_0)f_7A^2cos(2(ω_0t+θ))
    +(t-t_0)f_8A^2sin(2(ω_0t+θ)) 
       +f_9A^3cos(3(ω_0t+θ))+f_10A^3sin(3(ω_0t+θ))
    +f_11A^2+(t-t_0)f_12A^2,

where f_1a, f_1b, f_2a, f_2b, and f_k (k=3⋯ 12) are constants depending on Eqs. (<ref>), (<ref>), and (<ref>).
Note that the solution is constructed so that the secular terms 

vanish at t=t_0. 

Thus, collecting all of the terms, the approximate solution to Eq. (<ref>) up to the third order 
of ε reads

    x(t;t_0)=   x_0+ε Acos(ω_0 t+θ)
    +ε^2{(t-t_0)d_1Acos(ω_0 t+θ)
    +(t-t_0)d_2Asin(ω_0 t+θ) 
       
    +d_3A^2cos(2(ω_0t+θ))
    +d_4A^2sin(2(ω_0t+θ))
    +d_5A^2}
       +ε^3{ (t-t_0)(f_1aA^3+f_1bA)cos(ω_0t+θ)
    +(t-t_0)(f_2aA^3+f_2bA)sin(ω_0t+θ) 
       +(t-t_0)^2f_3Acos(ω_0t+θ)
    +(t-t_0)^2f_4Asin(ω_0t+θ)
    +f_5A^2cos(2(ω_0t+θ)) 
       +f_6A^2sin(2(ω_0t+θ))
    +(t-t_0)f_7A^2cos(2(ω_0t+θ))
    +(t-t_0)f_8A^2sin(2(ω_0t+θ)) 
       +f_9A^3cos(3(ω_0t+θ))+f_10A^3sin(3(ω_0t+θ))
    +f_11A^2+(t-t_0)f_12A^2} +o(ε^3).

Because Eq. (<ref>) contains the secular terms, 


this solution 

is valid only locally around t=t_0,
but 

it exhibits a divergent behavior as | t-t_0| 
goes infinity. 
In fact, this is a rather common behavior occurring in naïve perturbation expansions.

Next, we use a geometrical viewpoint to circumvent the disastrous situation 
following <cit.>.
The solution (<ref>) gives a family of curves with t_0 being the parameter specifying each curve 
in the t-x plane. Each curve gives a good approximate solution to the original equation in a local 
domain around t=t_0. The idea is that the envelope curve of the family of curves
hopefully gives an approximate but valid solution in the global domain including the arbitrary time t_0.
Indeed, this has rigorously been demonstrated to be the case
<cit.>. 
Now, the envelope curve can be constructed using the following envelope equation <cit.>:

    . 
    dx(t;t_0)/dt_0| _t_0=t
    =.
    ∂ x/∂ t_0| _t_0=t
    +.
    dA(t_0)/dt_0∂ x/∂ A| _t_0=t
    +.
    dθ (t_0)/dt_0∂ x/∂θ| _t_0=t 
    =0.

Note that we have 

taken into account the fact that
the integral constants A and θ depend on the `initial time' t=t_0, and (<ref>) actually gives the dynamical equations for these variables. As will be done shortly, the insertion of the solutions to the dynamic equation into (<ref>) gives 
an approximate but globally valid solution to the original equation.
Because the envelope equation (<ref>) takes a similar form 
as the RG equation in quantum field theory,
it is also called the RG equation, and the asymptotic/global
analysis based on this equation was named the RG method <cit.>.

Substituting Eq. (<ref>) into Eq. (<ref>), we have


    0=   ε{dA/dt-ε^2f_1aA^3-ε (d_1+ε f_1b)A
    }cos(ω_0 t+θ) 
       +ε A{
    -dθ/dt-ε^2 f_2aA^2-ε (d_2+ε f_2b)
    }sin(ω_0 t+θ) 
       +ε^2A{
    2(d_3+ε f_5)dA/dt+2(d_4+ε f_6)Adθ/dt-ε f_7A
    }cos(2(ω_0t+θ)) 
       +ε^2A{
    2(d_4+ε f_6)dA/dt-2(d_3+ε f_5)Adθ/dt-ε f_8A
    }sin(2(ω_0t+θ)) 
       +3ε^3A^2{
    f_9dA/dt+f_10Adθ/dt}cos(3(ω_0t+θ))
    +3ε^3A^2{
    f_10dA/dt-f_9Adθ/dt}sin(3(ω_0t+θ)) 
       +2ε^2(d_5+ε f_11)AdA/dt
    -ε^3 f_12A+o(ε^3).


Because dA/dt and dθ/dt are of order ε, the coefficients cos(2(ω_0t+θ)), 
sin(2(ω_0+θ)), 
cos(3(ω_0t+θ)), sin(3(ω_0t+θ)), and AdA/dt are of order ε^3 or 
higher. 
To make Eq. (<ref>) hold for any t, 

we only must ensure that 
the coefficients of the independent functions, namely cos(ω_0t+θ) and sin(ω_0t+θ),
vanish, and hence, we have 



    dA/dt=ε^2f_1aA^3+ε(d_1+ε f_1b)A+o(ε^2),  
    
       dθ/dt=-ε^2f_2aA^2-ε (d_2+ε f_2b)+o(ε^2),

which are the dynamic equations governing the 'integral constants' A and θ. We now see that
the integral constants have been lifted to dynamic variables 
through the RG/envelope equation.
The amplitude equation (<ref>) can be readily solved analytically. 
For instance, when f_1a < 0 and d_1+ε f_1b > 0, 
it yields 

    A(t)=A_0 A/√( A^2+(A_0^2- A^2)e^-2α t),

where 
α=ε(d_1+ε f_1b) and 



    A_0=√(-d_1+ε f_1b/ε f_1a),

with A being the initial amplitude.

Equation (<ref>) indicates that the amplitude approaches A_0 monotonically as 
t → ∞,
implying that A_0 is nothing but the amplitude of the limit cycle admitted 
in the original equation (<ref>). 
Furthermore, Eq. (<ref>) indicates that the angular frequency on the limit cycle reads


    ω=ω_0+(dθ/dt)|_A=A_0
    =
    ω_0-εd_2f_1a-d_1f_2a/f_1a,

which is constant.
The globally valid solution is given as the envelope of the 
family of curves, as previously stated. 
Thus, the solution on the limit cycle, which valid in a global
domain in the asymptotic regime, reads 


    x(t)=   x(t;t_0)|_t_0=t
    
    =   x_0+ε A_0cos(ω t+θ_0)
    +ε^2{ d_3A_0^2cos(2(ω t+θ_0))
    +d_4A_0^2sin(2(ω t+θ_0))+d_5A_0^2}
       +ε^3{ f_5A_0^2cos(2(ω t+θ_0))
    +f_6A_0^2sin(2(ω t+θ_0))
    +f_9A_0^3cos(3(ω t+θ_0)) 
       +f_10A_0^3sin(3(ω t+θ_0))
    +f_11A_0^2} +o(ε^3).


The RG method is a powerful method for obtaining a globally valid solution. 
This method can be applied to various models including discrete, stochastic, and partial differential equations, 

as given in <cit.>.



 §.§ A.2 Derivation of a time-evolution solution in a circadian rhythm model using the RG method


In this subsection, we applied the RG method to derive an approximate but globally valid solution of 
the 

circadian clock model given by Eq. (<ref>)-(<ref>), which is
 a system of first-order equations with three variables.
We first convert the system into a single equation with higher-order derivatives as


    d^3x_3/dt^3
    +s_1d^2x_3/dt^2
    +s_2dx_3/dt
    +s_3x_3=p_1p_2f(x_3)


where s_1=k_1+k_2+k_3, s_2=k_1k_2+k_2k_3+k_3k_1, s_3=k_1k_2k_3. 
To obtain the approximate solution, we set the transcriptional regulator f(x_3) as r/x_3^n. 
This model has a fixed point


    x_0=( s_3/p_1p_2r) ^-1/n+1,

which is destabilized through Hopf bifurcation for

    n> n_0≡s_4/s_3,

with s_4=(k_1+k_2)(k_2+k_3)(k_3+k_1).

Then, we expand the solution around t=t_0 as a series
of ε


    x_3(t;t_0) = x_0+ε u_1(t;t_0)+ε^2 u_2(t;t_0)+ε ^3u_3(t;t_0) +O(ε^3).


Substituting Eq. (<ref>) into Eq. (<ref>) and
equating the coefficients with the same powers of ε, we obtain


    O(ε): d^3u_1/dt^3+s_1d^2u_1/dt^2
    +s_2du_1/dt+s_1s_2u_1=0,
    
       O(ε^2): d^3u_2/dt^3+s_1d^2u_2/dt^2
    +s_2du_2/dt+s_1s_2u_2
    =-s_3u_1+B_1u_1^2,
    
       O(ε^3): d^3u_3/dt^3+s_1d^2u_3/dt^2
    +s_2du_3/dt+s_1s_2u_3
    =-s_3u_2+B_2u_1u_2+B_3u_1^2+B_4u_1^3,

where


    B_1=s_4s_2s_1/2s_3( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       B_2=s_4s_2s_1/s_3( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       B_3={
    s_4+s_3/2-( s_4s_3/2s_2s_1ln( s_3/p_1p_2r)
    )
    }( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       B_4=(s_2s_1+s_3)s_4s_1s_2/6s_3^2( 2s_3/p_1p_2r) ^s_3/s_1s_2.


Because Eq. (<ref>) has three eigenvalues, namely λ_1,2=± i√(s_2) and λ_3=-s_1, the general solution of Eq. (<ref>) reads


    u_1(t;t_0)=A(t_0)cos(ω_0t+θ (t_0))+c(t_0)e^-s_1t,      (ω_0:= √(s_2)),

where A, θ, and c are integral constants that depend on initial time t_0.


Considering the asymptotic regime in which the second term in Eq. (<ref>) is so small and negligible,

then the first-order solution can be written as


    u_1(t;t_0)=A(t_0)cos(ω_0t+θ (t_0)).



Substituting Eq. (<ref>) into Eq. (<ref>), we obtain


    d^3u_2/dt^3+s_1d^2u_2/dt^2
    +s_2du_2/dt+s_1s_2u_2
    =-s_3Acos(ω_0t+θ)+CA^2cos(2ω_0t+2θ)
    +CA^2,


where


    C=s_4s_1s_2/4s_3( s_3/p_1p_2r) ^s_3/s_1s_2.


Equation (<ref>) is an inhomogeneous equation that contains a zero mode of the linear operator, 
and the solution 

having a suitable form for applying the RG method is written as

    u_2(t;t_0)=
       D_1A(t-t_0)cos(ω_0t+θ)
    -D_2A(t-t_0)sin(ω_0t+θ) 
       -D_3A^2cos(2ω_0t+2θ)-D_4A^2sin(2ω_0t+2θ)
    +D_5A^2,


where the coefficients are


    D_1=s_3/2(s_1^2+s_2),
    
       D_2=s_1s_3/2√(s_2)(s_1^2+s_2),
    
       D_3=s_4s_1^2/12(s_1^2+4s_2)s_3( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       D_4=s_4s_1√(s_2)/6(s_1^2+4s_2)s_3( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       D_5=s_4/4s_3( s_3/p_1p_2r) ^s_3/s_1s_2.


Tentatively, after collecting all the obtained terms, we have the approximate solution in the second order as


    x_3(t;t_0)=   x_0+ε Acos(ω_0t+θ)
    +ε^2{ D_1A(t-t_0)cos(ω_0t+θ)
    -D_2A(t-t_0)sin(ω_0t+θ) 
       -D_3A^2cos(2ω_0t+2θ)-D_4A^2sin(2ω_0t+2θ)
    +D_5A^2} +o(ε^2).


Because it contains the secular terms, the solution diverges as |t-t_0| 

goes infinity.
To resum the would-be divergent terms, 
we apply the RG equation dx_3(t, t_0)/dt_0|_t_0=t, 
which leads to the equations governing the amplitude and phase as

    dA/dt=ε s_3A/2(s_1^2+s_2),

and 

    dθ /dt=ε s_1s_3/2√(s_2)(s_1^2+s_2),

nicely describing the slow motions of the amplitude and phase, respectively.
However, it fails to describe a transitional behavior approaching a limit cycle, as indicated by the present model.

Therefore, we analyzed the third-order equation, which might lead to 
a limit cycle solution. 
Substituting Eqs. (<ref>) and (<ref>) into Eq. (<ref>), we have


    d^3u_3/dt^3+s_1d^2u_3/dt^2
    +s_2du_3/dt+s_1s_2u_3
    =
    E_1A^3cos(ω_0t+θ)-E_2A^3sin(ω_0t+θ)
    
       -E_3A(t-t_0)cos(ω_0t+θ)
    +E_4A(t-t_0)sin(ω_0t+θ)
    +E_5A^2cos(2ω_0t+2θ)
    
       +E_6A^2sin(2ω_0t+2θ)
    +E_7A^2(t-t_0)cos(2ω_0t+2θ)
    -E_8A^2(t-t_0)sin(2ω_0t+2θ)
    
       -E_9A^3cos(3ω_0t+3θ)
    -E_10A^3sin(3ω_0t+3θ)
    +E_11A^2+E_12A^2(t-t_0).


where


    E_1=(s_1^2s_2-4s_1^2s_3+6s_1s_2^2-18s_2s_3)s_4s_1s_2/12(s_1^2+4s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       E_2=s_4^2s_1^2√(s_2^3)/12(s_1^2+4s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       E_3=s_3^2/2(s_1^2+s_2),
    
       E_4=s_1s_3^2/2(s_1^2+s_2)√(s_2),
    
       E_5={7s_1^3s_2-4s_1^2s_3+24s_1s_2^2-12s_2s_3/12(s_1^2+4s_2)
    -s_4s_3/4s_1s_2ln( s_3/p_1p_2r)
    }( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       E_6=s_4s_1√(s_2)/6(s_1^2+4s_2)( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       E_7=s_4s_1s_2/4(s_1^2+s_2)( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       E_8=s_4s_1^2√(s_2)/4(s_1^2+s_2)( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       E_9=((s_1^2+2s_2)s_1+2s_3)s_4s_1s_2^2/12(s_1^2+4s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       E_10=s_4^2s_1^2√(s_2^3)/12(s_1^2+4s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       E_11=1/4{
    s_1s_2
    -s_4s_3/s_1s_2ln( s_3/p_1p_2r)
    }( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       E_12=s_4s_1s_2/4(s_1^2+s_2)( s_3/p_1p_2r) ^s_3/s_1s_2.


The solution to Eq. (<ref>) 

is given by 

    u_3(t;t_0)=   -(F_1aA^3+F_1bA)(t-t_0)cos(ω_0t+θ)
    +(F_2aA^3+F_2bA)(t-t_0)sin(ω_0+θ)
    
       -F_3A(t-t_0)^2cos(ω_0t+θ)
    -F_4A(t-t_0)^2sin(ω_0t+θ)
    +F_5A^2cos(2ω_0t+2θ)
    
       +F_6A^2sin(2ω_0t+2θ)
    -F_7A^2(t-t_0)cos(2ω_0t+2θ)
    
       +F_8A^2(t-t_0)sin(2ω_0t+2θ)
    +F_9A^3cos(3ω_0t+3θ)
    +F_10A^3sin(3ω_0t+3θ)
    
       +F_11A^2+F_12A^2(t-t_0)

where

    F_1a=(2s_1s_2^2-(s_1^2+6s_2)s_3)s_4s_1s_2/8(s_1^2+4s_2)(s_1^2+s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       F_1b=s_1s_3^2/(s_1^2+s_2)^3,
    
       F_2a=((s_1^2+7s_2)s_1s_2-(4s_1^2+19s_2)s_3)s_4s_1^2√(s_2)/24(s_1^2+s_2)(s_1^2+4s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       F_2b=((s_1^2+6s_2)s_1^2-3s_2^2)s_3^2/8(s_1^2+s_2)^3√(s_2^3),
    
       F_3=(s_1^2-s_2)s_3^2/8(s_1^2+s_2)s_2,
    
       F_4=s_1s_3^2/4(s_1^2+s_2)√(s_2),
    
       F_5={s_4s_3/12(s_1^2+4s_2)s_2^2ln( s_3/p_1p_2r) .
    
       . -
    (3s_1^5+4s_1^3s_2+11s_1^2s_3+64s_1s_2^2-52s_2s_3)s_1/36(s_1^2+s_2)(s_1^2+4s_2)^2}( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       F_6={s_4s_3/6(s_1^2+4s_2)s_1√(s_2^3)ln( s_3/p_1p_2r) .
    
       . -
    7s_1^5s_2-s_1^4s_3-8s_1^3s_2^2+38s_1^2s_2s_3+48s_1s_2^3-24s_2^2s_3/36(s_1^2+s_2)(s_1^2+4s_2)^2√(s_2)}( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       F_7=s_4s_1^2/4(s_1^2+s_2)(s_1^2+4s_2)( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       F_8=(s_1^2-2s_2)s_4s_1/12(s_1^2+s_2)(s_1^2+4s_2)√(s_2)( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       F_9=(s_1(s_1^2-s_2)+5s_3)s_4s_1^2s_2/96(s_1^2+4s_2)(s_1^2+9s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       F_10=(2(2s_1^2+3s_2)s_1s_2+(6s_2-s_1^2)s_3)s_4s_1√(s_2)/96(s_1^2+4s_2)(s_1^2+9s_2)s_3^2( s_3/p_1p_2r) ^2s_3/s_1s_2,
    
       F_11={(s_1^3+s_3)/4(s_1^2+s_2)s_1
    -s_4s_3/4s_1^2s_2^2ln( s_3/p_1p_2r) 
    }( s_3/p_1p_2r) ^s_3/s_1s_2,
    
       F_12=s_4/4(s_1^2+s_2)( s_3/p_1p_2r) ^s_3/s_1s_2.


To obtain a globally valid solution using Eq. (<ref>), we apply

the RG method, which utilizes the RG equation

    . dx_3(t;t_0)/dt_0| _t_0=t
    =   . ∂ x_3(t;t_0)/∂ t_0| _t_0=t
    +. dA/dt_0∂ x_3(t;t_0)/∂ A| _t_0=t
    +. dθ/dt_0∂ x_3(t;t_0)/∂θ| _t_0=t
    
    =   {εdA/dt-ε D_1A+ε ^3(F_1aA^3+F_1bA) }cos(ω_0t+θ)
    
       +{ -ε Adθ/dt+ε ^2D_2A-ε ^3(F_2aA^3+F_2bA) }sin(ω_0t+θ)
    =0,

where we have neglected the higher-order terms o(ε^3).
For Eq. (<ref>) to hold for any t, the coefficients of the two independent functions should vanish. Thus, we obtain the dynamic equations for A and θ:

    dA/dt=ε D_1A-ε^2(F_1aA^3+F_1bA)
    +o(ε^2),
    
       dθ/dt=ε D_2-ε^2(F_2aA^2+F_2b)
    +o(ε^2).

The amplitude equation (<ref>) has a new fixed point 

    A_0=√(D_1-ε F_1b/ε F_1a)
    =√(4((s_1^2+s_2)^2-2ε s_1s_3)(s_1^2+4s_2)s_3^3/ε (2s_1s_2^2-(s_1^2+6s_2)s_3)(s_1^2+s_2)^2s_4s_1s_2)( p_1p_2r/s_3) ^s_3/s_1s_2,

which is 

nothing but the amplitude of the desired limit cycle.
The phase function θ(t) on the limit cycle is expressed as 

    θ (t)=(ε D_2-ε^2(F_2aA_0^2+F_2b)+o(ε^2))t+θ _0,

where θ _0 is the integral constant and it gives 
the initial phase at t=0. 
Substituting Eqs. (<ref>), (<ref>), (<ref>) and (<ref>) into (<ref>), θ (t) is reduced to

    θ (t)={ 
    -εs_1s_3s_4/6(2s_1s_2^2-(s_1^2+6s_2)s_3)√(s_2)+o(ε)
    }
    t +θ _0.

Thus, the solution describing the limit cycle, which is valid in a global domain in the asymptotic regime, is given by

    x_3(t)=   x_3(t;t_0)|_t_0=t
    
    =   x_0+ε A_0cos(ω t+θ _0)
    
       -ε^2{ D_3A_0^2cos(2(ω t+θ_0))
    +D_4A_0^2sin(2(ω t+θ_0))
    } +o(ε ^2).

In particular, if the initial phase is set to be θ_0 =-π /2, we have Eq. (<ref>).



 §.§ A.3 Numerical and 
 RG analyses of the Lotka-Volterra model


Equation (<ref>) and the numerical simulation indicate that non-sinusoidal power (NS) tends to be larger when the period
hardly changes and is stable in response to increases of the parameter values specifying the degradation rates, as presented in Fig. <ref>B. This suggests that the waveform becomes more distorted at higher temperatures when the circadian period is temperature-compensated. 
It was previously reported that the same conclusion holds for other oscillatory models, including a realistic mammalian circadian clock model, a post-translational model in cyanobacteria, and the van der Pol oscillator <cit.>. 
However, it is important to note that the findings using specific mathematical models might not be universally applicable to other models and actual organisms. 
Therefore, one should examine whether the waveform also plays a crucial role in the stability of the period in other oscillatory models. 

Thus, we conduct a numerical simulation to test the possible period-waveform correlation 
in the Lotka-Volterra model as done for the circadian clock model.
Needless to say, the Lotka-Volterra model is one of the most extensively studied mathematical models in biology <cit.>, and it effectively explains population dynamics in prey-predator systems.
The Lotka-Volterra model is given as a system with two variables as

    dx/dt=ax-ε xy,
    
       dy/dt=-by+ε 'xy,

where x(t) and y(t) are numbers of prey and predators. Parameter a is the growth rate of prey, b is the death rate of predators, ε is the death rate of prey attributable to predation, and ε ' is the growth rate of predators (Supplementary Fig. <ref>A). 




In the present numerical simulation, we first generated

100 reference parameter sets corresponding to the reference temperature.
The values of the model parameters a and b were generated randomly with a uniform distribution between 0 and 1, and similarly, the values of the other model parameters ε and ε ' 
were also randomly assigned values between 0 and 0.5. The period obtained for each parameter set with the initial condition (x_0,y_0)=(2,2) was denoted as τ_1. 
Next, a, b, ε, and ε ' 

were multiplied by a random factor within the range of 1.1-1.9 to simulate the increase in temperature, yielding 49 oscillatory parameter sets. Each resulting new period was denoted as τ_2, and thus, the ratio of the period =τ_2/τ_1, which is called the relative period,
was obtained.

The numerical simulations reveal a consistent positive correlation between NS and the relative period (Supplementary Fig. <ref>BC). 
A notable point is that the value of NS tends to increase along with the period when larger multiplicative factors are used in the simulation. 
This result again suggests that our findings that the waveform becomes more distorted when the period remains relatively stable in response to increased parameter values
is a rather universal phenomenon not restricted to the behavior observed in the circadian clock model (Fig. <ref>A).

The period of the Lotka-Volterra model, together with its approximate but globally valid solution, was previously derived analytically by one of the present authors <cit.> on the basis of the RG method; 
see also the pioneering work <cit.> based on a different method. 
Next, we will demonstrate that the expression explicitly reveals that the period of the Lotka-Volterra model almost linearly increases with the waveform distortion for the average of the prey-predator time series (i.e. NS). 
In the mathematical analysis, it proved convenient
to use new variables (ξ(t), η(t)) defined as

    x(t)=(b+εξ (t))/ε ' ,
       y(t)=a/ε +η (t).

The RG method performed in <cit.>
in the second order of ε leads to


    ξ(t)=   ( 1-ε^2a-b/4ab^2A^2/12) AsinΘ
    -ε^21/b√(ab)A^2/24cosΘ 
    -ε1/√(ab)A^2/6sin(2Θ) 
       -ε1/bA^2/3cos(2Θ)
    -ε^23a-b/4ab^2A^3/8sin(3Θ)
    +ε^21/b√(ab)A^2/8cos(3Θ) +o(ε^2),
    
    η(t)=   ε^21/b^2A^3/24sinΘ
    -√(ab)/b( 1+ε^2a-b/4ab^2A^2/12) AcosΘ 
    -ε√(ab)/b^2A^2/6sin(2Θ) 
       +ε1/bA^2/3cos(2Θ)
    +ε^21/b^2A^3/8sin(3Θ)
    +ε^2a-3b/4b^2√(ab)A^3/8cos(3Θ) 
    +o(ε^2),

where A and θ are the integral constants, which

are to be determined by the initial condition,

and Θ=ω̃ t+θ, with
ω̃ being the angular frequency given by


    ω̃=
    √(ab){ 1-ε^2A^2(a+b)/24ab^2},

from which we have the formula of the period of the system after some manipulation as

    τ=2π/ω̃=2π/√(ab){2/5( 1+5/2ε^2A^2(a+b)/24ab^2) +3/5} +o(ε^2).


From the waveforms for ξ(t) and η(t) given by (<ref>) and (<ref>), respectively, 
we can obtain the waveform distortion of each variable as follows:

    NS^(ξ) =1+ε^2A^2(4a+b)/24ab^2+o(ε^2), 
    
       NS^(η) =1+ε^2A^2(a+4b)/24ab^2+o(ε^2).


It is notable that the mean of NS^(ξ) and NS^(η) 
takes the form


    NS=1/2(NS^(ξ)+NS^(η))
    =1+5/2ε^2A^2(a+b)/24ab^2+o(ε^2).


Indeed, comparing (<ref>) and (<ref>), we arrive at


    τ=2π/√(ab)(2/5NS+3/5) 
    +o(ε^2),


which states that the period and the mean waveform distortion, namely NS, are linearly dependent on each other, and they tend to increase (or decrease) in a parallel manner. This is what we aimed to demonstrate for the
Lotka-Volterra model.




 §.§ A.4 
Numerical and RG analyses of the van der Pol model 


We consider the van del Pol model as follows:

    d^2x/dt^2+x=ε (1-x^2)dx/dt,


which is known as one of the fundamental non-linear oscillator models.
Previously, two of the authors derived the period formula τ =2π[∑_j=1^∞|a_j|^2j^2/∑_j=1^∞|a_j|^2]^1/2, meaning that the period of the model is also proportional to NS <cit.>. Then, we derive the approximate solution of the model using the RG method and confirm the proportionality between the period and the waveform distortion NS in detail. First, we represent the local solution around t=t_0 as a perturbation series

    x(t;t_0)=x_0(t;t_0)+ε x_1(t;t_0)+ε^2 x_2(t;t_0)+o(ε^2).


Then, substituting Eq. (<ref>) into Eq. (<ref>) and equating the terms with the same powers of ε, we obtain


    O(ε^0): d^2x_0/dt^2+x_0=0,
    
       O(ε^1): d^2x_1/dt^2+x_1
    =(1-x_0^2)dx_0/dt,
    
       O(ε^2): d^2x_2/dt^2+x_2
    =(1-x_0^2)dx_1/dt-2x_0x_1dx_0/dt.


The solution for the zeroth-order equation (<ref>) is


    x_0(t;t_0)=A(t_0)cos(t+θ(t_0)),

where A and θ are integral constants and they potentially depend on initial time t_0. Then, substituting Eq. (<ref>) into Eq. (<ref>), we have


    d^2x_1/dt^2+x_1
    =-A( 1-A^2/4) sin(t+θ)+A^3/4sin(3t+3θ).

By solving the solution of Eq. (<ref>) around t=t_0, the first-order solution is given by


    x_1(t;t_0)
    =A/2( 1-A^2/4) (t-t_0)cos(t+θ)
    -A^3/32sin(3t+3θ).

Similarly, by substituting zero-th and first-order solutions (<ref>) and (<ref>) into Eq. (<ref>), the second-order equation is

    d^2x_2/dt^2+x_2
    =   F_1(A)cos(t+θ)-F_2(A)(t-t_0)sin(t+θ)+F_3(A)cos(3t+3θ)
    
       -F_4(A)(t-t_0)sin(3t+3θ)+F_5(A)cos(5t+5θ),

where 

    F_1(A)=A/2( 13/64A^4-A^2+1 ) ,
    
       F_2(A)=A/2( 3/16A^4-A^2+1 ),
    
       F_3(A)=A^3/32( 5/2A^2-7 ),
    
       F_4(A)=3A^3/8( 1/4A^2-1 ),
    
       F_5(A)=5A^5/128.

The solution of the second-order equation (<ref>) is

    x_2(t;t_0)
    =   1/4(2F_1(A)-F_2(A))(t-t_0)sin(t+θ)
    +1/4F_2(A)(t-t_0)^2cos(t+θ) 
    
       +1/32(-4F_3(A)+3F_4(A))cos(3t+3θ)
    +1/8F_4(A)(t-t_0)sin(3t+3θ)
    
       -1/24F_5(A)cos(5t+5θ).

Therefore, the perturbative solution up to the second-order of ε is

    x(t;t_0)
    =   Acos(t+θ)
    +ε{A/2( 1-A^2/4) (t-t_0)cos(t+θ)
    -A^3/32sin(3t+3θ)
    }
       +ε^2{1/4(2F_1(A)-F_2(A))(t-t_0)sin(t+θ)
    +1/4F_2(A)(t-t_0)^2cos(t+θ) .
    
       +1/32(-4F_3(A)+3F_4(A))cos(3t+3θ)
    +1/8F_4(A)(t-t_0)sin(3t+3θ)
    
       .
    -1/24F_5(A)cos(5t+5θ).
    } +o(ε^2) .

To obtain the globally valid solution using Eq. (<ref>), we apply

the RG method, which utilizes the RG equation

    . dx/dt_0|_t_0=t=
       {dA/dt+ε1/2A( 1/4A^2-1 ) }cos(t+θ)
    
       +{ -Adθ/dt -ε^21/4(2F_1(A)-F_2(A)) }sin(t+θ)=0,

where we have neglected the higher-order terms o(ε^2).
For Eq. (<ref>) to hold for any t, the coefficients of the two independent functions should vanish. Thus, we end up with the 
dynamical equations for A and θ

    dA/dt=1/2ε A( 1-1/4A^2)
    
       dθ/dt=-ε^21/4A(2F_1(A)-F_2(A))
    =-ε^21/8( 7/32A^4-A^2+1 )


Equation (<ref>) has two fixed points, namely A=0 and 2. The amplitude A asymptotically approaches A=2, which is the limit cycle. 
Therefore, the dynamical behavior of θ (t) on the limit cycle reads

    θ(t) = -ε^21/16 t +θ_0,

where θ_0 is the initial phase at t=0. Thus, the globally valid solution on the limit cycle is given by


    x(t)=2cos(ω t+θ_0)
    -ε1/4sin(3ω t+3θ_0)
    -ε^23/32cos(3ω t+3θ_0)
    -ε^25/96cos(5ω t+5θ_0)
    +o(ε^2),


where the angular frequency ω up to the second order is expressed as

    ω = 1-ε^21/16+o(ε^2).


Using Eqs. (<ref>) and (<ref>), the period and NS of the van der Pol model read


    τ =2π/ω=2π/1-ε^2/16+o(ε^2)
    = 2π( 1+ε^21/16 +o(ε^2) ) ,
    
       NS=[ ∑_j=1^∞|a_j|^2j^2/∑_j=1^∞|a_j|^2]^1/2
    = [ 1+ε^2(9/64)+o(ε^2)/1+ε^2/64+o(ε^2)]^1/2
    = 1+ε^21/16+o(ε^2),


Accordingly, 


    τ =2π NS.


This expression clearly illustrates that the period τ and the waveform distortion NS in the van der Pol model tend to increase (or decrease) together in a proportional manner, which was also demonstrated using signal processing methods (Supplementary Fig. <ref>) <cit.>. 





§ SUPPLEMENTARY FIGURES






§ SUPPLEMENTARY TABLES
































































