Using Variability to Regulate Long Term Biological Rhythms

J. theor. Biol. (1999) 196, 455–471 Article No. jtbi.1998.0843, available online at http://www.idealibrary.com on

Using Variability to Regulate Long Term Biological Rhythms T S*†  K S‡ *Department of Chemical and Biochemical Engineering, Rutgers University, Piscataway, NJ 08854, U.S.A. and ‡Department of Neurobiology and Physiology, Northwestern University, Evanston, IL 60208, U.S.A. (Received on 8 April 1998, Accepted in revised form on 21 September 1998)

We present a model for the generation of precise, long term rhythms from a collection of imprecise, short term oscillators. The model uses variability between oscillators in conjunction with simple coupling rules to produce long term rhythms that are independent of rate equations (e.g. Arrhenius). The rhythms generated by the model are controlled by only two independent parameters and exhibit several physiologically interesting properties, including ready entrainment to external signals and splitting in response to strong constant signals. The model provides several predictions that can be tested in future experiments. 7 1999 Academic Press

1. Introduction Great progress has been made over the past decade towards understanding long term biological rhythms. Much of this progress results from the interaction between theoretical (Glass & Mackey, 1988; Winfree, 1990) and experimental studies. Experimental evidence suggests that circadian rhythms are generated by slow transcription/translation feedback loops involving one or only a very few proteins (Takahashi & Kornhauser, 1993; Rosbash, 1995; Loros, 1995; Kay & Millar, 1995). The fact that single gene mutations have profound effects on overt rhythmicity is consistent with this idea. However, such a mechanism cannot by itself account for several key observations in the field of mammalian circadian rhythms and theoretical models of coupled oscillators continue to make

†Author to whom correspondence should be addressed. E-mail: shinbrot.sol.rutgers.edu 0022–5193/99/004455 + 17 $30.00/0

significant contributions in the understanding of this physiological process. The master pacemaker in mammals resides in the suprachiasmatic nucleus (SCN), a neural structure of approximately 104 cells located within the hypothalamus (Moore-Ede et al., 1982; Klein et al., 1991). Many SCN neurons are spontaneously active in firing action potentials (Green & Gillette, 1982) yet experiments using tetrodotoxin to block Na+-dependent action potentials have proven conclusively that the clock’s timekeeping persists in the absence of neuronal firing (Schwartz et al., 1987; Welsh et al., 1995). Although there is increasingly convincing evidence that circadian rhythms exist within single cells (Welsh et al., 1995; Michel et al., 1993), experimental evidence has not provided clues as to how these slow oscillations may be coupled together to maintain synchronized phases and communicate temporal information to other brain structures. Thus, even if we knew exactly how cells generate a circadian 7 1999 Academic Press

456

.   . 

rhythm individually, the behavior of this collection of cells would still be largely unknown. This difficulty is compounded by the fact that neurons have their own dynamics, which there is good reason to believe are nonlinear (Wallenstein, 1993; Halloy et al., 1990; see also Klevecz & Bolen, 1993). Examples of nonlinear behavior by neurons include refractory periods and thresholds for firing. The model presented here makes use of the two entirely different types of oscillators known to be present in the mammalian SCN: high frequency nonlinear neuronal oscillations and slow intracellular biochemical oscillations. A second unexplained issue in the field of circadian rhythms is that transcription/ translation rhythms depend on Arrhenius rate equations—yet it has been well known for several decades that excepting transients during phase shifting, the periods of circadian rhythms are essentially temperature independent (Pittendrigh, 1954; Barrett & Takahashi, 1995). This result holds notwithstanding recent experiments which suggest that the amplitude of the circadian clock may vary as a function of temperature (Barrett & Takahashi, 1997). Third, models intended to account for a minimal set of factors involved in such a transcription/translation loop can include as many as 15 independent parameters (Goldbeter, 1995). If the system variability due to noise in a single parameter is 5%, for example, the total variability of such a system is expected to exceed 19%. By comparison, actual variability in mammalian circadian rhythms is typically two orders of magnitude smaller (Enright, 1980b). This key observation in the field of circadian rhythms will be considered here. Finally, certain physiological observations— such as the sharp onset of locomotor activity and the substantially less precise offset exhibited by hamsters and mice—have not yet been adequately explained in any model (Enright, 1980a). In the circadian context, it is known that despite irregularities that are observed in coupled nonlinear systems (Matthews et al., 1991; Strogatz, 1994), SCN neurons increase and decrease their firing rates regularly and reproducibly over the course of a day (Green & Gillette, 1982; Welsh et al., 1995). Moreover, in

nonlinear systems one typically sees a transfer of energy between short and long time-scales (Landau & Lifshitz, 1959)—hence it is not at all obvious how to entrain neuronal models which fire on a time-scale of tenths of seconds to a long term rhythm operating over a 24 hour period. Thus, modelers have some explaining to do. In this paper we present a candidate model with the following characteristics: 1. a robust and precise long term collective rhythm composed of much faster and more variable individual firing rates. This rhythm: 2. is independent in the long term of (Arrhenius) rate equations; 3. is controlled by only two independent parameters; 4. can exhibit very precise activity onset time and much less precise offset time; 5. is produced by nonlinear oscillators, yet readily entrains to low frequency external oscillations; 6. can exhibit ‘‘splitting’’ (Pittendrigh, 1974) in response to certain types of steady signals. In Section 2 we describe the model, in Section 3 we examine its endogenous behaviors, and in Section 4 we analyse interactions between the model oscillations and external rhythms. 2. Model  Recent mathematical and physiological research has demonstrated that noise can play a role in enhancing the sensitivity of biological detection mechanisms (Douglass et al., 1993; Longtin et al., 1991; Levin & Miller, 1996; Collins et al., 1996.) Since noise is ubiquitous in biological systems, one might ask whether noise can do more than enhance the capability of an existing system: one might ask whether noise can actually be used to control complex biological systems. In the context of circadian rhythms in mammals, two facts bear scrutiny. First, mammalian circadian rhythms are controlled by a locus of neurons. These neurons fire on a short time-scale, on the order of tenths of seconds, and biological rhythms on this time-scale are known 2.1.

      to exhibit non-stationary noise (Enright, 1980a; see also Mainen & Sejnowski, 1995; Arieli et al., 1996). The period of firing of these neurons varies by as much as a part in 10—yet the overt circadian period varies from cycle to cycle by only about a part in 500 (Enright, 1980b). Second, the SCN is uniquely characterized by a high level of local circuit synapses (van den Pol & Dudek, 1993; Sommers & Kopell, 1993) as well as by non-synaptic appositions between neuronal membranes (van den Pol, 1980) which may represent important functional coupling between oscillators (Bouskila & Dudek, 1993). Thus, the SCN can be viewed as a collection of noisy oscillators coupled to one another. In the following subsection we show that coupling multiple oscillators together in the presence of variability containing drift necessarily produces a new rhythm with a long time-scale (which qualitatively goes as the coupling strength over the amplitude of the drift component of the noise). Thus the troublesome aspects of noise can be turned to advantage to generate a long term, precise rhythm out of shorter term, imprecise components. We emphasize that coupling together noisy oscillators will produce a long term rhythm independent of any other biochemical rhythm that may exist. We hypothesize that this long term ‘‘neuronal’’ rhythm may be the workhorse for communicating circadian time to other brain areas regulating system-wide physiological rhythms. In our model, overt physiological rhythms depend on the ability of the circadian clock (arising from a transcription/translation feedback loop) to communicate phase information to a long term neuronal rhythm (arising from coupling of noisy oscillators). Circadian regulation of membrane conductance, as demonstrated in the pacemaker neurons of Bulla gouldiana (Michel et al., 1993), is one possible mechanism that a molecular clock could use to communicate phase information to *A notable exception may be found in the firefly, P. Malaccae; see Ermentrout, 1991; for other work on frequency variations see Matthews, 1990, 1991. †An interesting analog may be found in the dance of honeybees, where individual bees adjust their waggling frequency to match their neighbors’ (Stabentheiner, 1990). ‡The case f Q 0 is non-physical, for it leads to exponentially diverging frequencies.

457

a neuronal rhythm. Because of its relative independence from any 24 hour biochemical rhythm, this long term ‘‘neuronal’’ rhythm would also function to limit the effect of slight differences in phase among the multitude of individual circadian clocks in the SCN.      To gain insight into the mechanisms underlying biological rhythms, it has proven useful to model collections of neurons as coupled oscillators (Winfree, 1967; Pavlidis, 1969; Wilson & Cowan, 1972; Clay & DeHaan, 1979; Dye & Heligenberg, 1987). In most such models, one assumes that the oscillators have a fixed intrinsic frequency*, and one examines the effect of amplitude or phase modulation on these oscillators (Cohen et al., 1982; Ermentrout & Kopell, 1984, 1985; Kuramoto, 1975, 1984; Sherman & Kinzel, 1991, 1992; Matthews et al., 1990, 1991; Daido, 1993). We consider a collection of coupled limit cycle oscillators (van der Pol, 1926; Pedersen & Johnsson, 1994; Wever, 1984), with variable, non-stationary frequencies (cf. Matthews et al., 1991). Each oscillator is defined by the complex-valued state, Zi , obeying: 2.2.

Z j = n ·Zj (1−=Zj = 2) + ikj ·Zj+k(Zj − Zj ), (1) where n defines the level of dissipation of the limit cycles, kj is the frequency of the j-th oscillator, and Zj is the average state of all but the j-th oscillator. The coupling constant, k, is intended to account for short time, e.g. synaptic, neuronal coupling. This equation uses linear, all-to-all coupling and is among the simplest of coupled oscillator models. The states of the oscillators are coupled together, but in addition there is good reason to believe that a neuron will adjust its frequency in response to biochemical signals from its neighbors† (Traub et al., 1989). A simple mathematical description for adjusting the frequency in this way is: ktj = f ·(kj−kj ),

(2)

where f is a positive‡ constant, and again kj is the average frequency of all but the j-th oscillator. This equation provides that the j-th oscillator’s frequency will tend to the mean of its

.   . 

458

(4)

Here D is a constant† that gives rise to linear drift, and N0j (t) is a noise term containing zero-mean fluctuations. The drift takes on alternate signs with equal probability. This choice is adopted because there is no reason, a priori, to assume otherwise. This means that some neurons would tend, in isolation, to oscillate faster than average, while others would tend to oscillate slower. Thus, the addition of a simple form of biological variability to each oscillator divides the coupled oscillators into two subgroups. This is a crucial element of the model, for as we will see including biological variability in this manner leads directly to robust long time cyclicity. In separate simulations without drift, no low frequency rhythm was observed from the collection of coupled high frequency oscillators.

*Recall that k defines the individual oscillator frequency, i.e. 2–3 milliseconds for neurons, and k is determined by fast, e.g. synaptic, transmission. †This is the lowest order non-stationary form of variability. Other forms could also be investigated; for example it would be interesting to determine the effect of broadening the distribution of D values. We are grateful to an anonymous referee for bringing this point to our attention.

 In this section we have described a model which combines variability with frequency and phase coupling. The model contains two essential elements: a limit cycle oscillator [eqn (1)] and frequency adjustment [eqn (2)]. Drift plus zero-mean noise is added to the frequency

neighbors. Thus, this model couples both the phase (1) and the frequency (2) of its component oscillators. Equations (1) and (2) define a very simple model of coupled oscillators. As in most coupled oscillator models, this model has only fast* time scales, given here by 1/k and 1/k. In the presence of biological variability, however, a longer time-scale can emerge. The frequencies of the individual oscillators are, like all other biological parameters, subject to intrinsic biological variability. To incorporate this variability in our model, we add a term to eqn (2): ktj = f ·(kj−kj ) + Vj (t),

(3)

where we write the variability term, Vj (t), as: Vj (t) = 2D + N0j (t).

2.3.

Max

Min

Min

Amplitude

100

0

–100 0

2 f / 2D

4 f / 2D

6 f / 2D

8 f / 2D

10 f/ 2D

Time

F. 1. Collective behavior of 100 model oscillators obeying eqns (1), (3) and (4). Main figure: sum of real part of limit cycle oscillators is plotted as a function of time (in arbitrary units). All oscillators are initially at random states, but organize rapidly to exhibit regular collective rhythm of period t. Inset: detail shows large number of short term oscillations contained within the collective envelope.

      to simulate biological variability [eqns (3) and (4)]. We turn next to an examination of behaviors of this model. 3. Behaviors of Model   In Fig. 1, we show the sum of amplitudes of the real parts of 100 oscillators coupled according* to eqns (1), (3) and (4), where each oscillator is started at time t = 0 with random initial phase (and unit amplitude). The long time-scale, t, shown is many times the period of an individual oscillator (see inset). This rhythmic behavior is autonomous, emerges spontaneously†, and persists indefinitely‡, in the presence of noise, and despite the use of random initial phases. The zero-mean noise [N0j (t) in eqn (4)] produces fluctuations that are evident from scrutiny of the amplitudes in the main figure. N0j (t) in this figure was chosen to be uniform, with maximum amplitude of five times the value of D. The overt period, t, depends formally on five parameters: the linear coupling, k, the rate of frequency adaptation, f, the drift, D, the dissipation, n, and the mean frequency, ki . In practice, the asymptotic value of t does not depend perceptibly on changes in dissipation or mean frequency: the dissipation governs how 3.1

*In this plot, we display results from 4th order Runge-Kutta integration of eqns (1), (5) and (6) using a linear function F; similar results have been obtained for more complicated functions. The parameters used were as follows: k = 100, k = 1.25, f = 1.0, n = 100, D = 1/2p, and N(t) was a white noise process with maximum amplitude of 25/2p, applied every integration timestep Dt = 0.001. These values were chosen based on computational convenience and are not intended to be physiologically representative. †This rhythmic envelope emerges extremely quickly: note in Fig. 1 that the phases are initially randomized at t = 0. ‡In our numerical experiments, we integrated over more than 50 000 individual oscillations without change in the long term rhythm. §Oscillators which are disturbed from this state, due either to their initial condition or to noise, will be drawn to the bimodal state exponentially quickly, with characteristic time, 1/f. ¶The period given in eqn (6) applies to the weak coupling limit, k1. As k becomes substantial, its influence must be included, as we will discuss shortly.

459

rapidly the system will respond to external stimuli, and varying ki  numerically (discussed shortly) over two orders of magnitude produces a maximum change in t of under 2%. As for the remaining parameters, we can derive an approximate expression for t by averaging eqn (3) over time, obtaining the solution: ki = ki 2

D + oi e−ft, f

(5)

where oi is a constant determined by the initial state of the oscillators. Equation (5) tells us that oscillators coupled as we have described tend§ to a bimodal distribution of frequencies, centered around the two values k2 = k 2 D/f. Consequently the two sets of oscillators may be expected to beat on the time-scale,¶ t3

f . 2D

(6)

Two points bear mentioning. First, as we have indicated, this relation does not depend significantly on the mean oscillator period, kj . Thus if the mean period changes due, for example, to temperature fluctuations, this may affect the envelope rhythm transiently during the change, but should produce no effect on the ultimate period. Separate simulations confirm that changes in kj  by an order of magnitude in either direction have no significant effect on t. In this sense the t of this long term rhythm is ‘‘temperature’’ compensated, as is the circadian clock. From a physiological perspective it is perhaps surprising that the actual firing rate is much less important for the emergence of long term rhythmicity than is the model neuron’s spontaneous activity and noisiness. Second, eqn (6) is valid in the weak coupling limit: k1. For more strongly coupled systems, t is lengthened, because coupling tends to hold oscillators near the same state, thus resisting the change in state which leads to beating. For example, in Fig. 1, k = 1.25, and correspondingly t is somewhat larger than f/2D (t = 1.21 f/2D). Simulations performed using k = 0.01 gave t = 1.007 f/2D. So for small k, t increased linearly with f as expected from eqn (6). In Fig. 2(a), we plot variations in overt period, t, as a function of k, where the other parameters

.   . 

460 10

10

(a)

(b)

– f/2D

1

0.1

1

0.01

0.001 0.0

0.1 0.0

2.0

1.0

2.0

1.0 f

F. 2. Response of envelope period (t) to increases in either the phase coupling constant k or the frequency coupling constant f. Part (a) shows the response of t to increases in k when f = 1, part (b) shows the response of t to increases in f when k = 1.

are taken at the nominal values*, f = 1. t is very sensitive to changes in k (note the logarithmic scale). Figure 2(b) shows the dependence of t on f, when k is maintained constant. Clearly, t is highly sensitive to changes in f as well as k.

is of interest. In Fig. 3, we plot the effect of the number of components, M, on the precision of the envelope rhythm. In accord with the central limit theorem (Kreysig, 1979; see also Enright, 1980a), the precision increases approximately with 1/zM. Thus, there is a trade-off between increasing precision and increasing the number of components: one can equivalently study a small number of oscillators with a small amount of noise or a large number of components with a large amount of noise—the variability in t of these two systems are expected to be equivalent. For obvious computational reasons, it is desirable to limit the number of components. In the preceding results, we have included 100 oscillators in the model; in the results which follow we use only six oscillators. To ensure that

 Since this model explicitly depends on variability, the precision of its envelope rhythm 3.3.

*Except as noted, all subsequent simulations use the nominal paramater set: k = 100, k = 1.25, f = 1.0, n = 100, D = 1/2p, N0(t) = r[−5/2p,5/2p], where r[·] is randomly chosen on [·]. Simulations performed using 4th order Runge-Kutta integration of eqns (1), (3) and (4) using a time-step of 0.001 (confirmed to be an order of magnitude smaller than required to prevent anomalous numerical artefacts).

Variance in period,

2

1 2

0.1

=

(9±

3)·

M (–

1.0

±0

.1)

;R2

=

0.01

0.001 10

100

0.9

4

1000

10 000

Number of components, M

F. 3. The precision of the envelope period, t, depends on the number of coupled oscillators in the system. Variance in t is plotted against the number of components used in simulations. Each point represents the variance in envelope period calculated after completion of 10 long term cycles.

      spurious results* do not creep in, we have verified key dynamical features (e.g. the splitting transition reported in Section 4.2 following) using 100 oscillators. The simulation for Fig. 3 used N0j (t)$[−5,5], i.e. a zero-mean noise of up to five times the value of the drift, D. The resultant envelope period was comparatively precise, for a few hundred components the variance was less than 0.1 (in units where the envelope period, t = p) or less than 0.01 if more than 1000 individual oscillators were coupled. These results suggest that sufficient neurons do exist in the mammalian SCN to explain observed circadian precision through this mechanism, even assuming very large noise levels. This aspect of the model has interesting physiological ramifications. We do not know factors that might affect the relative amount of noise associated with neuronal oscillations, however one can imagine factors

*For example, it is conceivable that a ‘‘splay’’ state (Silber et al., 1993)—in which individual oscillator states repel one another—could be present in some simulations. Neither this nor any other anomalous state is observed in any of our simulations.

461

that may influence the number of functional oscillators within the SCN. Inherent in this model is the concept that circadian t should become less precise under any condition where the number of functional neuronal oscillators (or their coupling) is diminished. Indeed, lack of precision in circadian period occurs when transplanted SCN cells drive overt rhythmicity (Silver et al., 1990) and during normal aging (Scarbrough et al., 1997) where a decrease in the number of cells immunoreactive for SCN neuropeptides (Chee et al., 1988; Roozendaal et al., 1987) may be taken as evidence for a decrease in the number of functional oscillators.      This system has two control parameters, k and f. The envelope period is more sensitive to changes in k than f, but either parameter can be used to change this period. We will discuss ramifications of this fact shortly. For the time being, it is interesting to note that the qualitative shape of the envelope waveform differs for the cases of large and small k. That is, as indicated in Fig. 2, we can produce a specified period for multiple different pairs of k and f, but the shape 3.3.

Amplitude

100

= 0.1 f = 2.2

0

–100

Amplitude

100

= 2.0 f = 1.0

0

–100 0

2 f1/ 2D

4 f 1/ 2D

6 f 1/ 2D

8 f 1/ 2D

10 f 1/ 2D

Time

F. 4. Waveform of the envelope rhythm depends on choice of coupling constants. Top panel shows weak phase coupling, bottom panel shows strong phase coupling. Values of frequency coupling were chosen to yield the same t on the top and bottom panels.

462

.   . 

of the envelope will in general differ for different choices of these values. For example at the top of Fig. 4 we show a symmetric envelope waveform produced in the small phase coupling regime, k1, and for comparison at the bottom of the figure we show an asymmetric waveform produced by increasing k and decreasing f. This suggests that, to take a practical example, different organisms could produce the same overall period but could have different activity patterns. For instance, if we conjecture that an animal becomes active when the envelope rhythm grows beyond some threshold and rests when the rhythm falls below the same threshold (i.e. the dotted lines shown in Fig. 4) then for the symmetric, small phase-coupling case we would expect to see gradual and symmetric activity onset and offset times. By contrast, for the asymmetric, strongly phase-coupled case, we would expect to see a rather sudden and precise onset time and a more gradual and imprecise offset time. The underlying cause of the asymmetry in the envelope rhythm can be understood heuristically as follows. For small k, the envelope is simply defined by a symmetric sinusoid. This is because as k : 0, the envelope rhythm becomes insensitive to phase [cf. eqn (1)]. For larger values of k, however, this sinusoid is distorted because the coupling resists the separation between oscillator phases. To appreciate this point, it may be helpful to recall that the minimum collective amplitude occurs when the two sets of oscillators (one slightly faster; the other slightly slower) are p out of phase—hence their amplitudes sum to zero. This is an unstable state, and the coupling constant, k, is a measure of the magnitude of the propensity to return to the in-phase state: for large k, there is a large force pulling the oscillators back into phase with one another. Because of this large force, the oscillators are rapidly drawn to the maximum collective amplitude state where all oscillators are at the same phase. The converse is also true: once the oscillators are in phase and at the maximum amplitude, the coupling resists a return to the minimum, out of phase, state. This causes the envelope rhythm to rapidly move to the in-phase, maximum amplitude state, and to

slowly return to the out-of-phase, minimum state, as is observed in the lower plot of Fig. 4.  In this section we examined several qualitative as well as quantitative behaviors of the model defined earlier. We saw that the model produces a spontaneous long term rhythm (Fig. 1), and we suggested that this rhythm is analogous to beating [eqn (6)]. We found that increasing either the phase coupling constant, k, or the frequency coupling constant, f, causes a rapid lengthening in the period of this rhythm (Fig. 2). In addition, we found that the precision in this long term rhythm obeys the central limit theorem, as one might expect (Fig. 3). Finally, we described how the overt rhythm exhibits an abrupt onset for strong phase coupling. 3.4.

4. Interactions between the Coupled Oscillator Rhythm and Driving Rhythms We have so far discussed a mathematical model which uses variability between components and coupling to generate a long term and precise periodic rhythm. Both the model’s overt behavior and its response to simple parametric changes suggest some similarity to circadian behaviors. This model is, however, deficient in three respects. First, data strongly indicate that there exists a circadian rhythm in the SCN independent of neuronal firing (Schwartz et al., 1987; Earnest et al., 1991; Welsh et al., 1995). Second, as we have already mentioned there is now good evidence for the presence of a distinct long term transcription/translation rhythm that generates overt circadian rhythmicity in both animals and plants (Takahashi & Kornhauser, 1993; Kay & Millar, 1995; Loros, 1995; Rosbash, 1995). Third, the existence of circadian clocks in non-neuronal tissues and recent data regarding a circadian rhythm in melatonin secretion by cultured hamster retina (Tosini & Menaker, 1996) indicate that however the SCN regulates itself, other endogenous circadian rhythms surely exist. Thus, a study of interactions between the model rhythm and external stimuli—especially other periodic rhythms—is indicated. There are

      several ways of including the influence of other rhythms on this long term oscillation. We choose one method for coupling two rhythms and present the results of the interaction in a manner familiar to circadian biologists. We have seen that in the presence of drift, the coupled oscillator system described tends to self-assemble into an autonomous system of fast oscillators which beat over a much longer time-scale. Because this beating depends on variability between individual oscillators, the long time-scale is not highly sensitive to small fluctuations in oscillator dynamics. This characteristic distinguishes this model from other beating models.       Given a level of variability, the period of the beats depends on two parameters, f and k. The former governs the variation of frequency, while the latter directly affects the relative phase of the faster- and slower-sets of oscillators. In this sense, one parameter (f) can be viewed as a frequency-modulating term, and the second (k) as a phase-modulating term. This being the case, the most straightforward way to entrain the system to external cues is to vary f or k directly in response to the cues. Thus, our physiological correlate would be that the relevant effect on the clock of an external cue such as light, temperature, or pharmacological agents, would be to increase or decrease f or k. In our simulations, we have found that either parameter can be used to entrain the long term rhythm to external stimuli and the effects of varying f or k are nearly interchangeable. In the results which follow, we fix f and vary k. Physiologically speaking, there is no intrinsic reason to assume that a particular stimulus, e.g. light, should increase or decrease k, so we will explore both directions of change. In a later section we will discuss possible connections between these mathematical observations and experimental observations of the circadian system. 4.1.

*The Phase Transition Curve (PTC) may be more familiar to many theoretical biologists: the PRC differs from the PTC in that the PRC plots the difference between future and past phases as a function of time, while the PTC plots the future phase vs. the past phase.

463

First let us focus on the effect of a brief, pulsed, change in k on the state of the oscillators. In all of the simulations performed, we found t to be very rapidly restored after a pulse, so we plot only the change in phase of the oscillators. The pulse consists of increasing k briefly (the equivalent of 7 minutes in a 24 hour day) by a prescribed pulse amplitude and then returning k to its nominal value. The minima of the overt signal (cf. dashed lines in Fig. 1) are convenient, readily identifiable markers of clock phase, and to choose a standard measure, we plot the phase of the 5th minimum following a pulse, relative to the control case without a pulse. In Fig. 5 we plot the resulting change in phase in response to positive (open symbols) and negative (filled symbols) pulses. The axis units are arbitrary; for physiological comparisons we associate one complete cycle of the envelope shown in Fig. 1 with a 24 hour cycle, with 0:00 placed at the minimum of the envelope. Thus we construct what is known in the biological literature as the Phase Response Curve (PRC*) (DeCoursey, 1960; Pittendrigh, 1960; Winfree, 1990). Amplitudes of the pulses (in arbitrary units) are indicated in the inset. Our model ‘‘clock’’ exhibits a phase dependent difference in the amplitude of response to the pulse, qualitatively similar to that seen in the mammalian circadian system. That is, the amplitude of the phase shift depends on the time of day that the brief pulse was administered. From Fig. 5 we see that a positive stimulus before about 07:00 tends to delay the model clock, and vice versa for stimuli between 07:00 and about 23:00. Negative pulses of k advance the clock when administered between 0:00 and 07:00 and delay the clock when administered between 07:00 and 23:00. The clock responds with larger phase shifts during the period of time labeled 07:00 to 23:00. This has a fairly straightforward dynamical cause, which is linked to the asymmetry observed in the envelope of the collective oscillations shown in the lower plot of Fig. 4. Values of k and f that give rise to the symmetric envelope rhythm shown in the upper panel of Fig. 4 would be associated with symmetric PRCs, where the amplitude of the advance region would be very similar to the amplitude of the delay region.

.   . 

464

Moreover, from Fig. 5 we see that the circadian phase at which a positive pulse would yield advances is very different from the circadian phase at which a negative pulse would produce advances. The same is true for phase delays of the model clock. In addition, the relative magnitude of phase advances and phase delays differ according to the sign of the pulse. These aspects of the model may have a corollary with the two types of PRC observed in animals (Smith et al., 1992). One final observation from Fig. 5 is that large stimulus amplitudes give rise to phase changes in magnitude exceeding half the circadian cycle, but analysis of corresponding Phase Transition Curves (Winfree, 1990) reveal that this model nevertheless uniformly produces type 1 resetting. It is not known at this point whether this result is generic, or whether modifications of this model could produce other types of resetting. The abrupt change from advances to delays seen here in response to large pulses of k are also observed in rodent behavior in response to light pulses. The plot of the PRC shown here follows

convention in the circadian literature and no physiological significance is ascribed to the slopes of the curves.      Responses to discrete pulses of light are theorized to be the mechanism by which animals entrain to a light–dark cycle (Pittendrigh, 1965). To investigate the response of the model to a second periodic signal, we add the term A ·sin(2pt/t0), to the amplitude of the coupling constant, k. In our simulations we set A = 1. This second signal is intended to be analogous to circadian modulation of the phase coupling between the short term oscillators. Other methods of including periodic forcing are also possible. To distinguish semantically between the two rhythms, we call the rhythmic envelope of the coupled oscillators the ‘‘overt’’ rhythm, and the rhythm added to k the ‘‘signal’’ rhythm. This fits with the idea that there is an underlying rhythm of molecular origin which drives the neuronal rhythm, and that in animals this neuronal rhythm in turn drives overt functions 4.2.

+16:00

+12:00

Phase change

+8:00 +4:00

0:00 –4:00

–8:00

–12:00 –16:00 00:00

04:00

08:00

12:00

16:00

20:00

24:00

Pulse phase

F. 5. Phase response curves for oscillator model. We find two phases (here at 07:00 and 23:00) at which there is zero phase change in response to a pulse. Moreover, there are two distinct kinds of phase response curves, denoted in this figure by open vs. closed symbols. Similar results have been found in animals (Smith et al., 1992). Increased pulse amplitudes predictably accentuate the change in clock phase, and this change can exceed 12 hours. Nevertheless the phase resetting remains type 1. Pulse amplitudes: (w) + 10, (W) − 10, (q) + 20, (Q) − 20, (r) + 40, (R) − 40, (p) + 60, (P) − 60.

     

465

Signal rhythm

Overt rhythm

Period,

2 ˆ0

ˆ0

Baseline

0

ˆ0

2 ˆ0

Signal period,

0

F. 6. Variation in ‘‘overt’’ period, t, with changes in ‘‘signal’’ period, t0. The overt period entrains to the signal in a range of about 0.5tx0 Q t Q 1.5tx0, where tx0 is the overt period in the absence of a signal rhythm. For tx0 Q t, the overt period loses entrainment with the signal, while for tx0 q t, the overt period breaks into a period-two pattern, with alternating short and long cycles. Insets show relationship between signal rhythm (top of each inset) and overt rhythm (bottom).

such as rest/activity, feeding, or other metabolic cycles. In Fig. 6, we plot the response of the model to variations in t0, the period of the signal rhythm. Each open symbol represents the overt period on successive iterations (‘‘days’’) at the given signal period, and we include error bars at each signal period examined. During much of the range of entrainment the error is so small as to be occluded by the symbols. The ‘‘natural’’ period of the oscillator model in the absence of the signal rhythm is denoted tx0 in the figure, and the insets show the shape and qualitative behavior of the overt rhythm in response to various signal periods. There is a range in signal period, from about 0.5tx0 to 1.5tx0, for the parameter values studied here*, within which the overt oscillator period is entrained exactly to the signal period; this range varies with the chosen amplitude, A,

*For larger k, the range in entrainment grows (cf. Fig. 7 and footnote on p. 467).

of the signal rhythm. This result seems similar to the situation in animals where there is a definable range of entrainment to external driving (Moore-Ede et al., 1982). Outside of this entrainment range, more complicated behaviors are seen in the model (as well as in animals) with the period of the overt rhythm always averaging near tx0. For example, with signal periods slightly below 0.5tx0 the overt period becomes significantly less precise than at other signal periods; this precision is recovered for lower signal periods, but the shape of the overt rhythm changes qualitatively. This loss in precision of the overt rhythm may be similar to the phenomenon known as relative coordination (Moore-Ede et al., 1982), in which animal rhythms become unable to keep time with an external rhythm, yet are still influenced by it. On the other hand, in the rightmost inset for signal periods q1.5tx0 we see that the overt rhythm breaks into two alternating periods, one short and the next long. We know of no corollary reported to date for this phenomenon in animals.

466

.   . 

Thus for this form of coupling to a second signal, one can reliably vary the natural period of the overt rhythm by at most 50% in either direction. This range of entrainment for our ‘‘clock’’ is larger than observed in several animal species (Moore-Ede et al., 1982). Because the ability to phase shift determines the range of entrainment, it is interesting to note that the mammalian SCN in vitro consistently yields larger phase shifts in response to a variety of stimuli than intact animals (for example Prosser et al., 1989, 1990; Ding et al., 1994). Data like these suggest that the range of entrainment for isolated circadian clocks in vitro could be substantially larger than the range of entrainment observed in whole animals. However, testing this hypothesis is difficult because maintaining adult nervous tissue in vitro over the long term may be a formidable obstacle. Thus far, we have portrayed the model response to periodic signals as being the response of the long term neuronal firing rhythm to a second endogenous signal. Circadian rhythms are also known to be sensitive to external signals (e.g. light). Thus, one can alternatively interpret the results shown in Fig. 6 as the expected overt behavior in response to external signals. In this case, the external signals would presumably be filtered through one of the oscillatory systems: for example it appears likely that clock proteins are sensitive to light (Crosthwaite et al., 1995; Hunter-Ensor et al., 1996). In this scenario, the long term neuronal rhythm may respond to light only indirectly, after the light signal has been processed through the transcription/translation system. In addition, the abundance of some neurotransmitter peptides of the circadian clock are known to be modulated by light (Shinohara et al., 1993). However, it is difficult to predict the effect of such modulation on the long term neuronal rhythm modeled here. As increasingly sophisticated molecular techniques are used to analyse the detailed mechanisms behind the transcription/translation feedback loop, it may be possible in the future to model that system in tandem with the long term neuronal oscillation system discussed here. We anticipate that although Fig. 6 would not be reproducible in exact detail in animal

experiments, the qualitative changes indicated in the figure may be observed in vivo through careful experimentation. For example, the shape of the overt rhythm changes with increased signal period as indicated in the insets to Fig. 6. If we hypothesize that a rodent is inclined to be active whenever the envelope rhythm exceeds a given threshold, such changes in shape may have repercussions for the pattern of bouts of activity that are observed in rodent wheel-running records. As we have mentioned, at about t0 = 1.5tx0, entrainment is lost and the cycle breaks up into two qualitatively distinct periods, the first quite short and coinciding with the minimum in signal rhythm strength, and the second much longer and coinciding with the signal maximum. This is to some degree an artifact of the choice of parameters: as we have discussed, Fig. 6 is obtained using a comparatively large value of k, the phase coupling amplitude. This results in the previously mentioned asymmetry in the PRC, and has repercussions, as well, in asymmetries in the entrainment response of the model. Nevertheless, it is possible that some animals may exhibit these asymmetric behaviors, and it would be interesting to reinvestigate records from animals maintained under long T cycles.   The dynamics of the envelope rhythm become more complex when combined with periodic modulation of the coupling. As an example of this complexity, in Fig. 7 we plot the response to a periodic signal for various values of the mean coupling strength, k. The coupling strength is still modulated by the signal rhythm but in addition we consider the case where some additional stimulus also affects the mean k. This additional stimulus may be analogous to maintaining mammals (especially hamsters) in constant light. In the main plot, we show 16 successive periods (superimposed on error bars) of the overt rhythm after transients due to the choice of initial state have decayed away. As insets, we show several snapshots of the overt rhythm, and as before, we include the driving signal above each of these plots for phase comparison. The driving signal is chosen to have 4.3.

      period tx0 = 7/3 of the period* of the long term neuronal rhythm. At large values of mean k, the overt rhythm persists in synchrony with the driving rhythm, and the difference in phase between the maximum and the minimum of the overt rhythm diminishes with increasing mean k. For values of k between 1.5 and 3, the envelope shape of the overt rhythm changes, but the period of the rhythm remains constant to a high degree of accuracy. Below about k = 1.5, the shape of the overt rhythm becomes more complex, the mean period diminishes and the precision of the day-to-day t decreases dramati*The ratio 7/3 was chosen so that the overt and signal periods would be relatively prime, to prevent spurious resonances. Note that k ranges from 0 to 3 in Fig. 7, while in Fig. 6 k was fixed at k = 1. At large k, the overt rhythm is more strongly coupled to the signal rhythm, and consequently the maximum signal/overt period ratio is higher here than was possible in Fig. 6.

467

cally. At k of around 0.75, a precise ‘‘split’’ rhythm with period precisely half of the driving period emerges, and for smaller k the overt rhythm bifurcates into a period-2 cycle. These results lead to several predictions. First, it appears that the typical variability of the split rhythm (low k in Fig. 7) is considerably larger than the variability of the non-split rhythm (high k). The variability in the transition region (0.75 Q k Q 1.5) is even higher. Presumably this behavior should differ substantively from existing models (Pittendrigh & Daan, 1976; Daan & Berde, 1978; Kronauer et al., 1982) which have two 24 hour rhythms attracted to the ‘‘splay’’ state (Silber et al., 1993) in which each rhythm lies 12 hours out of phase from the other. Moreover, since the short-time oscillator component is conjectured to be due to neuronal firing rhythms, we would expect predictions specific to this model to fail in organisms in which the

Signal rhythm Overt rhythm

Period,

ˆ0

ˆ /2 0

0

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Mean coupling,

F. 7. Variation in ‘‘overt’’ period with changes in mean phase coupling strength, k. For strong coupling, the ‘‘overt’’ phase entrains exactly to the signal phase. As k decreases below k 3 1.5, the overt rhythm becomes dramatically more complex and less precise. As k drops below k 3 0.75, the overt rhythm splits to half the signal rhythm, and as k continues to decrease, this split rhythm period doubles into a sequence of short and long cycles. In the main plot, each open circle represents the ‘‘overt’’ t (2SE) during 16 successive iterations. tx0 is the initial period of the ‘‘overt’’ rhythm. Insets show relationship between signal rhythm (top of each inset) and overt rhythm (bottom).

468

.   . 

circadian cycle is controlled by another tissue such as the pineal gland. Second, the model suggests that raising vs. lowering mean k, the phase coupling constant, should have very different effects on the overt rhythm. For sufficiently strong external signals which lower mean k, the model exhibits ‘‘splitting’’ (where the overt period lies close to half the ‘‘natural’’ period) but this behavior is never seen for signals which raise k. Here again, this model provides an intriguing corollary to observations on circadian rhythms in hamsters. Splitting in hamsters occurs after a period of exposure to constant light but never in response to constant darkness (Pittendrigh, 1974; Earnest & Turek, 1982). In fact, constant darkness is known to cause the split locomotor rhythms to coalesce. The model also suggests that splitting may be associated with one type of stimuli (e.g. stimuli that produce a light-type PRC) but not in others (e.g. stimuli that produce a dark-type PRC) because these different types of PRC are produced depending upon whether k is raised or lowered (Fig. 5). To our knowledge, splitting has never been demonstrated in animals solely in response to a stimulus which has a dark-type PRC. A third prediction of the model that may be testable in the near future (as critical clock proteins are cloned) is that splitting first occurs in the long term neuronal rhythm we hypothesize to be the instrument directly driving most overt behavioral and metabolic rhythms. All elements downstream of the neuronal output of the SCN would be expected to be split but the putative transcription/ translation loop forming the molecular 24-hour clock need not split. In this model, the critical clock proteins making up the transcription/ translation feedback loop would exhibit a rhythm with t of approximately 24 hours, even in animals exhibiting split overt rhythms. Unfortunately, this hypothesis cannot be tested in the two systems, Drosophila and Neurospora, where clock proteins already have been identified. We are not aware of reports of splitting in either of these organisms, perhaps because constant light causes arhythmicity via its mediation of the destruction of tim in Drosophila (Hunter-Ensor et al., 1996) and via its effect

to constitutively raise frq in Neurospora (Crosthwaite et al., 1995). A fourth prediction derived from the model is that at external signal levels just short of the level needed to produce splitting, the envelope shape of the overt rhythm increases in complexity. As k is decreased, the mean overt period drops somewhat, and the variation in period from day to day increases dramatically (compare insets to Fig. 7 at k = 1.5 and at k = 1.0). To our knowledge, this has not been investigated before, and it would be particularly interesting to compare the complexities of animal activity records between animals which do and do not split under conditions of bright constant light. At still lower k, the model predicts that the split state of two 12 hour rhythms should bifurcate into two distinct rhythms totaling 24 hours. For example at k = 0.25, subjective ‘‘days’’ would alternate in length between about 11 and 13 hours. Although the dependence of this behavior on parameter values is not well understood at this point, and variability in the time of activity onset in animals may be great enough to partially mask this effect, it would be interesting to investigate whether this is observed in animal data. Finally, there appears from the model to be a subtle but real relation between the phase of the signal rhythm and the phase of the overt rhythm which seems to vary with changes in k. If we hypothesize that a stimulus like constant light decreases k directly, it may be possible to determine, for example, whether the phase relationship of the slow signal rhythm to the overt rhythm differs between animals kept under conditions of constant light and constant darkness.  In this section we studied dynamical changes in the model long term rhythm produced by coupled noisy high frequency oscillators. We plotted a Phase Response Curve, which exhibits asymmetric advances and delays (Fig. 5). The long term rhythm produced is Type 1 for all parameter variations studied. We studied the model’s response to periodic driving, as might occur as an interaction between a model neuronal rhythm and a second, e.g. 4.4.

      transcription/translation clock. We found that a limited range in period of this signal clock can entrain the model neuronal rhythm. Outside of this range, the coupled oscillators modeled here lose synchrony with the signal clock in a manner possibly similar to ‘‘relative coordination’’. We also studied the effect of changes in mean coupling strength between the model high frequency oscillators and how this affects their response to the signal rhythm hypothesized to come from a transcription/translation clock. Finally, we found that a phenomenon similar to ‘‘splitting’’ of the overt rhythm can be produced. This behavior is accompanied by a variety of rich dynamics, and we suggest that a search for these dynamics in animal experiments may prove worthwhile. 5. Conclusion The model that we presented exhibits several behaviors worth noting. Some of these behaviors—such as the transition to a split state, the sharp onset and slow offset of the overt rhythm, the independence of the overt rhythm from long term rate constants, and a precise long term collective rhythm in the presence of noisy short term oscillators—may have physiological relevance. Other aspects of the model—such as the existence of distinct neuronal and transcription/ translation rhythms, a specified variation in the shape of the overt rhythm with changes in entrainment period, an association between splitting and only one type of stimulus, and an increase of overt variability during splitting— seem ripe for experimental testing. Still other issues remain unresolved at this time. For example, we plan to investigate whether a ‘‘dead zone’’ in the PRC of the long term rhythm can be present for some parameter ranges in this model. We also plan to determine whether particular features seen in our simulations are generic or are limited to specific parameter ranges. It is our hope that future research involving collaborations between physiologists and dynamicists may better elucidate critical issues involved in basic biological questions. The authors wish to thank N. Kopell, S. H. Strogatz, P. M. Wise and members of the Turek and

469

Takahashi laboratories for helpful comments and criticisms.

REFERENCES A, A., S, A., G, A. & A, A. (1996). Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses. Science 273, 1868–1871. B, R. K. & T, J. S. (1995). Temperature compensation and temperature entrainment of the chick pineal cell circadian clock. J. Neurosci. 15, 5681–5692. B, R. K. & T, J. S. (1997). Lability of circadian pacemaker amplitude in chick pineal cells: a temperature-dependent process. J. Biol. Rhythms 12, 309–318. B, Y. & D, F. E. (1993). Neuronal synchronization without calcium-dependent synaptic transmission in the hypothalamus. Proc. Nat. Acad. Sci. U.S.A. 90, 3207–3210. C, C. A., R, B., S, D. F., G, E. & M, M. (1988). Vasoactive Intestinal Polypeptide neuron changes in the senile rat suprachiasmatic nucleus. Neurobiol. Aging 9, 307–312. C, J. R. & DH, R. L. (1979). Fluctuations in interbeat interval in rhythmic heart-cell clusters. Biophys. J. 28, 377–390. C, A. H., H, P. J. & R, R. H. (1982). The nature of the coupling between segmental oscillators of the lamprey spinal generator for locomotion: a mathematical model. J. Math. Biol. 13, 345–369. C, J. J., I, T. T. & G, P. (1996). Noiseenhance information transmission in rat SA1 cutaneous mechanoreceptors via aperiodic stochastic resonance. J. Neurophys. 76, 642–645. C, S. K., L, J. J. & D, J. C. (1995). Light-induced resetting of a circadian clock is mediated by a rapid increase in frequency transcript. Cell 81, 1003–1012. D, S. & B, C. (1978). Two coupled oscillators: simulations of the circadian pacemaker in mammalian activity rhythms. J. theor. Biol. 70, 297–313. D, H. (1993). A solvable model of coupled limit-cycle oscillators exhibiting partial perfect synchrony and novel frequency spectra. Physica D 69, 394–403. DC, P. J. (1960). Phase control of activity in a rodent. In: Cold Spring Harbor Symp. Quant. Biol. Cold Spring Harbor, Long Island, NY: The Biological Laboratory. pp. 49–55. D, J. M., C, D., W, E. T., F, L. E., R, M. A. & G, M. U. (1994). Resetting the biological clock: mediation of nocturnal circadian shifts by glutamate and NO. Science 266, 1713–1717. D, J. K., W, L., P, E. & M, F. (1993). Noise enhancement of information transfer in crayfish mechanoreceptors by stochastic resonance. Nature 365, 337–340. D, J. & H, W. (1987). Intracellular recording in the medullary pacemaker nucleus of the weakly electric fish, Apteronotus, during modulatory behaviors. J. Comp. Physiol. A 161, 187–200. E, D. J. & T, F. W. (1982). Splitting of the circadian rhythm of activity in hamsters: effects of

470

.   . 

exposure to constant darkness and subsequent reexposure to constant light. J. Comp. Physiol. 145, 405–411. E, D., DG, S. & S, C. (1991). Effects of tetrodotoxin on the circadian pacemaker mechanism in suprachiasmatic explants in vitro. Brain Res. Bull. 26, 677–682. E, J. T. (1980a). The Timing of Sleep and Wakefulness New York: Springer. E, J. T. (1980b). Temporal precision in circadian systems: a reliable neuronal clock from unreliable components? Science 209, 1542–1543. E, G. B. (1985). Synchronization in a pool of mutually coupled oscillators with random frequencies. J. Math. Biol. 22, 1–9. E, G. B. (1991). An adaptive model for synchrony in the firefly Pteroptyx malaccae. J. Math. Biol. 29, 571–585. E, G. B. & K, N. (1984). Frequency plateau in a chain of weakly coupled oscillators, I. SIAM J. Math. Anal. 15, 215–237. G, L. & M, M. C. (1988). From Clocks to Chaos: The Rhythms of Life. Princeton, NJ: Princeton University Press. G, A. (1995). A model for circadian oscillations in the Drosophila period protein (PER). Proc. R. Soc. Lond. B 261, 319–324. G, D. & G, R. (1982). Circadian rhythm of firing rate recorded from single cells in the rat suprachiasmatic brain slice. Brain Res. 245, 198–200. H, J., L, Y. X., M, J. L., W, B. & G, A. (1990). Coupling chaotic and periodic cells results in a period-doubling route to chaos in a model of cAMP oscillators in Dichtyostelium suspension. Phys. Lett. A 151, 33–36. H-E, M., O, A., & S, A. (1996). Regulation of the Drosophila protein timeless suggests a mechanism for resetting the circadian clock. Cell 84, 677–685. K, S. & M, A. (1995). New models in vogue for circadian clocks. Cell 83, 361–364. K, D. C., M, R. Y. & R, S. M. (1991). Suprachiasmatic Nucleus—The Mind’s Clock. New York: Oxford University Press. K, R. R. & B, J. L. (1993). Dynamic analysis of period mutants in Drosophila: how a precise circadian clock might emerge from a tissue composed of chaotic cellular oscillators. In: Molecular Genetics of Biological Rhythms, (Young, M. W., ed.) pp. 221–253. New York: Marcel Dekker. K, E. (1979). Advanced Engineering Mathematics p. 907. New York: Wiley. K, R. E., C, C. A., P, S. F., M-E, M. C. & W, E. D. (1982). Mathematical model of the human circadian system with two interacting oscillators. Am. J. Physiol. 242, R3–17. K, Y. (1975). Self-entrainment of a population of coupled nonlinear oscillators. In: Lecture Notes in Physics, 39th International Symposium on Mathematical Problems in Theoretical Physics (Araki, H., ed.) pp. 420–422. New York: Springer. K, Y. (1984). Chemical Oscillations, Waves, and Turbulence pp. 68–77. Berlin: Springer. L, L. D. & L, E. M. (1959). Fluid Mechanics p. 121. Reading, MA: Addison-Wesley.

L, J. E. & M, J. P. (1996). Broadband neural encoding in the cricket cercal sensory system enhanced by stochastic resonance. Nature 380, 165–168. L, A., B, A. & M, F. (1991). Time-interval sequences in bistable systems and the noise-induced transmission of information by sensory neurons. Phys. Rev. Lett. 67, 656–659. L, J. (1995). The molecular basis of the Neurospora clock. Sem. Neurosci. 7, 3–13. M, Z. F. & S, T. J. (1995). Reliability of spike timing in neocortical neurons. Science 268, 1503–1506. M, P. C. & S, S. H. (1990). Phase diagram for the collective behavior of limit-cycle oscillators. Phys. Rev. Lett. 65, 1701–1704. M, P. C., M, R. E. & S, S. H. (1991). Dynamics of a large system of coupled nonlinear oscillators. Physica D 52, 293–331. M, S., G, M. E., Z, J. J. & B, G. D. (1993). Circadian rhythm in membrane conductance expressed in isolated neurons. Science 259, 239–241. M-E, M. C., S, F. M. & F, C. A. (1982). The Clocks That Times Us. Cambridge: Harvard University Press. P, T. (1969). Populations of interacting oscillators and circadian rhythms. J. theor. Biol. 22, 418–436. P, M. & J, A. (1994). A study of the singularities in a mathematical model for circadian rhythms. Biosystems 33, 193–201. P, C. S. (1954). On temperature independence in the clock system controlling emergence in Drosophila. Proc. Nat. Acad. Sci. U.S.A. 40, 2697–2701. P, C. S. (1960). Circadian rhythms and the circadian organization of living systems. Cold Spring Harbor Symp. Quant. Biol. 25, 159–184. P, C. S. (1965). On the mechanism of the entrainment of a circadian rhythm by light cycles. In: Circadian Clocks (Aschoff, J., ed.) pp. 277–297. Amsterdam: North Holland. P, C. S. (1974). Circadian oscillations in cells and the circadian organization of multicellular systems. In: The Neurosciences Third Study Program (Schmitt, F. O. & Worden, F. G., eds.) pp. 437–458. Cambridge, MA: MIT Press. P, C. S. & D, S. (1976). A functional analysis of circadian pacemakers in nocturnal rodents. V. Pacemaker structure: a clock for all seasons. J. Comp. Physiol. 106, 333–355. P, R. A. & G, M. U. (1989). The mammalian circadian clock in the suprachiasmatic nuclei is reset in vitro by cAMP. J. Neurosci. 9, 1073–1081. P, R. A., M, J. D. & H, H. C. (1990). A serotonin agonist phase-shifts the circadian clock in the suprachiasmatic nuclei in vitro. Brain Res. 534, 336–339. R, B., V G, W. A., S, D. F., H, J. E. & M, M. (1987). Changes in vasopressin cells of the rat suprachiasmatic nucleus with aging. Brain Res 409, 259–264. R, M. (1995). Molecular control of circadian rhythms. Curr. Opin. Gen. & Dev. 5, 662–668. S, K., L-O, S., W, E. P. & T, F. W. (1997). Aging and photoperiod affect entrainment and quantitative aspects of locomotor behavior in Syrian hamsters. Am. J. Physiol. 272, R1219–R1225.

      S, W., G, R. A. & M, M. T. (1987). The suprachiasmatic nuclei contain a tetrodotoxin-resistant circadian pacemaker. Proc. Nat. Acad. Sci. U.S.A. 84, 1694–1698. S, A. & R, J. (1991). Model for synchronization of pancreatic b-cells by gap junction coupling. Biophys. J. 59, 547–559. S, A. & R, J. (1992). Rhythmogenic effects of weak electrotonic coupling in neuronal models. Proc. Nat. Acad. Sci. U.S.A. 89, 1–4. S, K., T, K., I, Y. & I, S.-I. T. (1993). Photic regulation of peptides located in the ventrolateral subdivision of the suprachiasmatic nucleus of the rat: daily variations of vasoactive intestinal peptide, gastrin-releasing peptide, and neuropeptide Y. J. Neurosci. 13, 793–800. S, M., F, L. & W, K. (1993). Stability results for in-phase and splay-states of solid-state laser arrays. J. Opt. Soc. Am. B 10, 1121–1129. S, R., L, M. N., G, M., G, W. R. & B, E. L. (1990). Dispersed cell suspensions of fetal SCN restore circadian rhythmicity in SCN-lesioned adult hamsters. Brain Res. 525, 45–58. S, R. D., T, F. W. & T, J. S. (1992). Two families of phase-response curves characterize the resetting of the hamster circadian clock. Am. J. Physiol. 262, R1149–R1153. S, D. & K, N. (1993). Rapid synchronization through fast threshold modulation. Biol. Cybern. 68, 393–407. S, A. & H¨, K. (1990). Sweet food means hot dancing in honeybees. Naturwissenschaften 78, 471–473. S, S. H. (1994). Nonlinear Dynamics and Chaos Reading, MA: Addison-Wesley.

471

T, J. S. & K, J. M. (1993). Molecular approaches to understanding circadian oscillations. Annu. Rev. Physiol. 55, 729–753. T, G. & M, M. (1996). Circadian rhythms in cultured mammalian retina. Science 272, 419–421. T, R. D., M, R. & W, R. K. S. (1989). Model of the origin of rhythmic population oscillations in the hippocampal slice. Science 243, 1319–1324.   P, A. N. & D, F. E. (1993). Cellular communication in the circadian clock, the suprachiasmatic nucleus. Neuroscience 56, 793–811.   P, A. N. (1980). The hypothalamic suprachiasmatic nucleus of rat: intrinsic anatomy. J. Comparative Neurol. 191, 661–702.   P, B. (1926). Relaxation-oscillations. Phil. Mag. 7, 978–992. W, G. V. (1993). Spatial, temporal, and global mode entropy in a thalamo-cortical network. Int. J. Bif. Chaos 3, 1487–1501. W, D., L, D., M, M. & R, S. (1995). Individual neurons dissociated from rat suprachiasmatic nucleus express independently phased circadian firing rhythms. Neuron 14, 697–706. W, R. (1984). Toward a mathematical model of circadian rhythmicity. In: Mathematical Models of the Circadian Sleep–Wake Cycle (Moore-Ede, M. C. & Czeisler, C. A., eds) pp. 17–79. New York: Raven. W, H. R. & C, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. W, A. T. (1967). Biological rhythms and the behavior of populations of coupled oscillators. J. theor. Biol. 16, 15–42. W, A. T. (1990). The Geometry of Biological Time. Berlin: Springer.