Modeling the Seasonal Adaptation of Circadian Clocks by Changes in the Network Structure of the Suprachiasmatic Nucleus

The circadian rhythm is a 24 h rhythm which can be found in many organisms ranging from cyanobacteria and fungi to mammals [1], [2], [3]. There is a huge interest in studying the circadian rhythm because of the well-known effects of jet-lag after traveling and diseases related to shift-work [4], [5]. In mammals the major pacemaker is the suprachiasmatic nucleus (SCN), which synchronizes all peripheral clocks in the body and controls the overall behavior [6], [7], [8]. It is a small region in the hypothalamus located below the third ventricle and directly above the optic chiasm. Light is the major entrainment factor of the SCN. An important dynamical property of the SCN is that it adapts to different photoperiods in summer and in winter [9], [10]. This means that the behavioral activity should be longer in summer days than in winter days, which is advantageous for the organism [11], [12].

The SCN is a symmetric structure, consisting of approximately 20 000 neurons, where each part is usually classified into a ventrolateral (VL) and dorsomedial (DM) region [13], [14], [15]. Cells in the VL region mainly express the neuropeptide, vasoactive intestinal peptide (VIP), comprising around 24% of the SCN neurons in rats [15], and receive the light information from photosensitive retinal ganglion cells [7], [16]. However, also other cell types have been identified that show light induced gene expression. For example in hamsters, cells expressing calbindin (CalB) are responsible for the mediation of light information [17]. Importantly, it has been found that these cells are nonrhythmic and uncoupled from each other [17], 18,19,20. In contrast to this cells in the DM region, expressing mainly vasopressin (AVP) are rhythmic but do not receive light input [17]. It has been shown that separating a part of the DM region leads to nonsynchronous rhythms in the individual cells [21], indicating that the coupling in the DM region is not sufficient to synchronize these cells. Indeed it has been shown that interactions between cells in this region are restricted to short-range connections [22]. On the other hand long range connections of neurons in the VL region to neurons in the DM region have been identified [13], [23], which seem to be important to synchronize cells in the DM and VL region to each other and ensure entrainment of DM cells to an external light-cycle. We are aware that the distinction into a VL and DM region is an oversimplification and not that clear in other species than rats since the distribution of neurotransmitters and retinal inputs shows a more complex spatial organization [17], [24], [25]. Nevertheless, since our study is a first attempt to model the influence of long-range couplings we aim at a simple and manageable model. Moreover, the separation into a “core" part, represented by the VL cells, and a “shell" part, represented by DM cells, is commonly accepted also for other rodents [17].

The importance of the neuronal network circuitry mediated by chemical synaptic interactions as opposed to simple global coupling has been demonstrated in several experimental studies [13], [21]. All these studies show that the functioning of the SCN and the regulation of the circadian rhythm in general is based on a very complex neuronal network consisting of short- and long-range connections. The question arises how the structural properties of this complex neuronal network are linked with the electrical activity of the SCN, which is known to crucially affect the circadian gene expression [26].

It is known that the seasonal adaption of the SCN is closely related to its electrical activity [9], [10]. In particular, single cells in the SCN have a peak in electrical activity (measured as firing activity) from 4–5 h, regardless of monitoring in winter or summer conditions [9], [27], [28], [29], [30], [31]. Therefore, it has been suggested that the phase distribution of electrical activity inside the SCN neurons, which is narrower in winter than in summer, leads to a shortened behavioral activity in winter [10]. Different phase distributions in winter and summer have been previously modeled by introducing delay times in the synaptic connections [32].

Here we hypothesize that the adaptation to seasonal changes in the phase distributions can also be related to changes in the structural properties, i.e., the topology and coupling of the complex neuronal network of the SCN. This is a reasonable hypothesis since, as already mentioned in the previous paragraph, the individual neurons do not change their electrical activity in summer and in winter. Many studies have already examined the network topology of the SCN in various manners [32], [33], [34], [35], . These studies use different single cell models that range from a generic van-der-Pol oscillator [32], [34] to more detailed biochemical models [36], [37], [38]. Moreover, they model coupling between the cells in different ways. Whereas some studies considered homogenous coupling between all cells [33], [35] others take into account the heterogeneity in the SCN cellular network [34], [36], [37], [38]. Furthermore, Vasalou et al. [37] have reported that a small-world network architecture of the SCN can firmly mimic the dynamical behavior of mean-field coupled models, but is on the other hand much more efficient in terms of connectivity cost. Also Hafner et al. [38] in a recent study analyze different network topologies with respect to rhythm output and jet-lag adaptation and find that coupling different network topologies leads to robustness of the overall rhythm with respect to perturbations. Nevertheless, all of the abovementioned studies were mainly focused on the synchronization and amplitude properties, entrainment and robustness of the SCN.

In our study, we focus on the role of long-range connections between the neurons, as they are known to lead to networks characterized by small-world properties [39], [40]. Moreover, in jet-lag experiments evidence was found for a connection between the ventral and dorsal part of the SCN [10], indicating a role for long-range connections between both parts. Using our model, defined by coupled ordinary differential equations, we show that the number of long-range connections between the cells in the VL and DM region is a fine-tunable parameter to adjust the phase distribution and consequently the length of behavioral activity. Our results thus indicate that the seasonal summer/winter dynamics of the circadian clocks can effectively be regulated by the plasticity of the SCN network structure.

We entrained our network model over 40 cycles with a 24 h rhythm and variable photoperiod lengths. In all numerical calculations the initial conditions were randomly distributed around x = 1 and y = 0, according to a normal distribution with a standard deviation 0.2. In all simulations, the final dynamical state was attained after just a few periods (see Figure S1 in Text S1). Therefore, we chose to discard an initial transient of 15 cycles. In Figure 2 we show the time courses of activity in the x-coordinate of 10 randomly chosen cells for and . This parameter affects the number of long-range connections in the SCN network. It should be emphasized that the coupling of the cells in the network leads to an approximately 10-fold increase in the amplitude of the activity of the single cells, as observed also experimentally [21]. Moreover, it can be seen that decreasing and thus the number of long-range connections leads to a broadening of the peak-phase distribution (see also Figure S2 in Text S1), whereas the activity peak of the individual cells always has a width of only about 4–5 h as observed experimentally [10].

To capture the effect of the broader phase distribution on the overall electrical activity of the SCN we calculated the mean activity over the entire SCN and shifted this activity into the positive quadrant (Figure 3A). It can be seen that the mean electrical activity is broader for small values of , which we therefore refer to as “summer-topology" and is shorter for large values of , which we therefore refer to as “winter-topology". Therefore, our model is able to explain the different experimentally observed phase-distributions in summer and in winter [10], by changes in the network topology. It has also been observed that after entrainment to different photoperiods the width of the phase distribution is preserved under free-running conditions as well. In Figure 3B we show the network's dynamical response without an entraining stimulus for the winter- and summer-topology. We can observe that even without entraining the phase distribution is also preserved under free-running conditions. Thus we can conclude that the width of the activity is mainly related to SCN plasticity [9]. Moreover, even in the absence of light input the non-rhythmic cells in the VL region become rhythmic due to the coupling to the oscillators with a self-sustained rhythm in the DM region as already observed in a previous study on synchronization induced rhythmicity in the SCN [36]. Another surprising observation is that a coupling strength of is sufficient to synchronize the oscillators, despite the significant heterogeneity in the oscillator periods. We attribute this to the weakness (small radial relaxation rate ) of the oscillators.

The sensitivity of the amplitude and width w (defined as the width of the peak at half the amplitude) of the mean electrical activity with respect to the parameter is shown in Figure 4. It can be seen that by changing the number of long-range connections the width of the mean electrical activity can be adjusted to winter and summer photoperiods. When moving from winter- to summer-topology also the mean activity amplitude decreases. To ensure that the width and amplitude of the mean activity is mainly determined by the network topology we also calculated them for short and long photoperiods. It can be seen that both only slightly depend on the length of the photoperiod (see Figure 4). Remarkably, the same conclusions can be drawn from results presented in Figure S3 in Text S1, where the mean-field of the summer () and winter topology () for different durations of the light input is shown.

To quantify the synchronization between the neuronal cells we calculated the correlation matrices for the summer- and winter-topology under the corresponding entrainment cycle (Figure 5A–B). The ij-th element of the matrix is defined as follows:

where is one of the Cartesian coordinates which signifies the neuronal activity and and are the mean value and the standard deviation of the time series of the i-th oscillator, respectively. The sum in Eq. (11) runs over the whole temporal series, whereby every 25th point of resulting from numerical integration is recorded. It can be seen that the cells receiving light input in the VL region are more synchronous to each other for both topologies. On the other hand, cells in the DM region in the summer-topology are asynchronous to each other reflecting the broader peak phase distribution, but are more synchronous to the cells in the VL region probably to the ones they have long-range connections to. For the winter-topology also the cells in the DM region show good synchrony to each other, which reflects the narrow phase distribution observed in short photoperiods [10]. The overall better synchrony with more long-range connections in the winter-topology also reflects itself in the average correlation coefficient (Figure 5C). Besides, the level of correlation between cells only slightly depends on the duration of the photoperiod, hence once again confirming that the intercellular network structure is the key agent governing the characteristics of the collective neuronal activity.

Furthermore, it has been experimentally observed that the PRC of the SCN shows a larger magnitude in winter than in summer [63]. We have numerically calculated the PRC for the winter- and summer-topology (Figure 6). In accordance with the experimental observations, the PRC for the winter-topology has a higher magnitude then the summer-topology. Notably, the extent of the phase shift under both conditions is regulated by the amplitude of the perturbation signal (see Figure S4 in Text S1). To further analyze the underlying reasons for the differences, we calculated the individual instantaneous PRCs of all light-receiving oscillators (Figure 7A–B), which when summed up represent the overall instantaneous PRC (Figure 7C–D). It can be seen that the individual iPRCs are also diminished in magnitude in the summer topology compared to the winter topology. Moreover, the individual iPRCs are more shifted to each other in the summer topology (correlation coefficient compared to ), most likely due to the lesser synchronization between the individual light-receiving cells. Both effects lead to a diminished overall PRC. The effect due to shifting of the curves is enhanced for even smaller values of : for (Figure S5 in Text S1).

To characterize the different network topologies in winter and in summer we calculated the network efficiency and clustering coefficient, two properties that allow the characterization of small-world networks (Figure 8). It has been shown that the product of the two measures is the largest when the network possesses small-world properties 57,64. Using this measure characterizes the winter-topology as a small-world network, which is efficient in terms of costs for the neuronal connections and the achieved synchrony between the individual cells. On the other hand the summer-topology is not able to achieve synchrony due to the reduced number of long-range connections, which lessens the small-world properties of the network.

To additionally characterize the synchronization behavior of the network we calculated the singular value spectrum of the network's Laplacian matrix L for 25 replicates of our network structure (see Section 5 and Figure S6 in Text S1). For a symmetric coupling it has been shown that large eigenvalues indicate groups of oscillator that synchronize fast, whereas near zero eigenvalues indicate a community structure inside the network, with groups of oscillators that do not synchronize their phases to each other [59], [60]. In Sections 2, 3 and 4 in Text S1 we generalize these numerical findings to a network of weakly coupled heterogeneous oscillators and show analytically that for a non-symmetric coupling small near-zero singular values lead to a large variance in the synchronized phase distribution. Although, the theory is established for a network of only self-sustained oscillators the numerical calculations indicate that its predicted effects are still applicable for the here considered network of damped and self-sustained oscillators. The singular value spectra for the summer topology possesses many near-zero eigenvalues indicating that indeed a community structure is present, where oscillators within a community are more connected to each other than to members from other communities (see Section 5 in Text S1). Our theoretical results now also explain how this effects the stable distribution of phases. In particular the near-zero eigenvalues present for the summer topology lead to a dramatic increase in the width of the phase distribution. To illustrate this even more we have also calculated the mean number of near zero singular values over different values of (Figure 9). It can be seen that for only the trivial zero singularvalue exists, whereas at lower values of small non-zero singularvalues indicate groups of oscillators that do not fully synchronize their phases to each other and therefore lead to a broader SCN activity.

To further characterize our model we also analyzed the entrainment properties of our model. In Section 5 in Text S1 we analyze the entrainment dynamics of a single, uncoupled oscillator in an analytic way and derive, several explicit formulae for the entrainment borders (see also Figure S7 in Text S1). These results shows that the individual weak, spiking single oscillator we consider shows an entrainment behavior that markedly differs from the entrainment of a rigid amplitude-phase oscillator. Especially, after a certain threshold of the entrainment amplitude b the external forcing overrides the internal dynamics and the oscillator is practically entrainable to every period . We therefore compared the entrainment region of the single uncoupled oscillator (Figure S7 in Text S1, dotted lines) with the entrainment region of the oscillator network in summer and winter (Figure 10A–B). Note that due to the high computational costs, we considered only one network realization per condition. Therefore, we expect that the borders of the entrainment region slightly change with each realization. Nevertheless, even the calculations for just one realization give valuable insights. We find that due to the coupling the oscillator network behaves more like a rigid oscillator as found also in a previous study [44]. Moreover, the entrainment amplitude b needs to be significantly higher compared to the single oscillator to achieve entrainment. This effect is related to the amplitude expansion of the coupled oscillators. This is confirmed by calculations of the entrainment region for an oscillator network with increased coupling strength (Figure 10C–D). The entrainment region becomes even smaller due to increased rigidity and amplitude of the oscillator network in line with the previous results [44].

Nevertheless, we observe that the entrainment region of the winter topology is larger than that of the summer topology. This result is in line with the observed reduced phase response curve of the summer topology (Figure 6). Another interesting observation is that for a small coupling strength the entrainment region is not symmetric for large entrainment amplitudes and especially, entrainment to smaller periods is not possible.

In order to verify that our findings are indeed qualitatively independent of the system size we calculated how the amplitude and the width of the mean electrical activity change with respect to the parameter delta for a three times larger system size (1800 neurons, 600 of them being located in the VL region). First, we computed at which values of delta small-world network properties are obtained. Results shown in Figures S8A–B in Text S1 reveal that those characteristics, which reflect winter conditions, are found around . Accordingly, the SCN network structure reflecting summer conditions was chosen to be at . Remarkably, the results in Figures S8C–D in Text S1, showing the changes of w and h as a function of delta, clearly indicate that qualitative identical results are obtained for larger system sizes as well (compare to Figure 4).

Finally, we examine how the SCN network behaves in the case that a more complex and biologically relevant model – the Goodwin oscillator – drives the dynamics of individual cells. Results presented in Figures 11A–B show the time evolution of the concentration of the clock gene u in summer () and winter () conditions. Obviously, the phase distribution is much broader in summer conditions than in winter conditions. This observation is additionally confirmed with the results presented in Figure S9 in Text S1, where the time evolution of the mean SCN activity over several entraining cycles is shown. Similar as in the case of the simple amplitude-phase oscillator (see Figure 3), the mean field signal of the Goodwin oscillator network is higher and narrower in winter. Furthermore, Figure 11C shows the average correlation coefficient as a function of the network parameter, whereby a greater extent of synchronization can be observed as is increased. Thus, the results obtained with the more complex model for circadian oscillations are qualitatively very similar to those obtained with the simple amplitude-phase oscillator.

The importance of network plasticity for the adaptation of the circadian rhythm to different photoperiods was put forth in several experimental studies [9], [10]. Using our model, we identified one possible mechanism able to explain the adaptation to different photoperiods: the introduction of long-range connections between cells in the VL (ventrolateral) and DM (dorsomedial) region of the SCN. In particular, the length of the behavioral activity can be regulated as follows: dense long-range connections during winter lead to a narrow activity phase, while rare long-range connections during summer lead to a broad activity phase. We found that this result is independent of the number of neurons. Moreover, similar results were obtained when replacing the model governing the dynamics of the individual neurons, namely the simple amplitude-phase oscillator, with the more biologically relevant Goodwin oscillator. To go from a wide phase distribution in summer to a more narrow and synchronized phase distribution in winter roughly the doubled number of long-range connections had to be introduced. In our model of 600 neurons, this amounts to around 450 connections in the summer topology and 900 connections in the winter topology. Taking into account that the fraction of added long range connections leading to a network with small world properties is inversely proportional to the number of oscillators [57], [58], we can extrapolate the number of required long-range connections in the real SCN network consisting of 20000 neurons. In this case, roughly 15000 connections would need to be introduced to change the network from the summer-topology to the winter-topology.

It should be emphasized that the photoperiod gradually changes from summer to winter and hence the introduction of the additional neuronal connections would occur in a time span over roughly half a year. In fact the physical connections between neurons do not need to be introduced but synaptic plasticity changing the responsiveness of individual neurons could lead to a weakening or strengthening of connections. To support those ideas we also calculated how the dynamics of the SCN network changes if the coupling strength of the long-range connections in a winter topology network withis gradually reduced. The simulation results indicate that reducing the coupling strength by at least 75% for half of the long-range connections (reflecting the summer topology with) leads to very similar effects as a complete deletion of links (Figure S10 in). Moreover, a simulation with a rapid transition from winter to summer topology shows that the entrainment transient is very short and therefore negligible for gradual changes over the year (Figure S11 in). Text S1 Text S1

Our model is also able to simulate the increased magnitude of the overall phase response curve (PRC) in winter compared to summer. This can be explained by the reduced magnitude of the individual instantaneous PRCs and their shifting relatively to each other in the summer topology, leading to an overall smaller phase response. Both effects are due to the reduced synchrony in the summer-topology shown by the correlation matrices in Figure 5. Since the cells are at different phases of their inner clock when not completely synchronized this leads to a smaller overall response, much like is the case for coupled but unsynchronized pendulums that are perturbed by a pulse. Of course, the individual PRCs do not contain a “dead" zone with no phase shift. The reason for this is the simple model for the individual oscillators, which cannot account for the complicated phase response behavior of more realistic models. However, the focus of this work was not on modeling the individual oscillators but analyzing the impact of network changes on the synchrony of the overall SCN and how photoperiod adaptation is related to this.

Along with the phase response, we also analyzed the dependence of the entrainment region of the oscillator network on the network topology and the coupling strength between the individual neurons. Our results agree with previous studies that found a decreased entrainment region for increased coupling between the oscillators, due to an increased in rigidity and amplitude of the oscillators [44]. Moreover, in line with the observed reduced PRC of the summer topology we also observe that the entrainment region of the summer topology is smaller. This result is counterintuitive, since intuitively one would expect that the oscillators become more rigid in a network with more long-range connections. In addition, the amplitude of the oscillator networks mean field is larger in the winter topology as compared to the summer topology. A common statement in the analysis of oscillator synchronization is that large amplitude oscillators are harder to entrain. Our results show that this statement cannot simply be carried over to oscillator networks. It seems that the effect of adding specific connections between cells is different from an increase in the overall coupling strength. Whereas, the latter leads to a decrease in the entrainment region, due to amplitude expansion and increase in rigidity, the former leads to a larger entrainment region due to a more synchronized phase response. These findings underline the importance of the network structure connecting the oscillators.

Another possibility to achieve different phase distributions in summer and winter are changes in the intrinsic oscillator characteristics or in the transmission of signals between the cells. For example the introduction of a delay distribution in the synaptic transmission also leads to a distribution of peak phases [32]. Changing the delay distribution could also allow the adaptation to different photoperiods. However, it is difficult to explain a delay of several hours physically since it has been observed that the transcriptional induction of Per upon a light pulse is within 5–15 minutes during subjective night [13]. Since the same pathway induces transcription in cells not sensing light but neuropeptide release, we can assume that the delay times are on the same order of magnitude also for intercellular coupling via neuropeptides.

Our results support the ideas pointed out by Meijer et al. [10] in that the SCN neuronal network plasticity crucially affects the activity phase distribution among SCN neurons and therefore contributes to the adaptation to changes in day length. Interestingly, our results indicate that the intercellular communication network in the SCN has features of a small-world network. Such complex topologies have been identified in numerous biological systems including functional as well as anatomical connections in the nervous system [65]. They indeed seem to be advantageous for various living organisms. A closer inspection of Figures 4A–B and 5 reveals, that below the region where the small-world properties of the network are well expressed (, see Figure 8) the curves are very steep. For the SCN it would thus be advantageous operating in the proximity of the optimal small-world configuration. Namely, rather small modifications in the extent of long-range connections enable the regulation of synchronization behavior and phase distributions of electrical activity in the neuronal population. In this manner the plasticity as well as the realization of different neuronal coupling mechanisms between the ventral and dorsal SCN have a large impact during the adjustment to seasonal changes.

The effect of the broad phase-distribution in the summer topology is explained by the introduction of well-connected communities separated by bottlenecks into the network. Local clusters of SCN oscillator synchronize fast and well to each other. However, the synchronization between these clusters is hindered due to the few connections between them. The community structure of a network can be quantified by the so called algebraic connectivity given by the second smallest eigenvalue of the networks Laplacian matrix (see Eq. 5) [66], [67]. If this eigenvalue is near zero the network can be easily separated into groups. It was shown numerically that near-zero eigenvalues lead to a more unsynchronized state of the oscillators [60], [68]. Our analytical results (Sections 2, 3 and 4 in Text S1) support these findings and show that the distribution variance of oscillator phases is strongly influenced by near-zero eigenvalues. Remarkably, these theoretical results do not depend on the underlying model governing the dynamics of individual cells and therefore generalize the proposed mechanism of photoperiod adaptation by controlling the number of long-range connections and consequently the community structure in the SCN network. Another advantage of the theoretical results is the insight into the relation between synchronized phase and period distribution on the one hand and network properties on the other hand (see for example Eq. 18 in Text S1). Future studies could use this relation to extend our model to account for other features of the phase distribution in the SCN, for example the observed bimodal phase distribution in long photoperiods [30], [31].

This research was supported by the Slovenian Research Agency (Grant J1-4055). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.