Early warning signals for desynchronization in periodically forced systems
EEarly warning signals for desynchronization inperiodically forced systems
Pablo Rodríguez-SánchezEgbert H. Van NesMarten Scheffer2020-02-29
Abstract
Conditions such as insomnia, cardiac arrhythmia and jet-lag share a common feature: theyare all related to the ability of biological systems to synchronize with the day-night cycle. Whenorganisms lose resilience, this ability of synchronizing can become weaker till they eventually becomedesynchronized in a state of malfunctioning or sickness. It would be useful to measure this loss ofresilience before the full desynchronization takes place. Several dynamical indicators of resilience(DIORs) have been proposed to account for the loss of resilience of a dynamical system. Theperformance of these indicators depends on the underlying mechanism of the critical transition,usually a saddle-node bifurcation. Before such bifurcation the recovery rate from perturbations ofthe system becomes slower, a mechanism known as critical slowing down. Here we show that, for awide class of biological systems, desynchronization happens through another bifurcation, namely thesaddle-node of cycles, for which critical slowing down cannot be directly detected. Such a bifurcationrepresents a system transitioning from synchronized (phase locked) to a desynchronized state, or viceversa. We show that after an appropriate transformation we can also detect this bifurcation usingdynamical indicators of resilience. We test this method with data generated by models of sleep-wakecycles.
The phenomenon of endogenous circadian rhythms, first observed by the French polymath Jean-Jacquesd’Ortous de Mairan in 1729 (d’Ortous de Mairan [1729]), has transcended science to become part of thepopular culture, often referred to as the inner clock . The evolutionary convenience of synchronizing such inner clocks with the external cues, usually provided by regular astronomical events such as day-nightperiods and seasons, is well established (Foster and Kreitzman [2017]). Synchronization, thus, provesuseful for living systems and a difficulty to synchronize (and sometimes also to desynchronize) can bean indicator of sickness or malfunctioning. Some synchronization-related conditions include insomnia,jet-lag, arrhythmia or epilepsy (Glass [2001]).The transition from a synchronized to a desynchronized regime is discontinuous. The system is eithersynchronized or not. Therefore it could be that synchronization is a special kind of critical transition(Scheffer [2009]). This is relevant as there have been developed ways to foresee whether at criticaltransition is likely to occur (Scheffer et al. [2009]). These dynamic indicators of resilience (DIORs)are based on the phenomenon of “critical slowing down” (Wissel [1984], Van Nes and Scheffer [2007]).According to this theory, the recovery rate from perturbations decreases if systems are close to a criticaltransition. In time series we can measure critical slowing down using different indicators, such as increasedautocorrelation and variance (Dakos et al. [2012]).In the present work we illustrate with simple models that some transitions from synchronized to desyn-chronized states indeed can be related to a special kind of critical transition, namely a saddle-nodebifurcation of cycles. We show that after an appropriate transformation of the data, we can still usecritical-slowing down indicators to see if one of these transitions is likely to happen.1 a r X i v : . [ q - b i o . Q M ] M a r Methods
Our goal is to develop generic indicators for the risk of desynchronization of biological cycles such asthe sleep-wake cycle. To understand the properties of this system, we analyze a generic model of suchperiodically forced cyclic systems. This minimal model consists of two oscillators: a master (representingthe external forcing, for instance of a diurnal rhythm) and a slave (representing the organism’s state, forinstance its sleep/awake status). We represent each oscillator by its most basic feature: phase ( θ (cid:12) for themaster and θ for the slave) The master’s frequency is constant (i.e. the phase grows steadily from 0 to2 π in 24 h), and it is not affected by the slave’s dynamics. The slave’s dynamics are more complex: inthe absence of coupling it has a natural frequency, and an increasing tendency to synchronize with themaster if the coupling gets more intense. These features are captured by model (1). ( dθdt = ω − k · f ( θ − θ (cid:12) ) dθ (cid:12) dt = ω (cid:12) (1)In model (1) each oscillator shows a natural frequency ( ω and ω (cid:12) ). The first oscillator shows a tendencyto slow-down if θ is ahead of θ (cid:12) , and to speed-up otherwise. The function f measures the differencebetween θ and θ (cid:12) . Note that f has to be a periodic function (in the sense of f ( x + 2 π ) = f ( x )). This isa consequence of the cyclic nature of phases: by definition phases θ and θ + 2 π represent the same pointin a cycle, and thus, the same physical reality. In most applications f is also continuous and smooth.The strength of the coupling is given by the positive constant parameter k . If the coupling is not strongenough (relative to the difference in natural frequencies), synchronization doesn’t happen.The system (1) becomes simpler (and even analytically tractable) if we use the phase difference φ ( t ) ≡ θ ( t ) − θ (cid:12) ( t ) as a new state variable. With this change of state variable, the system takes the form (2),where, for convenience, we made Ω ≡ ω − ω (cid:12) . dφdt = Ω − k · f ( φ ) (2)For the sake of clarity, we will use f ( φ ) = sin( φ ) in the rest of this work. As we discuss in the onlineappendix, we can do this without loss of generality. With this choice, our model becomes a simple subcaseof the classical Kuramoto model (see Kuramoto [1975], Strogatz [2000]). Equating (2) to zero, the stableand unstable equilibria of our system are easily found to be φ ∗ s = ∆ and φ ∗ u = π − ∆, where ∆ ≡ arcsin Ω k .It is important to note that, for those equilibria to exist, condition (3) should be satisfied. Intuitively,this means that our system can only synchronize cycles whose difference in natural frequencies (Ω) haveat most the same order of magnitude as the coupling term ( k ). | Ω k |≤ φ is constant, so both oscillators have the same frequency ( ω (cid:12) ) and are thus synchronized.If, on the contrary, condition (3) is not satisfied, the phase difference φ never stabilizes and consequentlysynchronization is not possible.But, what happens at the border between both cases, that is, when Ω /k approaches 1? In such a situation,the stable and unstable solutions collide and annihilate each other at φ ∗ = π (see first row of panel 1).This mechanism of losing stability is known as a saddle-node bifurcation (Kuznetsov [1998], Strogatz[1994]). For an extensive discussion about the choice of this system, please refer to the appendix.We’ll take advantage of the periodicity of our system by plotting its trajectories over a 2 π × π squarewith periodic boundary conditions (or, equivalently, on the surface of a torus). When the phase hitsany border, it reappears at the opposite side (just like in old-school video games such as Pac Man or Asteroids ). In figure 1 we see three different configurations of such a system. As our parameter approachesthe saddle-node bifurcation, both the stable and unstable cycles get closer. When we introduce additive2 k * Configuration A/ k = 0.40 k Configuration B/ k = 0.95 k Configuration C/ k = 1.10 Figure 1: Each of the columns corresponds to a different configuration for system (1), identified by thevalue of the bifurcation parameter Ω /k , which represents the synchronization capacity. From left to right,each column represents less coupling strength. Each of the rows corresponds to a different representationof the dynamics. In the first row we see the bifurcation diagram of the phase difference ( φ ). The redarrows in the first row represent the flow on the line. In the second and third rows, the continuous red linerepresents the stable branch, and the dotted one the unstable branch. Saddle-node bifurcations happen at Ω k = ±
1. If | Ω k | > σ = 0 . θ, θ (cid:12) ) as coordinates(phase space) and the third row uses ( t, φ ) coordinates (time series). Notice that in the second column,even if Ω k <
1, a noise induced transition may happen due to the proximity of the stable and unstablecycles.noise to the dynamics, transitions can happen before the bifurcation is reached if the noise is strongenough to make the state jump the gap between both cycles (figure 1, column B). Note that due tothe periodic boundaries the system is only momentarily desynchronized as it “collapses” back to thesynchronized dynamics.
From the previous subsection it should be clear that synchronization is related to a fold bifurcation thatoccurs in the phase difference of the internal clock of a system respective to that of the forcing. This phasedifference is usually not directly measurable. Instead, experimental data of periodic phenomena usuallygives us indirect information about the phase. The angle of the Sun respective to the local meridian,the height of the tide or even the subjective feeling of sleepiness or hunger along the day are obviouslyaffected by the phase of the cycle under study (see figure 2). But, can we use these indirect measurementsto robustly infer the phase?In order to answer this question, we will translate the ideas illustrated in the previous paragraph andfigure 2 to mathematical language. Particularly, we’ll assume, as a working hypothesis, that there is acertain functional relationship M between the phase of the cycle θ ( t ) and our observations y ( t ) (equation(4)). Due to the periodic nature of our problem, we expect M to have a period of 2 π .3 idnightnoonlow tidehigh tideasleepawake 0 24 48 72 t ( h )02 Figure 2: The first row shows the Sun’s angular height from a local horizon. Second row represents theheight of the tide. Third row shows a sleep wake cycle of a healthy individual. The fourth and last rowshows a common phase for the three above-mentioned phenomena (thus, the time series in all rows canbe expressed in the form given in equation (5)). All series have been plotted for three whole periods. y ( t ) = M [ θ ( t )] (4)We define a reference cycle y ref ( t ) based in our knowledge about the system under study. For instance, ifwe are studying sleep cycles and y ( t ) represents the asleep state, a reasonable choice for y ref ( t ) couldbe y ref ( t ) = 0 (awake) if t is between 8 and 24 h, and 1 (asleep) otherwise. Such a function representsthe idealized sleeping cycle of a healthy individual. We assume the reference cycle to be the result ofapplying the unknown function M to the phase of the external forcing θ (cid:12) (equation (5)). y ref ( t ) = M [ θ (cid:12) ( t )] (5)If our system is either synchronized or subject to slow variations in its external conditions, we can considerthe phase difference ( φ ( t ) ≡ θ ( t ) − θ (cid:12) ( t )) approximately constant over a given time span [ t a , t b ]. It canbe shown (see appendix section A.1) that under this circumstances we can expect that y ( t ) is just shiftedin time respective to y ref ( t ) by a certain time delay λ (equation (6)). y ( t ) = y ref ( t + λ ) (6)We find the time delay λ that best fits our data by minimizing the sum of squares between the time-shiftedreference cycle and our measurements (see figure 3 and equation (7)). This time delay λ min is proportionalto the phase difference. In section A.1 of the appendix, we show that, specifically, φ = ω (cid:12) λ min . D ( λ ) = b X i = a ( y i − y ref ( t i + λ )) (7)4 t a t e expected ( y ref ( t ))measured data ( y ( t i ))0 24 48 72time (h) s t a t e displaced ( y ref ( t + ))measured data ( y ( t i )) Figure 3: In both panels the red crosses represent the hypothetical activity of a human being experiencinga jet lag, measured every 60 minutes during 3 days. In the upper row we see the expected daily activity inblue (no activity while sleeping between 0 and 8 hours, and activity the rest of the day). In the lower rowwe see, in orange, the expected daily activity, but now displaced 6 h in the time axis. This displacementprovides the best fit for the data, and is calculated by minimizing the function D ( λ ) given in equation(7). 5igure 4: Schematic outline of our method.Applying the method described above to different time windows allows us to use data to estimate a timeseries of the phase differences ( φ ( t j )) even if the precise analytical form of M is unknown (see first row infigure 5 and second row in figure 6). See appendix section A.1 for details. Saddle-node bifurcations are often preceded by the phenomenon of critical slowing down (Scheffer et al.[2009]). Such a phenomenon can be directly observed in the time series even if the underlying dynamicsare unknown. In the present manuscript we used the above-mentioned minimization algorithm along amoving window of typically 1-day width to extract the phase difference. Then, we applied the methodsproposed by Dakos (Dakos et al. [2012]) to analyze this time series. Particularly, we first detrendedthe time series by simply subtracting the average value over non-intersecting windows of 1 day length.Afterwards, we calculated the standard deviation and autocorrelation of the residuals, using a rollingwindow with a length around 25-50% of the original time series’ length. The optimal parameters (windowlength, autocorrelation lag, etc.) depend on the time scale and characteristics of the data under study. Formore details about this method, see Dakos et al. [2012]. For an extended discussion about the limitationsof these methods, see Dakos et al. [2015].Figure 4 summarizes all the steps, inputs and outputs of our method.
We tested our method with two model-generated time series.The first time series was generated with the help of the sleep-wake model of Strogatz (Strogatz [1987]).We configured the system so the time series represents the sleep-wake dynamics of an individual thatis becoming progressively more prone to insomnia. The insomnia effect was simulated by allowing thecoupling term k to linearly decrease to zero along a period of 135 days. This makes the individual’sinner clock progressively less capable of coupling with the day-night cycle, and eventually completelyunable to do so. Strogatz’s model (Strogatz [1987]) can be understood as a Kuramoto oscillator followedby a postprocessing function M that transforms the inner clock’s phase θ into a sleep-wake time series.Particularly, M ( θ ) returns 1 (awake) if the inner clock’s phase θ is between 2 π/ π radians(corresponding to 8 and 24 hours in the inner clock), and 0 otherwise (asleep). We used the generatedtime series to estimate the phase difference. As this model contains an explicit phase, we can use it as acontrol, and compare it with our estimated phase as a verification of our method for extracting phasesfrom data (see first row of figure 5).To show the generality of our method, we applied it to a second time series generated with a more realisticmodel, the Phillips-Robinson model (Phillips and Robinson [2007]). The Phillips-Robinson model is adeterministic sleep-wake model based on neurological considerations, and it doesn’t contain an explicitphase. It describes the time evolution of three state variables: V v the activity of the ventrolateral preopticarea (prompting the body to stay asleep), V m the activity of the mono aminergic group (prompting thebody to stay awake) and H the homeostatic pressure (an auxiliary variable that quantifies the need forsleep). The dynamics of the model are given by the equations (8), where F ( V ) is a saturation functiongiven by (9) and C ( t ) (defined in equation (10)) is a time-dependent external forcing, representing theastronomical light/dark cycle. The remaining elements in equation (8), including the influence of theacetylcholine group ( V a ), are just constants. The parameters used are the same as in (Phillips andRobinson [2007]); see section A.3 in the online appendix. Additionally, this appendix section provides agraphical representation of the relationships in equation (8).6 τ v dV v dt = − V v − ν vm S ( V m ) + ν vh H − ν vc C ( t ) τ m dV m dt = − V m − ν mv S ( V v ) + ν ma S ( V a ) χ dHdt = − H + µS ( V m ) (8) S ( V ) = Q max e − V − θσ (9) C ( t ) = 12 (1 + cos( ωt + α )) (10)Once again, we simulated an individual whose sleep quality is slowly deteriorating. We achieved thiseffect by allowing the coupling parameter ν vc to decrease linearly from its normal value of 6 . mV to 0 mV along a period of three months. By doing this, the ability of the subject to synchronize his internalclock with the external time cues slowly disappears. The first episode of insomnia/desynchronizationhappens on the 83rd day (see first row in figure 6).In order to simulate the fluctuations expected in any biological system we added noise to the integrationof both our time series. In particular, we modeled our systems as Wiener processes. The deterministicterms have been described in the previous paragraphs. The stochastic terms ( dW = σdt ) for Strogatz’smodel where set to σ = 0 .
05 for the inner clock’s phase, and to 0 for the driver. For the Phillips-Robinsonmodel, the stochastic term was set to σ = 1 for the states V v and V m , and 0 for H . The integration wasperformed numerically with the Python package sdeint . We applied the minimization algorithm described in the methods section to the sleep-awake time seriesgenerated with Strogatz’s model (see methods). As a reference of a healthy sleep-wake cycle, we used asimple assumption: a healthy individual is awake between 8 in the morning and 24 at night, and asleepotherwise. We managed to reconstruct correctly the phase difference. The reconstructed phase differenceshows the classical signs of slowing down (namely, increase in standard deviation and autocorrelation)when the system is approaching the bifurcation (see figure 5).Our method was also applied with success to the time series generated with the Phillips-Robinson model(equation (8)). We focused our attention only in the time series corresponding to the state variable H . Weused the time series corresponding to day 1 as a reference to estimate the phase difference. Particularly,we built y ref ( t ) as a quadratic interpolator of the measurements corresponding to the first day. As wecan see in figure 6, the first episode of insomnia (83rd day) is preceded by an increase in both standarddeviation and autocorrelation of the estimated phase difference. In the current work we presented a way of deriving dynamic indicators of resilience (DIORs) for systemstransitioning from synchronized to desynchronized states through the family of bifurcations known assaddle-node of cycles. Our method is designed for time series, and doesn’t require detailed knowledgeof the deterministic dynamics of the system. This makes it particularly suitable for biological systemswhere a loss of synchronization may have an undesired effect (such as insomnia or arrhythmia (Glass[2001])) or may be an indicator of a loss of resilience (such as the disruption in daily activity patterns incows after calving (van Dixhoorn et al. [2018])).It may be argued that our method rests on the particular choice of the model given in equation (1). Aswe discuss in the appendix, equation (1) represents the simplest, albeit non-trivial representative of abroader family of synchronization dynamics. Different choices yield different geometries in the bifurcationdiagram (figure 1, first row), but the main characteristic, the fact that at least one saddle-node bifurcationexists, remains true. This, together with the method to extract phase differences for general timeseries of periodically forced systems, makes our approach valid under very general circumstances. Two7 .51.01.5 P h a s e d i ff Exact phase differenceEstimated phase difference0.00.10.2 S D A C Figure 5: The black dots in the upper panel represent the phase difference, as approximated by ourmethod. The exact phase difference is also shown (green line) as a reference of the method’s accuracy.The central and lower panel show the standard deviation (in blue) and the autocorrelation with a 24 hlag (in orange) of the estimated phase difference, both of them calculated for a window 50% the length ofthe data. In both panels, as the time increases, the resilience of our system gets weaker.8 S t a t e s VvVmH0246 P h a s e d i ff S D A C Figure 6: The upper panel shows a simulated time series obtained by integrating the Phillips-Robinsonmodel under the influence of stochastic noise, and with the parameter ν vc decreasing linearly in time tosimulate an increasing difficulty in synchronizing. The simulation was initialized with non-transient values,to ensure the first days represent a healthy sleep wake cycle. The second panel shows the estimated phasedifference (in hours) using the method described in our paper (sub-sampling the whole time series onceper day, and using day 1 as reference). Our reference time series was chosen to be the somnogen level H during the first day. The two lower panels show the standard deviation (in blue) and the autocorrelationwith a 48 h lag (in orange), calculated for a window of 45 days into the past, of the estimated phasedifferences. The red line the 83rd day marks the first episode of insomnia/desynchronization.9pplication examples of time series that were generated with two different sleep-wake models, Strogatz’sand Phillips-Robinson’s are analyzed.The method to extract phase differences requires an approximate reference time series. In the Strogatz’smodel application example we used the very simple assumption that a healthy individual sleeps from 0h at night to 8 h in the morning. Some problems may benefit from or even require more sophisticatedassumptions. In the absence of any detailed knowledge of the system under study, another approach couldbe using the dynamics of an arbitrary day as a reference. This is what we did in the Phillips-Robinsonapplication example.Additionally, our method requires high quality time series. Those time series should be long (as we needmany cycles to infer the indicators) and should have a high density of data points (typically of the orderof 10 points measured per cycle, depending on the shape of the time series). This makes our method lesssuitable to be applied with success in fields where the data is difficult and/or expensive to collect. Luckily,data-rich systems such as the ones provided by wearable devices are becoming increasingly popular inmedicine or veterinary sciences.Even after extracting the phase difference, the critical slowing down may be difficult to detect under somecircumstances. As already noted in Dakos et al. [2015], his method to forecast saddle-node bifurcationshas some fundamental limitations. For instance, it is required that the time-scale of the changes in theexternal forcing to be slower than the natural time-scale of the system (that is, the transitions shouldn’tbe too sudden). The role of noise is also a delicate issue. On one hand, noise is required in order toobserve the phenomenon of critical slowing down in the vicinity of a saddle-node bifurcation. On theother hand, it obscures the deterministic dynamics. Our analysis will prove weak for systems whosedynamics are strongly dominated by noise. Our method, based on Dakos’ indicators, shares this set oflimitations.Even with those applicability challenges, we consider our method to be a step in the direction of forecastingtransitions between synchronized and desynchronized states. The fact that the disruption in certainphysiological rhythms is associated with disease (Glass [2001]), together with the recent increase in theavailability of high-quality biometric time series, makes the analysis and potential forecasting of theserelevant kind of transitions a topic worth being explored. We thank Ingrid van Dixhoorn and Rudi de Mol for their useful comments and suggestions.This work was supported by funding from the European Union’s
Horizon 2020 research and innovationprogramme for the
ITN CRITICS under Grant Agreement Number 643073.10 ppendix
A.1 Detailed derivation of the phase extracting method
In order to justify the results of section 2.2, we will make use of equations (1), (4), (5) and (6). As a firststep, we will write both sides of equation (6) in terms of M . We can do this by directly applying (4) and(5). The result is shown in equation (A.1). ( y ( t ) = M [ θ ( t )] y ref ( t + λ ) = M [ θ (cid:12) ( t + λ )] (A.1)In order be able to compare the functions given by equations (4) and (5) we will change their coordinates.Using the phase difference ( φ ( t ) ≡ θ ( t ) − θ (cid:12) ( t )), we can rewrite the equation (4) in terms of θ (cid:12) and φ (see first line in equation (A.2)). Introducing a leftwards shift λ in the time coordinate of equation (5) ittakes the form y ref ( t + λ ) = M [ θ (cid:12) ( t + λ )], where θ (cid:12) ( t + λ ) can be evaluated exactly by its first orderTaylor expansion, that is, θ (cid:12) ( t + λ ) = θ (cid:12) ( t ) + ω (cid:12) λ (cf. second line of equation (1)). The results of bothcoordinate transformations appear in equation (A.2). ( y ( t ) = M [ θ (cid:12) ( t ) + φ ( t )] y ref ( t + λ ) = M [ θ (cid:12) ( t ) + ω (cid:12) λ ] (A.2)Note that if the system is synchronized and/or if the adiabatic approximation (that is, that the externalconditions vary slowly) holds in our region of interest ( t ∈ [ t a , t b ]), φ ( t ) can be approximated by a constant φ . We can estimate the value of φ by finding the shift λ min that minimizes the square distance betweenboth functions (equation (A.3)). By direct inspection of equation (A.2) we see that this optimal valuecorresponds to a phase difference of φ = ω (cid:12) λ min . D ( λ ) = Z t b t a ( y ( s ) − y ref ( s + λ )) ds (A.3)When faced with experimental data, we’ll have a collection of N measured values y i sampled at times t i (that is, y i = y ( t i )). The discrete equivalent of equation (A.3), representing the square distance betweenour measured and the expected points, is given in (7). By finding the value of the time displacement λ that minimizes D ( λ ), we find the time delay that better fits our data (see figure 3). A.2 Further generalization
By manipulating the parameter Ω /k in an equation like (2) with any non trivial continuous function f ( φ )we are sure of encountering at least one saddle-node bifurcation. Even more, the saddle-node is the onlykind of bifurcation that may happen.This can be proven graphically. As we discussed in the methods section, the coupling function f in themodel given by equation (2) can be any non-constant, continuous, smooth and periodic function, notnecessarily a sine. A function f satisfying these properties will have at least one local minimum and onelocal maximum per period. This is also true for the right-hand side of equation (2). The effect of theparameter Ω /k is to move up and down the curve defined by y ( φ ) = Ω − kf ( φ ), whose roots representthe equilibria. This rules out the pitchfork and the transcritical bifurcations, as those require a change inthe shape of the curve y ( φ ) (Strogatz [2003]). By manipulating the parameter Ω /k , the only possiblebifurcations are collisions of stable and unstable equilibria, that is, saddle-node bifurcations (Strogatz[2003]). Those bifurcations happen when a minimum or a maximum equals 0 (see figure A.1).Those readers familiar with analysis may prefer noticing that, in the vicinity of a minimum/maximum( φ ), the second order Taylor expansion of the right-hand side of equation (2) can be written as theequation of a parabola (A.4). 11 K f () / K = 0.5 / K = 1.0 / K = 1.5 Figure A.1: Here we plot the curve y ( φ ) = Ω − kf ( φ ) for a non-sinusoidal coupling function f . Thefunction is non-constant, continuous, smooth and periodic. We plotted it for three different values of thebifurcation parameter Ω k . The roots of each curve represent the equilibria (filled dots if stable, white ifunstable). Ω − kf ( φ ) ≈ Ω − kf ( φ ) + kf ( φ )2 ( φ − φ ) (A.4)By using the new variable x ≡ q kf ( φ )2 ( φ − φ ) (representing a shift and re-scale of the horizontalaxis), and renaming Ω − kf ( φ ) as r , the right-hand side of equation (A.4) adopts the canonical form ofsaddle-node bifurcation, i.e.: r + x (Kuznetsov [1998]).Due to the generality of the conditions requested to the coupling function f , we expect saddle-nodebifurcations in the phase difference to be a widespread mechanism of synchronization and desynchronization.Consequently, we expect those bifurcations to be susceptible of being detected by the method describedin this manuscript. A.3 Parameters for Phillips-Robinson model
The parameters used in equation (8) are the same as in (Phillips and Robinson [2007]), with the exceptionof ν vc , that decreases linearly from 6 . τ m / hτ v / hχ . hν vm . / mV · hν mv . / mV · hν vh . mV · nM − µ − nM · hν vc . − mVν ma S ( V a ) 1 mVQ max · h − θ mVσ mVω π/ h − α v V m H C ( t ) V a ν vc ν ma µν mv ν vm ν vh Figure A.2: Schematic summary of the dynamics of the Phillips-Robinson model. The light blue nodesrepresent the system’s states ( V v the activity of the ventrolateral preoptic area, V m the activity of themono aminergic group and H the homeostatic pressure). The pink nodes represent the external sources( C ( t ), the astronomical light/dark forcing, and V a , the acetylcholine group constant influence). Thepositive effects are coded as green arrows. Negative ones as red arrows. Blue arrows represent oscillatingeffects. References
Vasilis Dakos, Stephen R Carpenter, William A Brock, Aaron M Ellison, Vishwesha Guttal, Anthony RIves, Sonia Kéfi, Valerie Livina, David A Seekell, Egbert H van Nes, and Marten Scheffer. Methodsfor detecting early warnings of critical transitions in time series illustrated using simulated ecologicaldata.
PloS one , 7(7):e41010, jan 2012. ISSN 1932-6203. doi: 10.1371/journal.pone.0041010. URLhttp://journals.plos.org/plosone/article?id=10.1371/journal.pone.0041010.Vasilis Dakos, Stephen R. Carpenter, Egbert H. van Nes, and Maórten Scheffer. Resilience indicators:Prospects and limitations for early warnings of regime shifts.
Philosophical Transactions of the RoyalSociety B: Biological Sciences , 370(1659):1–10, jan 2015. ISSN 14712970. doi: 10.1098/rstb.2013.0263.Jean-Jacques d’Ortous de Mairan. Observation botanique.
Histoire de l’Academie Royale des Sciences ,31:35–36, 1729.Russell G. Foster and Leon Kreitzman.
Circadian rhythms. A very short introduction . Oxford UniversityPress, 2017. ISBN 9780198717683.Leon Glass. Synchronization and rhythmic processes in physiology.
Nature
InternationalSymposium on Mathematical Problems in Theoretical Physics , pages 420–422. 1975.Yuri A. Kuznetsov. Elements of Applied Bifurcation Theory, Second Edition.
Library , 1998. ISSN0066-5452. doi: 10.1007/b98848.A.J.K. Phillips and P.A. Robinson. A Quantitative Model of Sleep-Wake Dynamics Based on thePhysiology of the Brainstem Ascending Arousal System.
Journal of Biological Rhythms , 22(2):167–179,apr 2007. ISSN 0748-7304. doi: 10.1177/0748730406297512. URL http://journals.sagepub.com/doi/10.1177/0748730406297512.Marten Scheffer.
Critical Transitions in Nature and Society . 2009. ISBN 9780691122038. doi: 10.5860/CHOICE.47-1380. URL http://books.google.com/books?id=jYSZgaaxRv0C.13arten Scheffer, Jordi Bascompte, William a Brock, Victor Brovkin, Stephen R Carpenter, Vasilis Dakos,Hermann Held, Egbert H van Nes, Max Rietkerk, and George Sugihara. Early-warning signals forcritical transitions.
Nature , 461(7260):53–59, 2009. ISSN 0028-0836. doi: 10.1038/nature08227. URLhttp://dx.doi.org/10.1038/nature08227.Steven H. Strogatz. Human sleep and circadian rhythms: a simple model based on two coupled oscillators.
Journal of Mathematical Biology , 25(3):327–347, jul 1987. ISSN 0303-6812. doi: 10.1007/BF00276440.URL http://link.springer.com/10.1007/BF00276440.Steven H Strogatz.
Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, ChemistryAnd Engineering . 1994. ISBN 9788187169857.Steven H. Strogatz. From Kuramoto to Crawford: exploring the onset of synchronization in populationsof coupled oscillators.
Physica D: Nonlinear Phenomena , 143(1):1–20, 2000. ISSN 01672789. doi:10.1016/S0167-2789(00)00094-4. URL https://doi.org/10.1016/S0167-2789(00)00094-4.Steven H Strogatz.
SYNC: The Emerging Science of Spontaneous Order . Penguin Books, 2003. ISBN978-0-14-100763-2.I.D.E. van Dixhoorn, R.M. de Mol, J.T.N. van der Werf, S. van Mourik, and C.G. van Reenen. Indicatorsof resilience during the transition period in dairy cows: A case study.
Journal of Dairy Science
American Naturalist , 169(6):738–747, jun 2007. ISSN 00030147. doi:10.1086/516845.C. Wissel. A universal law of the characteristic return time near thresholds.