Quantitative Study of Dual Circadian Oscillator Models under Different Skeleton Photoperiods

Materials And Methods

To investigate the effects of photoperiod on the E-M model, we used 2 coupled nonlinear limit-cycle oscillators, governed by the Pittendrigh-Pavlidis equations. These equations have been previously used to model circadian rhythms in Drosophila (Pittendrigh, 1981; Pittendrigh et al., 1991) and the E and M oscillators in rodents (Oda and Friesen, 2002; Oda et al., 2000). Each oscillator is described by a set of 2 equations (1-4). Following previous conventions (Oda et al., 2000), the terms from E and M oscillators are identified by subscribed letters.

Evening oscillator (E):

d R E / d t = R E − c E S E − b E S E 2 + ( d E − L ) + K

(1)

d S E / d t = R E − a E S E + C M E S M

(2)

Morning oscillator (M):

d R M / d t = R M − c M S M − b M S M 2 + ( d M − L ) + K

(3)

d S M / d t = R M − a M S M + C E M S E

(4)

State variables R and S together describe the state of one oscillator at a given time point. The R variable is prevented from reaching negative values. Parameters a, b, c, and d are set to fixed values and collectively define an oscillator configuration, with intrinsic period and amplitude. K (Kyner) is a small nonlinear term that guarantees numerical smoothness (K = k₁/(1 + k₂R²), k₁ = 1, k₂ = 100; Oda et al., 2000). The L term simulates the light input (forcing variable). Its value is maintained at zero to simulate the dark hours. When lights are turned on, it is changed to a positive amplitude value in arbitrary units (a.u.), which simulates the light intensity. In one particular case indicated below, we used pulses of negative amplitude. The coupling term C controls the strength with which the E oscillator influences the morning oscillator (C_EM) and vice versa (C_ME). When coupling is symmetric, C_EM = C_ME.

In most cases, we used a standard configuration (a = 0.85, b = 0.3, c = 0.8, d = 0.5), previously employed to model circadian rhythms in rodents (Oda et al., 2000). Other configurations have been explored to investigate the dependence of phenomena on the period and amplitude of oscillators (see the supplemental materials). Computer simulations were performed in the CircadianDynamix extension of the Neurodynamix II software (Friesen and Friesen, 1994) using the Euler method for numerical integration, with 1000 integration steps per 24 h.

In our simulations, we portray the phases of the circadian oscillators (squares and triangles) and of the behavioral rest/activity rhythm (black bars) separately (Fig. 2), to emphasize the temporal interdependence between oscillator and behavioral activity. In Figure 2, we explicitly show how the phases corresponding to activity were transduced from the phases of the circadian oscillator. This was achieved by defining a threshold, which was set to 40% of the state variable amplitude. Whenever the S variable of either oscillator went above its threshold, the model organism was considered as resting (activity was inhibited); otherwise, it was considered as active. Thus, in the dual-oscillator model, the total activity band was determined by E and M together (Fig. 2), with activity set in antiphase to the peak amplitude of the circadian oscillators. This is equivalent to the phase relationship between behavioral activity and peak SCN electrical discharge in nocturnal mammals (Houben et al., 2009).

OPEN IN VIEWER

Figure 2.

Conversion of oscillators’ state variables to daily activity/rest rhythms in the simulated animal. The process is illustrated for a 48-h section of a simulation record. The oscillator is exposed to skeleton photoperiods indicated by “LD cycle” at the top, with white squares representing light pulses in a darkness (gray) background. The state variable S from the evening (S_E) and morning (S_M) oscillators is recorded continuously and compared with the activity thresholds. The peaks of the state variables are also recorded and represented as squares (evening) and triangles (morning). Whenever the state variables from both oscillators are below the threshold, the model animal is considered active; otherwise, it is set as inactive. This preliminary activity phase as well as the state variable peaks are represented by dashed symbols in the section “Activity (before 5h lag).” The pulses of the skeleton photoperiod are indicated as white squares. An addition step is made to better fit the phase to experimental data reported in the literature. We set a 5-h lag (arrows) to the composite phases of oscillators and activity, relative to the light-dark cycle, as indicated in “Activity (after 5h lag).” The data for consecutive days is plotted one day below the other in the “Actogram.” Each day represents 48 h of data, with days 1 to 2 in the first line, 2 to 3 in the second line, and so forth. In the actogram, state variable peaks are shown only every 3 days to make visualization easier. Color versions are available online.

Finally, an extra step was added in the simulations to better fit experimental data for nocturnal animals under skeleton photoperiods. The composite phases of oscillators and activity were delayed by 5 h throughout all simulations, establishing an interval lag relative to the light signal. Importantly, we maintained the same fixed phase relationship between oscillators and activity rhythm (Fig. 2). This step was necessary because although activity patterns in the simulations matched closely to experimental data, there was a discrepancy between the exact time of their occurrences. This discrepancy was resolved by setting a phase difference between the light pulses and the simulated oscillator.

The daily light/dark cycle was modeled as a skeleton photoperiod, with 2 light pulses every 24 h (Fig. 1B, top; Pittendrigh and Daan, 1976a; Spoelstra et al., 2014). We simulated 4 configurations of the model to represent different putative physiological mechanisms (Fig. 3, top). As a control, in model I, we first simulated a single oscillator subjected to the skeleton photoperiod (isolated E, with null coupling; Fig. 3A, top). All of the remaining configurations had both E and M oscillators, with symmetric coupling. In model II, each of the 2 pulses of the skeleton photoperiod was applied to only 1 of the oscillators (Fig. 3B, top). In model III, both pulses were applied to both oscillators (Fig. 3C, top). Finally, in model IV, the M oscillator received pulses with negative amplitude (Fig. 3D, top).

OPEN IN VIEWER

Figure 3.

Four configurations of the model and their entrainment to different photoperiods. (A-D, top) Schematic representation of model configurations depicting the skeleton photoperiods (gray-white bars), the oscillators (circular symbols), and the effects of the light pulses on the oscillators (arrows). When 2 oscillators are present (models II-IV), curved arrows indicate coupling between them. (A-D, bottom) Double-plotted actograms illustrating the activity rhythm of the model organism. Plotting conventions are the same used in Figure 2. (A) Model I, single oscillator responsive to both pulses of the skeleton photoperiod. (B) Model II, each oscillator is responsive to one of the pulses of the skeleton photoperiod. (C) Model III, both E-M oscillators respond to both pulses of the skeleton photoperiod. (D) Model IV, same as previous model, except that the M oscillator receives pulses of negative amplitude (−), while the E oscillator receives regular pulses of positive amplitude (+). The days during the psi-jump are highlighted with a vertical line on the right of each actogram. (E-L) Phase of entrainment and activity duration of the 4 model configurations, under different skeleton photoperiods. Measurements were taken on the last day of exposure to each photoperiod, with stable entrainment of the activity rhythm. (E-H) Phase of entrainment. The times of activity onset and offset (black-filled circles) are measured relative to the midpoint of the interpulse interval, in a circular scale (x-axis). Onset and offset points are connected by a black dashed line. Diagonal white lines indicate the times of the pulses, as the skeleton photoperiod changes in the y-axis. Peaks of the state variables are indicated for the evening (squares) and morning (triangles) oscillators in each photoperiod. (I-L) Activity duration (α) under different photoperiods. Entrainment was not stable (missing points) at 12-h “night” intervals in E, H, I, and L and at 10-h “night” intervals in F, G, J, and K. In all versions of the model, oscillators have the standard configuration: a = 0.85, b = 0.3, c = 0.8, d = 0.5. Coupling is symmetric and positive (C_EM = C_ME = 0.07). Pulse duration: 1 h; pulse amplitude: 1.1 or −1.1 (arbitrary units). Color versions are available online.

A model in which one oscillator receives light input and the other does not was not considered, because the light-insensitive oscillator would just passively follow the phase of the sensitive one, with negligible changes in the phase relationship.

The 4 configurations of the model were exposed to different skeleton photoperiods, with different intervals between the daily light pulses. Intervals were reduced from 22 to 2 h, in 2-h steps applied sequentially, each for 15 days. Importantly, this 15-day interval ensures the end of transients and that the entrainment phase achieved by the circadian oscillator is the steady-state phase for each different photoperiod. The lengths of the intervals were measured from the beginning of the first pulse to the beginning of the second pulse. On the last day of exposure to each photoperiod, we measured α and the phase of entrainment. α was measured as the interval between the daily onset and the daily end of activity. The phase of entrainment was measured as the interval between the midpoint of pulses and 2 reference phases of the rhythm, activity onset and offset, and then converted to the 0- to 24-h circular scale.

Simulation Results

We tested the photoperiodic response of the 4 putative models. When exposed to different skeleton photoperiods, all configurations of the model are successfully synchronized to 24 h. Representative actograms are shown in Figure 3A-D (bottom). For the single oscillator in model I (Fig. 3A), only the E oscillator is simulated. The actogram in Figure 3A (bottom) depicts a daily activity band in black, with the 2 pulses of the skeleton photoperiod in white and the oscillator variable peaks as empty squares. We assign “dusk” to the pulse that is locked to activity onset and “dawn” to the other pulse. When the “night” interval between pulses is decreased from 22 to 14 h, the onset of activity remains phase locked to the dusk pulse. This is summarized in Figure 3E, which shows the phase of entrainment quantified on the last day of exposure to each photoperiod. With a single oscillator, there is no noticeable change in α between photoperiods, as observed in the actograms (Fig. 3A, bottom) and quantifications (Fig. 3I). When the interval between pulses is decreased to 12 h, a psi-jump is observed (Fig. 3A, bottom), and from that point onward, the activity onset follows the dawn pulse of the skeleton photoperiod (Fig. 3A, bottom; Fig. 3E). In summary, with 1 oscillator, we can observe the psi-jump phenomenon, but α remains unchanged throughout different photoperiods.

Having observed the behavior of a single oscillator, we can then explore putative E-M mechanisms. In model II, the actogram depicts the activity band, together with squares and triangles representing E and M state-variable peaks, respectively (Fig. 3B, bottom). For simplicity, we chose to model a completely symmetric system in terms of oscillators and coupling, which leads to an overlap between E and M peaks in constant conditions (Fig. 3B, bottom), as reported previously (Oda and Friesen, 2002). When the system is exposed to different photoperiods, again, the activity onset is phase locked to the dusk pulse, up to the 12-h “night” interval (Fig. 3B, bottom; Fig. 3F). Furthermore, the E and M peak phases separate from each other upon initiation of the skeleton photoperiods (squares and triangles in Fig. 3B, bottom), and this becomes more pronounced as the “night” interval between pulses is decreased from 22 to 12 h (Fig. 3B, bottom). Consequently, there is a corresponding decrease in α (Fig. 3B, bottom; Fig. 3J). When the interval between pulses is further decreased to 10 and then 8 h, there is a psi-jump (Fig. 3B, bottom; Fig. 3F). After this, α is gradually increased (Fig. 3J), and the activity onset tracks the opposite, dawn pulse of the skeleton photoperiod (Fig. 3F). These results indicated that model II accounts for photoperiod-induced α adjustments.

One feature observed experimentally in α adjustments is the maintenance of α on the first days after an animal is released in constant darkness from a photoperiod protocol (Pittendrigh and Daan, 1976a). We released the oscillators of model II in constant conditions, after being held in different photoperiods, and observed that α was very similar on the first day. Throughout the days in constant conditions, α gradually lengthened to a steady-state value (Suppl. Fig. S1).

We then explored variations in the coupling strength between E and M in model II (Fig. 4). With very weak coupling, each oscillator tracks one of the pulses, following the pulse times almost passively (Fig. 4A). This results in a great modulation of α (Fig. 4A). Upon a slightly stronger coupling (Fig. 4B), the oscillators remain closer together in the more extreme photoperiods, with a more stable and larger α when the interval between pulses is too long. As the “night” interval between pulses is changed to intermediate photoperiods, E and M separate from each other, and α is modulated (Fig. 4B). With stronger coupling, higher than the standard value used in Figure 3B, E and M stick together throughout most photoperiods (Fig. 4C, D), initially following the dusk pulse and then shifting to the opposite pulse. Consequently, α remains more stable throughout photoperiods. Strikingly, when we compare the weaker coupling with the stronger ones, there is a change in the threshold photoperiod before the psi-jump. With weaker coupling, psi is reversed when the “night” interval between pulses is set to 10 h (Fig. 4A, B), whereas with higher coupling strengths, the psi-jump happens at “night” intervals of 8 and 6 h (Fig. 4C, D).

OPEN IN VIEWER

Figure 4.

Effects of the coupling strength between evening (E) and morning (M) on the phase of entrainment for model II (Fig. 3B, top). Graphs are plotted following the same conventions of Figure 3E-H. (A-D) Different (symmetric) coupling strengths, as indicated above the graphs. Other model parameters as in Figure 3. Coupling determines how independently the phases of the oscillators change with different photoperiods. Coupling also affects the threshold photoperiod before the psi-jump. Color versions are available online.

The E-M model above is responsive to photoperiod variations; however, it would imply that, in the SCN, distinct subregions corresponding to E and M responded separately to dawn or dusk inputs. Alternatively, it is possible that E and M subregions are equally responsive to dawn and dusk stimuli. To test this hypothesis, in model III, we again have a completely symmetric set of E and M oscillators, but both oscillators are directly synchronized by both pulses of the skeleton photoperiod (Fig. 3C, top). When exposed to the different photoperiods, the E and M phases remained completely overlapped, as seen from the squares and triangles (Fig. 3C, bottom). This overlap remained regardless of the interval between pulses (Fig. 3C, bottom). As a result, α remained basically unchanged (Fig. 3K). The psi-jump occurred just like in the previous versions of the model (Fig. 3G). In short, a completely symmetric E-M model behaves essentially as the single oscillator, when both E and M are equally responsive to the pulses of the skeleton photoperiod. We then hypothesized that if a phase difference between E and M existed in model III, it could generate different responses to light. As reported previously (Oda et al., 2000), a phase difference can generate separate responses to single light pulses in a DD background (Elliott and Tamarkin, 1994). To this end, simulations were performed with an asymmetric configuration. We used 2 oscillators with different intrinsic periods (Suppl. Fig. S2A, B), which results in a basal phase difference in constant conditions (Suppl Fig. S2C). Despite this basal phase difference, E and M phases still overlap in the presence of 2 light inputs in a model III configuration (Suppl. Fig. S2C, D). As shown in Supplemental Figure S2E, each individual oscillator, when isolated, presents very similar phases of entrainment to the different skeleton photoperiods. This explains why, in the coupled system, the phases of the 2 oscillators are brought together by the pulses. Thus, model III does not allow independent E and M shifts that could generate α changes.

The original E-M conceptual model suggested that the component oscillators should have different intrinsic properties, such as period and phase response to light (Pittendrigh and Daan, 1976a, 1976b). As a result, each component oscillator would phase lock to a different pulse of the skeleton photoperiod. As the interval between pulses changed with different photoperiods, the oscillators would tend to separate from each other and promote α adjustments. In the asymmetric E-M case described above, however, both oscillators tracked the dusk pulse, despite their asymmetries (Suppl. Fig. S2E). In this sense, we hypothesized that a fourth version of the model could be made responsive to photoperiodic variations if we added novel, specific asymmetries between E and M, so that each oscillator would track a different pulse of the skeleton photoperiod when isolated.

After several attempts, as described in the supplemental material (and Suppl. Figs. S3–8), we arrived at model IV (Fig. 3D, top). Model IV is a variation of model III, in which the M oscillator receives negative pulses. This guarantees that the isolated E and M go to opposite directions when the photoperiod is changed, because of a change in the phase response to the pulses (Suppl. Fig. S8). When isolated, E behaves as the single oscillator in Figure 3A, following the dusk pulse, while M follows the dawn pulse (Suppl. Fig. S8). In the coupled system (Fig. 3D, top), the free-running rhythm is again equivalent to E-M models II and III (Fig. 3D, bottom). After exposure to skeleton photoperiods, the E and M activities separate (Fig. 3D, bottom). Moreover, the phase relationship between E and M is adjusted with varying intervals between pulses (Fig. 3D, bottom). Hence, there is a change in α, especially in photoperiods close to the 12-h interval (Fig. 3L). As with all previous models, the compound system first tracks one pulse, then there is a psi-jump when the “night” interval is decreased beyond 10 h, and the system tracks the opposite pulse (Fig. 3H). Altogether, these results indicate that the addition of specific asymmetries can rescue a photoperiodic response in a coupled E-M model in which both oscillators are directly synchronized to both pulses of the skeleton photoperiod.

Model Limitations

Our dual-oscillator model, with only 2 emergent components, has been enough to coordinate our understanding of photoperiod-induced α changes and psi-jumps, which are phenomena observed at the behavioral level. However, several experiments at the neural level indicate that the photoperiodic adjustments in the rhythm of the compound nuclei result from the phase distribution of the rhythms from single SCN cells or subregions composed of groups of cells (Brown and Piggins, 2009; Buijink et al., 2016; Inagaki et al., 2007; Jagota et al., 2000; Naito et al., 2008; Rohling et al., 2006; Schaap et al., 2003; VanderLeest et al., 2007; Yoshikawa et al., 2017). It is likely that E and M correspond to SCN subregions, formed by coherent groups of cells. The particular properties of each oscillator may arise not only from variations in intracellular mechanisms (Nuesslein-Hildesheim et al., 2000) but also from differences in intercellular coupling within cell groups (Evans et al., 2013; Farajnia et al., 2014; Lucassen et al., 2012). A model that included changes in single-cell and oscillator network properties could certainly contribute to elucidating the explicit mechanisms and refined aspects of photoperiod encoding in multioscillatory organs. For example, an in silico study on evolving genetic networks suggests that a greater complexity (more oscillators and interactions) can be advantageous for a proper response to photoperiod when it is accompanied by noise fluctuations in environmental light (Troein et al., 2009).

Different photoperiods impose modifications within the fine structure of circadian oscillators, which persist under constant conditions (Mrugala et al., 2000). For instance, pineal glands from birds, explanted in vitro, express melatonin patterns that match the previous photoperiodic conditions experienced by the animals in vivo (Brandstätter et al., 2000). In mammals, it has been shown that, unlike the wild-type mice, vasoactive intestinal polypeptide (VIP) knockout mice fail to display aftereffects of entrainment to different photoperiods (Lucassen et al., 2012). Since VIP is a peptide that promotes interneuronal coupling in the SCN, this indicates that cell-to-cell communication is indispensable for endogenous encoding of photoperiodic information. These examples strongly imply that photoperiod coding has long-lasting effects, likely involving coupling reorganization within the cellular population that comprises the circadian clocks. Buijink et al. (2016) reported greater period instability in Per2::Luc rhythms of individual cells in the mouse SCN on long days versus short days, which further indicates a photoperiod-dependent change in coupling strength.

Coupling changes certainly entail oscillator amplitude changes (Oda et al., 2000), which in turn affect the clock phase-response curve to light (Lakin-Thomas et al., 1991). The ultimate mechanism behind photoperiod encoding could possibly lie on the seasonal modulation of the clock responses to light inputs (Evans et al., 2004; Pittendrigh et al., 1991). In this sense, the model proposed in the present work could be refined by splitting E and M into a distinct subpopulation of oscillators and manipulating individual coupling parameters in a complex neuronal network. This strategy was used by Taylor et al. (2017) to model different regions of the SCN, and they found a mild photoperiodic modulation in the phase difference between the regions. In their case, however, only one region received light inputs, which may explain the lack of pronounced photoperiodic effects.

Finally, our models do not consider the interaction of a homeostatic sleep drive and a circadian oscillator (Borbély et al., 2016; Deboer, 2018) in determining activity onsets and offsets under different photoperiods. Although prior activity is known to affect sleep propensity, several studies have indicated that the timing of sleep under different photoperiods is more strongly determined by the circadian pacemaker. In rats, sleep deprivation under different photoperiods is compensated for on the next day by different sleep compositions, with little change in sleep timing (Franken et al., 1995). In Borbély and Neuhaus (1978), the timing of sleep and wake onsets on the first day of abrupt photoperiod change was still shaped by the photoperiod of the previous day. Trachsel et al. (1986) showed that this also happens in rats released into constant darkness. We recognize that a sleep homeostat could potentially refine timing of activity. Thus, the model could be further complemented by the addition of a homeostatic oscillator. However, we would still consider a system composed of 2 coupled circadian oscillators to explain the complex dynamics of activity onsets and offsets in different light conditions.

Discussion

From an ecological viewpoint, it is tempting to interpret the duration of daily activity in different seasons as a direct response of organisms to the availability of favorable conditions for activity along the 24 h (Bennie et al., 2014). For instance, compression of daily activity in diurnal organisms during autumn/winter can be viewed simply as a reaction to fewer available light hours or to the concentration of warm temperatures around noon. However, laboratory studies with photoperiod manipulations have established that seasonal variation in α has an endogenous component, as shown by the establishment of nonlinear α and ψ variations when photoperiod is varied linearly (Fig. 1B, C; Elliott, 1976). Understanding the mechanisms of photoperiod dependence in α and ψ is still a challenge that requires mathematical modeling. It has been proposed that, in mammals, photoperiod variations of daily activity patterns involve (at least) 2 oscillators, E and M, which control activity onset and offset and are differentially entrained by dawn and dusk (Pittendrigh and Daan, 1976b).

In the mammalian circadian system, the SCN mediates photoperiodic adjustments in activity α; therefore, the nuclei have been strong candidates for the physical substrate of the E and M oscillators (Daan et al., 2001; Helfrich-Förster, 2009). The SCN functions as an overt oscillatory unit, expressing a 24-h profile in multiunit electrical activity (Houben et al., 2009). The duration of this high-electrical-activity state depends on the length of the photophase (reviewed in Coomans et al., 2015) and is mirrored in the daily α duration of overt activity rhythms (Houben et al., 2009) and melatonin profiles (Illnerová and Vaněček, 1982).

Our quantitative analysis of putative E-M models is an attempt to elucidate the requirements of the dual-oscillator system, which allows a response to different photoperiods. The general model predictions should be valid, irrespective of the precise tissue/cellular correspondence, as long as E and M behave as limit-cycle oscillators. In the current study, we searched for model configurations that presented the experimentally observed responses to (skeleton) photoperiods (Fig. 1B). After changes in the interval between pulses, a successful model should present (1) photoperiod-driven adjustments in α, with one oscillator tracking dawn and the other tracking dusk, and (2) a psi-jump after the interpulse interval reaches a threshold duration.

One emergent feature of all simulations is that when 2 light pulses are applied per day, the phase of the oscillator system tracks one of the pulses and remains this way upon systematic changes in the interval between pulses, until a certain threshold value. After this threshold, it tracks the other pulse, and the transition of these 2 regimes occurs abruptly, resulting in a psi-jump. This is true for a 1- or 2-oscillator system and independent of the configurations that were tested (Fig. 3E-H). The switch in the pulse association is typical of models with 2 zeitgeber inputs when the phase difference between them is systematically varied (Oda and Friesen, 2011), and it illustrates the contribution of mathematical modeling in coordinating our interpretation of complex dynamical phenomena. The 2-zeitgeber model configuration also predicted that bistability phenomena would occur when the 2 zeitgebers attained the vicinity of an antiphase relationship. This is in fact observed when skeleton photoperiods are systematically varied and the light pulses attain intervals around 12 h, when the “night interval” chosen by the nocturnal individual depends on the previous phase of activity (Pittendrigh and Daan, 1976a; Takahashi and Menaker, 1982).

While psi-jumps are easily replicated, photoperiod-driven adjustments in α appear to be more challenging. In the simplest configuration, model I, a single oscillator presented no considerable change in α. Model II, on the other hand, expressed photoperiod modulation of α. In this model, the coupling strength between oscillators determined how independently they moved under different photoperiods. In a loosely coupled system, the oscillators were almost independent from each other (Fig. 4A), passively following dawn and dusk. Therefore, model II was able to replicate photoperiod-induced changes in α, especially with low to mid coupling strengths (Fig. 4A, B). Interestingly, these α changes persisted when the simulated light pulses were turned off, mimicking photoperiod aftereffects on α (Suppl. Fig. S1). However, the aftereffects persisted for a short time, up to 8 days, depending on the initial value. The quick α decompression likely happened because our model did not include coupling changes in different photoperiods, which may enhance the aftereffects (Lucassen et al., 2012). But it is all the more striking to show that, even without these internal coupling changes, photoperiod aftereffects on α can be simulated in model II.

Despite the success of model II in replicating experimental results, to our knowledge, there has not been any evidence that specific SCN subpopulations receive dawn-only or dusk-only signals from its neuronal inputs. Signals from the retina reach the SCN through direct retinohypothalamic (Abrahamson and Moore, 2001) and indirect geniculohypothalamic (Jacob et al., 1999) pathways. In the former, light signals are conveyed by glutamatergic synapses to that VIP-containing neurons (Jones et al., 2018). In the indirect geniculohypothalamic pathway, neuropeptide Y (NPY) projections from the intergeniculate leaflet (IGL) input zeitgeber information to the SCN (Harrington et al., 1985). This pathway has been previously associated with photoperiodic responses. IGL-NPY signals have been implicated in entrainment of activity rhythms to fixed skeleton photoperiods (Edelstein and Amir, 1999), to complete photoperiods (Freeman et al., 2004; Kim and Harrington, 2008), and in the photoperiodic adjustments in the SCN (Menet et al., 2001). Moreover, NPY levels are crepuscular in the rat SCN (Shinohara et al., 1993), and the peptide has been shown to interfere with phase shifts in vivo (Weber and Rea, 1997) and in the SCN in vitro (Biello et al., 1997). Nevertheless, both afferent pathways input the SCN on its ventrolateral region. There remains, however, a possibility that dawn and dusk signals are equally received in the ventrolateral SCN, but there is a downstream switch that processes and transmits this information differently to the SCN network, depending on time of day.

In the dorsoventral plane of the mouse SCN, a photoperiod-dependent phase difference was recorded in PER2::LUC rhythms between the ventrolateral and dorsomedial regions (Evans et al., 2013). Photoperiod also modulates the spatiotemporal organization in clock gene rhythms in the anteroposterior plane (Inagaki et al., 2007; Yoshikawa et al., 2017). However, we find it hard to relate E-M to a binary subdivision of the SCN, since photoperiod affects the rhythms not only between the different subregions in a plane but also within each subregion. As reviewed in Coomas et al. (2015), the duration of nocturnal c-fos expression is longer in long nights versus short nights, within both the ventral and the dorsal SCN. Also, Buijink et al. (2016) recently showed a photoperiod-dependent change in the phase distribution within the anterior SCN.

Considering the difficulty in relating model II to known physiological mechanisms, a more plausible configuration was tried in model III, in which both E and M are sensitive to both dawn and dusk. However, in this configuration the 2 oscillators track the same pulse, and α does not change at all under different photoperiods, regardless of oscillator asymmetries. It essentially resembles the behavior of a single oscillator in model I. Model IV was an attempt to introduce a stereotyped asymmetry in the E and M responses by inverting the sign of the light pulses to the M oscillator. This resulted in an inversion of the shape of its phase-response curve (Johnson, 1999) relative to the E oscillator. Although some changes in α were observed, the resulting patterns do not match experimental ones.

In sum, model II best replicates experimental results, but it requires 2 SCN subpopulations, each exclusively responsive to either dawn or dusk. While it still seems challenging to assign any of these features to the SCN and their light inputs, this model has a potential physiological counterpart at another step of photoperiod transduction beyond the SCN. As mentioned before, the SCN are the first step in a series of complex neuroendocrine pathways that regulate seasonal physiology, specifically seasonal reproduction (Nishiwaki-Ohkawa and Yoshimura, 2016; Wood and Loudon, 2014). In short, photoperiod-shaped neural outputs from the SCN control the nocturnal release of the pineal hormone melatonin, such that the duration of the melatonin signal matches the length of the night (Bartness et al., 1993). Melatonin information is, in turn, decoded in the pars tuberalis (PT) of the pituitary. Studies in sheep and hamsters suggest that, in the PT, the processing of melatonin duration depends on 2 events, the daily times of melatonin rise and fall (Lincoln et al., 2003). Increasing melatonin levels at dusk control the phase of the clock gene Cryptochrome1 (Cry1), whereas the decreasing melatonin levels at dawn control the phase of the clock gene Period1 (Per1; Lincoln et al., 2002; Messager et al., 2000; West et al., 2013). It has been proposed that the phase relationship between the 2 clock genes transduces the photoperiodic information to downstream molecular pathways, which ultimately regulate reproductive physiology (Lincoln, 2006; Lincoln et al., 2003). This mechanism is analogous to the “internal coincidence” model postulated long ago for photoperiod decoding (Pittendrigh, 1972). The molecular events in the PT are compatible with our E-M model II, in which each oscillator is responsive to either dawn or dusk alone. Melatonin rise at dusk and melatonin fall at dawn are sensed by different oscillating molecules (Cry1 and Per1), which are likely mutually regulated (coupled) in a molecular network (Lincoln, 2006; Lincoln et al., 2003). Alternatively, however, other mechanisms have also been proposed to explain the processing of the melatonin signal in the PT (Dardente and Cermakian, 2007; Masumoto et al., 2010).

In summary, in mammals there are 2 physiological steps in the decoding of photoperiodic information, which depend on the phase relationship between dawn and dusk events. The different working versions of the mathematical model could represent the decoding of daily light/dark proportion within the SCN or the transduction of the nocturnal melatonin signal duration within the PT. We found it challenging to quantitatively simulate a plausible E-M model that could match known physiological properties in the SCN and, at the same time, respond to photoperiod by promoting adjustments in α and psi-jumping at critical interpulse intervals. The models at which we arrived are overall approximations, each with its pros and cons. We hope this first approach can provide coordinating steps toward other ongoing modeling efforts.

Supplemental Material

Supplemental Material

Supplemental Material

Supplemental Material

Supplemental Material

Supplemental Material

Supplemental Material

Supplemental Material

Supplemental Material