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

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 m m within groups and strength K n m between groups. In a periodic external field with period T and angular frequency ω F = 2 π / T , the local field phase ω F t + ϕ ( 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

d θ i ( m ) d t = ω i ( m ) + F ( m ) sin ( ω F t + ϕ ( m ) − θ i ( m ) ) .

(1)

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

d θ i ( m ) d t = ω i ( m ) + 1 N n ∑ n ∑ j ∈ G n K n m sin ( θ j ( n ) − θ i ( m ) ) + F ( m ) sin ( ω F t + ϕ ( m ) − θ i ( m ) ) .

(2)

The macroscopic state of each group is characterized by the synchronization index ρ m and group phase ψ m of the complex order parameter

z m = ρ m e i ψ m ≡ 1 N m ∑ j ∈ G m e i θ j ( m ) .

(3)

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 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

d ρ υ d t = − ρ υ Δ υ + 1 2 ρ υ K υ υ ( 1 − ρ υ 2 ) + 1 2 F ( 1 − ρ υ 2 ) cos ψ υ + 1 2 ρ d K d υ ( 1 − ρ υ 2 ) cos ( ψ d − ψ υ ) ,

(4)

d ψ υ d t = − Ω υ − 1 2 F ( 1 + ρ υ 2 ) ρ υ sin ψ υ + 1 2 ρ d K d υ ( 1 + ρ v 2 ) ρ v sin ( ψ d − ψ υ ) ,

(5)

d ρ d d t = − ρ d Δ d + 1 2 ρ d K d d ( 1 − ρ d 2 ) + 1 2 ρ υ K υ d ( 1 − ρ d 2 ) cos ( ψ υ − ψ d ) ,

(6)

d ψ d d t = − Ω d + 1 2 ρ υ K υ d ( 1 + ρ d 2 ) ρ d sin ( ψ υ − ψ d ) ,

(7)

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 υ d = K d υ = 0 , these equations are reduced to the explicit equations for the forced Kuramoto model derived in Childs and Strogatz (2008).

OPEN IN VIEWER

Figure 1.

Schematic representation of the core-shell model consisting of 2 groups of oscillators, the core (ventral part υ ) and shell (dorsal part d ). Parameters F and w F are the strength and frequency of an external light cue acting exclusively on the core. The core (shell) is characterized by a mean free-running frequency w ¯ u ( d ) , the statistical spread of free-running frequencies ∆ u ( d ) (Lorentzian half width at half maximum), and the intracoupling K u u ( d d ) between oscillators within the core (shell). Coupling between core and shell oscillators is determined by intercouplings K u d (core on shell) and K d u (shell on core).

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 n m / Δ υ → K n m and time t Δ υ → t . Dimensionless parameters are summarized in Table 1, along with their dimensional equivalents, and reference values are found in the literature.

OPEN IN VIEWER

Table 1.

List of parameters (parameters) and a summary of the corresponding values reported by Taylor et al. (2017) (reference), the core-shell model (model) of equations (4)-(7).

Parameters	Reference	Model
τ υ	25 . 1 h	25 . 1 h
τ d	23 . 9 h	23 . 3 h
σ υ	1 . 3 h	1 . 3 h
σ d	1 . 9 h	1 . 9 h
K υ υ	—	5.6
K d d	—	4.0
K υ d	—	1.1
K d υ	—	0.5
F	—	1.5
ω ¯ υ	—	19.3
ω ¯ d	—	20.8
∆ υ	—	1.0
∆ d	—	1.7

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 υ υ ( d d ) , 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: ω ¯ υ ( d ) = 2 π / τ υ ( d ) , 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).

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 = 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 . 1 h 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 ω ¯ υ ( d ) = 2 π / τ υ ( d ) . From the definition of standard deviation, and given τ υ ( d ) ≫ σ υ ( d ) , the standard deviation of the free-running frequencies is approximated as S D υ ( 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 S D υ ( d ) provide a measure of the statistical spread of free-running frequencies about the mean value, so we simply take ∆ υ ( d ) = S D υ ( 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 υ υ , K d d , 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 ( K υ υ , K d d , K υ d , K d υ ≫ 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 d d . 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 = K d υ = 0 and F = 0 ), intracoupling increases with the fraction of synchronized oscillators ρ υ ( d ) as follows:

K υ υ ( d d ) = 2 Δ υ ( d ) 1 − ρ υ ( d ) 2 or ρ υ ( d ) = 1 − 2 Δ υ ( d ) K υ υ ( d d ) ,

(8)

which can be seen from equations (4) and (6). For positive intracoupling, it follows that ρ υ > ρ d whenever K υ υ / K d d > ∆_v/∆_d, which is the case since K υ υ > K d d 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 K d d = 4 . 0 . The critical value of the core (shell) intracoupling K c υ υ = 2 ∆_v ( K c d d = 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 K c d d = 3 . 4 in the shell, where ( K d d − K c d d ) / K c d d = 0 . 18 , than to the critical value K c υ υ = 2 . 0 in the core, where ( K υ υ − K c υ υ ) / K c υ υ = 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 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 ( K υ d > K d υ ). 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 + K d υ > ω ¯ d − ω ¯ υ = 0 . 97 . Starting from these bounds and intracouplings K υ υ = 5 . 6 and K d d = 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 υ ( p ) − t d ( p ) . 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 υ ( p ) − t 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). Parameters τ d , K υ d , K d υ , and F were adjusted by inspection, until t υ ( p ) − t d ( p ) = 2 . 3 h in the entrained state ( T = 24 h ), for τ d = 23 . 3 h , K υ d = 1 . 1 , K d υ = 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 K d d 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 (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):

sin ( ψ d − ψ υ ) = 2 ρ d ( ω ¯ d − ω F ) ρ υ ( 1 + ρ d 2 ) K υ d .

(9)

When the synchronization index of the entrained SCN is large ( ρ υ ≈ ρ d ≈ 1 ), or more generally when 2 ρ d ≈ ρ υ ( 1 + ρ d 2 ) , it follows that

sin ( ψ d − ψ υ ) ≈ ω ¯ d − ω F K υ d .

(10)

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 :

K υ d > ω ¯ d − ω F ,

(11)

since | sin ( ψ d − ψ υ ) | ≤ 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 = ω F − ω ¯ d and of the intercoupling K υ 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 Ω d = 2 π / T − 2 π / τ d < 0 . In this case, the phase difference Δ ψ = ψ d − ψ υ is positive since K υ d is also positive ( K υ d = 1 . 1 ). 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 ω F = 2 π / T . Symbolically, Re( z υ ( d ) ) = ρ υ ( d ) cos ( ω F t + ψ υ ( d ) ) . 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( z υ ( d ) ) is maximum in the core (shell) at time t υ ( d ) ( p ) = 2 π n / ω F − ψ υ ( d ) / ω F , where n is an integer number. From here, it follows that the difference in peak activity time between the shell and the core is

t υ ( p ) − t d ( p ) = T ( ψ d − ψ υ ) 2 π .

(12)

Thus, the peak time difference t υ ( p ) − t d ( p ) is determined by the phase difference ψ d − ψ υ . 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 :

t d ( p ) = t υ ( p ) − T 2 π arcsin ( ω ¯ d − ω F K υ d ) .

(13)

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 ) > t 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).

OPEN IN VIEWER

Figure 2.

Results of numerical simulations of equations (4)-(7) with the parameters in Table 1. (a) Times t υ ( p ) (red circles) and t d ( p ) (navy squares) at which the real part of complex order parameters z υ and z d peak in the core and shell, respectively, and the core-shell difference between peak times t υ ( p ) − t d ( p ) (green diamonds) in the entrained state, as a function of the external light-dark period T . (b) Time dependence of the real part of the order parameter in the core ( Re( z υ ) , red line) and shell ( Re( z d ) , navy line) for light-dark cycles with period T = 24 h . The peak time in the shell precedes the peak time in the core by 2 . 3 h .

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 ) − t d ( p ) = 2 . 3 h corresponds to the phase difference ψ d − ψ υ = 0 . 607 rad ≈ 35 degrees . 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.

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: L L E < T < U L E , where L L E and U L E 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 L L E ∼ ∼ 22 h up to U L E ∼ ∼ 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( z υ ( d ) ) = ρ υ ( d ) cos ( ω F t + ψ υ ( 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.

OPEN IN VIEWER

Figure 3.

Results of numerical simulations of equations (4)-(7) with the parameters in Table 1. (a) Period T s r 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 s r . Black crosses are experimental data reported by Campuzano et al. (1998). The dashed horizontal line is the period of SCN oscillations τ D D = 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 ( s r , solid makers) components of SCN activity, in the core ( υ , circles) and shell ( d , squares). Vertical lines indicate the lower and upper bounds of the entrainment range, L L E = 23 . 26 h and U L E = 25 . 28 h , respectively. For comparison, the mean free-running periods of the shell and the core are τ d = 23 . 3 h and τ υ = 25 . 1 h , respectively. (c) Dependence of the entrainment range (gray horizontal bars) on the difference between the free-running periods of core and shell, τ υ − τ d . Here, τ υ is fixed as τ υ = 25 . 1 h (black dots) and τ d (blue dots) varies. Abbreviation: LLE = lower limit of entrainment;SCN = suprachiasmatic nucleus; ULE = upper limit of entrainment.

First, Figure 3a shows that a second rhythm ( s r ) with period T s r emerges outside the entrainment range. Above the upper bound ( U L E = 25 . 28 h ), T s r is smaller than T and decreases with increasing T . Below the lower bound ( L L E = 23 . 26 h ), T s r 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 s r = 25 h at T = U L E = 28 . 5 h (see Figure 3 by Usui et al., 2000), and T s r = 23 . 5 h at T = L L E = 23 h (see Figure 4 by Usui et al., 2000). The authors also report that T s r increases with decreasing T below L L E and decreases with increasing T above U L E , 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 % .

OPEN IN VIEWER

Figure 4.

Period τ of independent activity components in the core as a function of parameter B under constant light conditions (from numerical simulations of equations (4)-(7) with the parameters in Table 1). The period τ L L of the free-running component is indicated by circles, and the period T L L ( s r ) of a second (dissociated) component is indicated by navy stars. Negative B corresponds to the nocturnal mouse (see Table 1 for parameters) and positive B to a fictitious diurnal counterpart. The second rhythm appears at B c ≃ − 0 . 24 in the nocturnal case and B c ≃ 3 . 23 in the diurnal case.

Second, as shown in Figure 3b, the entrainment range is characterized by a single frequency component of intensity I F ( υ ) = I F ( d ) in the core and shell. This component corresponds to an entrained rhythm with period T between L L E = 23 . 26 h (lower bound) and U L E = 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 U L E and below L L E , as shown in Figure 3b. Below L L E , core and shell contributions to the second rhythm at period T s r increase, while contributions to the rhythm at period T decrease. However, when the LD period T F is lengthened above U L E , the core’s contribution is almost exclusively to the component at the LD period (near constant I F ( υ ) and vanishingly small I s r ( υ ) ), 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 > U L E (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 L L E , 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 L L E and U L E . The bifurcation at L L E is characterized by the emergence of a second rhythm with period T s r ∼ ∼ T , as shown in Figure 3a. The bifurcation at U L E , however, is characterized by a large drop in T s r relative to T = U L E . 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 L L E − T s r = 0 . 5 h , whereas T s r − U L E = 3 . 5 h for rats. For comparison, our model predicts L L E − T s r = 0 . 2 h and T s r − U L E = 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 U L E = 25 . 28 h is caused by the super Hopf bifurcation, while the saddle-node bifurcation occurs at L L E = 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 free-running period of the core τ υ = 25 . 1 h and shell τ d = 23 . 3 h for the model parameters in Table 1:

τ d ∼ ∼ L L E < T < U L E ∼ ∼ τ υ .

(14)

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 ) − t 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 τ D D across the whole SCN. Since τ D D , τ υ , 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 ω D D = 2 π / τ D D 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 ω F → ω D D , from which it follows that

ω D D = a ( a + b ) ω ¯ υ + b ( a + b ) ω ¯ d ,

(15)

where

a = ρ υ 2 K υ d ( 1 + ρ d 2 ) ,

(16)

b = ρ d 2 K d υ ( 1 + ρ υ 2 ) .

(17)

The synchronized SCN is characterized by core and shell synchronization indices ρ υ and ρ d . To find τ D D , 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 K d υ , the core-shell model of equations (4)-(7) predicts that ω ¯ υ < ω D D < ω ¯ d . Equivalently, the period τ D D = 2 π / ω D D must be in the range

τ d < τ D D < τ υ .

(18)

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

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

ω D D = ω ¯ υ − | b | ( ω ¯ d − ω ¯ υ ) ( a − | b | )

(19)

from which it follows that τ D D > τ υ for ω ¯ υ < ω ¯ d ( τ υ > τ d ), under the additional assumption that | K d υ | < 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 τ D D ∼ ∼ τ υ when | K d υ | ≪ K υ d . Conversely, if K d υ is instead positive and K υ d negative, it can be shown that the rhythm in the core becomes increasingly similar to the shell for | K d υ | ≪ K υ d . This follows from re-writing equation (19) for K υ d < 0 (and K d υ ≥ 0 ) based on the same experimental observations ( τ υ > τ d and K d υ < | K υ d | ), for which τ D D < τ d . The results of equations (18) and (19) show that experimental measurements of τ D D , τ υ , and τ d can indicate whether K d υ is positive or negative.

The relative strength of intercouplings K υ d and K d υ can also be estimated when the intracouplings K υ υ and K d d are large by comparison. Under these conditions, ρ υ ≈ ρ d ≈ 1 , so that equation (15) takes the form

ω D D = K υ d ( K υ d + K d υ ) ω ¯ υ + K d υ ( K υ d + K d υ ) ω ¯ d .

(20)

This equation can be rearranged to yield the relative intercoupling strength

K υ d K d υ = τ D D / τ d − 1 1 − τ D D / τ υ ,

(21)

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

A lower bound on the total intercoupling | K υ d + K d υ | can be established with reference to the phase difference ψ d − ψ υ under DD conditions. From equation (9) under LD conditions, the replacement ω F → ω D D yields

sin ( ψ d − ψ υ ) = 2 ρ d ( ω ¯ d − ω D D ) ρ υ K υ d ( 1 + ρ d 2 ) = 2 ρ υ ρ d ( ω ¯ d − ω ¯ υ ) ( a + b ) .

(22)

First, this shows that ψ d leads ψ υ under DD conditions, as experimental observations indicate that τ υ > τ d ( ω ¯ υ < ω ¯ d ) and K d υ < 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 K d υ becomes negative, their relative strength can cause the core to lead the shell under DD conditions, that is, if K υ d ≫ K d υ or K d υ ≫ 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 ∼ ∼ 2 K υ d and b ∼ ∼ 2 K d υ ; equation (22) then becomes

sin ( ψ d − ψ υ ) ≈ ω ¯ d − ω ¯ υ K υ d + K d υ .

(23)

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:

| K υ d + K d υ | ≥ | ω ¯ d − ω ¯ υ | ,

(24)

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 ( τ d ), and the SCN as a whole ( τ D D ), 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 τ υ > τ d , 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 | K υ d + K d υ | = | ω ¯ d − ω ¯ υ | , the SCN is no longer able to synchronize.

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 :

F sin ( ω F t + ϕ − θ i ( υ ) ) → B .

(25)

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

ω i ( υ ) ⇒ ω i ( υ ) + B ,

(26)

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

ω ¯ υ * = ω ¯ υ + B ,

(27)

and, equivalently, the renormalized mean period of oscillations is

τ υ * = 2 π ω ¯ υ * = τ υ 1 + B τ υ 2 π ,

(28)

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 ω ¯ υ < B < 0 , (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 τ L L 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 τ D D → τ L L . From equation (20) in particular, it then follows that the frequency ω L L of the steady circadian rhythm is

ω L L = a ( a + b ) ω ¯ υ * + b ( a + b ) ω ¯ d = a ( a + b ) ( ω ¯ υ + B ) + b ( a + b ) ω ¯ d .

(29)

Based on experimental observations of different animal species, Aschoff’s first rule predicts that τ L L < τ D D in diurnal animals and τ L L > τ D D in nocturnal animals, where τ D D and τ L L 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 ω L L − ω D D = a B ω ¯ υ / ( 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 τ L L < τ D D ( ω L L > ω D D ). Conversely, negative B yields τ L L > τ D D ( ω L L < ω D D ) 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 τ L L (increases ω L L ), whereas for nocturnal animals ( B < 0 ), increasing | B | increases τ L L (decreases ω L L ), 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., τ L L 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 L L ( s r ) at a critical brightness | B c | , with an order of magnitude difference between nocturnal and diurnal cases. In the nocturnal case, | B c | ≃ 0 . 24 , for which τ L L = 25 . 2 h , whereas in the diurnal case, | B c | ≃ 3 . 23 , for which τ L L = 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 τ L L = 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 L L ( s r ) , 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 τ L L , the dissociated rhythm T L L ( s r ) 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

sin ( ψ d − ψ υ ) = 2 ρ d ( ω ¯ d − ω L L ) ρ υ K υ d ( 1 + ρ d 2 ) = 2 ρ υ ρ υ ( ω ¯ d − ω ¯ υ − B ) ( a + b ) ,

(30)

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.