Are circadian amplitudes and periods correlated? A new twist in the story

Conservative linear and non-linear oscillators: Harmonic vs. Duffing oscillators

In classical mechanics, a mass-spring harmonic oscillator is a system of mass $m$ that, when displaced from its equilibrium position, experiences a restoring force $F$ proportional to the displacement $x$ , namely $F=-\mathit{kx}$ , where $k$ is the spring constant. When $F$ is the only restoring force, the system undergoes harmonic motion (sinusoidal oscillations) around its equilibrium point. In the absence of damping terms, the harmonic motion can be mathematically described by the following linear second order ordinary differential equation (ODE)

The solution to this differential equation is given by the function $x\left(t\right)=Acos\left(\frac{2\pi }{\tau }+\varphi \right)$ , ^{29 , 30} where $A$ represents the amplitude; $\varphi$ , the phase; and $\tau$ represents the period of the motion, $\tau =2\pi \sqrt{\frac{m}{k}}$ . Thus the oscillatory period is determined only by the mass $m$ and the spring constant $k$ . The amplitude $A$ , on the other hand, is determined solely by the starting conditions (by both initial displacement $x$ and velocity $v=\dot{x}$ ).

Introducing a non-linear term in the restoring force such that $F=-\mathit{kx}-\beta x^3$ allows the conversion of the simple harmonic oscillator into a Duffing oscillator. ³¹ Depending on the sign of $\beta$ , the coefficient that determines the strength of the non-linear term, the spring is termed hard or soft oscillator. The equation of the Duffing oscillator, in the absence of damping terms, reads

Due to the non-linear term introduced in equation 2, it is helpful to write this system in a form that can be easily treated via numerical integration. Considering the following change of variable $v=\dot{x}$ , the equation can be reformulated as a system of two first order ODEs:

Goodwin-like models

The three-variable Goodwin model is a minimal model based on a single negative feedback loop that can describe the emergence of oscillations in simple biochemical systems. All synthesis and degradation terms are linear, with the exception of the repression that $z$ exerts on $x$ which is modeled with a sigmoidal Hill curve. The equations that describe the dynamics read

The Goodwin model, however, requires a very large Hill exponent ( $n>8$ ) to produce self-sustained oscillations, ³² which biologists and modelers have often considered unrealistic. Gonze ³³ and others have shown that, by introducing additional nonlinearities in the system, the need for such high value of the Hill exponent can be reduced. In contrast to the linear degradation of variables of the original Goodwin model, the Gonze model ³³ describes degradation processes with Michaelis-Menten kinetics as follows

To mimic clock heterogeneity and evaluate the amplitude-period correlations among ensembles of Gonze or Goodwin oscillators, degradation rates were varied around $\pm 10\%$ their default value (Table 1).

Almeida model

The Almeida model ³⁶ is a protein model of the mammalian clockwork that includes 7 core clock proteins along with the PER:CRY complex and the regulation (activatory or repressive) that these exert on DNA binding sites known as clock-controlled elements to regulate circadian gene expression. The regulation at these clock-controlled elements, namely E-boxes, D-boxes and RORE elements, is described by the following terms:

The system of ODEs that describes the dynamics of the clock proteins in the Almeida model reads

All parameter descriptions along with their default values are given in Table 2.

To analyze the twist effects that arise from a population of heterogeneous clocks (i.e., parametric twist), all 18 model parameters were randomly varied around $\pm 20\%$ their default value (Table 2) one at a time. Since this model can lead to period-doubling effects upon changes of certain parameters, ³⁷ those oscillations whose period change resulted in oscillations with period-doubling were removed from the analysis. Moreover, if changing the default parameter in the ensemble resulted in a range of ratio of amplitude variation relative to the default amplitude $<0.1$ , then we considered that ensemble to have no twist for that particular control parameter. Periods and amplitudes were determined from the peaks and troughs of oscillations using the continuation software XPP-AUTO.

Poincaré model

The intrinsic dynamical properties from single oscillators and their interaction to external stimuli can be very conveniently described by means of a Poincaré model. ³⁸ We here propose a modification of its generic formulation that explicitly takes into account twist effects through the phase space twist parameter $\mathrm{ϵ}$ . The modified Poincaré model with twist reads

The first equation describes the rate of change of the radial coordinate $r\left(t\right)$ (i.e., the time-dependent distance from the origin), whereas the second equation determines the rate of change of the angular coordinate $\varphi \left(t\right)$ , where $\omega =\frac{2\pi }{\tau }$ . The parameters $\tau$ , $A$ , $\lambda$ and $\mathrm{ϵ}$ denote the free-running period (in units of time), amplitude (arbitrary units), amplitude relaxation rate (in units of time^-1) and phase space twist of the oscillator (in units of time^-1), respectively. In the absence of twist, namely $\mathrm{ϵ}=0$ , the phase changes constantly along the limit cycle at a rate $\frac{2\pi }{\tau }$ , independently of the radius. In the case of $\mathrm{ϵ}\ne 0$ , the phase changes at a constant rate only when $r=A$ ; if any perturbation is to modify $r$ such that $r\ne A$ , then the phase change will be accelerated or decelerated depending on the sign of $\mathrm{ϵ}$ and on whether $r>A$ or $r<A$ . The model parameters, unless otherwise specified in the figures or captions, are the following: $A=1$ a.u., $\lambda =0.05{\mathrm{h}}^{-1}$ , $\tau =24$ h and $\mathrm{ϵ}$ values of 0 or $\pm 0.1{\mathrm{h}}^{-1}$ .

To study the effects of twist on the oscillator’s response to pulse-like perturbations $\text{pert}\left(t\right)$ , periodic zeitgeber input $Z\left(t\right)$ or mean-field coupling $M$ , the individual Poincaré oscillators $i$ were converted into Cartesian coordinates and the respective terms were added in the equations of the $x_i$ variable as follows:

Computer simulations and data analysis

All numerical simulations were performed and analyzed in Python with the numpy, scipy, pandas and astropy libraries. The function odeint from scipy was used to numerically solve all ordinary differential equations. Bifurcation analyses were computed in XPP-AUTO ⁴⁰ using the parameters $N_{\mathit{tst}}=150,N_{\mathit{max}}=20000,{\mathit{Ds}}_{\mathit{min}}=0.0001$ and ${\mathit{Ds}}_{\mathit{max}}=0.0002$ .

Throughout our analyses, periods were determined by (i) normalizing the solutions to their mean, (ii) centering them around 0 (by subtracting one unit from the normalized solution), and (iii) by then computing the zeroes of the normalized rhythms. The period was defined as the distance between two consecutive zeros with a negative slope. Amplitudes were determined as the average peak-to-trough distance of the last (normalized) oscillations after removing transients.

Nonlinearities can result in amplitude-period correlations across oscillators

The simple harmonic oscillator (Figure 2A) is a classical model of a system that oscillates with a restoring force proportional to its displacement, namely $F=-\mathit{kx}$ . ⁴¹ The ordinary differential equation that describes the motion of a mass attached on a spring is linear (equation 1 in Materials and Methods) and the solution can be found analytically. ^{29 , 30} The period is determined by the size of the mass $m$ and the force constant $k$ (see Materials and Methods), while the amplitude and phase are determined by the starting position and the velocity. Thus, an ensemble of harmonic oscillators with different initial conditions will produce results that differ in amplitudes but whose periods are the same (Figure 2B, C). When plotting amplitudes against periods, no correlation or twist is observed: the period of a simple harmonic oscillator is independent of its amplitude (Figure 2D).

Figure 2. Nonlinearities introduce amplitude-period correlations in oscillator models.

The mass-spring harmonic oscillator (A) is a linear oscillator that, in the absence of friction, produces persistent undamped oscillations (B, in time series and C, in phase space). Different starting conditions produce oscillations that differ in amplitudes but whose period is independent of amplitude (D). The harmonic oscillator is converted to a Duffing oscillator (E, I) by introducing nonlinearities in the restoring force of the spring (the non-linear behavior is represented with deformations of the spring). The new restoring force results in slight changes in the initial conditions producing oscillations whose amplitude and period change (F, G, J, K) and become co-dependent (H, L) and where the correlations depend on the sign of the non-linear term. The terms soft and hard oscillators refer to, by convention, those with positive and negative twist, respectively. Period values in (F–H, J–L) are color-coded.

The simple harmonic oscillator can be converted into a Duffing oscillator ³¹ by including a cubic nonlinearity in its equation. In the Duffing oscillator (equations 2, 3 in Materials and Methods), the restoring force is no longer linear (see deformed springs in Figure 2) but instead described by $F=-\mathit{kx}-\beta x^3$ , where $\beta$ represents the coefficient of non-linear elasticity. These classical conservative oscillators instead are known to have an amplitude-dependent period. A network of Duffing clocks with different starting conditions will produce oscillations of different amplitudes as well as periods. Duffing oscillators with a negative cubic term have been termed soft oscillators and display periods that grow with amplitudes (Figure 2E–H), similar to Kepler’s Third Law of planetary motion, ⁴² where planets with larger distances to the Sun run at slower periods than those that are closer. On the other hand, Duffing oscillators with a positive $\beta$ term are known as hard oscillators and show negative twist (amplitude-period correlations) (Figure 2I–L). ^{43 , 44}

In short, nonlinearities in oscillator models can introduce twist effects among ensembles of oscillators with slight differences in their properties (initial conditions, parameters, etc.). Thus, models for the circadian clock, which are based on nonlinearities, are expected to show amplitude-period correlations.

Oscillator heterogeneity produces parametric twist effects in limit cycle clock models

Most circadian clock models generate stable limit cycle oscillations. Limit cycles are isolated closed periodic orbits with a given amplitude and period, where neighboring trajectories (e.g. perturbations applied to the cycle) spiral either towards or out of the limit cycle. ^{38 , 45} Stable limit cycles are examples of attractors: they imply self-sustained oscillations. The closed trajectory describes the perfect periodic behavior of the system, and any small perturbation from this trajectory causes the system to return to it, to be attracted back to it.

Limit cycles are inherently non-linear phenomena and they cannot occur in a linear system (i.e., a system in the form of $\dot{\overrightarrow{x}}=A\overrightarrow{x}$ , like the harmonic oscillator). From the previous section we have learned that non-linear terms can introduce amplitude-period co-dependencies. In this section we show that kinetic limit cycle models of the circadian clock also show twist effects. Instead of studying the amplitude-period correlation of an oscillator model with fixed parameters and changing initial conditions as in Figure 2 (as all initial conditions would be attracted to the same limit cycle), we study here the correlations that arise among different uncoupled oscillators due to oscillator heterogeneity (i.e., differences in biochemical parameters), as found experimentally. ^{14 , 26} We refer to amplitude-period correlations that arise due to oscillator heterogeneity as parametric twist.

Single negative feedback loop models

The Goodwin model is a simple kinetic oscillator model ⁴⁶ that is based on a delayed negative feedback loop, where the final product of a 3-step chain of reactions inhibits the production of the first component (equation 4 in Materials and Methods). In the context of circadian rhythms, ^{34 , 35} the model can be interpreted as a clock activator $x$ that produces a clock protein $y$ that, in turn, activates a transcriptional inhibitor $z$ that represses $x$ (Figure 3A). The Goodwin model has been extensively studied and fine-tuned by Gonze and others ³³ to study fundamental properties of circadian clocks ^{47 – 50} or synchronization and entrainment. ^{33 , 51 , 52} The Gonze model ³³ (equation 5) includes additional nonlinearities, where the degradation of all 3 variables is modeled with non-linear Michaelis Menten kinetics, to reduce the need of very large Hill exponents ( $n>8$ ), required in the original Goodwin model to generate self-sustained oscillations, ³² that have been questioned to be biologically meaningful. These Michaelian degradation processes can be interpreted as positive feedback loops which aid in the generation of oscillations ⁵³ (Figure 3B).

Figure 3. Parametric twist in an ensemble of 100 uncoupled Goodwin-like models: correlations are model- and parameter-dependent.

Scheme of the Goodwin (A) and Gonze (B) oscillator models. Changes in degradation parameters result in twist effects when varied randomly around $\pm$ 10% their default parameter value, mimicking oscillator heterogeneity of the Goodwin (C) and Gonze (D) models. Simultaneous co-variation of degradation parameters results in parametric twist: changes of $k_2$ and $k_4$ produce significant positive amplitude-period correlations in both models (E, F); co-variations of $k_6$ and $k_4$ result in amplitude-period correlations that are not significant in the Goodwin model (G), but significant positive parametric twist in the Gonze model (H). Shown in all panels is Spearman’s $R$ and the significance of the correlation (n. s.: not significant; ***: p value $<$ 0.001). Default model parameters are summarized in Table 1. Relative amplitudes are computed as the peak-to-trough distance of the $x$ variable. Oscillators from the ensemble whose amplitude was $<0.1$ have been removed from the plots; the total number of oscillators ( $n$ ) is indicated in the plots.

To study whether twist effects are present in an ensemble of 100 uncoupled Goodwin and Gonze oscillators with different parameter values, we randomly varied the degradation parameters of $x$ , $y$ or $z$ ( $k_2$ , $k_4$ and $k_6$ , respectively) in each oscillator around $\pm 10\%$ their default parameter value (given in Table 1) and analyzed the resulting amplitude-period correlation. The resulting parametric twist depends on the model and on the parameter being changed: in the Goodwin model, variations in the degradation rate of the transcriptional activator ( $k_2$ ) or of the clock protein ( $k_4$ ) result in positive twist effects (i.e., soft twist-control), while changes in the transcriptional repressor’s degradation rate ( $k_6$ ) result in a hard twist-control, i.e., negative amplitude-period correlation (Figure 3C). In the Gonze model, $k_2$ , $k_4$ and $k_6$ are all soft twist-control parameters, as individual changes of any of them all produce positive parametric twist effects (Figure 3D).

To mimic cell-to-cell variability in a more realistic manner, we introduced heterogeneity by changing combinations of the degradation parameters simultaneously around $\pm 10$ % their default parameter value and analyzing the resulting periods and amplitudes. We observed that the overall twist behavior depends on the particular influence that each parameter has individually. Random co-variations of $k_2$ and $k_4$ produce significant positive parametric twist effects in ensembles of both Goodwin (Figure 3E) or Gonze clocks (Figure 3F), consistent with the positive correlations when either of the parameters is changed individually (Figure 3C, D). Nevertheless, when $k_6$ , which has a negative twist effect in the Goodwin model, is changed at the same time as $k_4$ (or $k_2$ , data not shown) the correlation is not significant (Figure 3G). Random co-variation of $k_4$ and $k_6$ in an ensemble of Gonze clocks results in significant positive parametric twist effects (Figure 3H).

Changes in parameter values can result in significant alterations to a system’s long-term behavior, which can include differences in the number of steady-states, limit cycles, or their stability properties. Such qualitative changes in non-linear dynamics are known as bifurcations, with the corresponding parameter values at which they occur being referred to as bifurcation points. In oscillatory systems, Hopf bifurcations are an important type of bifurcation point. They occur when a limit cycle arises from a stable steady-state that loses its stability. A 1-dimensional bifurcation diagram illustrates how changes in a mathematical model’s control parameter affect its final states, for example the period or amplitude of oscillations. The Hopf bubble refers to the region in parameter space where the limit cycle exists, and it is commonly represented with the peaks and troughs of a measured variable in the $y$ axis, with the control parameter plotted on the $x$ axis (Figure 4A). The term “bubble” is used because of the shape of the curve, that resembles a bubble that grows or shrinks as the parameter is changed. Such bifurcation analyses can be used to predict the type of twist that a system shows: the amplitude will increase with parameter changes if the default parameter value is close to the opening of the Hopf bifurcation, or will decrease if the value is near the closing of the bubble.

Figure 4. The type of parametric twist depends on the parameter location within the Hopf bubble.

Bifurcation analysis of Gonze oscillators with changing $k_4$ , simulating heterogeneous oscillators: (A) Hopf bubble representing the peaks (maxima) and troughs (minima) of variable $x$ as a function of changes in the degradation rate of $y$ ( $k_4$ ); (B) monotonic period decrease for increasing $k_4$ . Represented in (C) is the combination of (A) and (B), namely how the difference between the maxima and minima of $x$ changes with the period: a non-monotonic behavior of parametric twist is observed, which depends on the values of $k_4$ . Note that, since $k_4$ results in a decrease in period length, the $k_4$ dimension is included implicitly in (C) and $k_4$ increases from the right to the left of the plot. Stars indicate the default parameter value, period or amplitude.

We observed that Gonze oscillators display self-sustained oscillations for values of $k_4$ between 0.2 and 0.43 (Figure 4A) and that, for increasing $k_4$ , periods decrease monotonically (Figure 4B). For oscillators with increasing $k_4$ , the parametric twist effects from the ensemble are first negative (amplitudes increase and periods decrease); however, when Gonze clocks have $k_4$ values that are part of the region of the bubble where amplitudes decrease, positive amplitude-period correlations appear in the ensemble (Figure 4C). In summary, the type of parametric twist depends on both the model and the parameter being studied but also on where the parameter is located within the Hopf bubble. Other more complex kinetic models of the mammalian circadian clockwork have shown that changes in some of the degradation parameters produce non-monotonic period changes, ⁵⁴ and thus the twist picture is expected to become even more complex.

Model with synergistic feedback loops

The Almeida model ³⁶ is a more detailed model of the core clock network in mammals that includes seven core clock proteins that exert their regulation at E-boxes, D-boxes and ROR binding elements (RORE) through multiple positive and negative feedback loops (Figure 5A, equation 6 in Materials and Methods). To study parametric twist in an ensemble of uncoupled Almeida oscillators, we randomly changed all 18 parameters individually around $\pm 20\%$ their default value (Table 2 in Materials and Methods) and computed the amplitude-period correlations. In this case, we calculated the ratio of amplitude variation after the parameter change relative to the default amplitude.

Figure 5. The type of parametric twist in models of the circadian clockwork with synergistic feedback loops depends on the parameter being studied and on the variable that is measured.

(A) Scheme of the Almeida model, ³⁶ illustrating activation (green arrows) and inhibition (red arrows) of different clock-controlled elements (at E-boxes, D-boxes and ROR elements) by the respective core clock proteins. Parametric twist was studied by modeling a heterogeneous uncoupled ensemble where model parameters were randomly and individually varied around $\pm 20\%$ of their default value (shown in Table 2) and assessing the period-amplitude correlation of the ensemble. (B–D) Almeida oscillators ( $n=40$ ) with heterogeneous values in the activation rate of D-boxes ( $V_D$ ): the twist is negative from the perspective of BMAL1, positive from the perspective of PER:CRY, but the amplitude of PER is not greatly affected by changing $V_D$ . (E–G) Almeida oscillators ( $n=40$ ) with heterogeneous values in the degradation rate of PER ( ${\gamma }_P$ ) show positive parametric twist effects for BMAL1, PER and PER:CRY. The twist in panels (B, E) is evaluated by comparing the relative amplitude of the respective protein/protein complex (computed as the average peak-to-trough distance) after parameter change to the default amplitude. Shown in (C, D, F, G) are representative oscillations with long and short periods: dashed lines represent rhythms of smallest amplitudes.

We found, interestingly, that the overall twist effects depend not only on the parameter being studied, but also on the variable which is measured. For example, changes in the rate of D-box activation parameter $V_D$ result in negative twist effects for BMAL1 (i.e., lower amplitude BMAL1 rhythms run slower than oscillators with higher amplitude BMAL1 rhythms), positive parametric twist for the PER:CRY complex but almost no parametric twist from the perspective of PER (Figure 5B–D). Changes in PER degradation ${\gamma }_P$ produce positive parametric twist for BMAL1, PER and PER:CRY but of different magnitudes (Figure 5E–G). Supplementary Table S1 (in ⁵⁷ ) summarizes the parametric twist effects for the additional parameters from the Almeida model, highlighting the complexity that arises with synergies of feedback loops and bringing us again to the question of defining what the relevant amplitude of a complex oscillator is.

Phase space twist in single oscillators: On amplitude-phase models, isochrones and perturbed trajectories

Up until now, our results have focused on studying amplitude-correlations among ensembles of self-sustained oscillators with differences in their intrinsic properties (e.g. biochemical parameters) in the absence of external cues. However, when a clock is exposed to external stimuli, its amplitude and period undergo adaptation. We refer to these amplitude-period correlations in individual clocks as they return to their steady-state oscillation after a stimulus as phase space twist. In contrast to parametric twist, which emerged as a result of a heterogenous population of oscillators, here we focus on individual clocks.

To further explore the correlation and interdependence between the frequency of an oscillator and its amplitude upon an external stimulus, more generalized models can be of use. The Poincaré oscillator model (equation 7 in Materials and Methods) is a simple conceptual oscillator model with only two variables, amplitude and phase, that has been widely used in chronobiology research. ^{20 , 22 , 24 , 38 , 39 , 55} This amplitude-phase model, regardless of molecular details, can capture the dynamics of an oscillating system and what happens when perturbations push the system away from the limit cycle. When a pulse is applied and an oscillating system is ‘kicked out’ from the limit cycle, the perturbed trajectory is attracted back to it at a rate $\lambda$ . Within the limit cycle, the dynamics are strictly periodic: if one takes a point in phase space and observes where the system returns to after exactly one period, the answer is trivial: to exactly the same spot (Figure 6, red dots). Nevertheless, outside the limit cycle, during the transient relaxation time (i.e., time between the perturbation and the moment that the trajectory reaches the limit cycle), the time between two consecutive peaks might be shorter or longer than the period of the limit cycle orbit. We will come back to this point in the next section.

Figure 6. Schematic of an isochrone in a Poincaré oscillator with (A) no twist, (B) positive twist or (C) negative twist.

To understand the concept of isochrones, a simple experiment is performed: if one considers a point in phase space (shown in red) within a limit cycle (black) and observes where the system returns to after exactly one period, the answer is trivial: to the same spot. If, however, now a perturbation is applied (red arrow) such that the system starts in a section of the state space which is outside of the stable periodic orbit, and one maps the points that the system ‘leaves as a footprint’ after exactly one period as it relaxes back to the limit cycle (blue trajectories), these points (shown in yellow) will outline an isochrone (grey line). The twist parameter $\mathrm{ϵ}$ affects the curvature of the isochrone, with implications in the response of the oscillating system to the perturbation. The oscillator with positive twist (B) arrives at a later phase than that with no twist (A), resulting in a phase delay with respect the clock with no twist, whereas the clock with negative twist (C) arrives to the limit cycle at an earlier phase (i.e., advanced with respect to the clock with no twist).

Arthur T. Winfree introduced a practical term, isochrones, to conceptualize timing relations in oscillators perturbed off their attracting cycles, ³⁹ which becomes important in physiological applications because often, biological oscillators are not on their attracting limit cycles. To understand what isochrones are, Winfree proposes in ³⁹ a simple experiment where a pulse-like perturbation applied in a system produces a ‘bump’ in the limit cycle. As the perturbation relaxes and spirals towards the attracting limit cycle, one records the position of the oscillator in time steps which equal the period of the unperturbed limit cycle. The sequence of points left as a footprint will converge to a fixed point on the cycle and will outline an isochrone (Figure 6). Isochrones are related to oscillator twist, and whereas in a Poincaré oscillator with no twist, the isochrones are straight and radial (Figure 6A), when twist is present, the isochrones become bent or skewed (Figure 6B, C).

In an individual oscillator with positive twist, a perturbation that increases the oscillator’s instantaneous amplitude will relax back and intersect the isochrone at an angle $<{360}^{\circ }$ (Figure 6B) than an oscillator with no twist in the time of one period, whereas an oscillator with negative twist will cover an angle $>{360}^{\circ }$ in the time of one period (Figure 6C). Thus, isochrones and consequently twist illustrate how a perturbation away from the limit cycle is decelerated or accelerated during the course of relaxation.

Mathematically, the ‘skewness’ of isochrones in an amplitude-phase oscillator is directly related to the twist parameter ϵ, that adds a radius dependency on the phase dynamics. The phase changes at a constant rate $\omega =\frac{2\pi }{\tau }$ in the limit cycle, but in conditions outside the limit cycle, the phase is made to be modulated by the radius through the twist parameter ϵ until $r=A$ (equation 7 in Materials and Methods). This phase space twist parameter ϵ now represents the amplitude-period correlations (acceleration/deceleration) as an individual oscillator returns to its steady-state oscillation.

Phase space twist in single oscillators affects their interaction with the environment

We have seen how, in an individual oscillator, the twist parameter $\mathrm{ϵ}$ acts to slow down (in case of positive twist) or to speed up (in case of negative twist) trajectories further away from the limit cycle compared to those closer to it. As a result, $\mathrm{ϵ}$ has a direct consequence on how the individual oscillator responds to pulse-like perturbations or periodic zeitgebers coming from the environment.

To analyze how phase space twist affects the response of an oscillator to a zeitgeber pulse, we applied a pulse-like perturbation to individual Poincaré oscillators with different phase space twist values $\mathrm{ϵ}$ . The pulse was applied at CT3 and was made to increase the instantaneous amplitude (such that $r>A$ ). If the oscillator has no twist ( $\mathrm{ϵ}=0$ , Figure 7A), the isochrones are straight and radial, but for a soft or a hard oscillator with positive or negative $\mathrm{ϵ}$ values, respectively, isochrones get skewed (Figure 7B, C, also Figure 6). In the case of a soft oscillator with positive $\mathrm{ϵ}$ (Figure 7B), the pulse at CT3 gets decelerated during the course of its relaxation, and as a consequence, the oscillator arrives back to the limit cycle at an later phase than the oscillator with $\mathrm{ϵ}=0$ (Figure 7D, E). The same can also be inferred mathematically: for this particular perturbation where $r>A$ , $A-r<0$ holds and hence the phase velocity outside the limit cycle is smaller ( $\dot{\varphi }<\omega$ ) than at the limit cycle for the oscillator with positive twist, since $\dot{\varphi }=\omega +\mathrm{ϵ}\left(A-r\right)$ (equation 7). The opposite holds for the oscillator with negative twist: this clock arrives at an earlier phase to the limit cycle (Figure 7C, F).

The positive and negative amplitude-period correlations characteristic of twist effects are evident in how, upon the pulse-like perturbation, the peak-to-peak distance changes (compared to the limit cycle period) as the perturbation returns to the periodic orbit. For a clock with $\mathrm{ϵ}=0$ , the peak-to-peak distance outside the limit cycle after a perturbation coincides with the 24 h peak-to-peak distance within the periodic orbit (Figure 7G). In the case of the Poincaré oscillator with positive twist, the peak-to-peak distance of the oscillator after the perturbation that increases the instantaneous amplitude is longer than the period of the steady-state rhythm (24 h), since $\dot{\varphi }<\omega$ . As the perturbed amplitude decreases and is attracted back to the stable cycle, the peak-to-peak distance approaches 24 hours (Figure 7H, consistent with the formulation ‘positive’ twist, which implies positive amplitude-period correlations). Conversely, the clock with negative twist exhibits a shorter peak-to-peak distance than the 24 h rhythm during the relaxation time $\left(\dot{\varphi }>\omega \right)$ , which subsequently increases as the system returns to the stable periodic orbit (Figure 7I). When considering the findings from the last two paragraphs as a whole, it becomes clear how the extent of the phase shift between a perturbed and an unperturbed oscillator (red and black dashed lines in Figure 7D–F, respectively) depends on $\mathrm{ϵ}$ and thus the shape of the phase response curve (PRC) and magnitude of the phase shift depend on this twist parameter (Figure 7J–L).

We then evaluated how twist affects the response of an individual oscillator to a periodic zeitgeber input. We observed how the range of entrainment becomes larger with larger values of absolute value $\mathrm{ϵ}$ , as seen by the wider Arnold tongues in Figure 7M–O. It is widely known that, when the frequency of an applied periodic force is equal or close to the natural frequency of the system on which it acts, the amplitude of the oscillator on which the driving force (zeitgeber) acts on increases due to resonance effects. Interestingly, also the resonance curve was affected by twist: for the oscillator with no twist, the maximum amplitude occurred as expected for a zeitgeber period of 24 h, matching the intrinsic oscillator’s frequency (Figure 7P). In oscillators with positive or negative twist, however, the maximum amplitude after zeitgeber forcing increases and decreases with zeitgeber period, respectively (Figure 7Q, R), resulting in skewed resonance curves. ⁵⁶ Interestingly, we observed that twist also affects the entrainment of an individual clock, and oscillators with high absolute values of twist cannot entrain to a periodic signal (Supplementary Figure S1 in ⁵⁷ ), a phenomenon that has previously been described as shear-induced chaos ^{27 , 58 , 59} and that arises because the isochrones become so skewed, that the trajectory of the relaxation gets stretched and folded.

Figure 7. Phase space twist in a single oscillator and its effect on responses to pulse-like perturbations and to entrainment to periodic zeitgebers.

(A–C) Poincaré oscillator models for different twist values $\mathrm{ϵ}$ , shown in phase space. The limit cycle is shown in black; perturbed trajectories are shown in red; isochrones are depicted in grey. (D–F) Time series corresponding to panels (A–C): the unperturbed limit cycle oscillations are shown with a dashed black line; perturbed trajectories are shown in red. (G–I) Peak-to-trough distance as a function of the peak-to-peak distance from the perturbed time series from (D–F) relaxing back to the limit cycle (the implicit time dimension is indicated in the panels). (J–L) Phase response curves: the shape of the PRC and the extent of the phase shifts (in hours) depends on the twist $\mathrm{ϵ}$ . (M–O) Arnold tongues illustrating entrainment ranges of an individual oscillator with different phase space twist values to a periodic sinusoidal zeitgeber input of different periods $T$ . The amplitude of the entrained clock is color-coded, with yellow colors corresponding to larger amplitudes. (P–R) The amplitude of an individual oscillator driven by a sinusoidal zeitgeber also depends on the twist of the individual oscillator. Shown are the resonance curves as a response to a zeitgeber input of strength $F=0.05$ and varying periods $T$ . Amplitudes are defined as the average peak-to-trough distances of the entrained signal.

Relaxation rate affects the response of oscillators to perturbations

Not only the twist parameter $\mathrm{ϵ}$ , but also the amplitude relaxation rate $\lambda$ affects the response of oscillators to perturbations and consequently the entrainment range. We have seen how, in the Poincaré model (equation 7), the twist parameter $\mathrm{ϵ}$ dictates how much trajectories are ‘slowed down’ or ‘sped up’ dependent on their distance from the limit cycle. The amplitude relaxation rate $\lambda$ describes the rate of attraction back to the limit cycle, which is independent of $\mathrm{ϵ}$ . Together, these parameters both dictate how ‘skewed’ the isochrones are in phase space (see analytical expression in Materials and Methods, equation 11). For Poincaré oscillators with twist, the isochrones become more radial and straight as the amplitude relaxation rate increases (Figure 8A–C), with the implications that this has on PRCs and entrainment that have been mentioned before. A more ‘plastic’ clock (with a lower $\lambda$ value) responds to an external pulse with larger phase shifts than a more rigid clock (with higher $\lambda$ ), resulting in phase response curves of larger amplitude (Figure 8D–F). This is also intuitive, as it is clear that for larger values of $\lambda$ , the twist has shorter effective time to act upon perturbed trajectories because the perturbation spends less time outside the limit cycle.

Figure 8. Role of amplitude relaxation rate on the skewing of isochrones and consequently on the response of oscillators to perturbations.

(A–C) Poincaré oscillator models with different amplitude relaxation values $\lambda$ , but otherwise identical (free running period $\tau =24$ h, amplitude $A=1$ and twist $\mathrm{ϵ}=0.05{\mathrm{h}}^{-1}$ ), shown in phase space. The limit cycle is shown in black; perturbed trajectories are shown in red; isochrones are depicted in grey. (D–F) Phase response curves: the shape of the PRC and the extent of the phase shifts depends on the relaxation rate $\lambda$ . See Materials and Methods for details on the analytical derivation of the equation of isochrones and roles of $\lambda$ and $\mathrm{ϵ}$ .

Single oscillator phase space twist affects coupled networks

Biological oscillators are rarely alone and uncoupled. Coupled oscillators, instead, are at the heart of a wide spectrum of living things: pacemaker cells in the heart, ³⁸ insulin- and glucagon-secreting cells in the pancreas ⁶⁰ or neural networks in the brain and spinal cord that control rhythmic behavior as breathing, running and chewing. ⁶¹ A number of studies have shown that networks of coupled oscillators behave in a fundamentally different way than ‘plain’ uncoupled oscillators. ^{33 , 62 , 63} This suggests that single oscillator twist might also affect how coupling synchronizes a network of oscillators.

To analyze the effect of twist on coupled oscillators, we study systems of Poincaré oscillators coupled through a mean-field (equation 8). We start by coupling two identical oscillators and analyzing the effect that different twist values have on the behavior of the coupled network. The numerical simulations show that increasing coupling strengths affect the period of the coupled system: oscillators with positive twist show longer periods as the coupling strength increases, whereas oscillators with negative twist show period-shortening (Figure 9A and Supplementary Figure S2A–C in ⁵⁷ ). Coupling results in an increase in amplitude but, for our default relaxation rate value ( $\lambda =0.05{\mathrm{h}}^{-1}$ ), no significant differences across oscillators with different twist values were found (Supplementary Figure S2 in ⁵⁷ ). Increasing relaxation rates (oscillators that are attracted faster back to the limit cycle upon a perturbation) nevertheless resulted in less coupling-induced period or amplitude changes (Supplementary Figure S3 in ⁵⁷ ).

Figure 9. Phase space twist $\mathrm{ϵ}$ affects the response of a network to coupling.

(A) Mean-field coupling of two identical Poincaré oscillators (sharing the same value of twist ϵ) and the effect on the period of the coupled system: whereas coupling among oscillators with positive twist results in longer periods of the coupled network, mean-field coupling among oscillators with negative twist results in the network running at a faster pace (period shortening). (B) Mean-field coupling can synchronize a network of oscillators with different twist values. Shown are three oscillators that only differ in their twist parameters (but else identical: amplitude, amplitude relaxation rate and period are the same for the three oscillators). All oscillations overlap during the first part of the time series, since twist has no effect within the limit cycle. At a certain point (indicated by the red arrow), mean-field coupling is turned on and it is observed how the oscillators respond differently to that mean-field coupling, ending up with different individual amplitudes as well as relative phases to the mean-field (depicted in dashed black line).

We then analyzed what happens when oscillators with different phase space twist values are coupled through their mean-field. For this purpose, we took a system of 3 oscillators (with twist values $\mathrm{ϵ}=0$ , $\mathrm{ϵ}<0$ and $\mathrm{ϵ}>0$ ) and turned on mean-field coupling after a certain time. We observed that all three oscillators synchronize to the mean-field for low values of the absolute value of twist (Figure 9B) although, due to the coupling-induced period differences, the relative phases of the individual oscillators to the mean-field depend on the specific twist values of the individual clocks: those with negative twist oscillate ‘ahead’ of the mean-field (with a phase advance) but clocks with positive twist cycle with a phase delay (Figure 9B). Very large absolute values of twist, interestingly, produce again complex chaotic dynamics (Supplementary Figure S4 in ⁵⁷ ). The oscillator with a large value of ϵ does not synchronize to the mean-field and, as a result of the inter-oscillator coupling, this desynchronized oscillator ‘pulls’ to desynchronization the other two oscillators which otherwise would have remained in sync.