Photic Desynchronization of Two Subgroups of Circadian Oscillators in a Network Model of the Suprachiasmatic Nucleus with Dispersed Coupling Strengths

The suprachiasmatic nucleus (SCN) is the master circadian clock in mammals and is composed of thousands of neuronal oscillators expressing different intrinsic periods. These oscillators form a coupled network with a free-running period around 24 h in constant darkness and entrainable to the external light-dark cycle (T cycle). Coupling plays an important role in setting the period of the network and its range of entrainment. Experiments in rats have shown that two subgroups of oscillators within the SCN, a ventrolateral (VL) subgroup that receives photic input and a dorsomedial (DM) subgroup that is coupled to VL, can be desynchronized under a short (22-h) T cycle, with VL entrained to the cycle and DM free-running. We use a modified Goodwin model to understand how entrainment of the subgroups to short (22-h) and long (26-h) T cycles is influenced by light intensity, the proportion of neurons that receives photic input, and coupling heterogeneity. We find that the model’s critical value for the proportion of photically-sensitive neurons is in accord with actual experimental estimates, while the model’s inclusion of dispersed coupling can account for the experimental observation that VL and DM desynchronize more readily under the 22-h than under the 26-h T cycle. Heterogeneous intercellular coupling within the SCN is likely central to the generation of complex behavioral patterns.


Introduction
Circadian (,24 h) rhythms in physiological and behavioral measures are universal in living things, reflecting the period of the earth's rotation. In mammals, circadian rhythms are regulated by a master clock in the suprachiasmatic nucleus (SCN) of the hypothalamus, composed of approximately 20,000 neuronal oscillators; SCN neurons are nonidentical, express different intrinsic periods, and are coupled together to form a network with a coherent output [1]. The period of the network's output signal is adaptable. Under constant darkness, the rhythm has a free-running period close to 24 h; whereas under an external lightdark cycle (T cycle), it is precisely entrained to a period identical to the external cycle.
The SCN network is heterogeneous [2,3,4]. It can be divided into distinct functional subgroups, including a ventrolateral part (VL), which receives photic input from the retina, and a dorsomedial part (DM), which is coupled to VL; both VL and DM contribute to the generation of overt circadian rhythms in physiological and behavioral measures. Peptide neurotransmitters differ between the VL and DM subdivisions, with neurons that express vasoactive intestinal polypeptide (VIP) in the VL and arginine vasopressin in the DM. Periods may vary in different regions of the SCN, with DM running faster than VL in tissue slices [5]. Gamma aminobutyric acid (GABA) neurons are present throughout the SCN and may play a role in coupling the two subdivisions [6]. It has been shown that the circadian oscillation between VL and DM can desynchronize with exposure to short T cycles [7] or after a phase shift of the light-dark cycle [6,8,9]; the VL appears to set the final phase of the SCN after the phase shift [6,8,9].
Much experimental [10,11,12] and theoretical [13,14,15] work has been motivated by a desire to understand how this heterogeneous SCN network is reliably entrained and able to generate a coherent output signal, and neuropeptidergic mechanisms appear to be necessary elements [16,17,18]. Modeling studies suggest that the circadian clock's free-running period is proportional to the average intercellular coupling strength [13] and that coupling governs the clock's range of entrainment to T cycles [15]. However, coupling strength between cells in the SCN network is unlikely to be uniform. The effects of heterogeneous coupling on network synchronization have been studied previously in multi-oscillator models. Daido considered the dispersion of coupling strengths in the Kuramoto model and studied the synchronization property of the network [19,20]; coupling strength between two oscillators was chosen from a normal distribution. Hong and Strogatz considered a heterogeneous network with excitatory (positive) and inhibitory (negative) coupling in the Kuramoto model to understand the relative contributions of excitatory and inhibitory properties on network synchronization [21]. Our recent work (C.G. and Z.L.) has demonstrated that the dispersion of coupling strengths between SCN cellular oscillators can influence the emergent free running period of the network [22]. To our knowledge, however, there has been no work on the relationship between coupling dispersion and network entrainment.
We examine this issue in the present work, inspired by an interesting experiment performed by de la Iglesia et al. [7] in which rats were exposed to an artificially short 22-h T cycle (11 h light alternating with 11 h darkness). Individual animals expressed two separate circadian motor activity rhythms, with one rhythm entrained by the light and oscillating with a period equal to the external cycle, while the other was not entrained and expressed a period greater than 24 h. Analyses of SCN gene expression suggested that these two motor activity rhythms reflected the stable forced desynchronization of VL and DM subdivisions, respectively. Here we model how entrainment of the subdivisions is influenced by coupling dispersion, as well as by the proportion of cellular oscillators that receive photic input (i.e., that are within VL) and the light intensity.
We use the Goodwin model, a network model of coupled oscillators that has been widely used to describe the mammalian circadian clock [13,22,23,24,25,26] (defined in Methods). An individual cellular oscillator of the Goodwin model has three variables: a clock gene mRNA, a clock protein, and a transcriptional inhibitor, all of which form a transcription-translation negative feedback loop. It is assumed that light induces the clock gene mRNA, that a neurotransmitter is increased by the clock gene mRNA, and that neurotransmitters from different neurons form a mean field that couples the neurons together. We consider that pN neurons receive photic input, where N is the total number of neurons in the SCN network and p is the ratio of the number of VL neurons to the total number of SCN neurons. We take T cycles of 22 h and 26 h as examples, i.e., symmetrically distant from 24 h. We chose mean field coupling for all of the neurons in the Goodwin model. The coupling strength g i of all the N neurons satisfies a normal distribution with mean valueSgT and deviation g.

T-cycle Entrainment of an SCN Network without Dispersion of Coupling Strengths
To determine the effects of p and light intensity, K f , on the entrainment of VL and DM to T cycles, we have numerically simulated the Goodwin model with no dispersion of coupling strengths, i.e., g~0. Figure 1 shows the mean field time series of VL and DM oscillations in the 22-h light-dark cycle. Similar to previous observations [13], we find that the time series show quasiperiodic behavior with low light intensity. In (A), the behavior of VL follows the 22-h cycle and sustains a stable phase relationship to it, while the behavior of DM loses its phase relationship to the cycle and runs with a period close to the intrinsic period of the network. This dissociation mimics the forced desynchronization of motor activity rhythms in rats under such a T cycle, as noted previously [7,27]. When p is increased, both VL and DM can be entrained, as in (B). Here the peak of the mean field time series of VL appears around the onset of darkness, whereas that of the DM is phase delayed. This change also can be implemented by K f . If K f is reduced, neither VL nor DM entrain to the 22-h cycle. If K f is increased, both VL and DM can be entrained, as in (B). In sum, both the number of neurons receiving light and the light intensity are important factors for entrainment of the entire SCN network to the T-cycle.
To understand the influence of the parameters p and K f on entrainment, we have calculated the phase diagram of the period of the mean fields of VL and DM in the p-K f plane under short and long T cycles, i.e., of 22 h and 26 h, withg~0 ( Figure 2). (A) and (C) show that the behavior of VL follows the T cycle for all values of p, provided that the light intensity is greater than a critical value, such as K f w0:02. When K f v0:02, the period of VL may not be entrained to the T cycle, depending on p; for example, in Fig. 2A, the period of VL can be 23 h or 24 h. For DM to follow the cycle, however, p also must be larger than some threshold; that is, there must be a sufficient number of lightreceiving neurons in VL in order to drive the neurons in DM. For a given level of K f , increasing p may allow the entire network to entrain to the driving T cycle. Surprisingly, when p is decreased under the 26-h T cycle (D), there is a threshold of p at which the period of DM suddenly jumps to a value of 20.8 h to form a locking ratio of 4:5 with the 26-h T cycle; with further decreased p, the period monotonically increases to reach a value of 24 h at p~0. A comparison of (B) and (D) shows that the threshold of p for entrainment of DM to the 26 h T cycle is greater than that for the 22 h cycle, suggesting that desynchronization between VL and DM might be more likely under long than under short T cycles. Experimentally, however, this appears not to be the case [28,29], prompting us to consider the influence of heterogeneous coupling strengths on the behavior of the SCN network. Figure 3 shows the phase diagram of the period of the mean fields of VL and DM in the p-K f plane using two values for g. Although qualitatively similar to the diagrams in Figure 2, there are quantitative differences when network coupling strengths are dispersed. In the case of the 22-h T cycle, entrainment is only modestly affected; in contrast, in the case of the 26-h T cycle, increased coupling dispersion significantly reduces the critical value of p for DM entrainment, suggesting that the network can be entrained to the long T cycle with a relatively lower K f .

T-cycle Entrainment of an SCN Network with Dispersion of Coupling Strengths
For weak K f and different values of g, we find that the critical p (p c ) for the 22-h and 26-h T cycles reaches approximately the values of 0.28 and 0.20, respectively. Figure 4 represents the variation in p c for different values of K f and g. In the case of the 22-h T cycle (A), there is little variation in p c for different K f , e.g., p c is between 0.14 and 0.24 for a K f~0 :04. However, p c does change significantly in the case of the 26-h T-cycle (B). Thus, dispersion of coupling strengths affects entrainment in an asymmetric way, with an influence that is larger for the long than for the short T cycle.
Instead of randomly assigning coupling values to the network, we also studied the network by selectively assigning coupling values. In two separate trials, we assigned the strongest coupling values to either VL or DM. The phase diagram of the period of the mean fields of VL and DM in the p-K f plane was similar to that previously reported.
We also simulated the network with dispersed oscillator periods, rather than dispersed coupling strengths, by selecting different values for the standard deviation of period (s) for each individual oscillator, such that the period distribution has a mean of 24 h with variability. Without dispersed coupling, we do not observe DM entrainment to the long T cycle at any p until s is increased to a value greater than 5 h. Since such a large non-identical intrinsic period is not realistic, the dispersion of coupling strengths is likely a crucial factor affecting the entrainment of the network to different T cycles.
Importantly, dispersion of coupling influences the mean field amplitude ( Figure 5). The amplitude of VL as a function of p changes modestly as the dispersion g is increased in the 22-h as well as the 26-h T cycle (A and C). On the other hand, for DM in the 22-h T cycle (B), amplitude decreases as g increases. In the 26- h T cycle (D), DM amplitude as a function of p changes dramatically, with relatively diminished amplitude as p is increased; dispersion g counteracts this effect. The enhancement of DM amplitude by increased g in the 26-h cycle could be due to enhanced phase synchronization of the oscillators in the network, increased amplitude of the individual oscillators, or both. To begin to distinguish among these possibilities, we studied the effect of dispersed coupling on the order parameter, a measure that represents phase synchronization of the network.

Effect of Coupling Dispersion on the Order Parameter of the Network
Order parameter characterizes the synchronization property of a network [30,31], and it is defined here by estimating the phases of the oscillators in VL and DM (see Methods). The order parameter will be unity if all oscillators in the network are perfectly synchronized and zero if they are completely uncorrelated. When their behavior is between these two extremes, the order parameter will be in (0, 1), i.e., representing a phase difference between VL and DM or desynchronization of individual oscillators within VL and/or DM.
We have studied the influence of g on the order parameter. Figure 6 shows the dependence of order parameter R on the parameters p and K f in the p{K f plane. To reveal the effect of coupling dispersion, we have considered two cases, one with g~0:15 and the other with g~0:0. Under the 22-h T cycle, coupling dispersion reduces R for larger p and K f values; whereas under the 26-h T cycle, coupling dispersion enhancesR. As p and K f increase from (0,0), the relationship between VL and DM changes; comparison of Figure 6 with Figures 2 and 3 visualizes the regions where Rv1, i.e., either when VL and DM express different periods or when VL and DM express the same period but with a large phase difference between them. Thus, for the 22-h cycle, although higher p and K f values enhance both VL and DM entrainment to the cycle, the reduction of R with dispersed coupling suggests that individual oscillators are not fully synchronized within the network, with a greater vulnerability to perturbations of light intensity. For the 26-h cycle, coupling dispersion synchronizes network oscillation for p.0.6; the gradual  increase of DM mean field amplitude as p increases further ( Figure 5D) is thus attributable to increased individual oscillator amplitude.
These considerations imply that VL and DM desynchronize more readily under the 22-h than under the 26-h T cycle and that dispersion of coupling strengths improves network robustness preferentially under the 26-h cycle.

Effect of Coupling Dispersion on the Network's Phase Response Curve to a Light Pulse
The network's capacity to generate phase advances or delays can be quantified as a phase-response curve (PRC), measured by plotting the phase shifts that occur in the rhythm when discrete light pulses are applied at different phase points across the circadian cycle [24,32,33,34]. Figure 7 represents the family of PRC's obtained to a 1-h light pulse of increasing intensities, showing that the phase response region (i.e., the area under the delay and advance zones) increases in magnitude with increasing K f . Notably, as the value of g increases, the area under the delay zone increases relatively more than that under the advance zone, as calculated in Table 1, where S represents the ratio of the area under the delay zone to the area under the advance zone.
Advances should correspond to the capacity of the network to follow a T cycle less than 24 h, while delays should correspond to its capacity to follow a T cycle greater than 24 h [35].

Discussion
Here we analyze the photic desynchronization of two subgroups of circadian oscillators in a network model of the suprachiasmatic nucleus. As also demonstrated in experiments with rats exposed to a short T cycle of low light intensity [7,36], a subgroup of oscillators receiving photic input (VL) can entrain to the external cycle while the other, coupled subgroup (DM) expresses an unentrained period greater than 24 h.
Granada et al. [37] have modeled this forced desynchronization of rat activity rhythms as a single oscillator with oscillatory interactions (modulation and superposition) between the external cycle and the internal clock, while Schwartz et al. [38] have modeled entrainment to the T cycle by two coupled oscillators forced by a Zeitgeber. Casiraghi et al. [39] have used a two oscillator model to analyze a chronic jet lag paradigm that leads to forced desynchrony, and they observed an asymmetry in its behavior similar to our findings reported here. We have taken the Goodwin model and extended it to include p, the proportion of all SCN cellular oscillators that receive photic input, and g, the dispersion of coupling strengths. We find, first, that network desynchronization (with an entrained VL but an unentrained DM) depends on light intensity and the value of p. Relatively higher light intensities protect the network from desynchronization, as reported experimentally [36]. Experiments estimate that the value of p for the rodent SCN ranges from 20% to 33%, based on molecular, electrophysiological, and computational studies [40,41]. Comparing these results to our simulations in Figure 4, we find that g~0:15 is a good parameter value to  fit the experiments. At this value, there is a critical value of p for network entrainment to short (22 h) and long (26 h) T cycles of 0.28 and 0.20, respectively, and the critical value appears fairly insensitive to K f . We predict that the rat SCN is likely to have this very large heterogeneity in coupling strengths, given that the critical p best matches the experimental estimate for the larger values of g.
Second, we find that the inclusion of dispersed coupling strengths affects network entrainment in a preferential manner, such that increased g significantly reduces the critical value of p for DM entrainment to the 26-h T cycle. A consequence of this g influence is that network robustness is superior under the 26-h cycle while desynchronization is favored under the 22-h cycle. In fact, such an asymmetry has been observed experimentally, with no obvious desynchronization of rat motor activity rhythms during exposure to long T cycles [28].
The basis for this asymmetry is unclear. Coupling dispersion appears to generally increase the effect of light on the system, and since the delay zone of the PRC is greater than the advance zone, a preferential effect on entrainment to the 26-h T cycle might be expected. But such an explanation does not account for the ushaped, rather than monotonic, p c function for the 22-h T cycle. Moreover, differences in entrainment to short and long T cycles may not be a general feature of the SCN; it may not be true for other species [42], and the possible effects of locomotor activity itself on SCN network behavior (e.g., in or out of a running wheel, diurnal or nocturnal activity pattern) needs further investigation.
Our findings may provide new elements to the theory of coupled oscillators, especially with regard to chimera states in which one group of the system is synchronized and the other is desynchronized [43,44,45,46,47]; in these studies, the oscillators in the two groups are identical and the chimera states are generally induced by the initial conditions. However, in our case, the discovery that the subgroup VL may be entrained to the T cycle while the group DM remains free-running is similar to the chimera state, but this phenomenon does not depend on the initial conditions. Thus, a novel oscillator theory is needed to explain this robustness to initial conditions and should be a topic for further studies.
Heterogeneous intercellular coupling within the SCN is likely central to the generation of complex behavioral patterns. Nonuniform SCN network architecture also has been implicated in the phase-''splitting'' of locomotor activity cycles seen in hamsters maintained in constant environmental light [48]. In the future, we hope to consider how topology influences entrainment, in contrast to the mean field used here.

Methods
We represent each mammalian cell of the network as a Goodwin oscillator. The Goodwin model is a widely used mathematical model to represent the behavior of the gene regulatory network in single cellular circadian oscillators [13]. The model represents the transcription-translation behavior of the single cell by using three variables that include a clock gene mRNA (x), a clock protein (y), and a transcriptional inhibitor (z). As our network model, we consider the mean field-modified Goodwin model proposed by [13] with a global coupling strength. The modified Goodwin model with N oscillators is represented as follows: Where the state variables x i , y i , z i represent the concentrations, respectively, of a clock gene mRNA, a clock protein and a transcriptional inhibitor in each clock cell i. Neurotransmitter V i is induced by the mRNA x i . The intercellular coupling is implemented through the neurotransmitter F which acts as a mean field of V i , the coupling strength g i represents the sensitivity of the individual oscillator to the neurotransmitter and is required to be a positive value here, and L i is the light term. We considered the parameters as in [13]: (a 1~0 :7nM=h, k 1~1 :0nM, n~4:0, a 2~0 :35nM=h, k 2~1 :0nM, k 3~0 :7=h, a 4~0 :35nM=h, k 4~1 :0nM, k 5~0 :7=h, a 6~0 :35nM=h, k 6~1 :0nM, k 7~0 :35=h, a 8~1 :0nM=h, k 8~1 :0nM, a c~0 :4nM=h, k c~1 nM): The coupling strength g i is different for different oscillators and assumes a value from a normal distribution with a mean 0.5 and a standard deviation g. When g=0, the network is heterogeneous with distribution of coupling.
In order to understand the dissociation behavior observed under a T cycle outside the range of entrainment, we modified the Goodwin model to include the fact that light acts directly on only a portion of the neurons in the network. Furthermore, the light term L i that is applied to a fraction pN neurons with p being less than one and positive, is considered to be located in the VL subdivision. Mathematically, the effect of light is represented as: Where T L is the period of the light-dark cycle and K f is the light intensity.
As pointed out in our previous paper [22], the dispersion of coupling strengths influences the free-running period of the SCN. In order to compare the influence of different coupling dispersions on entrainment of the SCN network to T cycles, it is necessary to make the free-running period the same for different dispersions.
To set the free-running period to 24 h, we multiply a rescaling factor s to the left hand of equation (1) except for the light term and coupling term, i.e., multiply the same s to the parameters a 1 , a 2 , k 3 , a 4 , k 5 , a 6 , k 7 , a 8 , a c for the deviation g. For example, s is equal to 1.26 for g~0:0, 1.22 for g~0:05, 1.16 for g~0:1, and 1.13 for g~0:15.
For simplicity, we refer to the network that is comprised of pN neurons as VL and the network comprised of the remaining (1{p)N neurons as DM. To understand the behavior of the VL and DM subdivisions, we define the mean field of VL and DM respectively as In addition, to understand the synchronization properties between VL and DM, we have estimated the phase of the individual neurons by using the Hilbert transform [49,50]. From the estimated phase of VL and DM, we introduce an order parameter as: where h is the estimated phase from the mean field output time series of VL or DM and ST denotes average over time. The average of dh dt is defined as the angular frequency and the period is obtained by T~2 p h : . To determine entrainment of VL or DM to the T L , we estimated the period of VL or DM and estimated its absolute difference from T L . If the absolute difference in period is less than 0.25 h, the corresponding subgroup (VL or DM) is considered to be entrained. To numerically calculate the equations, we use the fourth order Runga-Kutta method with time step of 0.1 h. Initial 20000 time steps are neglected to avoid the effect of transients. The number of oscillators isN~100, and the simulations are performed five times, with initial conditions selected randomly from a uniform distribution in the range (0-1) for x, y, and z. We have also calculated the case of N~400 and time step of 0.01 h. Two additional simulations are performed with selective assignment of coupling in which larger values are assigned to either VL or DM without changing the intrinsic distribution of g i . To obtain the phase-response curve, we applied 1-h light pulses at different phases to the model, with intensity K f and with different values for p §p c . The corresponding advance or delay is detected from the output of the network. Advance corresponds to the capacity of the SCN network to follow a light-dark cycle with period less than the free running period, and delay is the capacity of the network to follow a light-dark cycle with period greater than the free running period [35].