Generation and Disruption of Circadian Rhythms in the Suprachiasmatic Nucleus: A Core-Shell 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 Nm oscillators, coupled with strength Kmm within groups and strength Knm between groups. In a periodic external field with period T and angular frequency ωF=2π/T, the local field phase ωFt+ϕ(m) and strength F(m) can differ between groups. If interaction between the oscillators is absent, then the time evolution of individual Kuramoto oscillators is

In a system with interaction, the time evolution of the phaseof each Kuramoto oscillator is θ i ( ) m

The macroscopic state of each group is characterized by the synchronization indexand group phaseof the complex order parameter ρ m ψ m

The amplitude (synchronization index)characterizes the extent of phase alignment between oscillators in groupand varies betweenand. When, oscillators within groupare in an asynchronous state, whilecorresponds to a completely synchronized state. The group phaseis the average phase or direction of the oscillators. ρ m m 0 1 ρ m = 0 m ρ m = 1 ψ m

Here, we restrict our attention to a 2-group version of equation (2), with which we will model the core-shell 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

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 ω¯υ (cd) 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υd=Kdυ=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 indicesandcan be found analytically by solving-() 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. ρ υ ρ d equations (4) 7

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.

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 Knm/Δυ→Knm 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 et al. (2017). In the present core-shell model, 12-h LD cycles correspond to sequential half cycles of the field with period T=24h. 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.9h, respectively, and standard deviation συ=1.3h and σd=1.9h. The mean free-running periods can be directly converted to mean free-running frequencies ω¯υ(d)=2π/τυ(d). From the definition of standard deviation, and given τυ(d)≫συ(d), the standard deviation of the free-running frequencies is approximated as SDυ(d)≈2πσυ(d)/τυ(d)2. Although the present core-shell 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 ∆υ(d)=SDυ(d).

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υυ, Kdd, Kυd, and Kdυ 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 (Kυυ,Kdd,Kυd,Kdυ≫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 Kdd. 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 (Kυd=Kdυ=0 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 whenever Kυυ/Kdd>∆_v/∆_d, which is the case since Kυυ>Kdd and ∆_d > ∆_v. 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 Kdd=4.0. The critical value of the core (shell) intracoupling Kcυυ=2∆_v (Kcdd=2∆_d) at which synchronization emerges follows directly from equation (8) in the limit where ρυ(d)=0. Thus, the chosen intracoupling is considerably closer to the critical value Kcdd=3.4 in the shell, where (Kdd−Kcdd)/Kcdd=0.18, than to the critical value Kcυυ=2.0 in the core, where (Kυυ−Kcυυ)/Kcυυ=1.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 Kdυ, 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 (Kυd>Kdυ). 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 (11) shows that Kυd>ω¯d−ωF=0.15, and equation (24) shows that Kυd+Kdυ>ω¯d−ω¯υ=0.97. Starting from these bounds and intracouplings Kυυ=5.6 and Kdd=4.0 (chosen in the previous paragraph), parameters Kυd, Kdυ, 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υ(p)−td(p). As shown further below in equation (12), the difference between the peak activity time of the core tυ(p) and the shell td(p) is directly related to the phase difference ψd−ψυ that characterizes the state of the entrained SCN network. The difference tυ(p)−td(p) has been experimentally measured at 2.4±0.9h in the entrained SCN under symmetric LD conditions, based on peak expression of the Per2 gene (Taylor et al., 2017). Parameters τd, Kυd, Kdυ, and F were adjusted by inspection, until tυ(p)−td(p)=2.3h in the entrained state (T=24h), for τd=23.3h, Kυd=1.1, Kdυ=0.5, and F=1.5 (dimensionless units).

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 Kdd within the core and shell; the intercouplings Kυd and Kdυ 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.

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 (Silver, 2018; 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=2π/T. 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 (), or more generally when, it follows that ρ υ ρ d ≈ ≈ 1 2 1 ρ d ρ υ ≈ ( + ) ρ d 2

Under these conditions, positivemust be larger than the difference between the mean natural frequency of the shelland the frequency of the external cue: K υ d ω ¯ d ω F

since. Direct inspection ofalso shows that the sign of the shell-core phase difference is determined by the sign of the shell’s detuning parameterand of the intercoupling. Since the mean free-running period of the shellis smaller than the period of the external cue, it follows that. In this case, the phase differenceis positive sinceis also positive (). The extent to whichleadsis determined by the period of LD cycles. As hinted by, the phase differencetends to decrease asdecreases belowh, toward the shell’s free-running period, and tends to increase asincreases above. In general, the entrained phase difference is also dependent on how the synchronization indicesandvary with. | ( − ) | ≤ sin ψ d ψ υ 1 Ω d ω F ω ¯ d = − K υ d τ d = . 23 9 h T h = 24 Ω d τ d = / − / < 2 2 π T π Δ ψ = − ψ d ψ υ K υ d K υ d = . 1 1 ψ d ψ υ T π = / ω F − / ≤ ≤ / π Δ ψ π 2 2 T 24 τ d = . 23 3 h T 24 h ρ υ ρ d T equation (9) equation (9)

In the laboratory (non-rotating) frame, the real part of core (shell) order parameterinoscillates with frequency. Symbolically,. 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. z υ d ( ) ω F = / 2 π T Re( ) cos t z υ d ( ) ρ υ d ( ) ω F ψ υ d ( ) = ( + ) equation (3)

The observableis maximum in the core (shell) at time, whereis an integer number. From here, it follows that the difference in peak activity time between the shell and the core is Re( ) z υ d ( ) t υ d ( ) ( ) p = / − / 2 π n ω F ψ υ d ( ) ω F n

Thus, the peak time differenceis determined by the phase difference. Using, 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,, and the coupling: t υ ( ) p t d ( ) p − ψ d ψ υ − ω d K υ d equation (10)

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υ(p)>td(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.9h (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 shows that, for T=24 h, the lag tυ(p)−td(p)=2.3h corresponds to the phase difference ψd−ψυ=0.607rad≈35degrees. 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 self-sustained 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.

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∼∼22h up to ULE∼∼28.5h, 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(zυ(d))=ρυ(d)cos(ωFt+ψυ(d)) 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 Tsr emerges outside the entrainment range. Above the upper bound (ULE=25.28h), Tsr is smaller than T and decreases with increasing T. Below the lower bound (LLE=23.26h), Tsr 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 Tsr=25h at T=ULE=28.5h (see Figure 3 by Usui et al., 2000), and Tsr=23.5h at T=LLE=23h (see Figure 4 by Usui et al., 2000). The authors also report that Tsr 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 IF(υ)=IF(d) in the core and shell. This component corresponds to an entrained rhythm with period T between LLE=23.26h (lower bound) and ULE=25.28h (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 Tsr increase, while contributions to the rhythm at period T decrease. However, when the LD period TF is lengthened above ULE, the core’s contribution is almost exclusively to the component at the LD period (near constant IF(υ) and vanishingly small Isr(υ)), 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 Tsr∼∼T, as shown in Figure 3a. The bifurcation at ULE, however, is characterized by a large drop in Tsr 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−Tsr=0.5h, whereas Tsr−ULE=3.5h for rats. For comparison, our model predicts LLE−Tsr=0.2h and Tsr−ULE=1.7h 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.28h is caused by the super Hopf bifurcation, while the saddle-node bifurcation occurs at LLE=23.26h.

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-(), 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. equations (4) 7

Notably, the upper and lower bounds of the entrainment range are very close to the mean free-running period of the core τυ=25.1h and shell τd=23.3h 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 free-running 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υ(p)−td(p) (anticipation), decreases with increasing τd (i.e., decreasing ωd). Therefore, one can have a broad entrainment range but weak anticipation, and vice versa.

This section studies the behavior of the core and shell under DD conditions and how core-shell intercouplingsandare related toacross the whole SCN. Since,, andcan be measured experimentally under DD conditions, the aim is to identify any experimentally testable relationships that reveal information about the value and sign of intercouplingsand. K υ d K d υ τ D D τ D D τ υ τ d K υ d K d υ

The free-running frequencyof the synchronized SCN under DD conditions can be obtained self-consistently from the stationary state of-(), by settingand replacing, from which it follows that ω D D τ D D = / 2 π F = 0 ω F ω D D → equations (4) 7

where

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[zυ(d)(t)]=ρυ(d)cos(ψυ(d)), 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 Kdυ, the core-shell model of equations (4)-(7) predicts that ω¯υ<ωDD<ω¯d. Equivalently, the period τDD=2π/ωDD must be in the range

Using the parameters in Table 1, equation (15) tells us that τDD=24.8h, in good agreement with experimental observations in rats by Campuzano et al. (1998), who reported τDD=24.5h.

Under DD conditions, the analytical result ofpredicts thatis bounded betweenand. This is valid for positive inter- and intracouplings. However, in the case of negative couplings,is no longer bound betweenand. In the particular case when(and),becomes equation (18) equation (15) τ D D τ υ τ d τ D D τ υ τ d K d υ < 0 K υ d ≥ 0

from which it follows that τDD>τυ for ω¯υ<ω¯d (τυ>τd), under the additional assumption that |Kdυ|<Kυd as suggested by experimental observations (Varadarajan et al., 2018) for wild-type (WT) mice. Equation (19) also shows that strong asymmetry in intercoupling strength can cause the rhythm in the shell to become increasingly similar to the core, since τDD∼∼τυ when |Kdυ|≪Kυd. Conversely, if Kdυ is instead positive and Kυd negative, it can be shown that the rhythm in the core becomes increasingly similar to the shell for |Kdυ|≪Kυd. This follows from re-writing equation (19) for Kυd<0 (and Kdυ≥0) based on the same experimental observations (τυ>τd and Kdυ<|Kυd|), for which τDD<τd. The results of equations (18) and (19) show that experimental measurements of τDD, τυ, and τd can indicate whether Kdυ is positive or negative.

The relative strength of intercouplingsandcan also be estimated when the intracouplingsandare large by comparison. Under these conditions,, so thattakes the form K υ d K d υ K υ υ K d d ρ υ ρ d ≈ ≈ 1 equation (15)

This equation can be rearranged to yield the relative intercoupling strength

which can therefore be determined from experimental values of,, and. τ D D τ υ τ d

A lower bound on the total intercouplingcan be established with reference to the phase differenceunder DD conditions. Fromunder LD conditions, the replacementyields | + | K υ d K d υ ψ d ψ υ − ω F ω D D → equation (9)

First, this shows that ψd leads ψυ under DD conditions, as experimental observations indicate that τυ>τd (ω¯υ<ω¯d) and Kdυ<Kυd (Taylor et al., 2017; Varadarajan et al., 2018), assuming positive intercoupling. In this case, peak activity in the shell precedes peak activity in the core, similar to LD conditions (see equation (12)). However, if either Kυd or Kdυ becomes negative, their relative strength can cause the core to lead the shell under DD conditions, that is, if Kυd≫Kdυ or Kdυ≫Kυd. For the same to take place under LD conditions, Kυd must simply be negative. Second, when the intracouplings are sufficiently larger than the intercouplings, ρυ≈ρd≈1, so that a∼∼2Kυd and b∼∼2Kdυ; equation (22) then becomes

Thus, the total intercoupling must be larger than or equal to the difference between the mean natural frequencies of the shell and core, in absolute value:

since. | ( − ) | ≤ sin ψ d ψ υ 1

The results of this section show that the sign, relative strength, and the minimum total intercoupling can be determined from the free-running periods of the core (), shell (), and the SCN as a whole (), which can be measured experimentally under DD conditions. Moreover, the results show that, if all other parameters remain the same, changes in the relative strength and sign of the intercouplings can cause the core to lead the shell even when, and cause the synchronized rhythm of the core to approximate the endogenous (free-running) rhythm of the shell or vice versa. Below the minimum total intercoupling, the SCN is no longer able to synchronize. τ υ τ d τ D D τ υ τ d > | + | = | − | K υ d K d υ ω ¯ d ω ¯ υ

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 inwith a constant: equation (2) B

In our reduced Kuramoto model, parameteris phenomenological. It can be positive or negative. We assume that the magnitudecharacterizes light intensity and chooses the sign ofto satisfy Aschoff’s first rule. From this, it follows thatis negative for nocturnal animals and positive for diurnal animals, as will be proven below. This choice of parameterallows us to describe the basic properties of the SCN under LL conditions, in nocturnal and diurnal animals. B | | B B B B

As one can see in, introducing parameterrenormalizes the endogenous free-running frequencies of core oscillators equation (2) B

from which it follows that the renormalized mean frequency of the core is

and, equivalently, the renormalized mean period of oscillations is

whereis the endogenous mean free-running period of isolated core oscillators. This shows that under LL conditions, the positive parameter(diurnal animals) forces core oscillators to run faster than under DD conditions (). Conversely, if parameteris negative,, (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 (), 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. τ υ B > 0 ω ¯ υ * > ω ¯ υ B − < < 2 ω ¯ υ B ω ¯ υ * < ω ¯ υ τ d * = τ d

To determinein 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-() as DD conditions (), by renormalizing the mean period as in. This means that we can use the equations derived for DD conditions by making the simple replacementsand. Fromin particular, it then follows that the frequencyof the steady circadian rhythm is τ L L F = 0 τ υ → τ υ * τ D D τ L L → ω L L equations (4) 7 equation (28) equation (20)

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, 1979; Pittendrigh and Daan, 1976a; Carpenter and Grossberg, 1984). Within our model, by comparison with equation (20), it follows that ωLL−ωDD=aBω¯υ/(a+b), which satisfies Aschoff’s first rule when a/(a+b)>0 (refer to the discussion on the sign and strength of the intercouplings in the section on DD conditions). Assuming that B is positive for diurnal animals yields τLL<τDD (ωLL>ωDD). Conversely, negative B yields τLL>τDD (ωLL<ωDD) 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 TLL(sr) at a critical brightness |Bc|, with an order of magnitude difference between nocturnal and diurnal cases. In the nocturnal case, |Bc|≃0.24, for which τLL=25.2h, whereas in the diurnal case, |Bc|≃3.23, for which τLL=21.7h. 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.15h 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, TLL(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 TLL(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 replacementin equation (). For diurnal animals (), increasing the light intensity decreases the core-shell phase difference in animals whereand both intercouplings are positive. In nocturnal animals, increasing the light intensity increases the core-shell phase difference. In both cases, sufficiently large light intensity,, 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. ω υ → ω υ * 22 B > 0 τ d τ υ < | | B