What this is
- This research explores the relationships between circadian clock parameters: period, amplitude, and phase.
- It introduces the concepts of and to describe amplitude-period correlations.
- These correlations can influence how oscillators respond to environmental stimuli, affecting their entrainment and recovery.
Essence
- Circadian clock parameters—period and amplitude—are interdependent, influencing how clocks adapt to environmental cues. This study defines and , which describe these correlations and their implications for oscillator behavior.
Key takeaways
- describes amplitude-period correlations in heterogeneous populations of circadian clocks, showing that different parameters can result in either positive or negative correlations.
- refers to the amplitude-period correlations of an individual oscillator responding to external stimuli, affecting its recovery and adaptation to perturbations.
- The study reveals that the twist parameter influences how oscillators respond to environmental signals, which can affect their synchronization and overall behavior in coupled networks.
Caveats
- The research is theoretical, relying on models that may not fully capture the complexities of real biological systems.
- The definitions of twist effects may vary across different oscillator models, complicating the generalization of findings.
Definitions
- parametric twist: Amplitude-period correlations observed in populations of heterogeneous oscillators.
- phase space twist: Amplitude-period correlations within an individual oscillator's time series in response to external stimuli.
AI simplified
Introduction
Oscillations are happening all around us, from the vibrating atoms constituting matter to the beating of the animal heart or to circadian clocks present in all kingdoms of life. Circadian clocks are autonomous clocks that tick in the absence of external timing cues with a period of about 24 h and regulate our behavior, physiology and metabolism. A fundamental property of circadian clocks is that their phase and periodicity can be adjusted to external timing signals (zeitgebers) in a process known as entrainment. It is believed that natural selection has acted on the phase relationship between biological rhythms and the environmental cycle, and thus, this phase of entrainment is of central importance for the fitness of the organism, allowing it to anticipate to changes in the external world. 1
Three key properties of a circadian rhythm are its period, amplitude and phase. Chronobiological studies have usually focused on period because there are established tools that allow their direct measurement, including running wheels for mice or race tubes for fungi. Phase, which refers to the position of a point on the oscillation cycle relative to a reference point, can be measured similarly. Circadian rhythms can be entrained to various zeitgeber periods T as reviewed in, 2 but under the natural conditions of T=24 h, variations of intrinsic periods lead to different phases of entrainment that are the basis of chronotypes: faster running clocks (shorter endogenous periods) lead to early phases (‘morning larks’) and slower clocks (longer periods) correspond to later phases (‘night owls’). 3–7
Measuring amplitudes, however, is less straightforward. Some studies have considered activity recordings 8,9 or conidiation in race tubes, 10 but one might argue that these measures do not represent the clockwork’s amplitude, as they reflect outputs of the circadian system. Other studies have quantified gene expression profiles after careful normalization, 11 but from the ∼20 core clock genes that constitute the mammalian circadian oscillator, 12 it is not immediately evident which gene or protein is best representing the core clock amplitude. Actually, reporter signals monitoring expression of different clock genes and proteins have been used to quantify amplitudes. 4,13,14 Others have approached the amplitude challenge indirectly by measuring the response of an oscillator to zeitgeber pulses. 15–18 While small-amplitude clocks exhibit larger pulse-induced phase shifts and are easier to phase-reset, larger amplitude rhythms display smaller phase shifts, 17–22 with consequences in the size of the phase response curve 5 or in jet lag duration. 23 Amplitudes, together with periods, also govern entrainment 23,24 and seasonality. 22,23,25 There have been various theoretical and experimental studies showing, for example, how clocks with larger amplitudes display narrower ranges of entrainment than rhythms of lower amplitude, 20 and how the phase of entrainment is modulated by oscillator amplitude. 22,24
Taken together, these observations indicate, firstly, that the phase of entrainment is correlated with the intrinsic period, and secondly, that both phase of entrainment and phase changes in response to perturbations also correlate with oscillator amplitude. This leads to the question of whether amplitudes and periods are also co-modulated and what insights these interdependencies provide about the underlying oscillator. These questions are the focus of this paper. Do faster-running clocks have larger or smaller amplitudes than slower clocks? What are the implications? Experimental observations have provided evidence for both: in a human osteosarcoma cell line in culture, clocks with longer periods display larger amplitudes; 26 but in cells from the choroid plexus, the major producer of cerebrospinal fluid of the central nervous system, clocks with shorter periods are associated with larger amplitudes 14 (scheme in Figure 1). This dependence between periods and amplitudes is what we here refer to as twist, also known as shear in the literature. 27,28 By convention, negative twist describes oscillators in which amplitude increases are accompanied by a decreasing period (also termed hard oscillators) and vice versa for positive twist (soft oscillators).
In this paper we provide definitions for two important concepts in the context of limit cycle oscillations: parametric twist and phase space twist. Parametric twist refers to the amplitude-period correlations observed when analyzing a population of heterogeneous clocks; phase space twist describes the amplitude-period correlations within the time series of an individual oscillator, as influenced by its interaction with the environment.
We show, firstly, that nonlinearities can introduce amplitude-period correlations in simple oscillator models. Moreover, clock models of different complexity can reproduce the experimentally observed positive 26 and negative 14 twist effects, with the type of correlation depending on the model, on parameters, as well as on the variable being measured, illustrating the complexity in defining circadian amplitudes. Lastly, we show how twist effects can speed up or slow down zeitgeber-induced amplitude changes in a simple oscillator model. This helps the clock phase adapt and modulate entrainment of oscillators to a periodic signal or their response to coupling. Our results support the use of oscillator theory as a framework to understand emerging properties of circadian clocks with different twist. Moreover, they provide insights into how temporal or spatial phase patterning might arise in coupled networks, as well as how amplitude changes could help in stabilizing the circadian period in the face of temperature changes. Although we focus on circadian clocks, the presented theory can be applied to any other oscillatory system such as cardiac rhythms, flashing fireflies or voice production.
![Schematic of experimental observations of circadian amplitude-period correlations. (A) Positive correlations (i.e., positive twist) have been observed in U-2 OS cells kept in culture; (B) negative correlations have been observed in cells from the choroid plexus in the mouse brain. [26] [14]](https://europepmc.org/articles/PMC10526121/bin/f1000research-12-162736-g0000.jpg.jpg "Click to view full size")
Schematic of experimental observations of circadian amplitude-period correlations. (A) Positive correlations (i.e., positive twist) have been observed in U-2 OS cells kept in culture; (B) negative correlations have been observed in cells from the choroid plexus in the mouse brain. [26] [14]
Materials and Methods
Conservative linear and non-linear oscillators: Harmonic vs. Duffing oscillators
In classical mechanics, a mass-spring harmonic oscillator is a system of massthat, when displaced from its equilibrium position, experiences a restoring forceproportional to the displacement, namely, whereis the spring constant. Whenis 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) m F x F kx = − k F d 2 x d t 2 k m + = . x 0 (1)
The solution to this differential equation is given by the function xt=Acos2πτ+ϕ, 29,30 where A represents the amplitude; ϕ, the phase; and τ represents the period of the motion, τ=2πmk. 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=x˙).
Introducing a non-linear term in the restoring force such that F=−kx−βx3 allows the conversion of the simple harmonic oscillator into a Duffing oscillator. 31 Depending on the sign of β, the coefficient that determines the strength of the non-linear term, the spring is termed hard or soft oscillator. The nonlinearity introduces a dependency between the amplitude and period, meaning that varying initial conditions will produce oscillations with different, yet related, periods and amplitudes. The equation of the Duffing oscillator, in the absence of damping terms, reads d2xdt2+kmx+βmx3=0.(2)
Due to the non-linear term introduced in, it is helpful to write this system in a form that can be easily treated via numerical integration. Considering the following change of variable, the equation can be reformulated as a system of two first order ODEs: equation 2 v = x ˙ d x d t = v d v d t k m β m = − − . x x 3 (3)
Goodwin-like models
The Goodwin model is a minimal model that has been widely used to describe the emergence of oscillations in simple biochemical systems. It is based on a single negative feedback loop, where the final product of a 3-step chain of reactions inhibits the production of the first component. In the context of circadian rhythms, the model can be interpreted as a clock gene mRNA x that gets translated into a clock protein y, that then activates the repressor z, which ultimately inhibits the transcription of x. All synthesis and degradation terms are linear, with the exception of the repression that z exerts on x which is modeled with a Hill curve. The equations that describe the dynamics read dxdt=k1K1nK1n+zn−k2xdydt=k3x−k4ydzdt=k5y−k6z,(4)where k1, k3 and k5 represent the rates of synthesis of x, y and z, respectively; k2, k4 and k6, the degradation rates; and n, the Hill exponent.
The Goodwin model, however, requires a very large Hill exponent ( n>8) to produce self-sustained oscillations, 32 which biologists and modelers have often considered unrealistic. Gonze 33 Kurosawa 34 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 33 describes degradation processes with Michaelis-Menten kinetics as follows dxdt=k1K1nK1n+zn−k2xK2+xdydt=k3x−k4yK4+ydzdt=k5y−k6zK6+z.(5)
To mimic clock heterogeneity and evaluate the amplitude-period correlations among ensembles of Gonze or Goodwin oscillators, degradation rates were varied around ±10% their default value ( Table 1).
| Parameter | Biological interpretation | Value in Goodwin model[35],[36] | Value in Gonze model[33] |
|---|---|---|---|
| k 1 | transcription rate of core clock gene x | 1 nM | 0.7 nM h − 1 |
| k 2 | degradation rate of x | 0.2 h − 1 | 0.35 nM h − 1 |
| k 3 | translation rate of core clock protein y | 1 h − 1 | 0.7 h − 1 |
| k 4 | degradation rate of y | 0.15 h − 1 * | 0.35 nM h − 1 |
| k 5 | nuclear import rate of the repressor z | 1 h − 1 | 0.7 h − 1 |
| k 6 | degradation rate of z | 0.1 h − 1 | 0.35 nM h − 1 |
| K 1 | Michaelis constant of-mediated repression z | 1 nM | 1 nM |
| K 2 | Michaelis constant ofdegradation x | - | 1 nM |
| K 4 | Michaelis constant ofdegradation y | - | 1 nM |
| K 6 | Michaelis constant ofdegradation z | - | 1 nM |
| n | Hill exponent of cooperative inhibition ofon z x | 9.5 * | 4 |
Almeida model
The Almeida model 37 is a protein model of the mammalian clockwork that includes 7 core clock proteins along with the PER:CRY complex. Unlike prior more extensive models, 38,39 it does not consider any post-translational modifications of clock proteins or any nuclear import or export processes, but it takes into account the regulation (activatory or repressive) that these clock proteins 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: Ebox=VEBMAL1BMAL1+kE+kErBMAL1CRYRORE=VRRORkR+RORkRr2kRr2+REV2Dbox=VDDBPDBP+kDkDrkDr+E4BP4.
The system of ODEs that describes the dynamics of the clock proteins in the Almeida model reads dBMAL 1 d t = − RORE BMAL PERCRY γ BP 1 dROR d t = + − E box γ Ror RORE ROR dREV d t = + − 2 E box D box γ Rev REV dDBP d t = − E box γ Db DBP dE BP 4 4 d t = − 2 4 4 R box γ E 4 E BP dCRY d t = + − + − E box R box γ PC γ CP γ C 2 PERCRY PERCRY CRY dPER d t = + − + − E box D box γ PC γ CP γ P PERCRY PERCRY PER dPERCRY d t = − − . γ PC γ CP γ BP PER CRY PERCRY BMAL PERCRY 1 (6)
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 ±20% their default value ( Table 2) one at a time. Since this model can lead to period-doubling effects upon changes of certain parameters, 40 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.
| Parameter | Biological interpretation | Value | Parameter | Biological interpretation | Value |
|---|---|---|---|---|---|
| V R | rate of RORE activation | 44.4 h − 1 | γ Ror | ROR degradation rate | 2.55 h − 1 |
| k R | strength of RORE activation | 3.54 | γ Rev | REV degradation rate | 0.241 h − 1 |
| k Rr | strength of RORE inhibition | 80.1 | γ P | PER degradation rate | 0.844 h − 1 |
| V E | rate of E box activation | 30.3 h − 1 | γ C | CRY degradation rate | 2.34 h − 1 |
| k E | strength of E box activation | 214 | γ DB | DBP degradation rate | 0.156 h − 1 |
| k Er | strength of E box inhibition | 1.24 | γ E 4 | E4BP4 degradation rate | 0.295 h − 1 |
| V D | rate of D box activation | 202 h − 1 | γ PC | PER:CRY degradation rate | 0.191 h − 1 |
| k D | strength of D box activation | 5.32 | γ CP | PER:CRY formation rate | 0.141 h − 1 |
| k Dr | strength of D box inhibition | 94.7 | γ BP | BMAL1 nuclear export rate | 2.58 h − 1 |
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. 41 We here propose a modification of its generic formulation that can take into account phase space twist effects through the twist parameter ϵ, explicitly introduced in the equations. The modified Poincaré model with twist reads drdt=λrA−rdϕdt=ω+ϵA−r.(7)
The first equation describes the rate of change of the radial coordinate rt (i.e., the time-dependent distance from the origin), whereas the second equation determines the rate of change of the angular coordinate ϕt, where ω=2πτ. The parameters τ, A, λ and ϵ denote the free-running period (in units of time), amplitude (arbitrary units), amplitude relaxation rate (in units of time -1) and twist parameter of the oscillator (in units of time -1), respectively. In the absence of twist, namely ϵ=0, the phase changes constantly along the limit cycle at a rate 2πτ, independently of the radius. In the case of ϵ≠0, the phase changes at a constant rate only when r=A; if any perturbation is to modify r such that r≠A, then the phase change will be accelerated or decelerated depending on the sign of ϵ and on whether r>A or r<A. This can generate amplitude-period correlations within the time series of the oscillator during its relaxation time, i.e., what we denote as phase space twist. The model parameters, unless otherwise specified in the figures or captions, are the following: A=1 a.u., λ=0.05h−1, τ=24 h and ϵ values of 0 or ±0.1h−1.
To study phase space twist and the effects of differentvalues on the oscillator’s response to pulse-like perturbations, periodic zeitgeber inputor mean-field coupling, the individual Poincaré oscillatorswere converted into Cartesian coordinates and the respective terms were added in the equations of thevariable as follows:where. The zeitgeberis given by:whererepresents the zeitgeber period andthe strength (amplitude) of the zeitgeber input. The mean-fieldcoupling is given by:whererepresents the coupling strength andis the number of oscillators in the coupled ensemble. Lastly, the square-like perturbationis defined as follows:whereis the strength of the perturbation andis the time at which the perturbation starts. The perturbation lasts 1 h and is set toa.u. in all our simulations. ϵ pert t Z t M i x i r i = x i 2 y i 2 + Z t T F Z M K N pert t F P t start F P = 0.7 d x i d t = − + + + λ Z M x i y i A − r i ω ϵ + A − r i t t pert d y i d t = + , λ y i x i A − r i ω ϵ + A − r i (8) Z t 2 π T π 2 t + = cos , F Z M = , K N ∑ i = 1 N x i t pert t 0 , if t start t start ≤ ≤ + t 1 F P otherwise , =
It is important to differentiate between the twist parameter, denoted as, and the concept of phase space twist, as they are related but not equivalent terms. Whileserves as one of the characteristic parameters in the Poincaré model, phase space twist represents the observable outcome of this parameter's influence on how the oscillator's amplitude and period respond to environmental perturbations. Phase space twist is measurable and influenced by, but also by the amplitude relaxation rate λ. This will become clear in a further section. ϵ ϵ ϵ
Analytical derivation of isochrones
Arthur T. Winfree introduced the concept of isochrones as any set of dynamical states which oscillate with the same phase when they reach the limit cycle at time t→∞ . 42 That asymptotic phase at t→∞ is what Winfree termed latent phase Φ. Since the dynamical flow of the Poincaré model has polar symmetry, the isochrones must also have polar symmetry such that Φ=gϕr=!ϕ−fr.(9)
By definition, the latent phase velocitynecessarily increases at units of the angular velocityas the oscillator follows its kinetic equation: Φ ˙ ω = 2 π τ dΦ d t = ! ω (10)
Combining the two previous equations and with the use of the chain rule, we can calculate the radius dependency of the latent phase: Φ dΦ d t d ϕ d t d f r d r d r d t = 9 = 10 − ω
The termsandare defined in the Poincaré model (), so that the equation above can be rewritten as: d ϕ d t d r d t equation 7 ω ϵ λr ω + − A r − A r − d f r d r = 10
Solving for: d f r d r d f r d r ϵ A r − λr A r − ϵ λr = =
Next, the solution forcan be found through integration: f r f d r r A r = = ln − ln ϵ λ 1 r ϵ λ ϵ λ ∫ r A = r
If(like in all our simulations), the last term can be neglected such that A = 1 f r r = ln ϵ λ
Finally, the solution ofcan be inserted into end up with the equation for isochrones as a function of the radius: f r equation 9 Φ ϕ r = − ln ϵ λ (11)
Isochrones are thus loci of polar coordinatesin phase space with the same latent phase.shows how the shape of the isochrone (and, consequently, phase space twist effects) for the Poincaré model independs on the ratio of twistto relaxation rate. Specifically, highervalues (in modulus) and lower relaxation rateslead to more curved isochrones and greater phase space twist effects. In the Results section, we further elaborate on how to interpret the isochrones graphically. ϕ r Φ ϵ λ ϵ λ Equation 11 equation 7
To plot the isochrones, we simply reorderto plotas a function ofat isochrones with fixed values of latent phases: equation 11 ϕ r Φ = 0 π 4 π 2 3 π 4 π − 3 π 4 − π 2 − π 4 ϕ r Φ = ln − ϵ λ
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 43 using the parameters Ntst=150,Nmax=20000,Dsmin=0.0001 and Dsmax=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.
Results
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=−kx . 44 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).
The simple harmonic oscillator can be converted into a Duffing oscillator 31 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=−kx−βx3, where β 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, 45 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 β term are known as hard oscillators and show negative twist (amplitude-period correlations) ( Figure 2I–L). 46,47
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.
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.
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. 41,48 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 x→˙=Ax→, 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. But twist in limit cycle oscillators has to be studied differently, as different inidial conditions all return to the same cycle. Instead of studying the amplitude-period correlation of an oscillator model with fixed parameters and changing initial conditions as in Figure 2, 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 become evident in a population of heterogeneous oscillators as parametric twist.
Single negative feedback loop models
The Goodwin model is a simple kinetic oscillator model 49 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, 35,36 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, 33 Kurosawa 34 and others to study fundamental properties of circadian clocks 50–53 or synchronization and entrainment. 33,54,55 The Gonze model 33 ( 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, 32 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 56 ( Figure 3B).
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 ( k2, k4 and k6, respectively) in each oscillator around ±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 ( k2) or of the clock protein ( k4) result in positive twist effects (i.e., soft twist-control), while changes in the transcriptional repressor’s degradation rate ( k6) result in a hard twist-control, i.e., negative amplitude-period correlation ( Figure 3C). In the Gonze model, k2, k4 and k6 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 ±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 k2 and k4 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 k6, which has a negative twist effect in the Goodwin model, is changed at the same time as k4 (or k2, data not shown) the correlation is not significant ( Figure 3G). Random co-variation of k4 and k6 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.
We observed that Gonze oscillators display self-sustained oscillations for values of k4 between 0.2 and 0.43 ( Figure 4A) and that, for increasing k4, periods decrease monotonically ( Figure 4B). For oscillators with increasing k4, the parametric twist effects from the ensemble are first negative (amplitudes increase and periods decrease); however, when Gonze clocks have k4 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, 38 and thus the twist picture is expected to become even more complex.
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 around10% their default parameter value, mimicking oscillator heterogeneity of the Goodwin (C) and Gonze (D) models. Simultaneous co-variation of representative degradation parameters results in parametric twist: changes ofandproduce significant positive amplitude-period correlations in both models (E, F); co-variations ofandresult 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’sand the significance of the correlation (n. s.: not significant; ***: p value0.001). Default model parameters are summarized in. Relative amplitudes are computed as the peak-to-trough distance of thevariable, normalized to its mean. Oscillators from the ensemble whose amplitude washave been removed from the plots; the total number of oscillators () is indicated in the plots. ± k 2 k 4 k 6 k 4 R < x < 0.1 n Table 1
The type of parametric twist depends on the parameter location within the Hopf bubble. Bifurcation analysis of Gonze oscillators with changing, simulating heterogeneous oscillators: (A) Hopf bubble representing the peaks (maxima) and troughs (minima) of variable(normalized to its mean) as a function of changes in the degradation rate of(); (B) monotonic period decrease for increasing. Represented in (C) is the combination of (A) and (B), namely how the difference between the maxima and minima ofchanges with the period: a non-monotonic behavior of parametric twist is observed, which depends on the values of. Note that, sinceresults in a decrease in period length, thedimension is included implicitly in (C) andincreases from the right to the left of the plot. Stars indicate the default parameter value, period or amplitude. All panels display the numerical results from bifurcation analyses of Gonze oscillators obtained with XPP-AUTO. k 4 x y k 4 k 4 x k 4 k 4 k 4 k 4
Model with synergistic feedback loops
The Almeida model 37 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). It is a simple model in the sense that, unlike more complex models like the Relógio model, 38 it neglects nuclear import/export processes or post-translational modifications of the clock proteins; but complex enough to capture the synergistic positive and negative regulations and feedback processes that the 7 core clock proteins exert on the DNA clock-controlled elements, unlike older models. 39 To study parametric twist in an ensemble of uncoupled Almeida oscillators, we randomly changed all 18 parameters individually around ±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.
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 VD 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 γP produce positive parametric twist for BMAL1, PER and PER:CRY but of different magnitudes ( Figure 5E–G). Supplementary Table S1 (in 57) 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.
![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 aroundof their default value (shown in) and assessing the period-amplitude correlation of the ensemble. (B–D) Almeida oscillators () with heterogeneous values in the activation rate of D-boxes (): 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. (E–G) Almeida oscillators () with heterogeneous values in the degradation rate of PER () 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 normalized to its mean) 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. [37] ± % 20 n = 40 V D V D n = 40 γ P Table 2](https://europepmc.org/articles/PMC10526121/bin/f1000research-12-162736-g0004.jpg.jpg "Click to view full size")
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 aroundof their default value (shown in) and assessing the period-amplitude correlation of the ensemble. (B–D) Almeida oscillators () with heterogeneous values in the activation rate of D-boxes (): 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. (E–G) Almeida oscillators () with heterogeneous values in the degradation rate of PER () 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 normalized to its mean) 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. [37] ± % 20 n = 40 V D V D n = 40 γ P Table 2
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, another form of amplitude-period correlations can occur within an oscillator's time series, because 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.
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,41,42,57 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 λ. 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.
Arthur T. Winfree introduced a practical term, isochrones, to conceptualize timing relations in oscillators perturbed off their attracting cycles, 42 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 42 a simple experiment where a pulse-like perturbation applied in a system produces a ‘bump’ in the limit cycle. For stable limit cycles, the perturbation will relax and spiral towards the attracting cycle (blue spiral in Figure 6, until eventually the simulation returns to the limit cycle after enough cycles). If during this relaxation, 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 (yellow points in Figure 6) will converge to a fixed point on the cycle and will outline an isochrone ( Figure 6). All points in the same isochrone have the same latent phase Φ (see equation 9).
Mathematically, isochrones are related to the oscillator's twist parameter ϵ, and whereas in a Poincaré oscillator with no twist, the isochrones are straight and radial ( Figure 6A), when twist is present, the isochrones can become bent or skewed ( Figure 6B, C) because of the radius dependency that ϵ introduces on the phase dynamics. The phase changes at a constant rate ω=2πτ 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). 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∘ ( 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° in the time of one period ( Figure 6C). This parameter ϵ thus determines how perturbations away from the limit cycle are decelerated or accelerated during the course of relaxation. Ultimately, this influences the amplitude-period correlations within the oscillator's time series (i.e., phase space twist effects) as it returns to its steady-state oscillation following a perturbation.
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). The twist parameteraffects 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). ϵ
Phase space twist in single oscillators affects their interaction with the environment
We have seen how, in an individual oscillator, the twist parameter ϵ acts to slow down (in case of positive twist, Figure 6B) or to speed up (in case of negative twist, Figure 6C) trajectories further away from the limit cycle compared to those closer to it. As a result, ϵ has a direct consequence on how the individual oscillator responds to pulse-like perturbations or periodic zeitgebers coming from the environment.
To analyze how the acceleration or deceleration of perturbations induced by ϵ affects phase space twist and the response of an oscillator to a zeitgeber pulse, we applied a pulse-like perturbation to individual Poincaré oscillators with different twist values ϵ. The pulse was applied at CT3 and was made to increase the instantaneous amplitude (such that r>A). If the oscillator has no twist ( ϵ=0, Figure 7A), the isochrones are straight and radial, but for a soft or a hard oscillator with positive or negative ϵ values, respectively, isochrones get skewed ( Figure 7B, C, also Figure 6). In the case of a soft oscillator with positive ϵ ( 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 ϵ=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 ( ϕ˙<ω) than at the limit cycle for the oscillator with positive twist, since ϕ˙=ω+ϵA−r ( 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 phase space 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 ϵ=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 ϕ˙<ω. As the perturbed amplitude decreases and is attracted back to the stable cycle, the peak-to-peak distance approaches 24 hours. Consequently, a positive phase space twist correlation is observed ( Figure 7H). Conversely, the clock with negative twist exhibits a shorter peak-to-peak distance than the 24 h rhythm during the relaxation time ϕ˙>ω, which subsequently increases as the system returns to the stable periodic orbit, resulting in negative phase space twist effects ( 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 ϵ and thus the shape of the phase response curve (PRC) and magnitude of the phase shift depend on this twist parameter ( Figure 7J–L) and on the amplitude-period correlations (i.e., phase space twist) that it induces.
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 ϵ, 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 oscillator's intrinsic period ( 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. 58 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 Ref. 59), a phenomenon that has previously been described as shear-induced chaos 27,60,61 and that arises because the isochrones become so skewed, that the trajectory of the relaxation gets stretched and folded.
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, 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) Representation of phase space twist effects: 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. (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. 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 strengthand varying periods. Amplitudes are defined as the average peak-to-trough distances of the entrained signal. ϵ ϵ T F = 0.05 T
Relaxation rate also modulates phase space twist effects
Not only the twist parameter ϵ, but also the amplitude relaxation rate λ 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 ϵ dictates how much trajectories are ‘slowed down’ or ‘sped up’ dependent on their distance from the limit cycle. The amplitude relaxation rate λ describes the rate of attraction back to the limit cycle, which is independent of ϵ. Together, these parameters both dictate how ‘skewed’ the isochrones are in phase space (see analytical expression in Materials and Methods, equation 11) and modulate phase space twist. For Poincaré oscillators with twist, the isochrones become more radial and straight as the amplitude relaxation rate increases ( Figure 8A–C). The correlations between peak-to-peak periods and peak-to-trough amplitudes within the oscillator's time series also change ( Figure 8D–F), with the implications that this has on PRCs and entrainment that have been mentioned before. A more ‘plastic’ clock (with a lower λ value) responds to an external pulse with larger phase shifts than a more rigid clock (with higher λ), resulting in phase response curves of larger amplitude ( Figure 8G–I). This is also intuitive, as it is clear that for larger values of λ, the twist ϵ has shorter effective time to act upon perturbed trajectories because the perturbation spends less time outside the limit cycle. As a consequence, phase space twist effects can become less evident, as shown in Figure 8D–F.
Role of amplitude relaxation rate on the skewing of isochrones, phase space twist and consequently on the response of oscillators to perturbations. (A–C) Poincaré oscillator models with different amplitude relaxation values, but otherwise identical (free running periodh, amplitudeand twist), 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 space twist effects: amplitude relaxation rate λ modulates the correlations between period (defined as peak-to-peak distance) and amplitude (peak-to-trough distance). (G–I) Phase response curves: the shape of the PRC and the extent of the phase shifts depends on the relaxation rate. See Materials and Methods for details on the analytical derivation of the equation of isochrones and roles ofand. λ τ = 24 A = 1 ϵ = 0.05 h − 1 λ λ ϵ
Phase space twist in single oscillators 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, 41 insulin- and glucagon-secreting cells in the pancreas 62 or neural networks in the brain and spinal cord that control rhythmic behavior as breathing, running and chewing. 63 A number of studies have shown that networks of coupled oscillators behave in a fundamentally different way than ‘plain’ uncoupled oscillators. 33,64,65 Moreover, the findings in Figure 7 highlight that individual oscillators respond differently to environmental perturbations depending on the presence or absence of twist. This implies that the impact of single oscillator twist goes beyond the individual behavior, extending to how coupling synchronizes an entire 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 Ref. 59). Coupling results in an increase in amplitude due to resonance but, for our default relaxation rate value ( λ=0.05h−1), no significant differences across oscillators with different twist values were found (Supplementary Figure S2 in Ref. 59). 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 Ref. 59).
We then analyzed what happens when oscillators with different twist values are coupled through their mean-field. For this purpose, we took a system of 3 oscillators (with twist values ϵ=0, ϵ<0 and ϵ>0), we let them free run in the absence of coupling, and turned on mean-field coupling after a certain time (red arrow in Figure 9B). We observed that all three oscillators synchronize to the mean-field for low values of the absolute value of twist ( Figure 9B). Nevertheless, due to the period differences induced by coupling ( Figure 9A), 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 Ref. 59). 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.
Phase space twistaffects 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). ϵ
Discussion
This study aimed at characterizing period-amplitude correlations across circadian clock models of different complexity. Body clocks have to cope with cellular heterogeneity, what results in cellular clocks being heterogeneous across networks and tissues. We addressed whether the amplitude-period correlations that have been observed experimentally can be explained through heterogeneity in the cellular clocks, and what design principles are needed to produce such twist effects that we term parametric twist. Moreover, clocks live in constantly changing environments, that requires them to get adapted in the face of external changes. With the concept of phase space twist we refer to the co-modulations of period and amplitude that an individual oscillator experiences within its time series when it encounters an external stimulus. Ultimately, these twist effects tune the oscillator’s response to coupling, entrainment and pulse-like perturbations. We also retrieve the ‘old’ terminology of hard versus soft oscillators to refer to oscillations with negative and positive amplitude-period correlations, respectively.
The concept twist was historically coined on conservative, non-limit cycle oscillations. These type of oscillations, like the Duffing oscillator, 31 have an amplitude defined by the initial conditions, and hence any change in initial conditions results in a new cycle. In contrast, limit cycles oscillate with a characteristic period and amplitude regardless of the initial conditions. For this reason, twist cannot be studied in non-limit cycle oscillators in the same way it can be done for limit cycle oscillators. This is the reason why, in the context of the Duffing oscillators, twist is not termed neither parametric nor phase space. We have coined the terms parametric twist and phase space twist as the amplitude-period correlations that one can study in limit cycle oscillators.
Parametric twist effects become evident when analyzing populations of cells and require nonlinearities in oscillator models. The feedback loops needed to generate oscillations, which are commonly modeled with non-linear terms, 33,37,38,50 result in variations of parameters producing oscillations of different but correlated amplitudes and periods. We found that the type of parametric twist (positive or negative) depends on a number of factors: Firstly, on the biochemical parameter being affected, since different model parameters control the oscillation properties (amplitude and period) differently ( Figures 3 and 5). Second, on the region within the Hopf bubble where the clock’s parameter set is at: as seen in Figure 4, small changes in the default parameter value can lead to either an increase or decrease in amplitude, depending on the specific value of the default parameter. Lastly, the type of parametric twist also depends on the type of model and variable being measured: simple models with single negative feedback loops show the same parametric twist effects for all variables because the parameter of interest has the same effect on the amplitude of all variables. In complex models with synergies of loops, a change in one parameter might increase the amplitude of one variable but decrease the amplitude of a second variable ( Figure 5, also Ref. 66 and Supplementary Table S1 in Ref. 59). Additional clock models have found different types of parametric twist: for instance, Goldbeter found that, in a model of the Drosophila clock, the rate of PER degradation led to amplitude increasing and period decreasing 67 (i.e., negative parametric twist), unlike in the Almeida model ( Figure 5E). This highlights the challenge of defining circadian amplitude (and twist effects) to find a metric for the whole oscillating system. Our findings explain why some experimental studies have found positive twist, 26 others negative, 14 whereas some other works found very little correlations. 68
Bokka et al.66 recently conducted a comprehensive analysis of parametric twist in 11 ODE-based oscillator models, including the Goodwin model. In their analysis of the Goodwin model, they found that changes in k2 resulted in amplitude increases with a constant period (no twist). Changes in k4 led to increased amplitude and decreased period (negative twist), while changes in k6 resulted in decreases in both amplitude and period (positive twist). 66 However, our bifurcation analyses revealed different correlations. We observed positive twist when modifying k2 or k4, but a negative correlation with changes in k6 ( Figure 3). Such discrepancies may arise from different variables used for amplitude determination, from varying parameter sets, or from different amplitude definitions (we use relative amplitude, whereas Bokka et al. define it in absolute terms, i.e., as the distance from peak to trough). There have been several studies that have addressed the question of defining amplitude in circadian studies. One approach has been to propose a metric in which amplitude is defined as the geometric mean of all biomolecular species in the model of interest. 66,69 This metric, however, poses difficulties for experimental application. Defining which are the essential oscillatory species in a system is challenging, and even if agreed upon, measuring all biomolecular species within one experiment is often unfeasible.
The populations of oscillators considered in our and other 66 analyses to study parametric twist are uncoupled. But coupling across network of oscillators can narrow the distribution of oscillator periods, 64 thereby also reducing the extent of parametric twist effects. However, the decrease in parametric twist effects upon coupling should be reversible; if inter-oscillator coupling is disrupted, the original amplitude-period correlations are expected to return.
It is important to point out that our study has focused on kinetic oscillator models that use Hill functions to model transcriptional repression. 33,37,49 However, different repression mechanisms exist across different species, and such mechanisms might be modelled with different terms. For example, in Drosophila70 and in mammals, 71–75 transcriptional inhibition involves the formation of a 1:1 stoichiometric complex between activators and inhibitors – a process known as protein sequestration-based repression. 76,77 Here, repressors bind tightly to activators, creating an inactive complex and effectively 'sequestering' the activator. If the ratio greatly differs from 1:1, oscillations quickly dampen out. 76,77 The description of such protein sequestration processes has been described with terms different from Hill functions. 76,77 In contrast, species like Neurospora show a lower than 1:1 stoichiometry between activators and repressors in nuclei. 78 In this fungus, self-sustained oscillations require numerous phosphorylations in the activator. 78,79 As the number of mutated phosphorylation sites increases, the circadian rhythms become weaker and finally arrhythmic. 78 Such multistep cooperative reactions can be described with Hill curves. 80,81
We also addressed the concept of phase space twist. It is important to clarify that the parameter ϵ that we introduced in the modified Poincaré oscillator is referred to as 'twist' parameter throughout our manuscript (in line with prior studies 14,41,42). But phase space twist, strictly speaking, refers to the observable phenomenon: the amplitude period correlations of a perturbed oscillator as it returns to its steady state rhythm. In the Poincaré oscillator, phase space twist effects are a consequence of the ϵ parameter, but other models might also display such correlations without explicit parameters. Zeitgeber pulses, recurring zeitgeber inputs or coupling can all be regarded as ‘perturbations’, as any of these inputs modifies the natural clock’s limit cycle in phase space.
A feature of phase space twist is that it characterizes the adaptation of an oscillator to a perturbation and thus has significant implications for phase response curves (PRCs), entrainment and coupling. This new concept is therefore particularly interesting because, although it does not give direct insights into what the mechanistic or kinetic details of a system are, it does provide predictions on the dynamical properties of oscillators.
We start by discussing the role of phase space twist in PRCs. The twist parameter ϵ induces amplitude-period correlations within the oscillator's relaxation time ( Figure 7G–I) thus affecting the extent of the phase shift after a zeitgeber pulse ( Figure 7D–F): increasing the twist parameter in absolute value results in alterations in the phase response curves and consequently in the resetting properties of oscillators. In particular, increasing |ϵ| increases the amplitude of the PRC ( Figure 7J–L) until the PRC is converted from a type 1 PRC to a type 0 PRC. 19,82 Consistently, amplitudes 17,18,20,21,83 and periods 18 of oscillators have also been shown to modulate their resetting properties. In these examples, clocks with short periods 18 and small amplitudes (in mathematical models 17,18,20 or in experiments with reduced coupling 20,21 or mutations 83 that disrupt the normal rhythmicity) are easier to reset.
One can easily extrapolate these findings to responses to jet lag. If twist modulates the extent of a phase shift in response to a perturbation, it will also be critical in the adaptation to jet lag. This has indeed been found by Ananthasubramaniam et al.23 Interestingly, the authors found that the instantaneous amplitude effects induced by jet lag are also modulated by twist. In particular, positive twist aided recovery to jet lag in simple Goodwin-like models: reduced amplitudes were accompanied by faster clocks (i.e., shorter periods) upon a phase advance (e.g. when traveling eastwards), but larger amplitudes coincided with longer periods when the phase had to be delayed (e.g. when traveling westwards). Consistent with these theoretical observations are experiments performed in mice, where compromising coupling in the SCN (and thus decreasing amplitudes) reduce jet lag drastically, since resetting signals are much more efficient. 21
Transient amplitude-period correlations also affect the entrainment properties of an oscillator. Phase space twist modulates the range of zeitgeber periods to which the Poincaré oscillator can entrain ( ϵ ≠ 0 results in larger entrainment ranges, i.e., wider Arnold tongues, Figure 7M–O) and skews the resonance curves ( Figure 7P–R). Some studies have even found coexisting limit cycles for driven oscillators with twist. 31,84 Of note, however, is that the period of the Poincaré oscillator in all our simulations was set to 24 h. Other theoretical studies have shown that clocks with different intrinsic periods show different amplitude responses to entrainment, 23,85 thus combining phase space and parametric twist effects.
Single cells harbor self-sustained clocks, but they are coupled (see Ref. 86 for a review on coupling) to produce a coherent rhythm at the level of tissues and organisms. Coupled networks show different rhythmic properties than clocks in isolation: phases synchronize and ensemble amplitudes increase for over-critical values of coupling strengths. 64 Our simulations provide an additional level of regulation, by showing that phase space twist influences period length of individual clocks due to coupling. The presence of coupling-induced changes in the oscillation period or amplitude can provide insights into the underlying oscillator type and the presence of twist. For instance, the longer period that has been observed in dispersed (and presumably uncoupled) U-2 OS cells compared to high-density cultures, 87 but not in dispersed SCN neurons in culture, 88,89 could be explained by different implicit oscillator twist and amplitude relaxation rate values in different tissues. However, a study by Myung et al. suggests that the picture may be more intricate: The authors show how pharmacological or physical disruption of SCN coupling does not affect the period of Per2 rhythms, in contrast to Bmal1 oscillations, which display longer periods. 90
In addition, the twist-induced modulation of period also affects the phase relation with which the clock oscillates respect to the mean-field, consistent with chronotypes. It should be noted, however, that different phase relationships to the mean-field can also be obtained without twist, in networks of coupled oscillators with different intrinsic periods. 64 But here again, the slower running clocks will tend to be phase-delayed in comparison to the faster-running clocks which will be phase-advanced. This shows that twist plays an important role in how synchronization arises in network of coupled oscillators.
Our simulations assume that the coupling strength is constant across all oscillators. This assumption, however, might be questioned, considering that oscillators may 'communicate' with different strengths to each other, especially because individual clocks are spatially organized within tissues. For example, rhythms in cellular activity across SCN neurons do not occur simultaneously: depending on the developmental stage and the environmental conditions, the SCN can show complex spatio-temporal patterns such as phase waves or phase clustering. 64,90,97 Such phase organizations have been shown to depend on the previous light schedule that the SCN network was entrained to 98: while small phase differences of PER2::LUC rhythms are observed in cultured SCN slices after 12:12 light-dark schedules, entrainment to long days with 20 h of light and 4 h of darkness leads to phase clusters and region-specific phase differences of up to 12 h. 92,93 Moreover, it has been shown that long summertime days result in a reduced synchronization within the SCN, which might facilitate fast adaptations to a new light-dark cycle. 99 Experimental evidence suggests that these seasonal changes are regulated by coupling via neurotransmitters such as GABA 100,101 or VIP. 102 These results imply that adaptation to day length and thus to seasonal variations requires network re-organizations of the SCN. (For complete reviews on coupling in the SCN and seasonality, refer to Refs. 86, 98.) In light of these findings, a more plausible scenario may involve local coupling, where clocks couple to neighbor clocks with a strength proportional to the oscillator distance. Alternatively, differences in coupling may occur, with the same agent acting as a synchronizing or desynchronizing factor based on the network's state, as has been suggested for GABA. 93,100 Thus, twist might not only be critical in how temporal synchronization but also spatial patterning arises in coupled ensembles.
It is important to remark that also the relaxation rate λ (i.e., how rigid/plastic an oscillator is) modulates phase space twist effects ( Figure 8D–F) and affects the response of oscillators to perturbations ( Figure 8G–I), consistent with previous computational work. 20,57 It is in fact the ratio of ϵ to λ what determines the skewing of isochrones (see the analytical derivation in Materials and Methods) and the oscillator’s response to zeitgebers. Larger values of λ (more ‘rigid’ oscillators) imply that any perturbation is attracted back to the limit cycle at a higher rate and thus, perturbations ‘spend less time’ outside the limit cycle. Consequently, the twist parameter ϵ has a shorter effective time to act on the perturbed trajectory, and the isochrones become more straight and radial ( Figure 8A–C). Relaxation rate also has implications on coupling and entrainment. Rigid oscillators with ϵ≠0 display less coupling-induced amplitude expansions (Supplementary Figure S3B in Ref. 59, also Ref. 20) and less period variations (Supplementary Figure S3A in Ref. 59) in response to coupling than more ‘plastic’ oscillators with lower λ values. Moreover, rigid oscillators have smaller ranges of entrainment and narrower Arnold tongues. 20 Again, our results imply that any observation of resonant behavior and/or period changes might provide information on whether the underlying oscillator is a rigid versus plastic and hard versus soft clock.
We finish with an open question. Throughout our work we have claimed that period-amplitude correlations are widespread in both in vivo and in in silico clocks, and how they are critical to define how oscillators function in their environment. But, how does one integrate the circadian clock’s property of temperature compensation within this framework? The circadian clockwork is temperature-compensated, 18,52,69,103,104 which means that increasing temperatures do not speed up significantly the clock, which still runs at approximately 24 h. A proposed hypothesis for temperature compensation suggests that temperature-sensitivity of oscillation amplitude could stabilize the period. 18,69 In support of this supposition is evidence from Neurospora18 and Gonyaulax105 suggesting that amplitude increases in response to temperature, but decreases in Drosophila. 106 What are the twist effects behind this adaptation? Twist might play a role in how fast that temperature variation is sensed in the clockwork and modulate the adaptation to find the new steady-state rhythm of larger or smaller amplitude at the new temperature.
Conclusions
Despite the complexities in quantifying amplitude, our models stress the important role of circadian amplitudes and their correlation with oscillator period. Although amplitudes are known to be regulated by a number of internal and external cellular factors, including cellular biochemical rates, 107 light conditions, 108 coupling, 64 genetic and epigenetic factors 109 or ageing, 110 our findings stress that twist effects (i.e., co-modulations of amplitudes and periods) also ‘feed back’ and affect the interaction of oscillators with the environment, facilitating entrainment, fastening response to pulse-like perturbations or modifying the response of a system to coupling. The theory of our conceptual models can also be applied to other oscillating system such as cardiac rhythms, somite formation, central pattern generators or voice production.