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

Three parameters are important to characterize a circadian and in general any biological clock: period, phase and amplitude. While circadian periods have been shown to correlate with entrainment phases, and clock amplitude influences the phase response of an oscillator to pulse-like zeitgeber signals, the co-modulations of amplitude and periods, which we term twist, have not been studied in detail. In this paper we define two concepts: parametric twist refers to amplitude-period correlations arising in ensembles of self-sustained, limit cycle clocks in the absence of external inputs, and phase space twist refers to the co-modulation of an individual clock’s amplitude and period in response to external zeitgebers. Our findings show that twist influences the interaction of oscillators with the environment, facilitating entrainment, speeding upfastening recovery to pulse-like perturbations or modifying the response of an individual clock to coupling. This theoretical framework might be applied to understand the emerging properties of other oscillating systems.


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.[5][6][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,146][17][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][18][19][20][21][22] with consequences in the size of the phase response curve 5 or in jet lag duration. 23Amplitudes, together with periods, also govern entrainment 23,24 and seasonality. 22,23,25There 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,24ken 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,28By 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 zeitgeberinduced 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.

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 mass m that, when displaced from its equilibrium position, experiences a restoring force F proportional to the displacement x, namely F ¼ À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 29,30 where A represents the amplitude; ϕ, the phase; and τ represents the period of the motion, τ ¼ 2π ffiffi ffi m k p .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 À βx 3 allows the conversion of the simple harmonic oscillator into a Duffing oscillator. 31Depending on the sign of β, the coefficient that determines the strength of the nonlinear 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 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 ¼ _ x, the equation can be reformulated as a system of two first order ODEs: (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 where k 1 , k 3 and k 5 represent the rates of synthesis of x, y and z, respectively; k 2 , k 4 and k 6 , 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 To mimic clock heterogeneity and evaluate the amplitude-period correlations among ensembles of Gonze or Goodwin oscillators, degradation rates were varied around AE10% their default value (Table 1).

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: 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 AE20% 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.

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. 41We 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 The first equation describes the rate of change of the radial coordinate r t ð Þ (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 ϵ 6 ¼ 0, the phase changes at a constant rate only when r ¼ A; if any perturbation is to modify r such that r 6 ¼ 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:05 h À1 , τ ¼ 24 h and ϵ values of 0 or AE0:1 h À1 .
To study phase space twist and the effects of different ϵ values on the oscillator's response to pulse-like perturbations pert t ð Þ, periodic zeitgeber input Z t ð Þ 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: where The zeitgeber Z t ð Þ is given by: where T represents the zeitgeber period and F Z the strength (amplitude) of the zeitgeber input.The mean-field M coupling is given by: where K represents the coupling strength and N is the number of oscillators in the coupled ensemble.Lastly, the squarelike perturbation pert t ð Þ is defined as follows: where F P is the strength of the perturbation and t start is the time at which the perturbation starts.The perturbation lasts 1 h and is set to F P ¼ 0:7 a.u. in all our simulations.
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.While ϵ serves 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 By definition, the latent phase velocity _ Φ necessarily increases at units of the angular velocity ω ¼ 2π τ as the oscillator follows its kinetic equation: Combining the two previous equations and with the use of the chain rule, we can calculate the radius dependency of the latent phase Φ: The terms dϕ dt and dr dt are defined in the Poincaré model (equation 7), so that the equation above can be rewritten as: Next, the solution for f r ð Þ can be found through integration: like in all our simulations), the last term can be neglected such that Finally, the solution of f r ð Þ can be inserted in equation 9 to end up with the equation for isochrones as a function of the radius: Isochrones are thus loci of polar coordinates ϕ,r ð Þin phase space with the same latent phase Φ. Equation 11 shows how the shape of the isochrone (and, consequently, phase space twist effects) for the Poincaré model in equation 7 depends on the ratio of twist ϵ to relaxation rate λ.Specifically, higher ϵ values (in modulus) and lower relaxation rates λ lead to more curved isochrones and greater phase space twist effects.In the Results section, we further elaborate on how to interpret the isochrones graphically.
To plot the isochrones, we simply reorder equation 11 to plot ϕ as a function of r at isochrones with fixed values of latent phases 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 N tst ¼ 150, N max ¼ 20000, Ds min ¼ 0:0001 and Ds 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 ¼ Àkx. 44The 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,30he 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 À βx 3 , 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 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. 41,48Stable 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 , 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,26We 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 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.
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][51][52][53] or synchronization and entrainment. 33,54,55The 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 (k 2 , k 4 and k 6 , respectively) in each oscillator around AE10% 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 twistcontrol, 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 AE10% 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).  1. Relative amplitudes are computed as the peak-to-trough distance of the x variable, normalized to its mean.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.
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 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, 38 and thus the twist picture is expected to become even more complex.

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. 39o study parametric twist in an ensemble of uncoupled Almeida oscillators, we randomly changed all 18 parameters individually around AE20% 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 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 γ 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 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 (normalized to its mean) 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 nonmonotonic 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.All panels display the numerical results from bifurcation analyses of Gonze oscillators obtained with XPP-AUTO.
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, 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,57This 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).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.
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 steadystate oscillation following a perturbation.
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 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 parameter ϵ 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).
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 . 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. 58nterestingly, 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.
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.
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. 63A number of studies have shown that networks of coupled oscillators behave in a fundamentally different way than 'plain' uncoupled oscillators. 33,64,65Moreover, 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 meanfield (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:05 h À1 ),  .Phase space twist ϵ 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).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.

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. 68kka 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 k 2 resulted in amplitude increases with a constant period (no twist).Changes in k 4 led to increased amplitude and decreased period (negative twist), while changes in k 6 resulted in decreases in both amplitude and period (positive twist). 66However, our bifurcation analyses revealed different correlations.We observed positive twist when modifying k 2 or k 4 , but a negative correlation with changes in k 6 (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,69This 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 interoscillator 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,49However, different repression mechanisms exist across different species, and such mechanisms might be modelled with different terms.For example, in Drosophila 70 and in mammals, [71][72][73][74][75] transcriptional inhibition involves the formation of a 1:1 stoichiometric complex between activators and inhibitorsa process known as protein sequestration-based repression. 76,77Here, 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,77The description of such protein sequestration processes has been described with terms different from Hill functions. 76,77n contrast, species like Neurospora show a lower than 1:1 stoichiometry between activators and repressors in nuclei. 78In this fungus, self-sustained oscillations require numerous phosphorylations in the activator. 78,79As the number of mutated phosphorylation sites increases, the circadian rhythms become weaker and finally arrhythmic. 78Such multistep cooperative reactions can be described with Hill curves. 80,81 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 jϵj 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,82Consistently, 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. 21ansient 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 (ϵ 6 ¼ 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,84Of 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. 64Our 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 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 meanfield can also be obtained without twist, in networks of coupled oscillators with different intrinsic periods. 64But 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,97Such 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,93Moreover, 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. 99Experimental evidence suggests that these seasonal changes are regulated by coupling via neurotransmitters such as GABA 100,101 or VIP. 102These 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,100Thus, 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,57It 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 ϵ 6 ¼ 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. 20Again, 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,69In support of this supposition is evidence from Neurospora 18 and Gonyaulax 105 suggesting that amplitude increases in response to temperature, but decreases in Drosophila. 106What 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.

Stephanie Taylor
Department of Computer Science, Colby College, Waterville, Maine, USA In this article, the authors describe the phenomenon of twist (period-amplitude inter-dependence) using a range of differential equation oscillator models, ultimately with a view to understanding circadian clocks.The simplest models (harmonic and Duffing) are used to demonstrate the nonlinearity introduces a relationship between amplitude and period.The authors then switch to limitcycle oscillators, because limit-cycle oscillators are more appropriate for studying circadian clocks.
They use two approaches to studying twist.The first, referred to as "parametric twist", computes the periods and amplitudes of oscillators with varied parameter values, creating graphs that relates period to amplitude (and thus demonstrate twist -if the correlation is positive, there is positive twist, and negative correlation indicates negative twist).This parametric twist analysis is performed using two limit cycle models, one with one feedback loop and one with more (synergistic) loops.The second, referred to as "phase space twist", examines the response of a single oscillator to a temporary change in one parameter, mimicking signaling.A separate limitcycle model (a Poincare oscillator) is used for this analysis.Important to the discussion of phase space twist is the concept of isochrons, which are useful for understanding the effects of perturbations on the phase and speed of the clock (does it speed up or slow down the clock?) and it is easy to see how amplitude and phase response are related, when we look at how curved the isochrons are.The authors relate these analyses to phase response curves, recovery from jetlag, the range of entrainment, and temperature compensation and conclude, that although measuring amplitude and understanding relationship between amplitude and period are highly complicated, it is nonetheless important to be mindful of the role of amplitude and its correlation to period when studying clocks.
On the one hand, I really appreciated that the authors chose models that best illustrated the points and that the paper is organized to go from simple ideas to more complex ideas.It was helpful to have the initial explanation of twist use duffing oscillator (which has non-linearity) as a contrast to a harmonic oscillator (which doesn't).It provided a nice introduction to the cause of twist and helped me to develop an intuition.Then, it made sense to use a very simple (3-ODE Goodwin oscillator) as the first example of a limit cycle oscillator.Since it has so few parameters, it was possible to perturb them one at a time and in pairs to illustrate a good portion of the possibilities for how one might observe heterogeneity in the system.Then, moving to Almeida's oscillator demonstrated how quickly the analysis can become very complex.All of this provides a nice progression of complexity.Then, we pivot a bit to the Poincare oscillator, to study twist in a more dynamical way -instead of permanently choosing different parameters, the authors model light signaling (at least, initially) and study the effects of the model returning to its limit cycle after the perturbation is over.This model is particularly appropriate for studying the dynamics of recovery because it is a 2-ode model (so we can examine it in "phase space", which is helpful for visualizing the clock) and because there is a parameter that controls how related period and amplitude are.It is also straight-forward to adjust how curved the isochrons are, which makes illustrating twist very clear.Then, after studying the response of a single oscillator to light, they use the same model to study the effects of intercellular signaling and the effects of different twist on the properties of the ensemble oscillation (e.g.how long the ensemble period is in contrast to that of individual oscillators).
On the other hand, I would have appreciated a more explicit connection between the discussion of "parametric twist" and "phase space twist".They are both due to changes in parameter values causing changes in period and amplitude.In the one case, we measure the period and amplitude of oscillators with different (constant) parameter values.In the second case, we measure instantaneous changes in period and amplitude due to a temporary changes in parameter values (I am referring to Z(t) as a parameter here).In both cases, the shape of the isochrons in the unperturbed condition should help us to predict what will happen in the perturbed condition.There are two challenges with making that point clear with the Goodwin and Almeida models -1) the isochrons aren't lines, they are multi-dimensional manifolds, and 2) the isochrons in the unperturbed and perturbed systems are different.So I am asking for too much to relate the two discussions by using images of isochrons, but a more explicit comparison would help the reader understand how the analyses are related to each other.
In the analysis of the modified Poincare oscillator, the term "twist" is used differently than it is in the "parametric twist" analysis.In the parametric analysis, the twist captures the full effects of the perturbations (we plot the measured periods and amplitudes across the heterogeneous set of oscillators).In the phase space analysis, "twist" refers to a particular parameter \epsilon that, when non-zero, introduces a dependence of the speed and amplitude on each other.However, that twist parameter (\epsilon) is not the only parameter that affects the skewness of the isochrons.The relaxation rate \lambda also affects the skewness of the isochrons.And the skewness of the isochrons affects the changes in speed and amplitude when the oscillator is signaled (i.e. a parameter is temporarily perturbed).So, if we were to measure the effects on period and amplitude (as is done with figures 7G, H, and I), that reflects that values of both \lambda and \epsilon.When I was thinking about how to relate the "parametric twist" and "phase twist" analysis, I was looking for aspects of figures 3 (Gonze model), 5 (Almeida model), and 7 (poincare model) to relate and it seemed to me that that the subfigures that related the amplitude of the oscillator to its period (3C, 3D, 3E, 3F, 3G, 3H, 5B, 5B, 5E, 7G, 7H, 7I) should be the figures that I can connect.They are labeled as twist in figures 3 and 5.In Figure 7, they show an approximation of instantaneous amplitude and period (rather than steady-state amplitude and period), so it isn't a perfect 1-1 mapping, but it is still the closest that makes sense to me.And in all of these cases, the output could be labeled twist if \epsilon didn't get to claim that label.I think it would make the concept of twist clearer to the reader if it were consistent across both analyses that twist is the "effect" we can measure.My reading of the text as it is currently written is that twist is a phenomenon we can measure in "parametric twist" analyses and that twist is a parameter that causes one change in the isochrons (but not all of the changes that can happen) that ultimately leads to something we can measure.Also, labeling the parameter as twist in the Poincare oscillator makes it seem like we need to use the Poincare oscillator to perform a "phase space twist" analysis.But we aren't limited to the Poincare oscillator for such analyses.

Specific Points
It would be helpful to have a little more framing text that indicates the different models are used to make different points.
○ I think of the "phase space twist" analysis as analysis of recovery from temporary changes to parameters and "parametric twist" as analysis of permanent differences in parameters.If that is the way you intend readers to think about it, then I would suggest changing the names to reflect those ideas more directly.Right now, the term "phase space twist" stands out as a highly specific analysis -the model most be 2D so you can draw the isochrons.But I think that is limiting.I totally embrace idea that there is both heterogeneity among oscillators and temporary changes due to intercellular signaling, so I would love to see it characterized that way more explicitly.What about "twist due to heterogeneity" and "twist due to signaling"?They aren't quite zippy enough and rely on my characterization of twist as a measurable phenomenon, so I don't push those terms exactly.

○
On page 8, when you are introducing isochrons, it would be helpful to have a figure to refer to.The most natural figure is Figure 7, but that would force you to make an unnatural reordering of the figures.So maybe introduce a new figure?
○ On Page 11, (analyses shown in Figure 3), how did you choose which parameters to pair together?
○ On Page 11, there is a discussion of the Hopf bifurcation that is interesting (and highlights how sensitivity twist is to our choice of parameters).It occurs to me that this has implications for modeling SCN oscillators as they recover from TTX treatment.The individual oscillator amplitude crashes and the oscillators lose synchrony during TTX and both recover after TTX is washed out.Is there anything you could add to the discussion about this? ○ On Page 12, in the paragraph about how overall twist effects depend not only on the parameter being studies, but also on the variable that is measured, I was reminded of early modeling studies that showed not-immediately intuitive results.In his 1995 fly model paper (A model of circadian oscillations in the Drosophila period protein, Proc.R. Soc.Lond.B, 1995), Goldbeter showed that the rate of (doubly-phosphorylated) PER degradation affected both period and amplitude and that increasing this rate led to PER amplitude increasing (rather than decreasing, as one might expect).This can be reasoned through, but it is another illustration of the complexity of the non-linear feedback and the relationship between period and amplitude.Is this worth bringing up here or on the discussion?○ Figure 6.Why are there no yellow dots on the limit cycle for A and B? Also, maybe make the dot color on the limit cycle different from yellow (e.g. to blue), so that you can describe what is happening over time with even more clarity (i.e. that the simulation eventually reaches the "blue" dot after enough cycles).
○ Figure 8 and text referring to it.I appreciate that you relate twist to entrainment and PRCs.I think this is interesting and shows why twist is important to understand.
○ Page 16 first paragraph."This suggests that single oscillator twist might also affect how coupling synchronizes a network of oscillators" is too weak.Given that coupling (M) and the perturbation used in Figure 7 (Z(t)) both affect the model the same way (they both are added to the rate equation for the x variable), Figure 7 shows that the PRCs are affected by twist.Not only are they affected, they are affected dramatically -Figures 7K and 7L look like they are opposite each other (i.e.flipped sign).This means we are guaranteed that the oscillators will respond differently to intercellular signaling and that that the ensemble oscillation will be different.

○
Page 16, in the discussion of the ensemble period, we see that for this model that positive twist leads to longer periods (and negative twist to shorter).What is missing from this discussion is that the relative phase of the signal and of the PRC of the receiver determine the period of the ensemble.Positive twist leads to longer periods because the mean field is determined by the value of the x variable.If it were determined by a variable with the opposite phase, then we would have the opposite effect.This affects the related paragraph in the Discussion as well.
○ Page 17.The meaning of "For this reason we have not termed twist, in the context of conservative oscillations, parametric or phase space" wasn't clear to me.I think you mean you can't study twist in non-limit-cycle oscillators in the same way you can for limit cycle oscillators.Is that correct?Reviewer Expertise: I study the response of limit cycle oscillators to signals and sensitivity analysis involving isochrons.
I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Is the study design appropriate and is the work technically sound? Yes
Are sufficient details of methods and analysis provided to allow replication by others?Yes If applicable, is the statistical analysis and its interpretation appropriate?Yes Are all the source data underlying the results available to ensure full reproducibility?Yes

Are the conclusions drawn adequately supported by the results? Yes
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Circadian rhythms, computational modeling of circadian clocks, suprachiasmatic nucleus (SCN), bioluminescence reporter imaging I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
The benefits of publishing with F1000Research: Your article is published within days, with no editorial bias • You can publish traditional articles, null/negative results, case reports, data notes and more • The peer review process is transparent and collaborative • Your article is indexed in PubMed after passing peer review • Dedicated customer support at every stage • For pre-submission enquiries, contact research@f1000.com

Figure 1 .
Figure 1.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 culture26 ; (B) negative correlations have been observed in cells from the choroid plexus in the mouse brain.14

Figure 2 .
Figure2.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.

Figure 3 .
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 AE10% 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 of k 2 and k 4 produce significant positive amplitudeperiod 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 Table1.Relative amplitudes are computed as the peak-to-trough distance of the x variable, normalized to its mean.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.

Figure 5 .
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,37 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 AE20% of their default value (shown in Table2) 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 (γ 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 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.

Figure 6 .
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).The twist parameter ϵ 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).

Figure 7 .
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 ϵ, 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 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.

Figure 8 .
Figure 8. 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 period τ ¼ 24 h, amplitude A ¼ 1 and twist ϵ ¼ 0:05 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 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 of λ and ϵ.

Figure 9
Figure 9. Phase space twist ϵ 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).
GitHub: Source code to generate and analyze the data and to reproduce all figures.https://github.com/olmom/twistData and codes are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).ultrasensitivity and oscillations.Interface Focus.2022; 12 (3): 20210084 PubMed Abstract | Publisher Full Text 4. Kim JK: Protein sequestration versus Hill-type repression in circadian clock models.IET Syst Biol.2016; 10 (4): 125-35 PubMed Abstract | Publisher Full Text Is the work clearly and accurately presented and does it cite the current literature?Partly Is the study design appropriate and is the work technically sound?Yes Are sufficient details of methods and analysis provided to allow replication by others?Yes If applicable, is the statistical analysis and its interpretation appropriate?Yes Are all the source data underlying the results available to ensure full reproducibility?Yes Are the conclusions drawn adequately supported by the results?Yes Competing Interests: No competing interests were disclosed.Reviewer Expertise: Circadian rhythms, mathematical modeling I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.Reviewer Report 18 September 2023 https://doi.org/10.5256/f1000research.148657.r202731© 2023 Taylor S. This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

○
Is the work clearly and accurately presented and does it cite the current literature?YesIs the study design appropriate and is the work technically sound?YesAre sufficient details of methods and analysis provided to allow replication by others?YesIf applicable, is the statistical analysis and its interpretation appropriate?YesAre all the source data underlying the results available to ensure full reproducibility?YesAre the conclusions drawn adequately supported by the results?YesCompeting Interests: No competing interests were disclosed.

Table 1 .
Default parameters of the Goodwin and Gonze models used for our simulations.Their biological interpretation as well as the default values from publications are shown.

Table 1 .
Continued *Parameters have been adapted to obtain limit cycle oscillations, since the default degradation parameters from Ref. 36 produced weakly damped rhythms.

Table 2 .
37st of parameters of the Almeida model.Time units are given in hours and concentration units as arbitrary units.The default values are taken from the original reference, Almeida et al.37