Authors: Marta del Olmo, Christoph Schmal, Camillo Mizaikoff, Saskia Grabe, Christian Gabriel, Achim Kramer, Hanspeter Herzel
Categories: Research Article, amplitudes, circadian clocks, coupling, entrainment, heterogeneity, mathematical modeling, periods, twist, Articles
Source: F1000Research
Are circadian amplitudes and periods correlated? A new twist in the story
Authors: Marta del Olmo, Christoph Schmal, Camillo Mizaikoff, Saskia Grabe, Christian Gabriel, Achim Kramer, Hanspeter Herzel
Abstract
Three parameters are important to characterize a circadian and in general any biological 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 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. 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 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 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 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.
Materials and Methods
Conservative linear and non-linear 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)
d2xdt2+kmx=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 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:
dxdt=vdvdt=−kmx−βmx3. (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).
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
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
dBMAL1dt=RORE−γBPBMAL1 PERCRYdRORdt=Ebox+RORE−γRor RORdREVdt=2Ebox+Dbox−γRev REVdDBPdt=Ebox−γDbDBPdE4BP4dt=2Rbox−γE4E4BP4dCRYdt=Ebox+2Rbox−γPCPERCRY+γCPPERCRY−γC CRYdPERdt=Ebox+Dbox−γPCPERCRY+γCPPERCRY−γPPERdPERCRYdt=γPCPER CRY−γCPPERCRY−γBPBMAL1 PERCRY. (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.
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
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 different
ϵ values on the oscillator’s response to pulse-like perturbations
pertt , periodic zeitgeber input
Zt 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
xi variable as
dxidt=λxiA−ri−yiω+ϵA−ri+Zt+M+perttdyidt=λyiA−ri+xiω+ϵA−ri, (8) where
ri=xi2+yi2 . The zeitgeber
Zt is given
Zt=FZcos2πTt+π2, where
T represents the zeitgeber period and
FZ the strength (amplitude) of the zeitgeber input. The mean-field
M coupling is given
M=KN∑i=1Nxit, where
K represents the coupling strength and
N is the number of oscillators in the coupled ensemble. Lastly, the square-like perturbation
pertt is defined as
pertt=0,iftstart≤t≤tstart+1FPotherwise, where
FP is the strength of the perturbation and
tstart is the time at which the perturbation starts. The perturbation lasts 1 h and is set to
FP=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
Φ=gϕr=!ϕ−fr. (9)
By definition, the latent phase velocity
Φ˙ necessarily increases at units of the angular velocity
ω=2πτ as the oscillator follows its kinetic
dΦdt=!ω (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Φdt=9dϕdt−dfrdrdrdt=10ω
The terms
dϕdt and
drdt are defined in the Poincaré model ( equation 7), so that the equation above can be rewritten
ω+ϵA−r−dfrdrλrA−r=10ω
Solving for
dfrdr :
dfrdr=ϵA−rλrA−r=ϵλr
Next, the solution for
fr can be found through
fr=ϵλ∫r=Ar1rdr=ϵλlnr−ϵλlnA
If
A=1 (like in all our simulations), the last term can be neglected such that
fr=ϵλlnr
Finally, the solution of
fr can be inserted in equation 9 to end up with the equation for isochrones as a function of the
Φ=ϕ−ϵλlnr (11)
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
Φ=0π4π23π4π−3π4−π2−π4 :
ϕ=ϵλlnr−Φ
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 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.
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 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 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 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.
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.
Phase space twist in single 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 if one takes a point in phase space and observes where the system returns to after exactly one period, the answer is 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.
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 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 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 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.
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 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 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 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 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 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 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 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
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).
^
66
^ However, 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
^
^,^
^
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 Drosophila ^ 70 ^ and in mammals, ^ 71 ^ ^–^ ^ 75 ^ transcriptional inhibition involves the formation of a 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, 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 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 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 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 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 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 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 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 Neurospora ^ 18 ^ and Gonyaulax ^ 105 ^ 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.