Generation and Disruption of Circadian Rhythms in the Suprachiasmatic Nucleus: A Core-Shell Model

We focus our research on how the core-shell organization controls behavior of the suprachiasmatic nucleus (SCN), how the core and shell are synchronized to the environment, what impact they have on the behavior of the SCN under different lighting conditions, and what mechanisms disrupt synchronization. To this end, we use a reduced Kuramoto model, with parameters inferred from experimental observations and calibrated for mice, and perform a detailed comparison between the model and experimental data under light-dark (LD), dark-dark (DD), and light-light (LL) conditions. The operating limits of free-running and entrained SCN activity under symmetric LD cycles are analyzed, with particular focus on the phenomena of anticipation and dissociation. Results reveal that the core-shell organization of the SCN enables anticipation of future events over circadian cycles. The model predicts the emergence of a second (dissociated) rhythm for large and small LD periods. Our results are in good qualitative and quantitative agreement with experimental observations of circadian dissociation. We further describe SCN activity under LL conditions and show that our model satisfies Aschoff’s first rule, according to which the endogenous free-running circadian period observed under complete darkness will shorten in diurnal animals and lengthen in nocturnal animals under constant light. Our results strongly suggest that the Kuramoto model captures essential features of synchronization and entrainment in the SCN. Moreover, our approach is easily extendible to an arbitrary number of groups, with dynamics described by explicit equations for the group phase and synchronization index. Viewed together, the reduced Kuramoto model presents itself as a useful tool for exploring open problems in the study of circadian rhythms, one that can account for evolving views of the circadian system’s organization, including peripheral clocks and inter-hemispheric interaction, and can be translated to other nocturnal and diurnal animals, including humans.

Circadian rhythms are periodic physiological and behavioral changes that take place in living organisms over 24-h cycles. Following decades of experimental and theoretical work, it is now understood that circadian rhythms are orchestrated by the brain's suprachiasmatic nucleus (SCN) located in the hypothalamus (Evans and Gorman, 2016;Michel and Meijer, 2020;. These endogenous rhythms are entrained to environmental cues such as light, and their disruption has a significant and sometimes detrimental impact on bodily function, including jet lag (disrupted sleep resulting from rapid travel across time zones), dissociation (emergence of a second rhythm in the SCN beyond the entrainment range, that is, for sufficiently large or small light-dark [LD] periods), and splitting (antiphase rhythms within the SCN) (Pittendrigh and Daan, 1976c;de la Iglesia et al., 2000;Michel et al., 2013). Despite considerable progress in our understanding of the SCN, the role of its organization in the generation of circadian rhythms remains an open problem (Hastings et al., 2018). At the population level, the current paradigm of SCN organization is based on 2 interconnected parts, the ventral (core) and dorsal (shell) regions of the SCN (Leak et al., 1999;Abrahamson and Moore, 2001;Evans and Gorman, 2016;Michel and Meijer, 2020;Hastings et al., 2018) where only the core is retinorecipient (Meijer et al., 1986). At the single cell level, self-sustained circadian rhythms are known to originate in feedback loops involving the transcription of clock genes (Ko and Takahashi, 2006), which persist in the absence of exogenous cues.
One approach to modeling circadian rhythms is based on molecular models of individual neuron activity (Leloup and Goldbeter, 2003;Vasalou and Henson, 2011;Hafner et al., 2012;Taylor et al., 2017). However, a single hemisphere of the SCN comprises approximately 10,000 neurons so that molecular models require a large number of equations and parameters to characterize both individual cells and the network organization of the SCN. This makes such models computationally harder and less amenable to analytical analysis. Moreover, experimental evidence and simplified theoretical models show that population-level activity in the SCN emerges from the dynamics of individual neurons (Schibler et al., 2015;Gu et al., 2016Gu et al., , 2018. Among simplified models of SCN dynamics, the Kuramoto model is a biologically plausible approach since clock cells in the SCN demonstrate self-sustained, nearly sinusoidal oscillations. It accounts for the heterogeneity of free-running periods observed across individual cellular oscillators (neurons) in the SCN (Liu et al., 1997;Pauls et al., 2014;Welsh et al., 2010;Taylor et al., 2017). The Kuramoto model has been used to explain the eastwest asymmetry in jet lag recovery times (Lu et al., 2016), seasonal adaptation of the circadian clock (Gu et al., 2016;Meijer et al., 2010), and phase splitting within the SCN (Rohling and Meylahn, 2020). With few exceptions (Hannay et al., 2020), existing work either does not account for the core-shell organization of the SCN, which continues to guide considerable research into the functioning of the SCN (Evans et al., 2015;Evans and Gorman, 2016;Michel and Meijer, 2020;, or explicitly models the dynamics of individual neurons. Here, we mainly consider how the core-shell organization of the SCN controls SCN behavior, how the core and shell are synchronized to the environment, what impacts they have on the behavior of the SCN under different lighting conditions, and the mechanisms that disrupt synchronization. For this purpose, we use a reduced Kuramoto model, which explicitly describes the phase and synchronization index of the core and shell. The benefit of our model is that it uses only 9 biologically meaningful parameters, which can be inferred from experiments. Using this model, we study the operating limits of free-running and entrained SCN activity under symmetric LD, darkdark (DD), and light-light (LL) cycles. In particular, we consider how the core-shell organization of the SCN enables anticipation of regular events and the emergence of a dissociated rhythm under different lighting conditions. Model parameters are inferred from experimental observations for mice and can be calibrated for other nocturnal and diurnal mammals. Our numerical and analytical results demonstrate good qualitative and quantitative agreement with experimental results.

The Reduced Kuramoto Model
According to anatomical data, neurons in the SCN are densely interconnected (Güldner, 1976;Moore and Bernstein, 1989;Varadarajan et al., 2018). Each SCN neuron forms between 300 and 1000 connections with other SCN neurons. Such dense connectivity means that many important features of SCN dynamics can be described by a mean-field model with all-to-all interaction (Dorogovtsev et al., 2008). The SCN in the left and right hemispheres comprises 2 main groups of clock cells, namely, the core and shell. As we have noted in the Introduction, since clock cells in the SCN demonstrate self-sustained, nearly sinusoidal oscillations, one can model them as weakly nonlinear oscillators with a stable limit cycle (Liu et al., 1997). The dynamics of this kind of oscillators is described by the Kuramoto model. Let us consider a generalized Kuramoto model where N heterogeneous phase oscillators are divided into M groups (or communities). Each oscillator i is characterized by a phase θ i m ( ) and has its own natural angular frequency ω i m ( ) , where m is the group index. Each group comprises N m oscillators, coupled with strength K mm within groups and strength K nm between groups. In a periodic external field with period T and angular frequency ω π F T = 2 / , the local field phase ω φ F m t + ( ) and strength F m ( ) can differ between groups. If interaction between the oscillators is absent, then the time evolution of individual Kuramoto oscillators is (1) In a system with interaction, the time evolution of the phase θ i m ( ) of each Kuramoto oscillator is The macroscopic state of each group is characterized by the synchronization index ρ m and group phase ψ m of the complex order parameter The amplitude (synchronization index) ρ m characterizes the extent of phase alignment between oscillators in group m and varies between 0 and 1 . When ρ m = 0 , oscillators within group m are in an asynchronous state, while ρ m = 1 corresponds to a completely synchronized state. The group phase ψ m is the average phase or direction of the oscillators.
Here, we restrict our attention to a 2-group version of equation (2), with which we will model the coreshell organization and dynamics of the SCN. As shown elsewhere, equation (2) can be reduced to explicit equations for the phase and amplitude of each group (Ott and Antonsen, 2008;Yoon et al., 2021). Following this approach, the time dependence of the amplitude and phase in the core ( υ ) and shell ( d ) becomes (1 ) ( ) , coupling between core and shell oscillators is determined by intercouplings K ud (core on shell) and K du (shell on core).
in the frame rotating at the external field frequency ω F . The field phase is set to 0 since it provides a reference time for the external cue. Equations (4)-(7) describe the core-shell model schematically depicted in Figure 1, where the parameters characterizing the heterogeneity of dynamical properties among SCN oscillators are reduced to a set of parameters characterizing the mean properties of the core and shell. These include the characteristic mean ω υ ( c d ) and half width at half maximum ∆ υ (∆ d ) of the natural frequency distribution in the core (shell). At a given field frequency, the mean natural frequency determines the detuning Ω υ υ = ω ω F − in the core and Ω d F d = ω ω − in the shell. Note that a periodic external field is taken explicitly into account by equations (4)-(7), which allows to study nonlinear impacts of the external field on the SCN. At K K d d υ υ = =0 , these equations are reduced to the explicit equations for the forced Kuramoto model derived in Childs and Strogatz (2008).
In exchange for the simplification of equation (2) into equations (4)-(7), the natural frequencies or free-running periods of SCN oscillators are assumed to follow a Cauchy-Lorentz distribution. Despite the restriction, the Cauchy-Lorentz distribution is unimodal and symmetric, like the Gaussian distribution, which has been used to model the distribution of free-running periods in the core and shell (Taylor et al., 2017). From the physical point of view, the Kuramoto model demonstrates qualitatively similar behavior for both the Cauchy-Lorentz and Gaussian distributions of natural frequencies. The main difference is in their analytical properties. The Cauchy-Lorentz distribution has residuals in the complex plane. This fact simplifies analytical calculations and allows an explicit reduction of microscopic equation (2) to a set of equations (4)-(7) that describe the evolution of the macroscopic characteristics of the Kuramoto model.
Note that stationary values of the synchronization indices ρ υ and ρ d can be found analytically by solving equations (4)-(7) at zero time derivatives on the left hand side. Unfortunately, the strong non-linearity of these equations makes it difficult to find a simple analytical expression, while numerical solutions can be readily computed.
Recently, a model similar to our model in Figure 1 was considered by Hannay et al. (2020). Using the Ott-Antonsen ansatz, Hannay et al. (2020) derived equations similar to equations (4)-(7) to describe the evolution of the core and shell. While Hannay et al. (2020) focused on seasonal day length, after-effects of light entrainment, and seasonal variations in light sensitivity in the mammalian circadian clock, in this article, we mainly focused on (1) dynamical behavior of the core and shell under LD, DD, and LL conditions, including anticipation in nocturnal animals; (2) dependence of the entrainment range on the properties of core and shell neurons; (3) the mechanisms of circadian rhythm dissociation with increasing (or decreasing) LD cycle period above (or below) a critical threshold; and (4) dissociation under LL conditions in nocturnal and diurnal animals. Another important difference between the work of Hannay et al. (2020) and our own lies in how the external cue (i.e., light) is taken into account. In our approach, a periodic external optic cue is explicitly taken into account by equations (4)-(7). This enables to study nonlinear effects in the impact of the cue on the SCN. Hannay et al. considered only a collective phase response on an external cue, assuming that a light perturbation shifts mean-field phases.
In summary, our approach is based on the observation that clock cells in the SCN are self-sustained, nearly sinusoidal oscillators with a stable limit cycle (Liu et al., 1997). Individual dynamics of this kind of oscillators can be approximated by a simple differential equation (1). Accounting for interactions between oscillators, we obtain equation (2), which provides a complete description at the microscopic (cell) level. Finally, using explicit mathematical methods, we reduce the microscopic equation to equations (4)-(7), which describe the macroscopic dynamics of the core and shell. The assumption, that the SCN cells are self-sustained and nearly sinusoidal oscillators, may not hold true in general. Account of processes that disturb these conditions and lead to non-sinusoidal individual dynamics is an open problem. In this regard, one possible line of research, beyond the scope of the present work, is to compare models with distinct clock cell dynamics.

Selecting Model Parameters
This section discusses the choice of parameters used in the numerical study of equations (4)-(7). For convenience, numerical calculations were performed using dimensionless parameters, with reference to frequency unit ∆ υ . This is equivalent to dividing both sides of equations (4)-(7) by ∆ υ , rescaling all parameters K K nm nm / ∆ υ → and time t t ∆ υ → . Dimensionless parameters are summarized in Table 1, along with their dimensional equivalents, and reference values are found in the literature.
The choice of ω υ( ) d and ∆υ( ) d is based on a study of mice exposed to symmetric 12-h LD cycles by Taylor   Table 1

Parameters
Reference Model Time-dimensional parameters include σ υ( ) d , the standard deviation of free-running periods reported by Taylor et al. (2017); τ υ( ) d , the mean free-running period of the core (shell). Dimensionless parameters are obtained by multiplying time-dimensional parameters and dividing frequency-dimensional parameters by the frequency unit ∆ υ υ υ τ = / 2 πσ . Frequency-dimensionless parameters include K dd υυ( ) , the intracoupling in the core (shell); K d υ and K dυ , the core-shell intercouplings; F , the strength of the light cue acting exclusively on the core. The dimensionless equivalents of the model parameters: , the mean free-running frequency of core (shell) oscillators; and ∆υ( ) d , the Lorentzian spread (half width at half maximum) of free-running frequencies in the core (shell). et al. (2017). In the present core-shell model, 12-h LD cycles correspond to sequential half cycles of the field with period T = 24 h . Based on experimental observations, Taylor et al. (2017) fitted the free-running periods of the core ( υ ) and shell ( d ) to Gaussian distributions with mean τ υ = 25.1h and τ d = 23.9 h , respectively, and standard deviation σ υ = 1.3 h and σ d = 1.9 h . The mean free-running periods can be directly converted to mean free-running frequencies ω π From the definition of standard deviation, and given τ υ σ , the standard deviation of the free-running frequencies is approximated . Although the present coreshell model is based on a Lorentzian distribution of free-running frequencies with half width at half maximum ∆υ( ) d , Gaussian and Lorentzian distributions both describe a system where most neurons have a free-running period approximately equal to the system average, and the number of neurons with larger or smaller than average free-running frequencies is symmetric and monotonically decreasing. In both cases, ∆υ( ) d and SD d υ( ) provide a measure of the statistical spread of free-running frequencies about the mean value, so we simply take ∆υ The coupling mechanisms between SCN neurons are an ongoing topic of research (see, for example, Tokuda et al., 2018;Hastings et al., 2018; and references therein). Although direct measurements of intra-and intercouplings K υυ , K dd , K d υ , and K dυ are not available in the literature, these may be inferred from experimental observations by treating coupling as a proxy for neuron connectivity/communication. When inferring parameters, we may also restrict ourselves to values for which the system can sense and adjust to small changes in the external cue. At one extreme, excessively large coupling , 1 ) makes the system insensitive to the external cue. At the other extreme, an external cue of excessively large strength ( F 1 ) makes it impossible for the system to adjust its state through coupling changes, such as may be induced by neuronal plasticity (Rohr et al., 2019). Thus, realistic couplings and cue strength should be selected from an intermediate range where SCN plasticity is possible.
From a functional point of view, studies of SCN slices in the mouse by Yamaguchi et al. (2003) and rat by Noguchi et al. (2004) have shown that the ventral SCN remains synchronized when isolated from the dorsal SCN. These observations suggest that coupling within the isolated ventral SCN (core) is large enough to support synchronization. By comparison, the isolated shell was only found to remain synchronized in the rat (Noguchi et al., 2004). From a structural point of view, dendritic arbors within the shell are sparse, whereas those within the core are more extensive and larger (Moore, 1996). This suggests that there are more synaptic contacts between neurons in the core than in the shell. Viewed together, functional and structural observations suggest that intra-core coupling K υυ is larger than the critical value 2∆ υ required to synchronize the isolated core and larger than the intra-shell coupling K dd . For simplicity, the intracouplings are estimated for the isolated core and shell under constant darkness. Although the intracoupling may vary in response to photic input if the core and shell are coupled, consequent changes in the state of the system may be captured by other model parameters, as discussed in the next paragraph. In the isolated core (shell) under constant darkness and F = 0 ), intracoupling increases with the fraction of synchronized oscillators ρ υ( ) d as follows: which can be seen from equations (4) and (6). For positive intracoupling, it follows that ρ ρ υ > d when- In the absence of suitable experimental data with which to infer ρ υ( ) d , we take ρ ρ υ = 2 d , for example, and refer to computational models for a reasonable upper bound on ρ υ . For instance, the theoretical model used by Taylor et al. (2017) results in an average synchronization index of 0.9 in the entrained SCN. Considering a slightly smaller value of ρ υ = 0.8 in the isolated core, the corresponding shell value is ρ d = 0.4 . In dimensionless units, equation (8) then yields K υυ = 5.6 and K dd = 4.0 . The critical value of the core (shell) intracoupling K c υυ = 2 ∆ v (K c dd = 2 ∆ d ) at which synchronization emerges follows directly from equation (8) . "Near-critical" coupling in the shell is compatible with the idea that variability between species and/or individuals can result in coupling changes around the critical value and may explain why the isolated shell has been observed to synchronize in rats (Noguchi et al., 2004) but not in mice (Yamaguchi et al., 2003).
Regarding the intercouplings K d υ and K dυ , experimental evidence shows that there are far more contacts made by neurons expressing vasoactive intestinal polypeptide (VIP) in the core onto neurons expressing arginine vasopressin (AVP) in the shell than the converse (Varadarajan et al., 2018). In addition, functional studies of resynchronization dynamics in the SCN have revealed that the core entrains the shell (Taylor et al., 2017). Viewed together, such findings support the hypothesis that signaling from the core to the shell is stronger than the reverse ( However, different contributions to coupling by different neurotransmitters and neuropeptides may affect the sign of the coupling, which impacts entrainment boundaries, the emergence of unstable states, and splitting of activity rhythms into 2 synchronized bouts cycling in antiphase (Daan and Berde, 1978;Oda and Friesen, 2002;Evans et al., 2005;Evans and Gorman, 2016). For instance, the neurotransmitter gamma-aminobutyric acid (GABA) has been found to both inhibit and promote synchrony in the SCN network depending on the state of the SCN (Evans et al., 2013). In an antiphase configuration (ψ π ψ d = + υ ) following prolonged light exposure, GABA works with the neuropeptide VIP to promote synchrony. In the steady-state induced by symmetric LD conditions, GABA counters the synchronizing effect of the neuropeptide VIP. Since experimental observations demonstrate that the core entrains the shell under symmetric LD conditions (Taylor et al., 2017), the inhibiting contribution of GABA signaling is assumed smaller than the synchronizing effect of VIP ( K d υ > 0 ). Reasonable bounds on the intercouplings can also be established through analysis of equations (4)- (7) under symmetric LD conditions, as discussed in the next section, and constant darkness, as discussed in the section on DD conditions. In particular, equation . Starting from these bounds and intracouplings K υυ = 5.6 and K dd = 4.0 (chosen in the previous paragraph), parameters K d υ , K dυ , F , and τ d were adjusted to ensure agreement with experimental observations of the difference in peak activity time between the core and the shell t t . As shown further below in equation (12), the difference between the peak activity time of the core t p υ ( ) and the shell t d p ( ) is directly related to the phase difference ψ ψ d − υ that characterizes the state of the entrained SCN network. The difference t t p d p υ ( ) ( ) − has been experimentally measured at 2.4 0.9 ± h in the entrained SCN under symmetric LD conditions, based on peak expression of the Per2 gene (Taylor et al., 2017).
Thus, the advantage of our approach based on the Kuramoto model is that it reduces the large number of parameters that characterize tens of thousands of neurons and their couplings within the SCN to only 9 parameters, which capture the mean characteristics of core and shell neurons (see Table 1): the free-running periods τ υ and τ d of the core and shell; the standard deviation of the free-running periods, σ υ and σ d ; the intracouplings K υυ and K dd within the core and shell; the intercouplings K d υ and K dυ between the core and shell; and the strength of the optic cue F . These parameters were obtained from an analysis of experimental data (Taylor et al., 2017). The remaining 4 parameters in Table 1 ( ω υ , ω d , ∆ υ , ∆ d ) are dimensionless equivalents of the model parameters. As we will show in the next sections, these 9 parameters describe the dynamical properties of the SCN under 3 lighting conditions (LD, DD, and LL), anticipation in the shell, and circadian rhythm dissociation in response both to varying LD period and to varying light intensity under LL conditions, in qualitative and quantitative agreement with experimental observations for diurnal and nocturnal mammals.

core-Shell Phase Difference and anticipation In The Entrained State
The SCN can trigger behavior in anticipation of regular events. Examples include water-seeking in advance of sleep time and food-seeking in advance of meal times Gizowski et al., 2016). From an experimental point of view, existing evidence supports a model of the SCN where the core receives a photic cue and the shell is responsible for transmitting phase information to the wider circadian system. For example, anatomical evidence shows that vasopressin (VP) neurons in the shell provide output to downstream tissues (see Gizowski et al., 2016, and references therein), and functional evidence shows that the shell entrains the phase of cellular oscillators in downstream tissues (Evans et al., 2015).
This section investigates how the shell of the entrained SCN can provide a reference phase for anticipation. To this end, an analytical and numerical analysis of equations (4)- (7) is presented, using the parameters in Table 1. Entrainment corresponds to a stationary solution of equations (4)- (7), which describes entrainment in a frame rotating at the frequency of the external cue ω π F T = 2 / . In an entrained system, the phase difference ∆ψ ψ ψ ≡ − d υ between shell and core groups of oscillators is readily obtained from equation (7): When the synchronization index of the entrained SCN is large ( ρ ρ υ ≈ ≈ d 1 ), or more generally when 2 ( 1 ) Under these conditions, positive K d υ must be larger than the difference between the mean natural frequency of the shell ω d and the frequency of the external cue ω F : 1 . Direct inspection of equation (9) also shows that the sign of the shell-core phase difference is determined by the sign of the shell's detuning parameter Ω d Since the mean free-running period of the shell τ d = 23.9 h is smaller than the period of the external cue T = 24 h , it follows that The extent to which ψ d leads ψ υ is determined by the period of LD cycles T F = / π ω . As hinted by equation (9), the phase difference − ≤ ≤ π ψ π / 2 / 2 ∆ tends to decrease as T decreases below 24 h, toward the shell's free-running period τ d = 23.3 h , and tends to increase as T increases above 24 h . In general, the entrained phase difference is also dependent on how the synchronization indices ρ υ and ρ d vary with T.
In the laboratory (non-rotating) frame, the real part of core (shell) order parameter z d υ( ) in equation (3) oscillates with frequency ω π . This time-dependent observable acts as a proxy for SCN activity, encoding synchronization index and phase at a given cue frequency or period. Note that the peak in core activity is locked to the peak in the external cue because the photic cue in our core-shell model acts exclusively on the core.
The observable Re( , where n is an integer number. From here, it follows that the difference in peak activity time between the shell and the core is Thus, the peak time difference t t Using equation (10), we find how far ahead in time is the shell activity with respect to the core activity in the dependence on the frequency of the endogenous cycles in the shell, ω d , and the coupling K d υ : As discussed in the preceding paragraph, ψ d leads ψ υ for positive K d υ , since τ d T < and therefore ω ω d F − > 0 . From equation (12), it then follows that with increasing time, the core must peak after the shell, that is, t t p d p υ ( ) ( ) > (see Figure 2b). At that time, the peak of the core activity is locked to the peak of the external cue. This result agrees with the experimental observation that peak Per2 expression in the shell precedes peak Per2 expression in the core by 2.4 0.9 ± h (Taylor et al., 2017). Viewed together, the results presented in this section show that the core-shell organization of the SCN enables the shell to anticipate events by advancing its phase relative to the phase of the core and the external cue while entrained to a symmetric LD cycle of period T . This lead is determined by the core-shell phase difference. For example, the numerical solution of equations (4)-(7) with the parameters in Table 1  Moreover, the shell can also increase its lead over the core in response to increasing T , as shown in Figure 2a.
To outline the importance of the anticipation, which is formed by the shell, we want to note that the SCN establishes phase coherence between selfsustained and cell-autonomous oscillators in peripheral organs, for example, liver, muscle, pancreas, heart, adipose tissue, and other areas of the brain (e.g., pineal gland which produces a hormone [melatonin] and modulates sleep-wake cycles). It is important to note that the organs receive signals mainly from the shell. This gives them the possibility to be prepared in advance for a change of activity. These organs are also influenced by other cues such as feeding/fasting cycles.

Entrainment Range and Dissociation Under LD conditions
Experimental investigations have revealed that animals can only entrain to a limited range of symmetric LD periods, known as the entrainment range: LLE T ULE < < , where LLE and ULE are the lower and upper limit of entrainment, respectively (Pittendrigh and Daan, 1976a;Campuzano et al., 1998;Usui et al., 2000;de la Iglesia et al., 2004;Goldbeter and Leloup, 2021). In nocturnal rodents, the entrainment range spans from LLE ∼ ∼ 22 h up to ULE ∼ ∼ 28.5 h , varying between species and individuals (Pittendrigh and Daan, 1976a;Campuzano et al., 1998;Usui et al., 2000). Beyond the entrainment range, it has been observed that animal behavior and SCN activity follow an additional rhythm, which differs from that of the external light cue. This phenomenon is known as dissociation (Campuzano et al., 1998;Usui et al., 2000;de la Iglesia et al., 2004;Schwartz et al., 2009).
This section studies dissociation from the point of view of an observer in the laboratory frame. To this end, it is assumed that the observer measures the observable Re( or the corresponding distribution of frequency components S d υ( ) ( ) ω (spectral density), as a proxy for SCN activity under varying LD period T . Here, the core (shell) synchronization index ρ υ( ) d and phase ψ υ( ) d are determined by numerical solution of equations (4)-(7), using the parameters in Table 1. The main results concerning independent components of the spectral density S d υ( ) ( ) ω are presented in Figure 3. First, Figure 3a shows that a second rhythm ( sr ) with period T sr emerges outside the entrainment range. Above the upper bound ( ULE = 25.28 h ), T sr is smaller than T and decreases with increasing T . Below the lower bound ( LLE = 23.26 h ), T sr is larger than T and increases with decreasing T . These results are in qualitative agreement with experimental observations in rats (Campuzano et al., 1998;Usui et al., 2000). For instance, Usui et al. (2000) studied recurring drinking behavior in rats exposed to LD cycles with period T and report the emergence of a second recurring peak in drinking activity with period T sr = 25 h at T ULE = =28.5 h (see Figure 3 by Usui et al., 2000), and T sr = 23.5 h at T LLE = =23 h (see Figure 4 by Usui et al., 2000). The authors also report that T sr increases with decreasing T below LLE and decreases with increasing T above ULE , all in agreement with the numerical results of Figure 3a. Another study reports a similar trend in motor activity in rats under shortening LD cycles (Campuzano et al., 1998), as shown in Figure 3a (see black crosses). In this figure, we compared our results on the period of the second (dissociated) rhythm with the experimental data from Campuzano et al. (1998). Despite the clear separation between calculated and experimental curves, the trend is generally the same, and the relative difference between the numerical result and experiment is only 4% .
Second, as shown in Figure 3b, the entrainment range is characterized by a single frequency component of intensity I I in the core and shell. This component corresponds to an entrained rhythm with period T between LLE = 23.26 h (lower bound) and ULE = 25.28 h (upper bound), as shown in Figure 3a.
Third, analysis of the intensity of core and shell components reveals distinct contributions to the second rhythm above ULE and below LLE , as shown in Figure 3b. Below LLE , core and shell contributions to the second rhythm at period T sr increase, while contributions to the rhythm at period T decrease. However, when the LD period T F is lengthened above ULE , the core's contribution is almost exclusively to the component at the LD period (near constant I F ( ) υ and vanishingly small I sr ( ) υ ), in qualitative agreement with experimental observations made by Schwartz et al. (2009). This difference is related to the coexistence of dynamically and topologically distinct states in the core and shell at T ULE > (these states were discussed in Yoon et al. (2021)). The core is on average locked to the LD cue, while the shell is drifting at the second rhythm. Below LLE , both the core and shell drift at the second rhythm. In this particular case (set of model parameters), the emergence of distinct states in the SCN is a consequence of distinct types of bifurcations in equations (4)-(7) at LLE and ULE . The bifurcation at LLE is characterized by the emergence of a second rhythm with period T T sr ∼ ∼ , as shown in Figure 3a. The bifurcation at ULE , however, is characterized by a large drop in T sr relative to T ULE = . These characteristic features of distinct bifurcations are also consistent with the experimental data discussed in the previous paragraph. For example, the data reported by Usui et al. (2000) show that LLE T sr − = 0.5 h , whereas T ULE sr − = 3.5 h for rats. For comparison, our model predicts LLE T sr − = 0.2 h and T ULE sr − = 1.7 h for the model parameters in Table 1. To identify the type of bifurcation, we also carried out the stability analysis of fixed points of equations (4)-(7). This analysis revealed that the dissociation at ULE = 25.28 h is caused by the super Hopf bifurcation, while the saddle-node bifurcation occurs at LLE = 23.26 h .
The results presented in this section are in good qualitative agreement with experimental observations of dissociation in the SCN. In the core-shell model of equations (4)- (7), the emergence of a dissociated rhythm takes place when entrainment is disrupted. The dynamical features and bifurcation mechanisms of disrupted states in the core and shell describe key aspects of dissociation, namely, the inverse dependence of the dissociated period on the LD period.
Notably, the upper and lower bounds of the entrainment range are very close to the mean freerunning period of the core τ υ = 25.1h and shell τ d = 23.3 h for the model parameters in Table 1: To determine the dependence of the entrainment range on the endogenous core and shell rhythms, we performed additional numerical calculations of equations (4)-(7) and analyzed the dependence of the entrainment range on the difference between the freerunning periods, τ τ υ − d , at fixed τ υ and increasing τ d . The results of these numerical calculations are presented in Figure 3c and show that decreasing the difference τ τ υ − d increases the entrainment range. Simultaneously, equation (13) shows that the difference between peaks in the core and shell, t t p d p υ ( ) ( ) − (anticipation), decreases with increasing τ d (i.e., decreasing ω d ). Therefore, one can have a broad entrainment range but weak anticipation, and vice versa.

DD conditions
This section studies the behavior of the core and shell under DD conditions and how core-shell intercouplings K d υ and K dυ are related to τ DD across the whole SCN. Since τ DD , τ υ , and τ d can be measured experimentally under DD conditions, the aim is to identify any experimentally testable relationships that reveal information about the value and sign of intercouplings K d υ and K dυ . The free-running frequency ω π DD DD = 2 / τ of the synchronized SCN under DD conditions can be obtained self-consistently from the stationary state of equations (4)-(7), by setting F = 0 and replacing ω ω  Table 1.

(a) Period T sr of independent Scn activity components under external light-dark cycles of varying period T (in hours). The component at the field period T is indicated in red open circles within the entrainment range and in gray circles outside the entrainment range. Red stars indicate the period of Scn activity at a second (dissociated) rhythm T sr . Black crosses are experimental data reported by campuzano et al. (1998). The dashed horizontal line is the period of Scn oscillations τ DD = 24.84 h under dark-dark conditions (cue strength F = 0), calculated from the mean value of equations (5) and (7) in the synchronized state. (b) Intensity of light-dark (F , hollow markers) and second rhythm (sr , solid makers) components of Scn activity, in the core (υ, circles) and shell (d, squares). Vertical lines indicate the lower and upper bounds
The synchronized SCN is characterized by core and shell synchronization indices ρ υ and ρ d . To find τ DD , we solve numerically equations (4)-(7) at DD conditions with parameters in Table 1. Then we apply the Fourier analysis to the time dependence of real parts, Re[ , of the order parameters (equation (3)) of the core and shell. It gives the frequency of steady oscillations at DD conditions. Motivated by the experimental observation in Taylor et al. (2017), we assume that τ τ υ > d (see Table 1). In the case for positive intercouplings K d υ and K dυ , the core-shell model of equations (4)-(7) predicts that ω ω ω υ < < DD d . Equivalently, the period τ DD DD = 2 / π ω must be in the range Using the parameters in Table 1, equation (15) tells us that τ DD = 24.8 h , in good agreement with experimental observations in rats by Campuzano et al. (1998), who reported τ DD = 24.5 h .
Under DD conditions, the analytical result of equation (18) predicts that τ DD is bounded between τ υ and τ d . This is valid for positive inter-and intracouplings. However, in the case of negative couplings, τ DD is no longer bound between τ υ and τ d . In the particular case when K dυ < 0 (and K d υ ≥ 0 ), equation (15) becomes

LL conditions
This section studies the SCN under constant light or LL conditions for diurnal and nocturnal animals. Within our approach, LL conditions can be seen as the limit when LD cycles become so slow that they are effectively constant from the point of view of individual oscillators in the core. In turn, this motivates the replacement of the periodic forcing term in equation (2) with a constant B : In our reduced Kuramoto model, parameter B is phenomenological. It can be positive or negative. We assume that the magnitude | | B characterizes light intensity and chooses the sign of B to satisfy Aschoff's first rule. From this, it follows that B is negative for nocturnal animals and positive for diurnal animals, as will be proven below. This choice of parameter B allows us to describe the basic properties of the SCN under LL conditions, in nocturnal and diurnal animals.
As one can see in equation (2), introducing parameter B renormalizes the endogenous free-running frequencies of core oscillators from which it follows that the renormalized mean frequency of the core is ω ω and, equivalently, the renormalized mean period of oscillations is where τ υ is the endogenous mean free-running period of isolated core oscillators. This shows that under LL conditions, the positive parameter B > 0 (diurnal animals) forces core oscillators to run faster than under DD conditions ( ω ω υ υ * > ). Conversely, if parameter B is negative, −2 < < 0 ω υ B , (nocturnal animals), then the constant light cue forces core oscillators to run slower than under DD conditions ( ω ω υ υ * < ). The free-running period of isolated shell oscillators, on the other hand, remains unaffected ( τ τ d d * = ), since it is not directly exposed to the light cue. It would be interesting to check experimentally this effect of LL conditions on isolated clock cells in the core.
To determine τ LL in our core-shell model, we must account for the intercoupling between the core and shell. The dynamics of the core-shell system under LL conditions are described by the same equations (4)-(7) as DD conditions ( F = 0 ), by renormalizing the mean period as in equation (28). This means that we can use the equations derived for DD conditions by making the simple replacements τ τ υ υ → * and τ τ DD LL → . From equation (20) in particular, it then follows that the frequency ω LL of the steady circadian rhythm is Based on experimental observations of different animal species, Aschoff's first rule predicts that τ τ LL DD < in diurnal animals and τ τ LL DD > in nocturnal animals, where τ DD and τ LL are the free-running circadian periods under DD and LL conditions, respectively (Aschoff, 1965(Aschoff, , 1979Pittendrigh and Daan, 1976a;Carpenter and Grossberg, 1984). ) in nocturnal animals. Moreover, Aschoff's first rule also states that the effects of LL conditions are intensity-dependent (Aschoff, 1965;Pittendrigh and Daan, 1976a;Carpenter and Grossberg, 1984). The free-running period of circadian activity rhythms usually decreases with increasing light intensity: the brighter the constant light, the faster the animal's clock runs. In nocturnal animals, the converse is usually observed: the free-running period of the rhythm increases with increasing light intensity. Within the model of equation (25) and equations (4)-(7), the parameter | | B characterizes the intensity of the optic cue. For diurnal animals ( B > 0 ), increasing | | B decreases τ LL (increases ω LL ), whereas for nocturnal animals ( B < 0 ), increasing | | B increases τ LL (decreases ω LL ), in agreement with Aschoff's first rule.
Next, this section explores the possibility of dissociation under increasing light intensity. To this end, the spectral analysis of the section on dissociation under LD conditions was repeated after numerically solving the model of equation (25) and equations (4)-(7), using the parameters in Table 1. These parameters were determined in reference to experimental observations of mice, which are nocturnal ( B < 0 ). The case B > 0 in Figure 4 describes a "fictitious diurnal mouse". We use the notion "fictitious diurnal mouse" to outline that parameter B is positive for diurnal animals, although the remaining model parameters (Table 1) used throughout this article are for nocturnal mice. The result presented in Figure 4 shows that LL conditions slow down the circadian clock (i.e., τ LL is lengthened), shortens the daily active phase, and reduces total daily activity in nocturnal rodents. In diurnal animals, including humans, LL conditions generally have the opposite effect (Mistlberger and Rusak, 2005), although there are some exceptions across species (diurnal primates in particular).
Spectral analysis revealed the emergence of a second rhythm with period T LL sr ( ) at a critical brightness | | B c , with an order of magnitude difference between nocturnal and diurnal cases. In the nocturnal case, | | 0.24 B c , for which τ LL = 25.2 h , whereas in the diurnal case, | | 3.23 B c , for which τ LL = 21.7 h . These results are in qualitative agreement with the experimental observation of rhythm dissociation at constant dim illumination of approximately 4.5 lux in rats for which τ LL = 25.15 h following 3 months of exposure (Albers et al., 1981). Unfortunately, due to a lack of detailed experimental data, we are unable to present a detailed numerical comparison as a function of illumination. However, we expect that the dependence of dissociated rhythm period, T LL sr ( ) , on light intensity | | B for diurnal and nocturnal animals can be checked experimentally. In addition, Figure 4 shows that with increasing the light intensity, | | B , the period of the dissociated rhythm decreases a little in the nocturnal case but slightly increases in the diurnal case. Unlike the free-running period τ LL , the dissociated rhythm T LL sr ( ) does not obey Aschoff's first rule. Finally, the model of equation (25) and equations (4)-(7) predicts that the light intensity, | | B , can also alter the core-shell phase difference under LL conditions compared with DD conditions. Under LL conditions, the phase difference ψ ψ d − υ can be determined from which follows from the replacement ω ω υ υ → * in equation ( 22 ). For diurnal animals ( B > 0 ), increasing the light intensity decreases the core-shell phase difference in animals where τ τ υ d < and both intercouplings are positive. In nocturnal animals, increasing the light intensity increases the core-shell phase difference. In both cases, sufficiently large light intensity, | | B , will invert the core-shell phase relationship, causing the core to lead the shell. Under these conditions, the shell cannot provide a reference phase for anticipation in peripheral clocks.

DIScUSSIOn
In this work, we studied how the core-shell organization controls the behavior of the SCN, the entrainment of the core and shell to the environmental cue, the mechanism of disruption of the synchronized state, and the impact of different light conditions on dynamics of the SCN. For this purpose, we used a core-shell model of the SCN based on reduced dynamical equations for the forced Kuramoto model. Our approach is based on the observation that clock cells in the SCN are self-sustained, nearly sinusoidal oscillators with a stable limit cycle. The first benefit of our model is that it only has 9 biologically meaningful parameters, instead of thousands of parameters for an equally large number of neurons. Second, the proposed dynamical equations for core and shell rhythms describe known SCN activity under 3 lighting conditions (LD, DD, and LL), including the anticipation and the dissociation of the circadian rhythms in the dependence on the period of external optic cue and light intensity. Third, the model can be calibrated for diurnal and nocturnal mammals. In this article, we calibrated the model by using parameters for mice. Because these nocturnal animals are good model animals, there are numerous experimental studies on the behavior of their circadian rhythms  Table 1). The period τ LL of the free-running component is indicated by circles, and the period T LL sr ( ) of a second (dissociated) component is indicated by navy stars. negative B corresponds to the nocturnal mouse (see Table 1  under different light conditions, and reliable quantitative and qualitative results were obtained. Fourth, since our model is based on explicit equations, it allows us to perform both an analytical analysis and detailed numerical comparison between the experiment and theory concerning the temporal behavior of the SCN as a whole and its principal modules, the core and shell, at different light conditions, in the entrained state and the states with dissociated rhythms. As far as we know, this detailed comparison was not performed within the existing models.
An important difference between the core and shell is that neuron populations within each subdivision have distinct mean free-running periods τ υ and τ d , respectively. The numerical results present in Figure 3 suggest that τ υ and τ d are important parameters that determine the entrainment range and the core-shell phase difference under LD conditions. From an evolutionary point of view, this suggests that the differentiation between endogenous rhythms in the core and shell is an adaptation to environmental variations in LD cycle length.
Under DD conditions, the analytical results of equation (18) predict that the SCN's free-running period τ DD is bounded between τ υ and τ d . Based on experimental evidence that τ d < 24 h and τ υ > 24 h (Taylor et al., 2017), equation (18) predicts that the core-shell SCN can maintain near-24-h periods of activity in the absence of an external light cue. Our numerical result, τ DD = 24.8 h , obtained for the model parameters calibrated for mice, agrees with Campuzano et al. (1998), who reported τ DD = 24.5 h for rats. From a biological point of view, the ability to retain circadian activity in constant darkness ensures minimum disruption of circadian activity. However, this ability is also dependent on the relative strength of core-shell communication, parametrized by intercouplings K d υ and K dυ . As discussed in the section on DD conditions, a large asymmetry in the strength of the intercouplings, which parametrize core-shell communication, is expected to shift the free-running period under constant darkness closer to the endogenous rhythm of the core ( The period of an external cue (LD cycles) as the control parameter can be easily adjusted in laboratory conditions as in many experimental investigations. Here, we demonstrated appearance of dissociation when the period of the LD cycle is larger (or smaller) than a critical value. There are also other important parameters, such as couplings between oscillators and mean free-running periods in the core and shell. Changes of these parameters also may result in the dissociation of circadian rhythms. The impact of abnormal changes of these parameters on the SCN dynamics is an important problem in experimental research. External factors, such as temperature or continuous intake of chemical substances, for example, melatonin or antidepressants, may change these parameters. In particular, the couplings among neurons are mediated by neurotransmitters, such as VIP and GABA, which can be influenced by drugs. The appearance of dissociated rhythms is a common phenomenon in the SCN of mammals. We believe that our results may also be translated to humans being diurnal animals. However, these interesting problems need a detailed literature search and further investigations that go out of the scope of the present article.
Our model shows that the SCN's functional separation into a core and shell enables anticipation, the ability of an organism to trigger physiological changes and behavior in anticipation of regular events. As argued in the section on anticipation in the entrained state, anticipation is supported by the coreshell phase difference ψ ψ d − υ , which is proportional to the difference in peak activity time between the shell and the core. Experimental evidence shows that the shell coordinates the phase of tissue clocks throughout the brain and body (Evans et al., 2015;Coomans et al., 2015;Gizowski et al., 2016), and peak activity in the shell precedes peak activity in the core (Taylor et al., 2017) under LD conditions. This corresponds to a situation where the shell phase ψ d leads the core phase, ψ υ , which is locked to the cue phase and acts as a reference cue for anticipation in peripheral clocks. Under LD cycles with angular frequency ω π F T = 2 / , the core-shell phase difference is determined by the relative magnitude and sign of intercoupling K d υ and the shell detuning ω ω d F − . In this case, the shell leads the core provided that communication from the core has a synchronizing effect ( K d υ > 0 ). For a given LD cycle with period T , the extent of anticipation (phase difference) varies inversely with the strength of K d υ . Under DD conditions, anticipation becomes dependent on the relative magnitude and sign of the difference τ τ υ − d and the total intercoupling K K d d υ υ + , as shown in the section on DD conditions. In this case, the shell phase can also precede the phase in the core, even when one of the intercouplings has a weak desynchronizing effect, and the extent of anticipation varies inversely with the total intercoupling.
We also considered the behavior of the SCN under LL condition by introducing a phenomenological parameter B which can be positive or negative. Its magnitude | | B characterizes the light intensity. By choosing the sign of B to satisfy Aschoff's first rule, negative B for nocturnal animals and positive B for diurnal animals, we described the basic properties of the SCN under LL conditions for nocturnal and diurnal animals in qualitative and quantitative agreements with observations. Compared to DD conditions, the analytical results for LL conditions predict that a constant light cue with intensity | | B decreases the core-shell phase difference in diurnal animals ( B > 0 ) but has the opposite effect in nocturnal animals ( B < 0 ). However, sufficiently large light intensity causes the shell to lag the core and is therefore expected to make anticipation impossible in both nocturnal and diurnal animals.
At the present time, there is surprisingly little research aimed specifically at determining the physiological, anatomical, or molecular mechanisms underlying Aschoff's rules. Genetic components of Aschoff's first rule were discussed by Muñoz et al. (2005). These authors found that LL lengthens the circadian period by inhibiting the normal dark-induced degradation of mPER2, and constitutively elevated levels of mPER2 act to enhance the phase-delaying limb of the molecular oscillator. Unfortunately, it is unclear how this inhibition-enhancement mechanism may relate to the sign of the phenomenological parameter B . It is an open problem for further experimental and theoretical investigations.
In this work, we also analyzed the difference in the collective behavior of the SCN oscillators under long and short photoperiod conditions. Operating limits on free-running and entrained SCN activity are partly determined by the same parameters as the core-shell phase difference ψ ψ d − υ . During steady-state activity, 0 <| | ψ ψ π d − ≤ υ , as can be inferred from the parameters in Table 1 and equations (9), (22), and (30). The numerical results of Figures 3 (LD conditions) and 4 (LL conditions) predict rhythm dissociation under varying LD period or the light intensity under LL. Under LD conditions, entrainment is ensured for any LD period T within the entrainment range LLE T ULE < < . Outside this range, a second rhythm appears, and its period T sr varies inversely with the LD period. At ULE , T sr is significantly smaller than T and decreases with increasing T . At LLE , T sr is infinitesimally larger than T and increases with decreasing T . These results are in good qualitative agreement with experimental observations of dissociation in rats (Campuzano et al., 1998;Usui et al., 2000). Interestingly, the spectral analysis presented in Figure 3b predicts that the dissociated rhythm at ULE is largely produced by the shell, unlike at LLE , where the core and shell contribute almost equally to the dissociated rhythm. Under LL conditions, dissociation appears under increasing light intensity | | B in both diurnal ( B > 0 ) and nocturnal animals ( B < 0 ), as shown in Figure 4. The period T LL sr ( ) of the dissociated rhythm is predicted to increase with light intensity in diurnal animals, but to decrease in nocturnal animals. In addition, the critical light intensity for dissociation was found to be an order of magnitude smaller in the nocturnal case, in agreement with the experimental observation of rhythm dissociation in rats under LL conditions (Albers et al., 1981).
In conclusion, the simplified model presented in this work captures important functional aspects of the SCN and its core-shell organization. Using this model, we extend existing work on dissociation by studying the period and strength of the dissociated rhythm. Numerical results were found to be in good qualitative agreement with experimental observations, in particular regarding the inverse relationship between the period of the LD cycle and the period of the dissociated rhythm. Within the model, the state of the core and shell is characterized by a synchronization index and phase, and determined by parameters which can either be measured experimentally or be inferred from experiments, as discussed in the section on the choice of parameter values. Among these parameters, the core-shell intercouplings and mean free-running periods constrain the core-shell phase difference, which determines the shell's ability to anticipate the core under different symmetric lighting conditions and intensities. In the particular case of constant lighting conditions, the dependence on light intensity can be modeled explicitly through the introduction of a phenomenological constant, as shown in the corresponding section. This simple modification to the model reproduces Aschoff's first rule for diurnal and nocturnal animals and reveals rhythm dissociation under increasing light intensity. Moreover, the introduction of the constant opens up the possibility to explore other open problems and experimental conditions in the study of circadian rhythms. Potential applications include asymmetric LD cycles, which are used to mimic seasonal variations in day length (Meijer et al., 2010;Myung et al., 2015, and the problem of antiphase rhythm splitting under constant light (Pittendrigh, 1960;Pittendrigh and Daan, 1976b;Ohta et al., 2005). From a theoretical point of view, another interesting problem is to extend the proposed model with other features of the SCN's network organization, namely, degree distribution (Gu et al., 2021), the existence of local groups of neurons that form intermediate structures in the SCN (Yoshikawa et al., 2021).
The results presented in this work strongly suggest that the Kuramoto model captures essential features of synchronization and entrainment in the SCN. Moreover, the reduced model is easily extendible to an arbitrary number of groups, with dynamics described by explicit equations for the group phase and synchronization index. Viewed together, the reduced Kuramoto model presents itself as a useful tool for exploring open problems in the study of circadian rhythms, one that can account for evolving views of the circadian system's organization, including inter-hemispheric interaction and peripheral clocks (e.g., in liver, muscle, pancreas, heart, adipose tissue). We believe that our model may be translated to the other animals, both nocturnal and diurnal animals, including humans.

acKnOwLEDGMEnTS
This work is funded by national funds (OE), through Portugal's FCT Fundação para a Ciência e Tecnologia, I.P., within the scope of the framework contract foreseen in paragraphs 4, 5, and 6 of article 23, of Decree-Law 57/2016, of 29 August, and amended by Law 57/2017, of 19 July. E.A.P. Wright acknowledges the financial support provided by FCT under PhD grant SFRH/ BD/121331/2016.