Kuramoto model with a bimodal frequency distribution

A system of N coupled phase oscillators whose dynamical evolution is given by

$$\begin }_i(t) = \omega _i + \sum _^N K_(\theta _j-\theta _i) \end$$

(1)

is known as a Kuramoto model; see also Eq. (14) of the Main text. In the section Photoperiodic encoding through network re-organizations, we assume a functional separation of the N oscillators into two groups or communities representing, e.g., the core and shell part of the mammalian core pacemaker, the suprachiasmatic nucleus (SCN). In general, the communities can be of different size as described by the fractions \(p \in [0,1]\) and \((1-p)\), respectively, with different mean frequencies \(\omega _1\) and \(\omega _2\) as well as different frequency spreads (or scale factors) \(\gamma _1\) and \(\gamma _2\) as given by a bimodal Cauchy–Lorentz distribution

$$\begin g(\omega ) = \frac + \frac \quad . \end$$

(2)

Using the Ott–Antonsen approach (Ott and Antonsen 2008, 2009), the temporal evolution of the communities order parameters \(R_1(t)\) and \(R_2(t)\) as well as the phase difference between the clusters \(\Delta \varphi (t)=\varphi _2(t)-\varphi _1(t)\) under the assumption of identical intra- and inter-community coupling strength \(K_=K\) reads as

$$\begin }_1&= - \gamma _1 R_1 + \frac (1-R_1^2)(pR_1+qR_2\cos (\Delta \varphi )) \end$$

(3)

$$\begin }_2&= - \gamma _2 R_2 + \frac (1-R_2^2)(qR_2+pR_1\cos (\Delta \varphi )) \end$$

(4)

$$\begin \Delta }&= \omega _2 - \omega _1 - \frac \sin (\Delta \varphi ) \frac.\end$$

(5)

For the sake of simplicity, we assume in Fig. 5 that both communities are of equal size \(p=1-p=0.5\) and have an identical frequency spread \(\gamma _1=\gamma _2=\gamma\). Under such conditions, one can focus to solutions satisfying the symmetry condition \(R(t):=R_1(t)=R_2(t)\), such that above equations facilitate to

$$\begin \dot&= \frac R \left( 1 - \frac - R^2 + (1-R^2) \cos (\Delta \psi )\right) \end$$

(6)

$$\begin }&= \Delta \omega - \frac (1+R^2) \sin (\Delta \psi ) \end$$

(7)

with \(\Delta \omega = \omega _2 - \omega _1\), being in line with results of Martens et al. (2009).

Numerical solutions

Numerical simulations underlying Figs. 4b and 5f have been obtained by integrating Eq. (14) via the odeint function of the SciPy package. Bifurcation diagrams based on Eqs. (6, 7) as depicted in Fig. 5e are obtained via the XPP-AUTO package as previously described (Schmal et al. 2014), using the parameters \(Ntst = 150\), \(N\max = 20,000\), \(Dsmin = 0.0001\), and \(Dsmax = 0.0005\).

Seasonal entrainmentConceptual models explain complex data

Konopka and Benzer were the first who discovered a single-gene mutation that affects circadian free-running rhythms in Drosophila melanogaster (Konopka and Benzer 1971), leading to a new era of molecular genetics in chronobiology that eventually revealed the molecular constituents of circadian clocks across various organisms, such as cynobacteria, Neurpospora crassa, Arabidopsis thaliana, Drosophila melanogaster, as well as mammals (Bell-Pedersen et al. 2005). Even long before the molecular cogs and levers of the regulatory feedback loops underlying circadian rhythm generation have been found, conceptual oscillator models have been developed and used to understand circadian behavior and photoperiodic responses (Wever 1964; Pavlidis 1967; Winfree 1967). Such conceptual or generic oscillator models do not consider molecular details specific to certain organisms, tissues, or cell types but rather focus on general oscillator properties and their potential to explain observed experimental data (Roenneberg et al. 2008).

Phase oscillator models

One of the most abstract and simple conceptual models are phase oscillators. The only variable used to describe the circadian clock dynamics in this class of models is the phase of its oscillation \(\theta (t)\) essentially evolving between 0 and \(2\pi\). By this, we tacitly assume that the clock self-sustains its oscillation with a robust period \(\tau\) or angular velocity \(\omega =\frac\). Yoshiki Kuramoto introduced an intuitive way to describe the interaction between a given oscillator \(\theta (t)\) and a second oscillator \(\varphi (t)\) by means of a sinusoidal coupling term

$$\begin \frac = \omega + z \sin (\varphi (t) - \theta (t)), \end$$

(8)

such that oscillator \(\theta (t)\) slows down in case its phase advances the second oscillator \(\varphi (t)\) and speeds up in case it is delayed compared to \(\varphi (t)\) as the term \(\sin (\varphi (t) - \theta (t))\) in (8) becomes negative or positive, respectively (Kuramoto 1975, 2003).

Fig. 1

Weak zeitgebers or strong clocks lead to a large phase variability. a Arnold tongue based on phase oscillator model (9) in the parameter plane spanned by the internal free-running period \(\tau\) and the amplitude or strength z of the external zeitgeber signal. Color-coded values depict the phase of entrainment \(\psi\) as given by Eq. (11). b Experimentally obtained entrainment phases in dependence of the intrinsic free-running period \(\tau\) for ruin lizards subject to temperature cycles of different amplitude, i.e., zeitgeber strength. Data have been extracted from Fig. 5 of (Hoffmann 1969) via the WebPlotDigitizer software (Rohatgi 2022). c Arnold tongue in the parameter plane spanned by the period T and amplitude or strength z of the external zeitgeber signal. Color-coded values depict the phase of entrainment \(\psi\) as given by Eq. (11). d Experimentally obtained entrainment phases \(\psi\) for different species subject to entrainment cycles of different external zeitgeber period T. Species have been categorized into vertebrates (purple lines) as well as plants and unicellular species (brown lines). Please refer to the original publication (Aschoff and Pohl 1978) for the detailed description of the investigated animals and entrainment properties. Data have been extracted from Fig. 2 of (Aschoff and Pohl 1978) via the WebPlotDigitizer software (Rohatgi 2022)

Even such a simple model allows to explain a variety of experimental results and can help to better understand properties of the circadian clock. Assuming that a circadian clock with period \(\tau\) and described by phase variable \(\theta (t)\) is driven by an external zeitgeber \(\varphi (t)\) of period T, the dynamical evolution of the phase difference \(\psi (t)=\varphi (t)-\theta (t)\) between the zeitgeber and clock phase is governed by the well-known Adler equation

$$\begin \frac=\Delta \omega - z \sin (\psi (t)), \end$$

(9)

where \(\Delta \omega = \frac - \frac\) is the difference of the angular velocity of the zeitgeber and internal clock and z is the effective zeitgeber strength. Here, we tacitly assumed that there is no feedback from the clock to the zeitgeber, such that the entrainment cue can be described via \(\frac=\frac\), with T being the zeitgeber period. From (9), it follows that the circadian clock is only able to entrain to the zeitgeber signal for small enough frequency detunings (\(\Delta \omega\)) or high enough zeitgeber strengths (z) given by the condition

$$\begin \left\| \frac \right\| < 1. \end$$

(10)

The range of periods for which the internal clock entrains to the external zeitgeber is termed entrainment range. For a given zeitgeber period, e.g., \(T=24\)h, the entrainment range generally increases for increasing zeitgeber strength z, leading to a wedge shaped entrainment region in the \(\tau\)-z parameter plane, known as the Arnold tongue (Fig. 1a). For combinations of free-running periods \(\tau\) and zeitgeber strength z that lie within the Arnold tongue, the internal clock and external zeitgeber signal oscillate with a common period (frequency-locking) and adopt a stable phase relationship \(\psi\) (phase-locking). This phase of entrainment \(\psi\) is of fundamental importance for the proper alignment of physiological processes around the solar day and thus under evolutionary selection.

Our conceptual phase oscillator model (9) predicts that the dependency of the phase of entrainment \(\psi\) on zeitgeber strength z and the frequency detuning \(\Delta \omega\) is given by

$$\begin \psi ^=\arcsin \left( \frac\right) . \end$$

(11)

Thus, for any given zeitgeber intensity z, the phase of entrainment can vary only by 180° with respect to changes of \(\tau\) or T. Rütger Wever, a pioneer in mathematical modeling of the circadian system, described such 180° rule already in 1964 for a conceptual oscillator model adopted from electrical engineering, the Van der Pol oscillator, that has been originally developed to describe oscillations in electrical circuits employing vacuum tubes (Wever 1964). The 180-degree-rule predicts a smaller phase variability with respect to variations of the intrinsic free-running period \(\tau\) for increasing zeitgeber strength due to larger entrainment ranges, i.e., phases are less compressed (see color-coded area in Fig. 1a. This is in line with results from early entrainment experiments by Klaus Hoffmann, showing that entrainment ranges increase, while the phase variability decreases with an increasing zeitgeber amplitude (strength) for the ruin lizard Lacerta sicula subject to temperature cycles of \(T=24\)h period (Hoffmann 1969), compare Fig. 1b and corresponding arrows in Fig. 1a. Positions of arrows in Fig. 1a, i.e., the zeitgeber strength z estimated from the experimental data, have been determined by linear regressions (dashed lines) in Fig. 1b.

The resistance of a self-sustained oscillator to entrain to a certain zeitgeber signal can be used to define strong and weak clocks. While strong clocks are characterized by relatively small entrainment ranges and, thus, large phase variabilities with respect to changes in T, weak clocks exhibit broad entrainment ranges and small phase variabilities, corresponding to small and large zeitgeber strengths z in Fig. 1c. Along these lines, entrainment experiments can be used to categorize circadian pacemakers into weak and strong clocks and to infer internal oscillator properties of a given organism, tissue, or cell type of interest. Aschoff and Pohl summarized the entrainment behavior of 19 different species subject to entrainment cues of different zeitgeber period T (Aschoff and Pohl 1978). By comparing the dependency of the entrainment phase \(\psi\) to changes in T (i.e., the slope of curves in Fig. 1d) with results from our conceptual model (Fig. 1c), it turns out that vertebrate clocks rather behave like relatively strong clocks (Fig. 1c, purple arrow), while clocks of plants and unicellular species behave more like weak oscillators (Fig. 1c, brown arrow). Along these lines, above-described \(180^\circ\) rule has been used to show that strong oscillators like the vertebrate clock with a high phase variability are able to translate a narrow distribution of internal free-running periods \(\tau\) in a population with standard deviations of as little as \(\sigma =0.2\hbox \) for humans (Duffy et al. 2011) into the experimentally found large spread of human chronotypes which can be related to a large spread in the distribution of entrainment phases \(\psi\) (Roenneberg et al. 2004; Granada et al. 2013; Schmal et al. 2020). A similar reasoning has been used to argue that the weak circadian clocks as observed for organisms living at high latitudes such as certain Drosophila strains (Beauchamp et al. 2018) or rain deer (van Oort et al. 2005) could be an adaptive advantage as weak oscillators are able to entrain better under extreme photoperiodic conditions such as long summer days or long winter nights in comparison to strong clocks (Vaze and Helfrich-Förster 2016). Analogously to a weaker circadian clock, entrainment can be also facilitated by increasing an organisms light sensitivity as suggested in a comparative study of a northern and southern line of the parasitoid wasp Nasonia vitripennis (Floessner et al. 2023).

Broad applicability of phase oscillator models

The conceptual phase oscillator approach described in the previous section solely relies on the assumption that oscillators exhibit self-sustained oscillations and that interactions between clocks are weak in a way that amplitude effects can be neglected and the overall system dynamics can be adequately described by its phase of oscillation. Due to the general validity of these assumptions among many systems, the phase oscillator approach has been applied to a plethora of physical, chemical, and biological systems, such as synchronizing fireflies, frog choruses, or the crowd synchronization of pedestrians on London’s Millennium bridge (Ermentrout and Rinzel 1984; Strogatz et al. 2005; Ota et al. 2020), to name a few.

Entrainment under varying photoperiods

So far, we discussed general principles of entrainment under the assumption of symmetric zeitgeber cues with equal durations of day and night. Due to the tilt of the Earth’s rotation axis with respect to its orbit around the Sun, properties of zeitgeber signals such as the photoperiod of light–dark cycles depend on latitude and season. In Schmal et al. (2015), the concept of Arnold tongues (Fig. 1a, c was extended to account for photoperiodic entrainment, i.e., to zeitgeber cycles of varying daylengths. Since pure phase descriptions as given by Eqs. (8, 9, 10, 11) are unable to directly account for amplitude-dependent effects on entrainment and phase resetting (Lakin-Thomas et al. 1991; Ananthasubramaniam et al. 2020), we use a conceptual amplitude-phase oscillator model

$$\begin \frac&= \lambda r (A-r) \end$$

(12)

$$\begin \frac&= \frac \end$$

(13)

also known as Poincaré oscillator (Glass and Mackey 1988), instead. In radial coordinates as given by Eqs. (12 and 13), variables r and \(\phi\) denote the time-dependent (instantaneous) amplitude and phase of the internal clock, respectively, while parameters A, \(\tau\) and \(\lambda\) conveniently describe properties of the internal clock such as the steady-state amplitude, period, and radial relaxation rate which can differ and be related to specific organisms, tissues, or cell types. The resulting entrainment regions in the photoperiod and zeitgeber period parameter plane have their largest entrainment range at the equinox and taper toward the internal clocks free-running period \(\tau\) under constant darkness and constant light (Fig. 2a). The tilt of this Arnold Onion is given by Aschoff’s rule, i.e., the difference between \(\tau\) under constant darkness and constant light, with the internal period under constant light being typically shorter or longer compared to the period under constant darkness in day-active animals and plants or night-active animals, respectively (Aschoff 1960; Pittendrigh 1960). A complementary theoretical treatise to explain the emergence and properties of Arnold onions using a pure phase oscillator description as given by Hoveijn (2016) connects these results with the mathematical approach of the previous section.

Fig. 2

Arnold onions capture essential features of seasonal entrainment. a Entrainment regions adopt an onion-shaped geometry in the photoperiod-zeitgeber period parameter plane. The tilt of the Arnold onion can be explained by Aschoff’s rule, i.e., the difference between the internal free-running period under constant darkness (photoperiod of \(0\%\)) and constant light (photoperiod of \(100\%\)), depicted by vertical dashed lines. Phases of entrainment \(\psi\) are color-coded within the region of entrainment. b Experimentally obtained entrainment phase \(\psi\) in dependence of the zeitgeber period T for the golden hamster Mesocricetus auratus subject to light–dark cycles with equinoctial (blue) and extremely short (orange) photoperiods. Data have been extracted from Fig. 3 of (Aschoff and Pohl 1978) via the WebPlotDigitizer software (Rohatgi 2022)

Again, such a straightforward conceptual model is able to explain a variety of experimental results on photoperiodic entrainment. For the model assumptions underlying Fig. 2a, the 180° rule holds true within the entrainment range at a given fixed photoperiod, analogous to the observation for pure oscillators in Fig. 1. From this, it follows that the phase variability with respect to changes in zeitgeber period T is lowest under equinoctial conditions and increases with increasing or decreasing photoperiods as experimentally observed for golden hamster entrained to light–dark cycles of different photoperiods; compare Fig. 2b. Another prediction from Fig. 2a is that a large tilt of the Arnold onion as given by Aschoff’s rule can lead to a situation where the internal clock might be able to entrain to zeitgeber signals of short but not of long photoperiods or vice versa. This phenomenon can explain why the drinking behavior of squirrel monkeys (Saimiri sciureus) synchronizes to 24 h light–dark schedules of extremely short photoperiods but not to those longer than 21 h (Schmal et al. 2015; Sulzman et al. 1982).

While Fig. 2a shows entrainment ranges and phases for zeitgeber signals with a square-wave-like waveform as used in the laboratory, entrainment to light–dark cycles similar to those observed under natural conditions have been studied in (Schmal et al. 2020).

Intrinsic oscillator properties affect seasonal entrainment

Entrainment characteristics of the circadian system do not rely only on properties of the zeitgeber signal and the internal period \(\tau\) as discussed in the previous paragraphs but also on other intrinsic properties of the circadian clock such as the amplitude, radial relaxation rate, waveform, or twist (i.e., the dependence of the internal period on amplitude). Modeling approaches have been used to show that increasing amplitudes and radial relaxation rates make an oscillator more resistant toward entrainment in comparison to clocks with relatively small amplitudes and relaxation rates (Lakin-Thomas et al. 1991; Abraham et al. 2010). The finding that collective amplitudes and relaxation rates increase due to resonance effects in ensembles of interacting clocks (Abraham et al. 2010; Bordyugov et al. 2011; Schmal et al. 2018) has been used to explain why strongly coupled systems, such as the mammalian core pacemaker, the suprachiasmatic nucleus (SCN), have a narrow entrainment range, while putatively weakly coupled systems like lung or heart tissue rather behave like a weak clock and entrain to more extreme zeitgeber periods (Abraham et al. 2010). This interpretation is further strengthened by the fact that pharmacological decoupling of SCN neurons by MDL or TTX leads to a better entrainability of cultured SCN slices subject to temperature cycles (Abraham et al. 2010) as well as the observation that a faster recovery from jet-lag is observed for mice lacking receptors for the coupling agent arginine vasopressin (AVP) (Yamaguchi et al. 2013). Along these lines, it has been proposed that genetic redundancy within the molecular regulatory network underlying the mammalian circadian rhythm generation strengthens the clock and, thus, leads to narrow entrainment ranges (Erzberger et al. 2013).

Bifurcations affect seasonal entrainment

Bifurcations are defined by qualitative changes of systems dynamics due to variations of an internal or external parameters. Many of such qualitative changes in the systems dynamics upon parameter variations have been described for circadian clocks of different organisms. For example in mammals, the dissociation of a single activity band into two bands, termed splitting or frequency doubling, has been observed as a response to changes in zeitgeber properties such as an increasing light intensity under constant conditions (Pittendrigh and Daan 1976b). A transition from self-sustained to damped oscillations has been reported for circadian KaiC rhythms in cyanobacteria after reducing the ambient temperature below 18.6 °C (Murayama et al. 2017).

Fig. 3

Intrinsic oscillator properties govern seasonal entrainment characteristics. a The Goodwin oscillator is considered a blueprint for models of molecular negative feedback loops. We assume that the (square-wave) zeitgeber signal affects the negative feedback loop by an additive term to the \(X_1\)-variable. b Example oscillations under free-running conditions for a parameter set that leads to self-sustained oscillations. Same parameters as those underlying Fig. 8 of (Ananthasubramaniam et al. 2020) have been used. c For an increasing constant zeitgeber strength (e.g., constant light of increasing intensity), the system eventually changes its qualitative behavior through a Hopf bifurcation and looses its ability to self-sustain the oscillations. d For large zeitgeber intensities, the bifurcation shown in panel (c) translates into a broad entrainment range under long photoperiods similar to the behavior of damped oscillators. Panels c, d are adapted from Fig. S5a and Fig. 8a of (Ananthasubramaniam et al. 2020), respectively (under CC BY 4.0 license)

Such changes in qualitative behavior of circadian oscillator properties will have an impact on the entrainment characteristics as recently reported in a mathematical study using the Goodwin model (Ananthasubramaniam et al. 2020). The Goodwin oscillator is a generic model of a delayed negative feedback loop where the final product \(X_3\) of a three-component activatory chain inhibits the production of the first component \(X_1\) (Goodwin 1965); see Fig. 3a, b. It fulfills all necessary requirements to produce self-sustained oscillations, such as a negative feedback, non-linearity, as well as delay, and has a long tradition of being applied in modeling circadian clocks (Ruoff et al. 2001; Gonze and Ruoff 2021). Assuming that light enters the model as an additive term to the \(X_1\) variable, one observes that increasing constant light finally drives the system to a dampened regime through a Hopf bifurcation; see Fig. 3c. Damped oscillators can be entrained much more easily to rhythmic zeitgeber signals in comparison to a strong self-sustained clock (Bain et al. 2004; Gonze et al. 2005). Thus, forcing the Goodwin oscillator with a zeitgeber intensity that corresponds to a value that would drive the system to a damped regime under constant conditions (e.g., gray dotted line in Fig. 3c) leads to broadening of the entrainment region under long photoperiods (Fig. 3d). A similar behavior that relies on the bifurcation structure of the underlying molecular feedback loop has been observed in Neurospora crassa subject to temperature entrainment, where increasingly damped oscillations for decreasing temperatures (Liu et al. 1997b) translate into an experimentally observed broader entrainment range under short thermoperiods (Burt et al.