Seasonality and light phase-resetting in the mammalian circadian rhythm

We study the impact of light on the mammalian circadian system using the theory of phase response curves. Using a recently developed ansatz we derive a low-dimensional macroscopic model for the core circadian clock in mammals. Significantly, the variables and parameters in our model have physiological interpretations and may be compared with experimental results. We focus on the effect of four key factors which help shape the mammalian phase response to light: heterogeneity in the population of oscillators, the structure of the typical light phase response curve, the fraction of oscillators which receive direct light input and changes in the coupling strengths associated with seasonal day-lengths. We find these factors can explain several experimental results and provide insight into the processing of light information in the mammalian circadian system. In particular, we find that the sensitivity of the circadian system to light may be modulated by changes in the relative coupling forces between the light sensing and non-sensing populations. Finally, we show how seasonal day-length, after-effects to light entrainment and seasonal variations in light sensitivity in the mammalian circadian clock are interrelated.

Daily or circadian cycles in behavior and metabolism can be observed for virtually all forms of life. The utility of circadian rhythms relies on the proper timing of these cycles relative to external environmental oscillations. Thus, a defining property of circadian rhythms is their ability to be entrained to external time cues or zeitgebers. The principal zeitgeber for the mammalian circadian clock is light 1 . Therefore, a crucial component to understanding mammalian circadian rhythms is an improved understanding of the impact of light on the circadian cycle. A first step in this endeavor is understanding the response of the circadian circuit to a brief light pulse.
The theory of phase response curves (PRC) provides a natural language for studying the effects of external stimuli on endogenous rhythms with a rich history of application to circadian biology 2,3 . Phase response curves characterize the phase shift induced by the application of the stimulus at different phases of the oscillation. For instance, the amplitude of the PRC gives the entrainment range of the system to a weak resetting signal and the zeros of the PRC determine the entrainment angle 4 .
Phase response curve theory figured prominently in early investigations of circadian rhythms, where organisms were exposed to a sensory stimulus at a sampling of points across the daily cycle and phase shifts were measured relative to some behavioral or physiological marker 3,5 . However, the discovery of the location of the master circadian clock in a small region of the hypothalamus known as the suprachaismatic nucleus (SCN) provided a neurological basis for circadian rhythms in mammals 6,7 . The SCN was found to contain thousands of coupled clock neurons which each contain a biochemical oscillator with a period of approximately 24 h 8 . Daily activity cycles are driven by this large ensemble of coupled oscillators acting collectively to produce a reliable circadian oscillation.
Thus, a light stimulus applied to the mammalian circadian rhythm does not act by shifting a single limit cycle oscillator, but rather acts by shifting the oscillations of individual clock neurons which in turn induce a shift in the collective rhythms produced by the ensemble. The recognition of this distinction helped motivate the development of the theory of collective phase response curves which describe the collective phase shift of an entire population of coupled oscillators subjected to a stimulus [9][10][11][12] . Collective phase resetting is especially important for the mammalian circadian response to light, because only a fraction of the clock cells are phase shifted in response to the stimulus 13 . This may induce non-trivial transient dynamics on the system following a light perturbation and forms a major focus of this work 14 .
In general, the phase-shifting behavior of a coupled ensemble of oscillators differs from the behavior of a single autonomous oscillator. In this work we study the transformation between the response of a single circadian cell to a light-pulse (microscopic PRC) and the collective phase response described by the shift in the mean-phase of the population of oscillators. A growing literature on collective phase resetting has revealed that coupling, oscillator heterogeneity and network structure can all lead to significant differences between the microscopic and collective PRCs in networks of coupled oscillators [9][10][11][12] . Within the circadian literature, examinations of collective phase resetting have led to the formation of a rule of thumb that increasing phase dispersion in the oscillator population leads to a monotonic increase in the amplitude of the collective phase response 15,16 and thus the entrainment range of the collective oscillator 3 .
In addition to the role of light input in ensuring the circadian clock is synchronized to the outside environment, the SCN is also responsible for storing seasonal day-length information [17][18][19] . The ability of seasonal day-lengths to alter the core circadian clock was established in early circadian studies, where it was noticed that entrainment of mammals to long/short day-lengths caused lasting changes in the endogenous circadian period when organisms were transferred to a dark environment 20 . These effects are known as seasonal after-effects and have been described for many mammal species 5,20 .
Recently, significant progress has been made in characterizing the physiological changes in the SCN underlying seasonal day-lengths changes [21][22][23][24][25] . The physiological changes in the SCN which encode the seasons have been shown to affect the phase response to brief light pulses. When organisms are entrained to long (summer) days the phase shift caused by a brief light pulse is seen to decrease 21,26 . This seemingly contradicts the rule of thumb for collective phase resetting, because experimental evidence has also shown that phase dispersion in the SCN increases in longer day-lengths 21 . A primary goal of this work is to provide a unified theory of collective phase resetting to light in mammals, consistent with seasonal changes in SCN physiology.
In order to study phase resetting to light we make use of the recently developed m 2 ansatz to derive a threedimensional model for the core circadian clock 27 . Significantly, the m 2 ansatz is supported by experimental evidence and the resulting model gives variables and parameters which may be interpreted physiologically in the SCN 27,28 . Our work utilizes the framework for studying collective phase resetting developed in Levnajic and Pikovsky 12 , and we apply this theory specifically to light induced phase resetting in mammals. We extend the results in Hannay et al. 29 to consider oscillator networks where only a fraction of the population receives the phase-shifting stimulus. Additionally, we develop a perturbation technique which can be applied generally to characterize the effects of coupling on collective phase resetting.
The principal biological impacts of our work are threefold. First, we provide a general theory for the effect of the collective amplitude on the phase-shifting capacity of the circadian clock. Our analysis reveals that the rule of thumb that lower amplitude rhythms with larger phase dispersion give larger phase shifts in response to a stimulus is incomplete, and more detailed analysis is required for many real-world phase response curves. Secondly, our analysis reveals the reduction in light-shifting capacity observed for organisms entrained to long daylengths may be explained by an adjustment of the coupling strengths between clock neurons as a function of the seasonal day-length. Finally, we find this adjustment of coupling strengths is consistent with current theories for seasonal day-length encoding and is required to explain seasonal after-effects to light entrainment in mammals.

Results
Formulation of the model. Circadian model. The suprachaismatic nucleus (SCN) is a collection of about 10,000 neurons which form the core circadian pacemaker in mammals. Individual neurons in the SCN contain a biochemical feedback loop which cycles with a period of approximately 24 h. While the SCN produces a variety of spatiotemporal activity patterns it can be functionally and physiologically broken into the ventral (core) and dorsal (shell) populations 30 .
In mammals, light input comes in through the eyes and is channeled to the SCN along the retinohypothamic tract (RHT) 13 . It is remarkable that the core circadian clock receives light information through a direct pathway from the eyes, which underscores the importance of light in entraining mammalian circadian rhythms. However, only a fraction of the clock cells in the SCN receive light input directly, with the majority of cells receiving input in the ventral region 13 . Therefore, for the purposes of our model we split the SCN into ventral and dorsal phase clusters and allow light-input into only the ventral population (Fig. 1).
Coupling between clock neurons in the SCN is mediated by a large suite of neurotransmitters 31 . In this work we focus on the functional coupling between the regions, although it may assist the reader to give an  www.nature.com/scientificreports/ interpretation of the coupling in reference to two predominant neurotransmitters in mammalian circadian rhythms: vasoactive intestinal polypeptide (VIP) and γ-aminobutyric acid (GABA) whose properties have been characterized experimentally 23,32,33 . Perhaps, the best understood coupling agent in the SCN is VIP which is released by the ventral population and is received by all the oscillators 32 . VIP is known to be a synchronizing force in the SCN (phase attractive) 33 . Recent experimental results 23,34 and detailed mathematical modeling 35 suggest that GABA mediated coupling is more subtle. GABA is released and received by all or nearly all clock neurons in the SCN and has been identified as both a synchronizing 36 and desynchronizing 37,38 agent among clock cells, although recent evidence suggests these properties vary spatially in the SCN 23,35 . Evidence suggests that both VIP and GABA are involved in the communication of phase shifts between the ventral and dorsal SCN as well as the storage of seasonal day-length information in the SCN 23,35,39,40 . Here we assume the combined action of VIP and GABA act to modulate the strength of the coupling between the ventral and dorsal phase clusters in the SCN. This conceptual model, summarized in Fig. 1, may be translated into a coupled phase oscillator system, where Q(φ) gives the microscopic phase response curve of the ventral oscillators to light and η v,d k defines a white noise process i.e. ( �η k (t)� = 0 and �η k (t)η l (t ′ )� = 2δ kl δ(t − t ′ ) ). The ε factor scales the microscopic phase response curve and is taken to be a small parameter. The coupling strengths are given as K from,to and we let M v,d indicate the total number of oscillators which fall into the ventral and dorsal phase clusters. Finally, we define q = M v /(M v + M d ) to be the fraction of ventral (sensing) oscillators in the population and p = 1 − q.
We allow the oscillators within the ventral and dorsal regions to be heterogeneous in their intrinsic frequencies and assume each cluster has a Cauchy (Lorentzian) distribution of frequencies, with the mean frequency ω v,d 0 and dispersion parameter γ. In addition to daily timekeeping, the SCN is also responsible for storing seasonal day-length information 17,18,41 . It has been shown by several experimental groups that the phase difference between the dorsal and ventral clusters grows with the seasonal day-length, making this a leading hypothesis for how seasonal information is encoded in the SCN 23,24 . Additionally, it has been suggested that the physiological root of this seasonal variation in the phase difference is alterations in the coupling forces in the SCN 23,25,35 . Thus, we incorporate seasonal effects (daylengths) into our model by allowing the coupling strengths K vd and K dv to vary with the seasonal day-length, as these coupling terms can control the phase difference between the ventral and dorsal populations.
Macroscopic model. The model for the mammalian SCN as given in Eq. 1 gives a high-dimensional representation of the dynamical state of the circadian rhythm as M v + M d = O (10 4 ) . This high dimensional representation of the system makes analytical analysis of the light-response difficult. Therefore, we make use of the recently developed m 2 ansatz to derive a low-dimensional macroscopic model for the ventral and dorsal phase clusters 27 . Mathematically, we have demonstrated that the m 2 ansatz provides an accurate approximation to noisy heterogeneous systems as seen in Eq. 1 27 . In the absence of noise (D = 0) and a Cauchy distribution of frequencies g(ω) the Ott-Antonsen ansatz provides an exact reduction 42 . Crucially, for this application the m 2 ansatz may be justified through comparison with experimental data on the phase distribution of cellular oscillators in the mammalian SCN 27,43 .
First, we define the Daido order parameters 44,45 for the ventral and dorsal phase clusters as: where n ∈ Z . The special case of n = 1 gives the classical Kuramoto order parameter Z v,d where R is known as the phase coherence and gives a measure of the overall synchrony in the population. When R 1 = 1 the population is phase locked in perfect synchrony and R 1 = 0 if the population is completely desynchronized and spread uniformly around the unit circle. Additionally, ψ 1 gives the mean phase of the population. Since our macroscopic model will be written in terms of the Kuramoto order parameters we will drop the subscript for those terms i.e. Z v,d Using these order parameter definitions we may rewrite Eq. 1 as, where the bar indicates the complex conjugate. Inserting the Fourier series representation into the continuity equation gives a system for the Fourier coefficients with the bar representing the complex conjugate. In the continuum limit the Daido order parameters Z v,d n are given by, using that all the terms in the Fourier series integrate to zero except the n = k term. Further, since we approximate the natural frequency distribution as a Cauchy/Lorentzian distribution (Eq. 2) we may evaluate the integral (Eq. 8) under the assumption that A k (ω, t) may be analytically continued into the complex ω plane 42 . Thus we have that, This substitution into Eq. 7 gives, Finally, we consider the system with k = 1 and apply the m 2 ansatz ( 27 . Applying the m 2 ansatz and separation into the real and imaginary parts of the expressions gives a four dimension system describing the phase coherence R v,d and mean phase ψ v,d of each cluster. However, with a change of variables θ = ψ d − ψ v ("phase gap") and letting �ω = ω d 0 − ω d 0 and γ =γ + D we arrive at a three dimensional system of equations: By setting ω = qω v 0 + pω d 0 , we can define � = qψ v + pψ d as the collective frequency of the system in a synchronous state, www.nature.com/scientificreports/ In order to facilitate our analysis we define a default parameter set for this model given in Table 1. Under this parameter set Eq. 11 evolve to a fixed point (R * v , R * d , θ * ) with a collective frequency * -we will use starred quantities refer to fixed-point solutions. Finally, we note that experimental evidence has shown that the phase gap between the dorsal and ventral populations is typically a small variable θ * ∈ [0, 0.5] radians for photoperiods in the range of 6-18 h of light 23 , although it has been shown in mice to grow considerably when they are kept in 20 h or more of light each day 24 .
Collective phase response curves. Components of the collective phase response curve. The collective phase response to a stimulus may be defined by the shift in the mean-phase ψ = Arg(Z) induced by the light perturbation. For a brief (Dirac δ function) stimulus we may break the collective phase shift into two components 12 : 1. The prompt phase shift ( 0 ) induced at t =t the instant the stimulus is applied. 2. The relaxation phase shift ( R ) which results from any phase shifts induced as the system relaxes to its asymptotic state.
The collective phase shift ( ∞ ) is then given by, where we define Z 0 as the order parameter just prior to the perturbation and Z A as the order parameter just after the perturbation. A subscript of A on quantities will refer to the quantity just after the perturbation is applied. It is also useful to define the amplitude resetting curve as a measure of the perturbations transient effect on the amplitude of the collective rhythm, Given the assumed stability of the limit cycle, perturbations of the amplitude R are expected to decay, thus the amplitude resetting is defined in terms of the initial amplitude reduction imposed on the system.
Single population case. We first consider the collective phase response for a single population of oscillators, that is, we consider the case where all oscillators in the population receive the light stimulus. These results will aid our consideration of the two population case, as they can be used to describe the phase shift in the ventral oscillator population. For a general phase response curve Q(φ) we have previously derived an asymptotic formula for the collective phase response of a single population 29 of Kuramoto-Sakaguchi oscillators making use of the Ott-Antonsen formalism. In this section we adapt those results to study phase shifts in a population which follows the m 2 ansatz as has been found in experimental measurements of the SCN phase distribution 27 . For times close to the perturbation t ≈t we may approximate the single population continuity equation as, From these expressions we can see the principal effect of the phase distribution on the shape of the PRC is to re-weight the Fourier harmonics according to f k (R) . As R → 1 we have that f k → 1 and the collective and microscopic prompt phase response curves coincide. However, when R < 1 the first harmonic of the microscopic phase response curve is amplified like R 3 + 1 R while the higher harmonics are damped out by higher powers of R. The imaginary part of Q (ψ) can be expressed as, In this case we see that the modulation term g k (R) goes to zero as R → 1 meaning the amplitude is unaffected by the stimulus in this limit. Additionally we observe that, A n e inφ +Ā n e −inφ = A 0 2 + ∞ n=1 a n sin(nφ) + b n cos(nφ). www.nature.com/scientificreports/ so we expect the amplitude shifts to be greatest around the zeros of the microscopic phase response curve, with transient increases in R around stable points and decreases around unstable zeros.
Application to light PRCs. The mammalian phase response curve to light has been characterized over a large variety of species and conditions 2,46 . Generally, the circadian rhythm shows small sensitivity to light during the subjective day, phase delays during the early subjective night and phase advances in the late subjective night 47 .
Although we expect our results to hold more generally we will focus our attention on phase response curves with this general shape. The single population results derived in the last section indicate that for microscopic PRC's dominated by their first harmonic, as are commonly assumed in the circadian literature 15,16 , a general amplification in the phase response is expected as the oscillators are more dispersed in phase or have a reduced amplitude (either through weaker coupling or greater frequency heterogeneity in the population). This expectation has surfaced in the circadian literature under a variety of guises in the context of both phenomenological and biochemically motivated models 16 .
These observations have led to the formation of a rule of thumb in the circadian literature that the phase response of a limit cycle oscillator should increase when the collective amplitude/phase coherence decreases 16 . However, our analysis shows this prediction will only hold when the underlying individual/microscopic phase response curve is dominated by its first harmonic-otherwise a overall reduction in the amplitude may be observed 29 (Eq. 22). An example of this effect is shown in Fig. 2 for a first harmonic phase response curve and a light-like PRC shape. The first harmonic curve shows a uniform increase in amplitude as the phase coherence of the sensing population decreases, whereas the light-like PRC shows a initial decrease in amplitude as the higher harmonics are dissipated.
The general shape of the mammalian phase response curve to light has significant power at higher harmonics 46 . Therefore, the rule of thumb does not necessarily apply in this case. For example, the phase coherence of the ventral (sensor) population is known to decrease with increasing day-length 48 which leads to the expectation of increasing amplitude in the response to light. However, as previously noted the opposite trend has been observed experimentally where organisms entrained to longer day-lengths show decreased sensitivity to light-pulses 21 .
This reduction in amplitude of the collective phase response, in spite of a reduction in the phase coherence of the sensing population, may be at least partially explained when the higher harmonics in the microscopic phase response to light are taken into consideration. This effect has been noted in the course of simulations 21 and is readily explained by the theory given here using the m 2 ansatz and detailed previously for cases adhering to the Ott-Antonsen ansatz 29 .
In the following section we investigate the effects of only having a fraction of the total population shift in response to a light pulse. This analysis reveals an additional effect which allows the amplitude of the collective Two population case. We now consider the collective phase response of the circadian model Eq. 11 to brief light pulses using that ∞ = 0 + R . In the first section we compute the prompt phase shifting behavior 0 for the circadian model with a subset of sensor cells and observe the effects on the initial phase shifting behavior ( 0 ). In the next subsection, we present a perturbation technique to determine the relaxation phase shift ( R ).
In the figures for this section we consider a microscopic phase response curve Q(ψ) which is fit to experimental measurements of the human phase response curve to brief light pulses 46 . Prompt resetting 0 . We begin by studying the prompt phase shifting curve 0 for the circadian model. By applying Eq. 20 and using that the dorsal population is unaffected by the perturbation, we find the order parameter Z = qZ v + pZ d just after the perturbation Z A = qZ v (1 + iεQ(ψ v )) + pZ d . Therefore the prompt phase shift 0 for the SCN model can be derived as: where η = p/q is the ratio of dorsal to ventral (sensors) in the population. We may now expand Eq. 24d, Which gives an analytical expression for the prompt resetting in our system using our expressions for v 0 and v for a single population of oscillators (Eq. 21). We note that Eq. 25 has the expected limits: As η → 0 , 0 → v 0 and the system converges to the behavior of a single population of oscillators, in addition as the dorsal population grows η → ∞ we see that 0 → 0 causing the system to become unresponsive to perturbations. This analytical approximation gives an accurate approximation for the prompt resetting curve when compared with numerical simulations (Fig. 3).
It is interesting to note the differences between our system and a single population of oscillators which all shift in response to the stimulus. In the two population system the damping of higher harmonics in the microscopic www.nature.com/scientificreports/ PRC is also observed, however we additionally see a decrease in the initial shift as a function of the fraction of oscillators which receive the light pulse. We note that under the assumption that θ ≈ 0 we may approximate Eq. 25a as, giving the intuitive result that the overall amplitude of the initial response scales with the fraction of ventral sensor cells in the population q.
In addition, unlike the single population case we see the prompt phase response curve depends on the amplitude response function v . This dependence leads to a slight change in the zeros (entrainment points) of the prompt PRC when compared to the microscopic PRC since v will be largest about the zeros of Q. Moreover, this effect is dependent on having a non-zero phase gap between the two populations ( θ = 0).
Relaxation phase shift R . We now consider the relaxation shift R which describes the phase shift induced during the return of the system to equilibrium following a perturbation. For the single population case this relaxation phase shift was directly computable from the Ott-Antonsen equations describing the collective dynamics 12 . This computation relied on the relatively simple spiral isochrons of the Kuramoto-Sakagucki coupling scheme 4,12,29 .
For the circadian model presented here this quantity can no longer be easily computed by direct integration. However, the relaxation phase shift occurs during the transient decay of the system back to the dynamical fixed point (R * v , R * d , θ * ) of the macroscopic model. The collective frequency of the unperturbed system is given by . Therefore, the phase shift induced as the system relaxes back to the equilibrium state is given as, Thus, we may calculate the relaxation phase shift by integrating the frequency mismatch between the perturbed system and the steady state system along the trajectory of the system as it returns to equilibrium. The relaxation trajectory may be approximated by a perturbation about the fixed point under the assumption the light-pulse does not induce a large deviation from (R * v , R * n , θ * ) . We set, where σ is a small parameter and with initial conditions R v (0) = R * v + �R v , R d (0) = R * d , θ(0) = θ * + �θ using that the dorsal population is initially unaffected by the light stimulus. The initial changes in R v and θ can be written in terms of the prompt phase and amplitude response curves for the ventral population: is the ventral phase just after the perturbation. Therefore, the leading order terms in σ for the relaxation phase shift is given by, with the (A,B) constants determined by the model parameters. In practice we find the leading order term in σ is sufficient to provide a good approximation to the numerical solutions (see Fig. 4) although higher order terms may be taken in the perturbation series (Eq. 28) if additional accuracy is required. We note that (A, B) are a measure of the sensitivity of the collective frequency of the system to perturbations in the amplitude ( R v ) and phase gap θ respectively, and are weighted by the stiffness of the system to perturbations in those directions.
To gain intuition of how the circadian model parameters will affect the relaxation phase shifts, we solve for the relaxation terms analytically for a simplified system (see supplementary information). If the amplitude of the ventral and dorsal populations are fixed (without loss of generality let R v,d = 1 ) we find that, From this simplification we can see that B ∈ [−p, q] when both coupling terms are positive. Moreover the change occurs at α = q p from a positive to negative value. Note, for symmetric coupling ( α = 1 ) between the regions and q = p = 1 2 we have B = 0 . In Fig. 5 we show the variation of B with α using both the perturbation approach and the simplified formula Eq. 30.
Collective phase response curve. When the prompt phase shift 0 is combined with the relaxation shift R we find the collective phase response curve ∞ , Scientific Reports | (2020) 10:19506 | https://doi.org/10.1038/s41598-020-74002-2 www.nature.com/scientificreports/ with the constants (A, B, C, D) as defined in the previous sections. In general, we find this approximate formula provides a good approximation to the numerically determined collective phase response (Fig. 6). To improve our understanding of the role of light-input to the mammalian circadian system we analyze these results further.
Of particular importance is the amplitude of the collective phase response as this determines the entrainment range for weakly forced systems 3 .
(31)  www.nature.com/scientificreports/ We first note that the collective PRC for the circadian system carries over many of the trends of the single population model. Namely, we expect that higher harmonics in the microscopic phase response curve will be damped with the first harmonic amplified like R 3 v + 1 R v . This transformation in the shape leads to an overall smoothing effect and is tied to the disorder in the underlying population.
Additionally, we note that the term proportional to the amplitude response curve (D − A)(1 − � v ) is expected to be of comparatively small magnitude as it is proportional to 1 Moreover the amplitude response curve reaches its maximum values around the zeros of the individual phase response curve. Thus, the amplitude of the collective phase response curve is largely determined by the (C − B)� v 0 term, while a shift in the entrainment points is determined by the From Eq. 31 we see the amplitude of the collective PRC is influenced by the sign of the B constant as determined by the relaxation dynamics. Positive values of B indicate the relaxation shift acts to decrease the initial phase shift thus providing a resistance to the phase shift. Negative values of B indicate the initial shift is reinforced/increased during the transient relaxation, while a value near zero indicates the relaxation shift has a small effect on the collective phase response.
The simplified expression in Eq. 30 allows for intuition on the scale and sign of B. In particular we can see the value of B scales with the ratio of the feedforward coupling strength K vd to the feedback strength K dv . For a balanced system K vd = K dv , q = 0.5 we observe that B = 0 , although by varying this ratio of coupling strengths the system can toggle between a resistant/reinforcing behavior to the initial phase shifts.
For a pure feedforward network, where K dv = 0 , the B term is positive and the effect of having a fractional sensor population on the amplitude of the collective phase response disappears. In this limit the phase shift in the ventral population is imposed on the dorsal population over time. Thus, when K dv = 0 the non-sensing dorsal population can act as a feedback on the phase shifts and integrate the current phase shift against the past history.
Therefore, we see that the entrainment range of a two-population system is crucially dependent on the ratio of the coupling strengths between the sensing (ventral) and non-sensing (dorsal) populations (Fig. 7). By adjusting the ratio of these coupling strengths the size of the light response may be modulated.

Seasonal effects on light resetting.
Our study of the collective phase response in a two-population model identified the relative strengths of the intra-population coupling forces as an important factor in determining the amplitude of the phase response. Experimental evidence has found that the amplitude of phase shifts induced by a light-stimulus decreases in mammals entrained to long day-lengths 21,26 . This could be explained in our framework by an increase in the network resistance to the phase shift (B), in animals exposed to long day-lengths.
To evaluate this hypothesis we may check for consistency against two other seasonal light effects on mammalian circadian rhythms: Seasonal encoding and light entrainment after-effects. First, we consider seasonal encoding. Experimental evidence has indicated that the phase difference between the dorsal and ventral phase clusters grows with the day-length 23,24 . Within our model this corresponds to the variable |θ * | growing with the entrained day-length. We make use of the hypothesis that these changes in θ occur through an adjustment of the coupling strengths rather than a change in the intrinsic periods in the ventral/dorsal SCN 23,35 . Thus, we consider the intrinsic periods to be constant while allowing the coupling strengths to vary.
Let χ ∈ [0, 1] be the photoperiod, or the fraction of the circadian day in which the organism is exposed to light. Considering Eq. 11 we see the steady state phase gap θ * is given by, www.nature.com/scientificreports/ From this we can see that we must have G ′ (χ) < 0 to allow the absolute value of the phase gap to increase with the photoperiod length. Additionally, we observe that θ * will only show significant variation for a small range G(χ) = O(�ω) and will asymptote to a small value for G(χ) outside this range. This nonlinear dependence of the phase gap θ * (χ) with asymptotic values for short photoperiods has been observed in experiments using both Per2 and Bmal circadian phase markers 23,24 .
Moreover, the macroscopic model also identifies a fundamental trade-off as the photoperiod is lengthened. As G(χ) decreases towards �ω the phase gap increases quickly. However, a further increase in χ will cause the system to undergo a bifurcation where the ventral and dorsal regions decouple from one another. Thus, we predict that organisms which show robust seasonal adjustment necessarily must approach a bifurcation to desynchrony/ large phase gaps at long photoperiods. Therefore, organisms which show robust decreases in G(χ) , thereby showing larger changes in the phase gap variable θ * as the photoperiod lengthens, will also approach a bifurcation to desynchrony more closely and display rhythm abnormalities.
In fact, several species of mammals show desynchrony or large phase gaps when perturbed outside their normal photoperiodic range by unnatural lighting conditions (e.g. constant light) 24,49,50 . Moreover, it has been observed that hamsters which show robust seasonal adjustments to short day lengths have a higher propensity for rhythm abnormalities under constant lighting conditions-illustrating the trade-off identified by our analysis 51 .
In order to relate these properties to reduced light phase-resetting responses at long photoperiods we now recall a predominant circadian after-effect to light-entrainment in mammals 1,23 : Mammals entrained under short day-lengths show a transient increase in period when moved to a dark environment. The long day length after-effect works in an opposite direction by inducing a short period with the magnitude of the period change increasing with the entrained day length. In the context of our model this implies that d� dχ > 0 in Eq. 12. Expanding this condition gives, Applying the assumption that θ * is a small variable we may simplify this condition to give, and we consider the case that �ω > 0 as indicated by experimental evidence 22,23,52 . Under the approximation that R v ≈ R d we have that 2H(χ) ≈ qK dv (χ) − pK vd (χ) and G(χ) ≈ K dv (χ) + K vd (χ) . This simplification allows us to express our approximation for the network resistance to phase shifts B (Eq. 30), in terms of H and G, Thus, taking the derivative with respect to the photoperiod χ, we see that the network resistance to phase shifts B will increase with the photoperiod directly from the aftereffect condition (Eq. 34). This gives the surprising result that seasonal entrainment after-effects and the reduced sensitivity to light-pulses at long photoperiods are intimately related to one another. In fact, the presence of one implies the other in our model. Additionally, we see that the seasonal adjustment condition G ′ (χ) < 0 is consistent with the after-effect and increasing phase shift resistance condition (Eq. 36). Therefore, we find the adjustment of the coupling strengths required to explain three predominant light mediated circadian effects are all mutually consistent within our model.

Discussion
In this work we focus on phase resetting to light in mammalian rhythms making use of the m 2 ansatz to derive a simplified model of the central clock. The reduced model holds the advantage that the collective variables (R v , R d , θ) all have physiological interpretations and are measurable in experimental treatments. We have focused on the effects of heterogeneity in oscillator frequencies, the shape of the microscopic phase response curve to light, the effects of only a fraction of populations receiving direct light input and variation of the coupling strengths between regions on the phase resetting response.
Similar to previous work on this subject we find heterogeneity of the population changes the shape, amplitude and zeros of the collective phase response curve 12,29 . Moreover, we note these alterations occur through a re-weighting of the Fourier components of the collective phase response with the first harmonic amplified and higher harmonics being damped out as the oscillators spread out in phase. We find examination of only the first harmonic terms may give misleading results when considering phase resetting to light in mammals.
The effect of having a fraction of the population receive light-input can lead to a reduction in the overall amplitude of the collective phase response. However, we find this effect is dependent on the coupling between the oscillators in the SCN. In a feedforward network, where the sensing ventral cells project more strongly on www.nature.com/scientificreports/ the non-sensing dorsal cells than the feedback connection, the initial shift induced on the sensing population is largely imposed on the total population over time. Thus, for pure feedforward network architectures the effect of a fractional sensing population is to induce a time delay on the shift of the population mean phase. However, when feedback coupling from the dorsal cells to the sensing ventral populations is significant we see the reduction in the initial shift caused by the fractional sensing is retained and even reinforced over time. This leads to the conclusion that the relative coupling strengths between the ventral and dorsal oscillators allows the system to weight phase shifts induced by light differently. The re-weighting of the coupling strengths between the sub-populations with seasonal changes in day length may act like the aperture on a camera by allowing the sensitivity of the clock to light to vary with the total amount of light input received. In short days the feed-forward connection is weighted more strongly to allow for larger responses when light input is more scarce. In contrast, relatively weighting the feedback connection more strongly in long days results in a reduction in light sensitivity when more light is received over the course of the day.
These results may also have applications to the study of aging of the circadian clock. In a similar manner to long-day lengths, older animals tend to show reduced phase coherence and the clock neurons are thought to be more weakly coupled as the animals age [53][54][55] . Moreover, these aged animals also show reduced phase-shifts in response to light stimuli and slower entrainment to shifted light schedules [56][57][58][59] . Our analysis reveals this reduction in phase shifting capacity may be explained in terms of a reduction in the ratio of the feed-forward to feedback coupling between the ventral and dorsal populations.
Our results are consistent with current ideas of how seasonal information is encoded the mammalian circadian clock 23,60 and provide an explanation for mammals showing reduced phase shifts to light when entrained to long day lengths 21,26 . Furthermore, our model reveals an intimate connection between seasonal day-length encoding, seasonal entrainment after-effects and the amplitude of the phase-response to light. Additionally, we find the change in coupling strengths with the day-length required to explain each of these phenomena are mutually consistent within our model.
Recently, it was found that two feedback loops which make up the cellular circadian clock show a differing light response 61 . Mathematical models of these two feedback loops as weakly coupled phase oscillators bear a strong resemblance to the network level models we consider in this work 62 . This work on the cellular scale also identified the inter-oscillator coupling as an important parameter in the dynamics 62 . Moreover, the interaction between the two biochemical feedback loops appears to vary across tissues, allowing for differential responses to be encoded across the body 63 . This possible convergence on a common mechanism across scales suggests these findings could be a unifying principle in circadian biology.
The model reduction technique ( m 2 ansatz) was especially chosen because of its agreement with experimental results for mammalian circadian rhythms 27,43 . However, we note that all of the results presented in this work could be expressed using an alternative moment closure-for instance the Ott-Antonsen technique 27,42 . The choice of moment closure enters our calculations at only two points: to close the Daido order parameter sequence in the macroscopic model (Eq. 10) and again in calculating the prompt resetting curve (Eq. 19). Finally, we note that changing between the Ott-Antonsen and m 2 ansatz will not qualitatively change the conclusions of this work.
These results remain to be strengthened both from a biological and theoretical standpoint. In order to derive these results we have assumed an all-to-all connectivity between clock cells in the SCN and simple sinusoidal coupling between the oscillators. However, the connectivity in the SCN is known to be much more complex 43 . An important extension of these results would be to consider phase resetting in a general circadian network building on previous results 64, 65 .
From a biological standpoint it remains to be tested whether an increasing resistance to phase shifts under long-day lengths underlies the decreased sensitivity to light pulses for organisms entrained to long day-lengths. Although, this prediction seems testable by measuring initial phase shifts to light and comparing these with the asymptotic phase shifts obtained over long-times. This would be particularly interesting if the effect was seen to vary with the entrained day-length and could provide evidence for a variation in the relative coupling strengths between the ventral and dorsal SCN with seasonal day-length. Recent experimental evidence suggests that the coupling strength, as determined globally in the SCN, decreases with increasing day-lengths which could provide indirect evidence for this hypothesis 25 .