Modeling Two-Oscillator Circadian Systems Entrained by Two Environmental Cycles

Several experimental studies have altered the phase relationship between photic and non-photic environmental, 24 h cycles (zeitgebers) in order to assess their role in the synchronization of circadian rhythms. To assist in the interpretation of the complex activity patterns that emerge from these “conflicting zeitgeber” protocols, we present computer simulations of coupled circadian oscillators forced by two independent zeitgebers. This circadian system configuration was first employed by Pittendrigh and Bruce (1959), to model their studies of the light and temperature entrainment of the eclosion oscillator in Drosophila. Whereas most of the recent experiments have restricted conflicting zeitgeber experiments to two experimental conditions, by comparing circadian oscillator phases under two distinct phase relationships between zeitgebers (usually 0 and 12 h), Pittendrigh and Bruce compared eclosion phase under 12 distinct phase relationships, spanning the 24 h interval. Our simulations using non-linear differential equations replicated complex non-linear phenomena, such as “phase jumps” and sudden switches in zeitgeber preferences, which had previously been difficult to interpret. Our simulations reveal that these phenomena generally arise when inter-oscillator coupling is high in relation to the zeitgeber strength. Manipulations in the structural symmetry of the model indicated that these results can be expected to apply to a wide range of system configurations. Finally, our studies recommend the use of the complete protocol employed by Pittendrigh and Bruce, because different system configurations can generate similar results when a “conflicting zeitgeber experiment” incorporates only two phase relationships between zeitgebers.


Introduction
Synchronization of the physiology and behavior of organisms to the earth's periodic environment is achieved in part by the entrainment of circadian oscillators to species-specific combinations of daily photic and non-photic environmental cycles, known as ''zeitgebers'' [1][2][3][4].
The dominant zeitgeber that entrains circadian oscillators appears to be the light/dark cycle; however, the importance of food [5,6] and temperature in circadian entrainment [7,8] has received increasing recognition. In particular, temperature entrainment in the circadian organization of both ectotherms and endotherms [9][10][11][12], together with light effects, uncovers the necessity of considering the simultaneous action of two zeitgebers. The complexity of circadian systems comprising multiple oscillators entrained by two, or more zeitgebers can be understood through theoretical studies, where the contributions of zeitgebers and the internal circadian structure can be dissected. Such studies can also guide experiments and provide interpretations of the complex activity patterns of organisms.
Several experimental studies have altered the phase relationship between zeitgebers that occurs in nature, in order to assess their role in the entrainment of circadian rhythms [13][14][15][16][17][18][19][20][21][22][23][24][25]. The common experimental procedure is to artificially generate a phase difference of 12 h, in effect subjecting organisms to conflicting environmental time cues. In such ''conflicting zeitgeber'' experiments, it is commonly assumed that the circadian oscillator is phase-locked to the strongest zeitgeber. Intuitively, one might expect that some oscillators would follow the phase-displaced zeitgeber, whereas others would remain phase-locked to the unaltered zeitgeber.
To the best of our knowledge, Pittendrigh and Bruce (1959) performed the most complete set of conflicting zeitgeber experiments, revealing complex dynamics in the phase of the overt rhythm. Their experiments were intended to assess the relative strengths of light/dark and temperature cycles in the entrainment of the circadian oscillators controlling adult eclosion in Drosophila pseudoobscura. Populations of flies were raised under 24 h light-dark and temperature cycles that are zeitgebers in this species [26], the latter being successively phase-shifted by 2 h steps relative to dawn (Fig. 1). The eclosion rhythm tracks the phase of temperature cycle during the first 6 conditions, but tracks the phase of light cycle in the subsequent 6 conditions. This observed switch in the phase association of eclosion rhythm, first to the temperature and then to the light cycle, has made the identification of the stronger zeitgeber for the control of eclosion in Drosophila uncertain until now.
In the present work, we describe our numerical simulations of limit-cycle oscillator models [27] that explain the main features of the complex behavioral results of Pittendrigh and Bruce (1959). More importantly, our simulations illuminate the spectrum of dynamics generated by conflicting zeitgeber experiments by revealing how coupled oscillators respond to progressive phase displacements between zeitgebers. Our simulations show how conflicting zeitgeber experiments disentangle the complex interactions between oscillators and zeitgebers if it is ''complete''; that is, when the phase relationship between zeitgebers is progressively increased in small steps. We argue that a single phase displacement between two zeitgebers may generate misleading models of the circadian system.

Methods
We performed numerical simulations inspired by the model of Pittendrigh et al. [28] and the experimental protocol of Pittendrigh and Bruce [13] with a system of coupled limit-cycle oscillators (A and B), with each oscillator affected by one zeitgeber, L or T ( Figure 1B).

Oscillator equations
The A and B oscillators were simulated by coupled Pittendrigh-Pavlidis equations (1-4), where R and S are state variables, and a, b, c and d, are parameters.

A=B
Parameters C AB and C BA set the coupling strengths of oscillator A to B and of oscillator B to A, respectively. Zeitgebers L, for oscillator A, and T, for oscillator B, are square-wave functions with a 24 h period. These equations differ from the Pavlidis equations [27] by a variable K (Kyner), which is a small, nonlinear term that ensures numerical smoothness [29]. The R variables are explicitly Figure 1. Experimental data and model of Pittendrigh and Bruce (1959 constrained to be positive. These equations were developed originally for studies of the oscillator controlling the eclosion rhythm in Drosophila [27] and were employed in our modeling of the general properties of mammalian circadian oscillators [29][30][31]. As in former applications of this model, for the sake of simplicity and to better evaluate the effects of varying inter-oscillator coupling or zeitgeber strength, we assume that the oscillators are identical by fixing the parameters such that a A~aB~0 :85, b A~bB~0 :3, c A~cB~0 :8 and d A~dB~0 :5. This parameter set generates an oscillator with intrinsic period < 24 hr. Short, 1 h pulses (single pulse T-cycles with T~24 h [32]) were used in the simulations because they are the simplest mathematical models of daily zeitgebers. Zeitgeber amplitude (L~T~2) and coupling strengths (C~0:01 to 0.18) were chosen in a range that allowed a single zeitgeber to entrain the weakly coupled oscillators when alone. This choice of default values enabled exploration of coupling values that had the same effective magnitude as the zeitgeber strength.
Phases are defined with respect to the 24-hour day as follows (Fig. 2): Q A = acrophase of oscillator variable S A , phase at which S A takes its maximum value.
Q B = acrophase of oscillator variable S B . Q AB~QB {Q A = phase difference between coupled oscillators A and B.
W L = phase of zeitgeber L pulse onset. W T = phase of zeitgeber T pulse onset. W LT~WT {W L = phase difference between zeitgebers L and T. All phase values are given with respect to W L , which is assigned a value of 12 h. Thus, for example, a temperature onset phase of 8 h means that the temperature parameter T was set to 2, 4 h before light onset.
Simulations were performed with the CircadianDynamix software, which was developed to explore problems related to coupled and forced oscillators in chronobiology. It is an extension of Neurodynamix II [33,34]. We used the Euler method for numerical integration, with 1000 integration steps per 24-hour day.

Simulation Protocol
The phases of each oscillator (Q A and Q B ) were evaluated under a series of entrainment conditions by successively increasing W T , in one-hour steps, from +12 to +24 h, and then from 0 to +12 h. The final state of each entrained condition was used as the initial state of the subsequent condition. The reverse sequence, from +24 to +12 h and from +12 h back to 0 h was also employed in order to test for dependence on initial conditions. Furthermore, we focused on how Q A and Q B , at each W T , are affected by changes in the strength and symmetry of the interoscillator coupling and zeitgeber strength. The reference system was completely symmetrical, i.e. the oscillators, zeitgebers and interoscillator coupling were identical; asymmetry was added by incremental changes in relative strengths of inter-oscillator couplings or zeitgeber amplitude.

Symmetric System: identical oscillators, zeitgebers and coupling
To learn how the phase difference between zeitgebers (W LT ) affects the phases of the A (Q A ) and B (Q B ) oscillators, as well as their phase relationship (Q AB ), we set Pittendrigh-Pavlidis oscillator parameters to the simplest configuration: two identical oscillators with equal bidirectional coupling (C AB~CBA~C ). This system was subjected two identical, independent zeitgebers, whose relative phases were stepped from 0 to 24 h.
We first investigated the effects of coupling strength between the two oscillators. The output of the symmetric system is shown in Figure 3, where the inter-oscillator coupling strength was increased from C~0   oscillator is phase-locked to its zeitgeber and the phase relationship between zeitgebers and oscillator acrophases remains constant (Fig. 3A). We found that even very weak coupling between oscillators was sufficient to modulate Q A and Q B , as revealed by the changing relationships between the zeitgebers and the oscillators (Fig. 3B, C). These effects become increasingly prominent with increasing inter-oscillator coupling strength (while the zeitgeber strength was maintained); eventually, the oscillators appear more strongly influenced by each other than by their zeitgebers (Fig. 3D,  E).
A particularly interesting phenomenon is observed when zeitgebers T and L are near antiphase (W LT &12h). When interoscillator coupling is weak relative to the zeitgeber strength ( Fig. 3A-C), oscillators attain large phase differences, up to 12 h, by nonlinear, but smooth Q A and Q B changes. These phase differences are shown more clearly in Figure 4, where simulation data for the symmetric systems are replotted as Q AB as a function of W LT . These curves are sigmoidal, but with a sharp transition zone (for C~0:07 or greater). The maximum coupling value for which the antiphasic relationship between oscillators occurs (Q AB~1 2h when W LT~1 2h) is C~0:07 (Fig. 4A). Below this coupling value, Q AB is a nearly linear function of W LT , but above it the oscillators become tightly phase-coupled, acting more as a single system. Under this tight coupling, large values of Q AB are not allowed ( Fig. 3D and E, Fig. 4B).
Moderately abrupt changes in Q AB occur at W LT &12h, corresponding to the inflection point of the Q AB vs W LT curve, henceforth called ''inflection phase''. In this vicinity, small changes in W LT result in relatively large changes in Q AB , hereafter described as a ''phase jump''. The stronger the coupling, the steeper the phase jump, such as shown for C~0:15.
For coupling values larger than C~0:15, the results look rather different. In Figure 4B, (C~0:18, corresponding to data in Fig. 3E), as W LT increases from 0 to 24 h, Q AB remains near 0 h even beyond W LT~1 2h, revealing an unexpected asymmetry in the graph. Then, as W LT increases further, at D 2 (Fig. 4B) the value of Q AB jumps suddenly from just above 0 h to just below 24 h. With greater W LT values, Q AB remains below 24 h and a nearly linear function of W LT , attaining 0 h = 24 h for W LT~0~2 4h.
A complete picture is revealed by considering the converse sequence of zeitgeber phase displacements, from W LT~2 4h to W LT~0 h. Although most Q AB values are identical to those obtained before, the phase jump occurs at W LT v12h (at D 1 , Fig. 4B). Thus Q AB attains two steady-states in the interval D 1 vW LT vD 2 , depending on whether the previous W LT was larger or smaller than 12 h. The complete graph, generated by increasing and decreasing W LT (indicated by arrows, Fig. 4B), restores the symmetry that was apparently missing; the picture now comprises bistability and hysteresis [35].
Bistability implies the existence of two stable steady states Q AB , for the same W LT values located in the ''bistability zone'' D 1 vW LT vD 2 , as is also seen in the explicit Q A and Q B values (lines 12,13 and 14 of Fig. 3E) . The phases attained by the system depend on initial conditions and this dependence on the path of parameter change is a hallmark of hysteresis (arrows, Fig. 4B). Preliminary exploration of alternative parameter sets of Pittendrigh-Pavlidis equations have shown, however, that the bistability zone is reduced as the free-running periods of oscillators deviate from 24 h.

Asymmetric Systems a) Asymmetry in inter oscillator coupling strengths
We next examine the role of asymmetry by retaining equal zeitgeber strengths, but introducing asymmetry in the relative interoscillator coupling strengths. Using Pittendrigh-Pavlidis equations and departing from the symmetric case of a weakly coupled configuration C AB~CBA~0 :07, C BA was reduced from C AB to zero.
When Q AB is plotted against W LT (Fig. 5A), the main qualitative feature of the symmetric system (Fig. 4) is recognized; namely, the sigmoidal shape of the curve. However, as asymmetry is increased by having C AB wC BA , a progressive shift of the inflection phase occurs from W LT~1 2h to W LT v12h (Fig. 5A). Maximum departure of the inflection phase from W LT~1 2h occurs at the extreme, C BA~0 case, which corresponds to unidirectional coupling; that is, a ''master-slave'' configuration.
In the case of unidirectional coupling from A to B, Q A is always phase-locked to zeitgeber L, whereas Q B is modulated by both zeitgeber T and by L (via inputs from A; Fig. 5). Therefore, the dynamics of this asymmetrical system is now strongly dependent on the relative interaction strengths between zeitgeber T and C AB , which are pulling the slave B oscillator in opposite phase directions. Having fixed zeitgeber strengths (L~T~2), the master-slave system was now simulated for different coupling strengths C AB (Fig. 5B,C).
For a representative value of weak coupling (C AB~0 :07), Q B is a smooth function of W T until W T~6 h, where a phase-jump occurs (7 th line, Fig.5B), and with the oscillator following the zeitgeber T smoothly thereafter. The sudden increase in Q B at W LT~6 h, not at W LT~1 2h, corresponds to the shift in the inflection phase of the sigmoidal curve.
The overall picture is similar for stronger coupling values (C AB~0 :11, 0.15 and 0.18), except that T is not sufficiently strong to phase-lock Q B at some W T (9 h and 10 h, in this case). Thus, there is no stable entrainment at these W T , and oscillator B is in relative coordination [36] (Fig. 5C, lines 10 and 11, where Q B was omitted). As coupling strength increases, the range of W LT that yields relative coordination enlarges.

b) Asymmetry in zeitgeber strengths
Asymmetry in zeitgeber strength, with symmetric inter-oscillator coupling, yields results similar to asymmetrically coupled oscillator systems; namely, a sigmoidal Q AB vs . W LT curve, with a shift in the phase of inflection. The inflection phase is greater for greater asymmetry between the two zeitgeber strengths (Fig. 6A). Phase jumps again occur at values that differ from W LT~1 2h. For comparison, oscillator phases are shown for a zeitgeber L which is 4 times stronger than T (Fig. 6B) and conversely, for a zeitgeber T which is 4 times stronger than L (Fig. 6C), resulting in a shift of inflection phase at W LT smaller or greater than 12 h, respectively.

Interpreting Pittendrigh and Bruces (1959) Drosophila eclosion data
The following associations were made, in order to apply our model simulations to the Drosophila data.The phase relationship between zeitgebers was assigned the value W LT~0 when the phase of temperature minima (W T~0 ) occurred at dawn (W L~0 ), as represented in the first horizontal bar in Fig. 1A. Other complementary experiments described in [13] have shown that the phase of eclosion is determined by a temperature dependent oscillator B, which corresponds in our simulations to Q B . Furthermore, parallel simulations, not presented in this manuscript, have shown that the main dynamical features of periodic single-pulse zeitgebers are replicated by other cyclic wave forms, if appropriate amplitude adjustments are made.
We now focus on the following main features of the phase dynamics in Pittendrigh and Bruces (1959) data (Fig. 1A): 1. The phase shifts of the eclosion rhythm tracked the successive phase shifts of the temperature cycle in the 0vW LT v12h interval (horizontal bars 1 to 6). 2. An abrupt phase jump of the eclosion peak was observed when zeitgebers attained maximum conflicting phase differences; i.e., near W LT~1 2h (horizontal bar 7). 3. The eclosion phase remained nearly unaltered thereafter, independent of the phase of the temperature cycle along the 12vW LT v24h interval (horizontal bars 8 to 12).
First, our simulations support the two-oscillator model of Pittendrigh and Bruce because an alternative model, comprising a single oscillator with temperature and light inputs is equivalent to a master-slave configuration system (Fig. 5) that does not replicate the above features.
Our simulations of the two-oscillator symmetric system with intermediate coupling relative to zeitgeber strength (Fig. 3D), qualitatively reproduces the findings in Drosophila. During the interval 0vW LT v12h, the phase-shifts of oscillator B track the phase shifts of the temperature zeitgeber. When W LT &12h, a ''phase jump'' of Q B occurs; thereupon Q B tracks the phase of the light zeitgeber for 12vW LT v24h. Our simulations have thus shown that the switch at W LT &12h of the preferential phase association of eclosion rhythm, first to the temperature and then to the light cycle can arise even in the most symmetric configuration, with identical oscillators coupled to identical zeitgebers (Fig. 3). By using a simple symmetric configuration, we show that additional complexity (unequal parameter values), even if it exists, is unnecessary for generating the observed complex phenomena. In other words, the results of the Pittendrigh and Bruce experiment do not indicate any zeitgeber dominance in Drosophila, neither of light/dark nor of temperature.
Most of the experimental studies that altered the natural phase relationship between photic and non-photic zeitgebers have restricted it to two W LT values, usually 0 and 12 h. Furthermore, they have also assumed linear relationships between the phases of zeitgebers and oscillators [18][19][20][21][22][23]25,37]. Such assumptions lead to the proposition that the phase of an output rhythm tracks the phase displacement of the stronger zeitgeber and thus two W LT conditions are sufficient for revealing the zeitgeber hierarchy. While this is correct when there is no coupling between oscillators (Fig. 3A), inter-oscillator coupling effects (Fig. 3B-E) add more complexity to the conflicting zeitgeber experiment, requiring a systematic change in W LT conditions to avoid ambiguous models.
The complex picture that emerged from Pittendrigh and Bruces experiment, comprising 12 W LT relations (Fig. 1A) is best replicated by the features of Fig. 3D, where light is as strong as temperature in the entrainment of the eclosion oscillator. However, if their experiment were restricted to only two phase relations between zeitgebers, for instance, W LT~0 and W LT~1 2h (horizontal bars 1 and 7 of Fig. 1), the more limited set of phase results could be replicated by several other configurations.Furthermore, the problem persists if a different W LT pair is chosen for a conflicting zeitgeber experiment.
Some predictions and guidelines for a complete conflicting zeitgeber experiment arise from our simulations. The hypothesized general system is composed of two oscillators, differentially entrained by two zeitgebers. If the phase of the output rhythm (or of the oscillator itself) is plotted against the phase difference between zeitgebers Z 1 and Z 2 , W Z1Z2 a curve with different characteristics is expected.
N If the resulting graph is linear, instead of sigmoidal, there is no coupling between oscillators. Alternatively, the strength of one zeitgeber may not have been sufficiently strong to entrain its corresponding oscillator or the observed output is controlled by a master oscillator without feed-back from the unobserved slave (as it would occur in Fig. 5B, if the phase of zeitgeber L were shifted, instead of T).
N For a sigmoidal graph, the steeper the curve, with more pronounced phase jumps, the stronger is the inter-oscillator coupling with respect to the zeitgeber strengths (e.g., Fig. 3).Phase-jumps occurring at W Z1Z2 =12h are indicative of either an asymmetric system, with unequal coupling/zeitgeber strength between the oscillators, or a bistable system, with two potential entrained oscillator phases for an interval of W Z1Z2. N Dependence on initial conditions, with two possible entrained phases, is predicted to occur around the phase jump region if the inter-oscillator coupling is strong compared to the zeitgebers. This can be tested by pre-entraining with either of two W Z1Z2 : in the natural phase relation or in W Z1Z2 greater than the phase of inflection.
N Relative coordination is expected to occur for some W Z1Z2 values, when coupling is strong with respect to zeitgeber strength.
Some other experiments have reported our predicted shift in the inflection phase, typical of asymmetric systems. This was observed by Bruce [14] and Pittendrigh [40], when they assayed, respectively, the circadian phototaxis rhythm in Euglena viridis and activity in cockroaches under light/dark and temperature cycles. Recently, Watari and Tanaka [24] verified that phase jumps occur at a W LT w12h in their conflicting zeitgeber experiments on the eclosion of onion fly Delia antiqua. Since only fixed initial conditions were used in these experiments, phase-jumps at W LT =12h could alternatively reflect a single branch of a bistable system (Fig. 3E), which would be testable by manipulating initial entrainment conditions.
There is much to be understood about the roles and interactions between non-photic and photic daily cues in the circadian organization of different species, in the context of different environments [15,[41][42][43][44][45]. Theoretical studies of multiple oscillator models explained phenomena that could not be accounted for by a single oscillator [32,[46][47] and widened opportunities for creative biological experiments [11,12,31,[48][49][50][51]. We are now facing the challenge of including multiple zeitgebers in this scenario, but in these early steps, our modeling study is appropriately limited to two general cases: symmetrical systems and those with simple asymmetries. Despite its simplicity, our system sufficed to give a glimpse of the rich dynamics behind multiple zeitgeber phenomena.