Amplitude effects allow short jetlags and large seasonal phase shifts in minimal clock models

Mathematical models of varying complexity have helped shed light on different aspects of circadian clock function. In this work, we question whether minimal clock models (Goodwin models) are sufficient to reproduce essential phenotypes of the clock: a small phase response curve (PRC), fast jetlag and seasonal phase shifts. Instead of building a single best model, we take an approach where we study the properties of a set of models satisfying certain constraints; here a PRC to a one-hour pulse of three hours and clock periods between 22h and 26h. Surprisingly, almost all these randomly parameterized models showed a phase shift of about four hours between long and short days and jetlag durations of three to seven days in advance and delay. Moreover, intrinsic clock period influenced jetlag duration and entrainment amplitude and phase. Fast jetlag was realized in this model by means of a novel amplitude effect: the association between clock amplitude and clock period termed ‘twist’. This twist allows amplitude changes to speed up and slow down clocks enabling faster shifts. These findings were robust to the addition of additional positive feedback to the model. In summary, the known design principles of rhythm generation – negative feedback, long delay and switch-like inhibition (we review these in detail) – are sufficient to reproduce the essential clock phenotypes. Furthermore, amplitudes play a role in determining clock properties and must be always considered, although they are difficult to measure.

In this paper, we step back and address the question whether or not simple models can reproduce phenotypic features, such as phase response curves (PRCs), entrainment, jetlag, and seasonality. We askwhat can we learn from a systematic analysis of generic models about the underlying mechanisms? How oscillator properties, such as period and amplitude, govern phenotypic features including jetlag and seasonality? We find that ensembles of quite basic models reproduce these phenotypic features surprisingly well. Moreover, it turns out that amplitudes have a strong effect on jetlag duration and entrainment phase.
As a starting point we summarize now major phenomenological features of circadian clocks: PRCs, jetlag, and seasonality. Then we introduce simple models of delayed negative feedback loops with switchlike inhibition. seasonal phase shift of the model in Figure 2A is about four hours ( Figure 2C). The phase of entrainment approximately tracks the midpoint of the dark phase. We also subjected the models to advancing and delaying jetlag of six hours under a 24h (LD 12:12) Zeitgeber input ( Figure 2D,E). An abrupt change in Zeitgeber phase (jetlag) shifts the phase of the model until a stable phase relative to the new Zeitgeber phase is achieved. The course of this change is the jetlag transient. The jetlag shifts simulated here represent flying from Berlin to New York (delay) or New York to Berlin (advance). The model shifted within 4-5 days to within 15 min of the final stable phase for both advancing and delaying jetlags. Moreover, the delaying jetlag shift was shorter than the advancing jetlag shift in this model.

Short jetlag and large seasonal shifts are common
The Monte Carlo approach outlined in the previous section generated an ensemble of models (N = 1186). All the models were self-sustained oscillators and had a PRC to a one-hour light pulse with a range of about three hours by design ( Figure S1); we allowed a tolerance of 0.25h in order to speed up optimization of the Zeitgeber strength. We present here ( Figure 3) the landscape of the duration of jetlag and magnitude of seasonal phase shifts within this ensemble.
The duration of the jetlag transient depends on the phase of the Zeitgeber at the instant of the jetlag shift in addition to the magnitude and direction of the jetlag shift. Therefore, jetlag was measured as the median jetlag duration starting at four phases 6h apart in order to capture the dependence on the starting phase. Since advancing and delaying jetlags were different (Figure 2), we present results for 6h advancing and 6h delaying jetlags separately.
The whole ensemble exhibited short jetlag shifts (about 3-7 days in length) and seasonal phase shifts of about four hours that matched the observed phenotypes described earlier. We expected specific optimized parameter sets alone to produce the necessary phenotypes [35]. Contrary to our expectation, randomly-chosen parameter sets with the same PRC yielded a highly clustered combination of reasonable jetlag duration and seasonal phase shift.
Advancing jetlag shifts were slower than delaying jetlag shifts ( Figure 3). This might be anticipated from the larger delay region compared to a advance region in the PRC ( Figure 2B). In fact, PRCs of all models in the ensemble had approximately a maximum advance of one hour and a maximum delay of two hours ( Figure S1). We can roughly estimate (using the PRC) that an advance of one hour or delay of two hours requires a day. If shifts on subsequent days can be carefully adjusted, the simulated 6h jetlags would need six days for the advance or three days for the delay. Surprisingly, the randomly-parameterized models in the ensemble shifted as fast or sometimes faster than the rough estimates. We examine, in the remaining sections, the ensemble properties to better understand the mechanisms that result in the landscape in Figure  3.

Entrainment phases correlate with clock period
The theory of entrainment describes the relationship between the Zeitgeber strength, Zeitgeber period, intrinsic clock period and intrinsic clock amplitude in a certain clock (model) needed for entrainment. Entrainment has occurred once the clock achieves a fixed phase (called phase of entrainment ψ) with respect to the Zeitgeber. Instead of observing a single oscillator under a range of Zeitgeber periods, we observed the ensemble (comprised of models with different intrinsic clock periods) after entrainment to a 24h (LD 12:12) Zeitgeber. The heterogeneity of human circadian clocks was similarly used to study the entrainment properties of humans [40,41].
Generally, short period clocks had earlier ψ and long period clocks had later ψ ( Figure 4A). The slope of ψ versus clock period τ was small (slope ∼ 0.5) across the ensemble. Unicellular organisms and plants under entrainment show such small slopes [42]. Moreover, entrainment theory predicts ∼12h change in the phase of entrainment across the entrainment range of an oscillator [2,30]. This suggests a large entrainment range is necessary to achieve this small slope. We shall return to the range of entrainment in Figure 8.
The entrained amplitudes of the models in the ensemble, on the other hand, showed a small but significant decreasing trend with intrinsic clock period τ ( Figures 4B-D); such an amplitude trend has been previously reported [43]. However, the intrinsic (unentrained) amplitude was uncorrelated with intrinsic clock period (ρ = 0.03, p-value=0.21). The observed trend arises from the phenomenon of 'twist' in these models; twist is an association between the amplitude and period of an oscillator. The representative model, for example, has a positive association between period and amplitude ( Figure 6C). Therefore, when the Zeitgeber period is shorter than the intrinsic clock period, the entrained amplitude decreases (e.g., τ = 26h). The corresponding argument for shorter τ explains the decreasing trend in data. We will return to this amplitude effect in Figure 6.

Intrinsic periods affect length of jetlag
Both advancing and delaying jetlag were generally short throughout the ensemble. The PRC yielded reasonable rough estimates of jetlag duration. We could use the ensemble of models with their heterogeneous PRCs to test whether the PRC indeed can predict the jetlag duration. If this hypothesis were true, models with larger advance regions should have shorter advance jetlags and similarly models with larger delay regions should have a shorter delay jetlags. Models with larger delay regions (and smaller advance regions, since delay + advance region was chosen to be about three hours) indeed had shorter delay jetlag transients ( Figure S3A). Intriguingly, models with larger delay regions (and smaller advance regions) also had shorter advance jetlag transients ( Figure S3B). Therefore, the PRC alone is insufficient to describe the model's jetlag behavior. To gain further insight into the models, we explored the dependence of the duration of the jetlag on model parameter and properties.
Short intrinsic clock periods favor shorter advancing jetlags (ρ spearman = 0.25, p-value < 0.001, Figure  5A). Similarly, longer intrinsic periods favor shorter delaying jetlags, although the effect was less pronounced (ρ spearman = -0.10, p-value = 0.008, Figure 5B). During an advancing jetlag transient, the entrained clock needs to speed up (reduce its period/increase its phase velocity) to achieve the new Zeitgeber phase. Thus, the shorter intrinsic period (faster clock) naturally aids the jetlag shift. Likewise, response to a delaying jetlag shift requires slowing the clock, which is favored in long intrinsic period clocks.
The Goodwin model requires a Hill coefficient of at least eight in order to oscillate. The Hill coefficient controls the degree of nonlinearity in the negative feedback loop, a necessary component of an oscillator. Models with stronger nonlinearity shifted faster in response to both advance (ρ spearman = -0.56, p-value < 0.001, Figure 5C) and delay jetlag shifts (ρ spearman = -0.75, p-value < 0.001, Figures 5D). The Goodwin model with higher Hill coefficients is associated with higher amplitude relaxation rates and faster entrainment ( Figure S2) [44]. Amplitude relaxation rate measures the time taken for the effect of a perturbation to disappear [12,44]. We suspect the faster time to entrainment is closely related to the duration of jetlag transient.

Amplitude changes facilitate short jetlag
We earlier observed that both the amplitude and the jetlag duration varied with the instrinsic clock period (Figures 4 and 5). We wondered whether the amplitude of the oscillator played a role in aiding the short jetlags. Therefore, we observed the behavior of the representative model from Figure 2 during simulated advancing and delaying jetlags.
The amplitude of the clock changed significantly over the course of the jetlag transient ( Figure 6). The clock amplitude decreased from the steady-state entrained amplitude during a 6h advance ( Figure 6A), but returned to the amplitude prior to jetlag after the phase shift had been achieved. On the other hand, the clock amplitudes increased relative to the steady-state amplitudes during a 6h delay. This leads one to question how these amplitude changes aid in shifting the clock phase during jetlag.
To this end, we characterized the twist or amplitude-period association of this model without Zeitgeber input. We observed how the amplitude and period of the clock varied near the steady-state rhythm (limit cycle). The instantaneous amplitude and period of the clock as the perturbed (unentrained) clock returns to the limit cycle is used to measure twist ( Figure 6C). It is clear that there is an approximately linear relationship between amplitude and period. Smaller periods are associated with smaller amplitudes and larger periods with larger amplitudes. When the clock phase needs to be advanced, during a jetlag advance, a reduced amplitude is accompanied by a faster clock (i.e., smaller period), which facilitates the advance. Similarly, during a jetlag delay, where phase needs to be delayed, an increased amplitude coincides with a slower clock (i.e., larger period) that aids the delay. This interpretation is consistent with amplitude changes seen in Figures 6A,B. The amplitude of even the minimal models used here is difficult to define uniquely. The amplitude of each variable in the model can be measured, but it is unclear how these amplitudes can be reduced to a single metric for the whole model. Fortunately, consistent changes in amplitude occur across all three model variables during jetlag and the definition of amplitude does not affect our conclusions.

Observations are robust to addition of positive feedback
We have previously observed that many negative feedback oscillators have auxiliary positive feedback loops [39]. Moreover, the core clock network does include many positive feedback loops in the form of Michaelis-Menten (MM) degradation, complex formation and nuclear import-export. We therefore investigated the effect of a positive feedback in the form of MM degradation of variable Y on our observations. A representative model with positive feedback oscillates ( Figure S4A) and has a PRC with about an one hour advance and a two hour delay ( Figure S4B), similar to the Goodwin model. The transient after a 6h phase advance or 6h phase delay lasts about four to six days ( Figure S4D,E) and the seasonal phase shift between long and short day is also around four hours ( Figure S4C).
At the ensemble level, the landscape of jetlag duration and seasonal phase shift ( Figure 7A) was tightly clustered like the landscape of the Goodwin model ( Figure 3). However, the jetlag duration for both delay and advance was slightly longer with positive feedback. Delay jetlag remained shorter than advance jetlag even with positive feedback. The seasonal phase shift was also consistently shorter by about 0.5h with positive feedback. Nevertheless, the qualitative features of short jetlag and seasonal phase shift was robust to these model changes.
The phase of entrainment ψ varied by almost the same slope with positive feedback ( Figure 7B). This suggests that the entrainment properties of the model are conserved. When we examined the dependence of short jetlag on the model parameters, the jetlag duration was very weakly affected by intrinsic clock period ( Figure 7C). The advancing jetlag transients were unaffected by the circadian period (ρ spearman = -0.05, p-value = 0.18), while delaying jetlag transients possessed a weaker reversed trend with respect to instrinsic clock period (ρ spearman = 0.11, p-value = 0.002). Higher Hill coefficients consistently favored shorter jetlag duration for both advance (ρ spearman = -0.54, p-value < 0.001) and delay (ρ spearman = -0.52, p-value < 0.001, Figure 7D). Positive feedback allows oscillations even for smaller Hill coefficients (below eight) [45], which can be seen from the greater spread of Hill coefficients in the ensemble. Furthermore, short jetlag was obtained consistently along a range of Hill coefficients. Amplitude changes were also observed in the representative model with positive feedback ( Figure  S4F,G). As with the Goodwin model, jetlag phase advances were accompanied by amplitude reduction and jetlag phase delay by amplitude expansion. This favors short jetlag transients, since the model even with positive feedback has a positive amplitude-period association or twist ( Figure S4H) [46]. This allows speeding up the clock with amplitude reduction and slowing down the clock with amplitude expansion that can be used to efficiently overcome jetlag.

Discussion
Circadian rhythms are self-sustained oscillations that can be entrained by light-dark cycles. Mathematical models of these rhythms exhibit limit cycles driven by pulses or periodic inputs from a Zeitgeber. Any useful model should reproduce known basic phenotypic features, such as the PRC, phase of entrainment and jetlag duration.
In our paper, we studied ensembles of generic feedback models to explore systematically oscillator properties that allow short jetlags and large seasonal phase shifts. These generic models are minimal and contain just six parameters. The Zeitgeber strength was adjusted to get a PRC with a range of three hours for a one-hour light pulse [27]. The remaining parameters were randomly sampled in reasonable ranges. Models exhibiting self-sustained oscillations with an intrinsic period between 22h and 26h were analyzed in detail. The resulting set of hundreds of models could be exploited to study connections between oscillator properties and phenotypic features. Figure 2 demonstrates that a minimal model with a reasonable PRC has large seasonal flexibility and short jetlags. We show in Figure 3 that these features are reproduced by the whole ensemble of models. Furthermore, we analyzed how periods, amplitudes and entrainment phases are related ( Figure 4) and studied dependencies of jetlag duration on intrinsic periods and strengths of nonlinearities ( Figure 5). Figure 6 illustrates an interesting mechanism: amplitude changes can affect jetlag duration. It has been described earlier in experimental studies that small amplitude oscillators are easy to reset [12,[47][48][49][50]. In these examples mutations and compromised coupling improved resetting and shortened jetlag. In Figure  6C, we point to another amplitude effect: correlations of amplitudes and periods (termed "twist") imply that Zeitgeber-induced amplitude changes can speed up or slow down oscillations. Consequently, transient amplitude changes as shown in Figures 6 and S4 can tune jetlag duration.
In order to evaluate the robustness of our result, we generated another ensemble of models with an auxiliary positive feedback. Figures 7 and S4 demonstrated that parameter dependencies and reproduced features do not change drastically. Many of the described results can be quantified by extensive bifurcation analysis of the models [19]. Figure 8 shows an example of a two-dimensional bifurcation diagram for the representative model: the entrainment range, phase ( Figure 8A) and amplitude ( Figure 8B) are presented for varying Zeitgeber period and photoperiod. It turns out that entrainment phases and amplitudes vary widely. Selected model timecourses ( Figure 8C) demonstrate again that amplitudes and entrainment phases co-vary.
Another interesting observation is the abrupt widening of the entrainment range for long days. Bifurcation analysis ( Figure S5) reveals that this finding results from systematic changes due to increased light on long days. As found previously in another model [19,51], increasing input can transform self-sustained oscillations to damped oscillations (a reverse "Hopf bifurcation"). Damped oscillations can be driven by any extrinsic periodic force and hence the "entrainment range" is huge. However, the amplitudes decrease drastically for large Zeitgeber periods ( Figure 8B).
In the following sections, we discuss general aspects of mathematical models of circadian clocks. We emphasize that models are useful on different levels, that long days and switches are necessary for selfsustained oscillations, and that amplitudes are quite important to understand jetlag duration and seasonal flexibility.

Modeling can be applied on different levels
We analyzed generic models of a negative feedback loop. The required delay was achieved by a chain of reactions and slow degradation rates. The switch-like inhibition was modeled by a high Hill coefficient. Such a Goodwin model was motivated by self-inhibition of the Period and Cryptochrome genes and the Bmal1-Rev-erba loop. However, in contrast to previous models [15,22,52,53], we did not fit the models to specific expression data. This implies that our models, while not quantitative in detail, can nevertheless be applicable in a wider context including in organisms other than mammals.
We emphasize that mathematical models of circadian rhythms can be applied on the single cell level describing gene-regulatory feedback loops. Other models might be appropiate on the tissue level, or on the organismic level. The available experimental data are quite heterogeneous. On the genetic level, the phenotypes of many mutants and downregulations are available [54,55]. On the organ level, many expression profiles have been measured [56,57]. On the organismic level many data on periods, entrainment phases and PRCs are documented [24,42]. Consequently, also a great diversity of mathematical models have been developed including phenomenological amplitude-phase models [58], oscillator networks [51], explicit delay models [52,59] and detailed kinetic gene regulation models [60]. Complex models can reproduce multiple mutants [20,21] and pharmacological interventions [60].
Here we studied basic Goodwin models [61] that reproduced PRCs, jetlag, and seasonal adaptation surprisingly well. It has been shown earlier that mechanisms of temperature compensation could be also addressed with a simple Goodwin model [13]. Such generic models focus on the most important features of self-sustained oscillations: delayed negative feedback loops and switch-like inhibition.

Circadian rhythm generation requires a 6 hour delay
Under quite general conditions mathematical theory predicts that self-sustained oscillations with a period of 24 hours require a delay of at least 6 hours [52,62,63]. Most transcriptional-translational feedback loops introduce delays of about an hour and thus the corresponding periods are in the range of a few hours [64][65][66]. In order to reach periods of 24 hours, specific extra delays are required. In Drosophila, nuclear import of PER and TIM proteins is delayed by about six hours [67]. In mammals and Neurospora crassa multiple phosphorylations of the intrinsically disordered proteins PER2 and FRQ contribute to long delays [68][69][70][71]. Furthermore, cytoplasmatic and nuclear complex formation [72] and epigenetic regulations [73] can induce delays.

Possible mechanisms of switch-like inhibition
In our simple models, a high Hill coefficient is needed to get self-sustained oscillations [45]. Extended models have more realistic exponents due to longer reaction chains [74], positive feedback [39], protein sequestration [75] and protein inactivation [76]. In any case, switch-like inhibitions are required to get selfsustained rhythms. Theoretical studies have shown that switches can be obtained by cooperativity [77,78], positive feedbacks [79], competition of antagonistic enzymes [80], and by sequestration [81]. In mammalian clocks, the most essential nonlinearities ("switches") are not precisely known. It is likely, that multiple phosphorylations [82], complex formations [72], and nuclear translocation [69] play a role. Presumably, the competition of histone acetylation and deacetylation plays an important role since several HATs and HDACs influence clock properties [83]. Moreover, histone modifications are modulated by PER and CRY binding [72]. Similar cellular mechanisms have been described in the Neurospora clock [84].

Conclusions
Our model simulations show that amplitudes have pronounced effects on jetlag duration and entrainment phase. Typically, periods are the main focus of chronobiological studies since there are established devices, such as running wheels or race tubes, to measure periods. The measurement of amplitudes is less straight forward. Quantification of activity records [85] or conidation in race tubes [86] reflects just specific outputs. Gene expression profiles have to be carefully normalized for amplitude quantification [52]. Furthermore it is not immediately evident which gene is representative of the core clock amplitude. In Neurospora, FRQ and WCC amplitudes have been considered [87,88]. In mammals, PER2 [89] and Rev-erba levels [90] have been used to quantify amplitudes. Measuring the response to pulses provides an indirect quantification of limit cycle amplitudes [89,91,92]. Small amplitude limit cycles exhibit larger pulse induced phase shifts than large amplitude limit cycles [25]. Despite the complexities in quantifying amplitude, we suggest studies take into consideration amplitude in addition to the standard phase-based metrics, such PRCs.
Our comprehensive analysis of generic feedback oscillators provides insight into the design principles of circadian clocks. The essential role of fine-tuned delays and molecular switches is emphasized by mathematical modeling. Reasonable jetlag durations and remarkably large seasonal variations of the entrainment phase can be reproduced by quite simple models. Our models stress the important role of circadian amplitudes. Simulations show that amplitudes contribute significantly to phenotypic features such as jetlag duration and seasonal adaptation.

Model structure -Goodwin model
We chose the Goodwin model to represent the core circadian clock feedback loop like many earlier works [13,51,93] (Figure 1B). Theoretically, the three variables in this model are the minimum number required to produce self-sustained oscillations. We further reduced the number of parameters in the standard model by using the fact that amplitudes can be rescaled; this is often called nondimensionalization. This yielded a model with five parameters (d 1 , d 2 , d 3 , K, h): If the model oscillates for a particular choice of parameters, the period of oscillation τ is some complex function of the parameters. Scaling the equations produces a model with a desired period of oscillation τ : Finally, the input to the clock from the Zeitgeber was modeled as an addition to the rate of production of variable X. The Zeitgeber input consisted of the Zeitgeber strength L times I(t), where I(t) = 1 during the day (light phase) and I(t) = 0 during the night (dark phase). The model used for studying entrainment, jetlag simulation and seasonal phase shift was:

Monte Carlo simulations
In order to study the generic behavior of the Goodwin model, we resorted to the Monte Carlo approach. In this approach, model parameters are chosen randomly to obtain an ensemble of models that can be studied for a variety of phenotypes. We now outline the schema of our simulations.
We randomly generated values for the five model parameters. In addition, we also chose a random intrinsic clock period for each model in order to have a uniform distribution of intrinsic periods in the ensemble. In this manner, 2000 different sets of six parameters were generated using latin hypercube sampling (LHS). LHS ensures that parameter sets explore the entire chosen range of each parameter. The degradation parameters could vary between 0.1 and 0.3 to get oscillation periods around 24h. The Hill coefficient could lie between five and twelve. The Goodwin model requires a Hill coefficient of at least eight [45], but the model with positive feedback can oscillate also with smaller Hill coefficients [39]. The half-maximum concentration of the Hill term was chosen between 0.25 and 1 in line with common values taken by the variable Z (see (2)). Finally, the instrinsic period was restricted to 22h to 26h based on the range of periods seen in mammalian clocks.
For each random choice of parameter values, we tested whether the model oscillated robustly. If it did, then we optimized the Zeitgeber strength L such that the PRC to a one hour pulse had a maximum phase range (largest advance -largest delay) of 3 ± 0.25h. Bear in mind that scaling the model to reduce parameters makes the oscillator amplitudes arbitrary. The entrainment behavior of the model depends on the ratio between the Zeitgeber strength and the oscillator amplitude. We used the constraint on the PRC to keep effect of the Zeitgeber consistent across models in the ensemble.
With the optimized Zeitgeber strength L * , we ensured that the model can entrain to a 24h (LD 12:12) Zeitgeber. Upon entrainment, the phase of entrainment ψ is computed as the acrophase of variable Z relative to the lights-on (0 to L transition) of the Zeitgeber. Only when the model oscillates and is entrainable, we computed phenotypes of the model.
For each model, the seasonal phase shift is the difference between the phase of entrainment ψ for long day (24h LD 16:8 Zeitgeber) and short day (24h LD 8:16 Zeitgeber). Furthermore, for each model, we simulated jetlags with both a 6h phase advance and a 6h phase delay. This ensured also that our comparisons are well defined and data for advance and delay are always paired. Since the jetlag duration depends on the Zeitgeber phase when the phase shift occurs, we computed the median of the jetlag durations for 6h shifts starting from four different initial Zeitgeber phases (0, 6h, 12h, 18h). We measured jetlag duration from time of the Zeitgeber shift to the time point after which the phase of entrainment ψ of the clock stayed within 15 min of the final (shifted) phase.

Addition of positive feedback
We have previously defined the basic Goodwin model with different possible auxiliary positive feedback loops [39] and their significance in the mammalian circadian clock. Here, we used a variation of Goodwin model with Michaelis-Menten (MM) degradation of the variable Y ( Figure 1C): where c controls the strength of the positive feedback. In reality, for each simulated standard Goodwin model, we computed the phenotypes for the related positive feedback model (4) by randomly choosing a c value between 0.5 and 5. Thus, the comparisons between standard Goodwin and Goodwin with positive feedback models are well-defined, since parameter sets/models are paired. In practice, LHS generated sets of seven variables (five model parameters + instrinsic period + positive feedback strength) and these values were used to simulate the Goodwin model and the corresponding positive feedback model. The computation of model phenotypes was exactly as previously described.

Arnold Onion
The Arnold Onion visualizes the entrainment behavior of the representative Goodwin model (Figure 2) for different combinations of Zeitgeber period and photoperiod. For each combination of Zeitgeber period T and photoperiod, we numerically solved (2) for a total integration time of 2020 entrainment cycles T at a step size of ∆t = 0.01h. Transient dynamics within the first 2000 cycles were neglected. Next, we saved the complete state of the system at the beginning of the 2001st entrainment cycle and determined whether the system returns to the close neighborhood of this state, given by an -ball of radius 0.01 as described in [19]. Stable recurrence times equal to the Zeitgeber period T suggests an entrained state of the Goodwin oscillator. This is cross-validated by demanding a minimal of oscillation peak values of the Z variable. Phases within the entrained region were determined as a the distance between dawn, i.e., lights-on, and the peak time of the oscillations of the Z variable. Amplitudes are determined as half of the difference between the peak and trough values.

Floquet Multipiers
The speed at which the model returns to the limit cycle after a perturbation is measured by the amplitude relaxation rate. The average relaxation rate over one cycle can be calculated using the continuation software AUTO07p [94]. The software computes the so-called Floquet multipliers for any model defined using ordinary differential equations. These multipliers then are related to the amplitude relaxation rate.

Bifurcation Analysis
The bifurcation analysis of the model in Figure 2 was performed using XPPAUT with a constant Zeitgeber input, whose intensity was varied ( Figure S5). Our approach is outlined in [95].