A statistical model of the human core-temperature circadian rhythm

Emery N. Brown1, Yong Choe1, Harry Luithardt1, and Charles A. Czeisler2

1 Statistics Research Laboratory, Department of Anesthesia and Critical Care, Massachusetts General Hospital, and Division of Health Sciences and Technology, Harvard Medical School-Massachusetts Institute of Technology, Boston 02114-2696; and 2 Circadian, Neuroendocrine and Sleep Disorders Section, Division of Endocrinology, Department of Medicine, Brigham and Women's Hospital, Harvard Medical School, Boston, Massachusetts 02115


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

We formulate a statistical model of the human core-temperature circadian rhythm in which the circadian signal is modeled as a van der Pol oscillator, the thermoregulatory response is represented as a first-order autoregressive process, and the evoked effect of activity is modeled with a function specific for each circadian protocol. The new model directly links differential equation-based simulation models and harmonic regression analysis methods and permits statistical analysis of both static and dynamical properties of the circadian pacemaker from experimental data. We estimate the model parameters by using numerically efficient maximum likelihood algorithms and analyze human core-temperature data from forced desynchrony, free-run, and constant-routine protocols. By representing explicitly the dynamical effects of ambient light input to the human circadian pacemaker, the new model can estimate with high precision the correct intrinsic period of this oscillator (~24 h) from both free-run and forced desynchrony studies. Although the van der Pol model approximates well the dynamical features of the circadian pacemaker, the optimal dynamical model of the human biological clock may have a harmonic structure different from that of the van der Pol oscillator.

biological clock; dynamical system; harmonic regression; perturbation expansion; thermoregulation; van der Pol equation


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

CIRCADIAN RHYTHMS ARE BIOLOGICAL rhythms generated by an organism or group of organisms that have an intrinsic period of ~24 h. The term circadian was coined by Franz Halberg from the Latin circa, meaning about, and dies, meaning day (40). In humans, the site of the circadian pacemaker or biological clock is the suprachiasmatic nucleus of the hypothalamus (39, 47). This ~24-h oscillator helps ensure correct timing among the body's physiological processes and between those processes and events in the outside world. The human circadian pacemaker is studied by measuring marker rhythms whose behaviors are known to be tightly coupled to the clock. The three principal marker rhythms used in human circadian studies are core temperature, plasma cortisol levels, and plasma melatonin levels. Of these three markers, core temperature and melatonin are believed to be the most reliable (18).

An advantage of using core temperature instead of melatonin to study the human circadian system is that it can be monitored continuously without the taking of blood samples. A drawback to the use of core temperature is that other physiological factors, such as posture, level of activity, and the dynamics of the body's thermoregulatory system, contribute significantly to its observed structure. Despite these potential confounding factors, core temperature is the most widely used marker rhythm in studies of the human circadian pacemaker (16, 19, 20, 44).

Mathematical models developed for analyses of core-temperature rhythms play an important role in research on the human circadian oscillator (5, 7-10, 12, 21-24, 26, 27, 29, 30-35, 37, 38, 48-50). These models divide into two categories: those designed to study the dynamical properties of the pacemaker and using differential equation-based simulations, and those designed to study static properties of the experimental data with statistical models. The simulation models are developed by specifying a set of dynamical behaviors that the model should describe, formulating a differential equation or family of equations capable of describing those behaviors, and determining the parameters for the model by comparing simulation findings with experimental observations. In these studies the core-temperature series is treated as the direct output of the oscillator. Once the model parameters have been determined, the model is used to test how well it represents the initially identified behaviors or to predict the behavior of the pacemaker under conditions not observed in the original experiments.

An important dynamical property of the human circadian system to describe is its limit cycle behavior. This is the ability of the circadian pacemaker to maintain a stable oscillation without an external stimulus and to return to this stable oscillation after perturbation (32). The limit cycle behavior is essential for describing fundamental properties of the pacemaker, such as its interaction with the sleep-wake cycle and its response to different artificial and natural light conditions (32-34). The limit cycle properties of the human circadian pacemaker are consistent with a weakly nonlinear oscillator (34). The van der Pol oscillator is the most commonly used weakly nonlinear differential equation model for core-temperature analyses, and it is defined as
<A><AC>s</AC><AC>¨</AC></A>(t)−&egr; <FR><NU>2&pgr;</NU><DE>&tgr;</DE></FR> <FENCE>1−<FR><NU>4</NU><DE>&ggr;<SUP>2</SUP></DE></FR> s<SUP>2</SUP>(t)</FENCE><A><AC>s</AC><AC>˙</AC></A>(t)+<FENCE><FR><NU>2&pgr;</NU><DE>&tgr;</DE></FR></FENCE><SUP>2</SUP>s(t)=0 (1)
where s(t) is the signal from the circadian pacemaker, tau  is the approximate period of the oscillation, gamma  is the approximate limiting amplitude, and epsilon  is the internal stiffness parameter. The internal stiffness parameter is a dimensionless quantity that governs the rate at which the solution approaches its limit cycle. The smaller the value of epsilon , the slower the approach of the oscillator to its limit cycle and the more sinusoidal its oscillations become. Weakly nonlinear means 0 < epsilon  1. If epsilon  = 0, the model in Eq. 1 becomes a sine wave or simple harmonic oscillator. The van der Pol model has a long history of application in both the physical and biological sciences (42). Because it is not precisely known how the human body produces the circadian oscillations in core temperature, the terms in the van der Pol equation cannot be given any further anatomic or physiological interpretations beyond those in the definitions of tau , gamma , and epsilon . Hence, the van der Pol equation is used primarily as a model of the circadian pacemaker, because it is a simple, analytically tractable system capable of representing the self-sustaining, weakly nonlinear oscillations characteristic of the human biological clock. Many authors have used various forms of the van der Pol oscillator in simulation studies of the human circadian pacemaker (23, 24, 29-35, 48, 50).

Statistical models of the human core-temperature rhythm have used harmonic regression methods because of the obvious sinusoidal structure of these data series. The rhythm is assumed to have a stable oscillation during the study interval, and the features of the rhythm typically characterized are static properties such as the rhythm's phases, its amplitude, and its period (8, 26, 49). The observed core-temperature rhythm also has correlated noise, which represents primarily the normal homeostatic response of the body's thermoregulatory system as well as effects due to a person's rest-activity cycle (8). The inclusion of correlated noise and the rest-activity components improve significantly goodness-of-fit and parameter standard error estimates. Nonharmonic regression techniques for core-temperature analysis include periodogram methods (21, 22), minimum variance methods (15), and cross-correlation methods (37, 38). The first two methods have been used primarily to estimate period, whereas the latter has been used to assess the magnitude of phase shifts after an intervention.

Analyses based on differential equation models have provided important insight into the dynamical properties of the human circadian pacemaker. With a few exceptions (5, 7, 9, 10, 27), formal statistical methods have received only limited use in these analyses. Therefore, the uncertainty in the inferences based on these models and their sensitivity to model specification and parameter error cannot be evaluated. In addition, the differential equation models generally make no attempt to account for other factors known to affect core-temperature measurements, such as thermoregulation and the rest-activity cycle. On the other hand, statistical models have been used solely to assess the static rather than the dynamical properties of the circadian pacemaker. The ideal model for core-temperature analysis would use in a formal statistical framework a dynamic representation of the circadian pacemaker based on a differential equation model (5, 7, 9, 10).

We present a statistical model of the human core-temperature circadian rhythm in which the circadian signal is represented as a van der Pol oscillator, the thermoregulatory process is modeled as a continuous first-order autoregressive process, and the effect of activity on the observed temperature rhythm is modeled with a protocol-specific function. We apply the model to the analysis of human core-temperature data from forced desynchrony, free-run, and constant-routine protocols.


    PROTOCOLS, MODEL FORMULATION, AND PARAMETER ESTIMATION
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Free-run protocol. Many of the early data on the human circadian system were collected on subjects studied in isolated environments free of time cues. Under these conditions, the subject's circadian pacemaker was believed to oscillate at its intrinsic period of ~25 h (50). This condition is called free-running. The free-run protocol has been, until recently, the "gold standard" for assessing the intrinsic properties of the circadian system. It is now appreciated that light, even at the level of ordinary room lighting, can shift the circadian pacemaker and that the period estimated under free-run conditions is not the intrinsic period of the human circadian pacemaker (18, 20). This is because the human circadian pacemaker, like those of lower animals, has a well-defined phase response curve (PRC) to light (20). Under free-run conditions, a subject self-selects his/her light-dark cycle and, hence, the phase at which the light is self-administered. Because of the asymmetric structure of the human PRC, this self-selection results more frequently in light administered during phases that favor delays of the human pacemaker, rather than advances, and therefore in an observed period that is longer than the pacemaker's intrinsic period.

Forced desynchrony protocol. Under the forced desynchrony protocol, a subject is monitored for an extended time in an isolated environment whose light-dark cycle is maintained outside the interval between 23 and 27 h (18). This 4-h interval, termed the range of entrainment, defines the set of light-dark cycle periods to which the human circadian pacemaker may be synchronized (50). During two-thirds of the light-dark cycle, the lights are maintained continuously on at a fixed intensity, whereas during the remaining one-third of the cycle, the subject is maintained in total darkness. In the recently designed forced desynchrony studies of Czeisler et al. (18), a light-dark cycle of either 20 or 28 h is typically chosen, the light intensity level is set at 20 lux or lower, and the behavior of the clock is monitored by recording marker rhythms such as core-temperature, plasma cortisol, and plasma melatonin levels (18). Because the light levels during forced desynchrony are low, and because the clock cannot synchronize to a light-dark cycle whose period is outside the range of entrainment, the circadian pacemaker oscillates at its intrinsic period.

Constant-routine protocol. The constant-routine protocol controls environmental conditions and a subject's level of activity to identify accurately the circadian component of an observed biological rhythm (16, 17, 19, 20, 38). This protocol was developed as an alternative to longer free-run studies, as a means of assessing the properties of an individual's circadian system in a shorter study interval. The protocol calls for subjects to remain continuously awake for 30-60 h in a semirecumbent posture exposed to constant indoor light (~150 lux) and to consume their daily caloric intake in evenly divided snacks at approximate hourly intervals (19, 20). Physiological and behavioral circadian variables are recorded, and the properties of the circadian pacemaker are studied by characterizing the phases and amplitudes of these marker rhythms such as core temperature, plasma melatonin levels, and plasma cortisol levels (16, 20). An advantage of this protocol is that it minimizes the effects of environmental conditions and level of activity on a subject's observed circadian rhythm. From a subject's constant-routine core-temperature data, it is possible to estimate the amplitude and phase of his/her circadian pacemaker. The pacemaker period cannot be determined, because, at best, only 1.5-2.5 cycles of the oscillation are observed on this protocol.

The core-temperature model. We assume that for a given circadian protocol, core-temperature data yt1,...,ytn are collected in the interval [0,T], where tn = nDelta t, n = 1,2,...,N, and NDelta t = T. For each protocol, the core-temperature measurement may be expressed as
y<SUB>t<SUB>n</SUB></SUB>=&mgr;+s<SUB>t<SUB>n</SUB></SUB>+x<SUB>t<SUB>n</SUB></SUB>+v<SUB>t<SUB>n</SUB></SUB>, (2)
where µ is the core-temperature mean, stn is the circadian oscillation, xtnis the evoked effect on core-temperature of the subject's activity level, and vtn is the fluctuation in core-temperature measurements due to the body's thermoregulatory response. The variable stn is represented as the solution to the modified van der Pol equation, defined in Ref. 31 as
  <FR><NU>d<IT>s</IT>(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR><IT>=z</IT>(<IT>t</IT>)<IT>+&egr;</IT><FENCE><FR><NU><IT>2&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE><FENCE><IT>s</IT>(<IT>t</IT>)<IT>−</IT><FR><NU><IT>4s</IT>(<IT>t</IT>)<SUP><IT>3</IT></SUP></NU><DE><IT>3&ggr;<SUP>2</SUP></IT></DE></FR></FENCE>

+<FENCE><FR><NU>2&pgr;</NU><DE>&tgr;</DE></FR></FENCE>&ggr;(1−ms(t))CI(t)<SUP>1/3</SUP> (3)

<FR><NU>d<IT>z</IT>(<IT>t</IT>)</NU><DE>d<IT>t</IT></DE></FR><IT>=</IT>−<FENCE><FR><NU><IT>2&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE><SUP><IT>2</IT></SUP><IT>s</IT>(<IT>t</IT>)<IT>+</IT><FR><NU><IT>2&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR><IT> z</IT>(<IT>t</IT>)<FENCE><FR><NU><IT>1−ms</IT>(<IT>t</IT>)</NU><DE><IT>3</IT></DE></FR></FENCE><IT>CI</IT>(<IT>t</IT>)<SUP><IT>1/3</IT></SUP>
where I(t) is the physical intensity of the ambient light at time t, m is the circadian light modulation index, C is a constant of proportionality, and the parameters epsilon , gamma , and tau  are as defined in Eq. 1. Equation 3 makes explicit the direct effect of ambient light on the circadian oscillator. By setting I(t) = 0 and applying the Liénard transformation (45), the van der Pol model in Eq. 3 is equivalent to Eq. 1. For the harmonic regression model, stn is represented as
s<SUB>t<SUB>n</SUB></SUB>=<LIM><OP>∑</OP><LL>r=1</LL><UL>d</UL></LIM> A<SUB>r</SUB> cos <FENCE><FR><NU><IT>2&pgr;rt<SUB>n</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE><IT>+B<SUB>r</SUB> </IT>sin <FENCE><FR><NU><IT>2&pgr;rt<SUB>n</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE> (4)
where d, the number of harmonics, is either 2 or 3. The choice of d = 2 follows from the harmonic regression analysis of core-temperature data under the constant-routine protocol (8). In the next section we show that the choice of d = 3 makes the harmonic regression model equivalent to the asymptotic series representation of the van der Pol oscillator.

The form of xtn depends on the circadian protocol. For the forced desynchrony protocol with a 28-h light-dark cycle, the regular, square-wave shape of xtn is well described by the harmonic regression model
x<SUB>t<SUB>n</SUB></SUB>=<LIM><OP>∑</OP><LL>r=1</LL><UL>4</UL></LIM> C<SUB>r</SUB> cos <FENCE><FR><NU><IT>2&pgr;rt<SUB>n</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE><IT>+D<SUB>r</SUB> </IT>sin <FENCE><FR><NU><IT>2&pgr;rt<SUB>n</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE> (5)
(12, 18). For the self-selected timing of activity during the free-run protocol, xtn is
x<SUB>t<SUB>n</SUB></SUB>=<FENCE><AR><R><C>&bgr;X(t<SUB>n</SUB>)</C><C>where<IT> X</IT>(<IT>t<SUB>n</SUB></IT>)<IT>=1, </IT>if<IT> I</IT>(<IT>t<SUB>n</SUB></IT>)<IT>>0</IT></C></R><R><C><IT>0</IT></C><C>   otherwise</C></R></AR></FENCE> (6)
This term models the effect of activity during free run as an increase in mean core temperature by an amount beta  for beta  > 0. Under the constant-routine protocol, the subject's activity is kept to a minimum so that xtn = 0 for all tn.

The random variable vtn is a discrete sample from a continuous first-order autoregressive [AR(1)] process and is defined as
v<SUB>t<SUB>n</SUB></SUB>=e<SUP>−&agr;&Dgr;t</SUP>v<SUB>t<SUB>n−1</SUB></SUB>+&eegr;<SUB>t<SUB>n</SUB></SUB> (7)
where alpha -1 is the time constant for the thermoregulatory response, and the eta tn are independent, Gaussian random variables with zero mean and variance sigma eta 2. Thermoregulation is controlled primarily by the anterior, posterior, and preoptic nuclei in the hypothalamus (28). A current hypothesis regarding the relation between the biological clock and the thermoregulatory system states that the circadian pacemaker creates a dynamic set point that the thermoregulatory system tracks (25). Our present model does not consider the interaction between the circadian pacemaker and the thermoregulatory system.

Relation between the van der Pol and harmonic regression models. The relation between the van der Pol model with 0 < epsilon   1 and the harmonic regression models can be made explicit by expanding Eq. 1 in an asymptotic series in powers of epsilon  (41). The nonlimit cycle perturbation solution of Eq. 1 to O(epsilon 2) is
s(t)=&ggr;(t) cos <FENCE><FR><NU><IT>2&pgr;t</IT></NU><DE><IT>&tgr;</IT></DE></FR><IT>+&psgr;<SUB>0</SUB></IT></FENCE>

<IT>−&egr; </IT><FR><NU><IT>&ggr;</IT>(<IT>t</IT>)<SUP><IT>3</IT></SUP></NU><DE><IT>8&ggr;<SUP>2</SUP></IT></DE></FR> sin <FENCE><FR><NU>6&pgr;t</NU><DE>&tgr;</DE></FR>+3&psgr;<SUB>0</SUB></FENCE>+O(&egr;<SUP>2</SUP>) (8)
where
&ggr;(t)=<FR><NU>2C<SUP>1/2</SUP><SUB>0</SUB> exp<FENCE><FR><NU><IT>&egr;&pgr;t</IT></NU><DE><IT>2&tgr;</IT></DE></FR></FENCE></NU><DE><FENCE><IT>1+</IT><FR><NU><IT>4C<SUB>0</SUB></IT></NU><DE><IT>&ggr;<SUP>2</SUP></IT></DE></FR> exp<FENCE><FR><NU><IT>&egr;&pgr;t</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></FENCE><SUP><IT>1/2</IT></SUP></DE></FR> (9)
and C0 is a constant of integration. Equations 8 and 9 are derived in APPENDIX A. Because gamma  = tright-arrow infinity limgamma (t), the asymptotic or limit cycle solution of Eq. 8 to O(epsilon 2) is
s(t)=&ggr; cos <FENCE><FR><NU><IT>2&pgr;t</IT></NU><DE><IT>&tgr;</IT></DE></FR><IT>+&psgr;<SUB>0</SUB></IT></FENCE><IT>−&egr; </IT><FR><NU><IT>&ggr;</IT></NU><DE><IT>8</IT></DE></FR> sin <FENCE><FR><NU><IT>6&pgr;t</IT></NU><DE><IT>&tgr;</IT></DE></FR><IT>+3&psgr;<SUB>0</SUB></IT></FENCE><IT>+O</IT>(<IT>&egr;<SUP>2</SUP></IT>) (10)
The asymptotic solution of the van der Pol equation when the oscillator is close to or on its limit cycle is equivalent to a harmonic regression model comprised of the first two odd harmonics of the frequency 2pi /tau . The asymptotic expansions can be carried out in principle to any desired order; however, second or third order suffices for most problems. Equation 10 suggests that if the human circadian pacemaker truly behaves like a van der Pol oscillator, then use of a three-harmonic regression model would best describe the stable oscillating conditions of the constant-routine protocol with low ambient light conditions.

Model estimation, parameter standard errors, and goodness-of fit. The number of model parameters to be estimated from the core-temperature data is large. For example, for the forced desynchrony model it is 18, whereas it is 11 for the free-run model. It is imperative, therefore, to implement efficient algorithms to fit the models to the core-temperature data. To do this, we have developed a Newton's method maximum likelihood algorithm based on Corradi's theorem of separable nonlinear optimization (8, 14). As we show below, we separate the maximization of the likelihood into the linear and nonlinear parameters. Given values of the nonlinear parameters, the optimal linear parameters can be computed explicitly as weighted least squares estimates in terms of the nonlinear parameter values. Corradi's theorem shows that this algorithm leads to the same maximum of the likelihood that would be obtained if the linear parameters were treated the same as the nonlinear parameters. The algorithm is computationally efficient, because the dimension of the parameter vector in the Newton's step corresponds only to the number of nonlinear parameters in the model.

Under the Gaussian assumption for the vtn, the log likelihood for each model has the form
L=−<FR><NU><IT>1</IT></NU><DE><IT>2</IT></DE></FR> log<IT> ‖&Ggr;</IT>(<IT>&agr;</IT>)<IT>‖−</IT><FR><NU><IT>N</IT></NU><DE><IT>2</IT></DE></FR> log (<IT>&sfgr;</IT><SUP><IT>2</IT></SUP><SUB><IT>&eegr;</IT></SUB>)<IT>−</IT><FR><NU><IT>1</IT></NU><DE><IT>2</IT></DE></FR> <FR><NU><IT>S<SUB>N</SUB></IT></NU><DE><IT>&sfgr;</IT><SUP><IT>2</IT></SUP><SUB><IT>&eegr;</IT></SUB></DE></FR> (11)
where constants not depending on the model parameters or the data have been omitted,  Gamma (alpha )  is the determinant of the covariance matrix of the AR(1) process, and SN is a quadratic form that depends on the data, the model, and the circadian protocol. Because <A><AC>&sfgr;</AC><AC>ˆ</AC></A>eta 2 = SN/N is the maximum likelihood estimate of sigma eta 2, the concentrated log likelihoods are
L′=−<FR><NU><IT>1</IT></NU><DE><IT>2</IT></DE></FR> log<IT> ‖&Ggr;</IT>(<IT>&agr;</IT>)<IT>‖−</IT><FR><NU><IT>N</IT></NU><DE><IT>2</IT></DE></FR> log (<IT>S<SUB>N</SUB></IT>) (12)
Let y = (yt1,...,ytn), and for the van der Pol model let theta  = (st1, ·t1,tau ,gamma ,epsilon ,m,C), and s(theta ) = [st1(theta ), ...,stn(theta )]. For the free-run protocol, define the regression coefficient xi FR = (µ,beta ) and the N × 2 design matrix
Z<SUB>FR</SUB>=<FENCE><AR><R><C>1</C><C>X(t<SUB>1</SUB>)</C></R><R><C>&vtdot;</C><C>&vtdot;</C></R><R><C>1</C><C>X(t<SUB>N</SUB>)</C></R></AR></FENCE> (13)
where X(t) is given in Eq. 6. For the forced desynchrony protocol, define the regression coefficient xi FD = (µ,C1,D1,...,C4,D4), and the N × 9 design matrix


Z<SUB>FD</SUB>=<FENCE><AR><R><C>1</C><C>cos <FENCE><FR><NU><IT>2&pgr;t<SUB>1</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C><C>sin <FENCE><FR><NU><IT>2&pgr;t<SUB>1</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C><C><IT>…</IT></C><C>cos <FENCE><FR><NU><IT>8&pgr;t<SUB>1</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C><C>sin <FENCE><FR><NU><IT>8&pgr;t<SUB>1</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C></R><R><C><IT>&vtdot;</IT></C><C></C><C></C><C></C><C></C><C><IT>&vtdot;</IT></C></R><R><C><IT>1</IT></C><C>cos <FENCE><FR><NU><IT>2&pgr;t<SUB>N</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C><C>sin <FENCE><FR><NU><IT>2&pgr;t<SUB>N</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C><C><IT>…</IT></C><C>cos <FENCE><FR><NU><IT>8&pgr;t<SUB>N</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C><C>sin <FENCE><FR><NU><IT>8&pgr;t<SUB>N</SUB></IT></NU><DE><IT>28</IT></DE></FR></FENCE></C></R></AR></FENCE> (14)

Then SN has the form
S<SUB>N</SUB>=(y−s(&thgr;)−Z&xgr;)<SUP>t</SUP>&Ggr;(&agr;)<SUP>−1</SUP>(y−s(&thgr;)−Z&xgr;) (15)
where xi  is either xi FR or xi FD, and Z is either ZFR or ZFD, depending on whether the protocol is free run or forced desynchrony, respectively. For a given value of theta , s(theta ) can be computed by a Runge-Kutta algorithm (1), whereas for a given value of alpha , Gamma (alpha )-1 and  Gamma (alpha )  are computed using the Kalman filter (11). The Kalman filter algorithm exploits the Markov structure in the correlated noise model to invert Gamma (alpha ) and compute its determinant efficiently. An important advantage of this algorithm is that it can also be used when observations are missing. Given values of theta  and alpha , the maximum likelihood estimate of xi  is
<A><AC>&xgr;</AC><AC>ˆ</AC></A>(&thgr;, &agr;)=[Z<SUP>t</SUP>&Ggr;(&agr;)<SUP>−1</SUP>Z]<SUP>−1</SUP>Z<SUP>t</SUP>&Ggr;(&agr;)<SUP>−1</SUP>[y−s(&thgr;)] (16)
where ^ denotes the estimate. Therefore, the maximum likelihood estimates of theta  and alpha  can be computed using Newton's procedure to maximize Eq. 12 with respect to theta  and alpha  at each iteration and computing <A><AC>&xgr;</AC><AC>ˆ</AC></A>(theta ,alpha ) from Eq. 16. The dimension of the parameter vector in the Newton's method is 8, i.e., the dimension of theta  plus alpha .

For the harmonic regression model, we define for the free-run protocol the regression coefficient xi HR = (A1,B1,A2,B2)t and the N × 4 design matrix
Z<SUB>HR</SUB>=<FENCE><AR><R><C>cos <FENCE><FR><NU><IT>2&pgr;t<SUB>1</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C><C>sin <FENCE><FR><NU><IT>2&pgr;t<SUB>1</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C><C>cos <FENCE><FR><NU><IT>4&pgr;t<SUB>1</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C><C>cos <FENCE><FR><NU><IT>4&pgr;t<SUB>1</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C></R><R><C><IT>&vtdot;</IT></C><C><IT>&vtdot;</IT></C><C><IT>&vtdot;</IT></C><C><IT>&vtdot;</IT></C></R><R><C>cos <FENCE><FR><NU><IT>2&pgr;t<SUB>N</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C><C>sin <FENCE><FR><NU><IT>2&pgr;t<SUB>N</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C><C>cos <FENCE><FR><NU><IT>4&pgr;t<SUB>N</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C><C>sin <FENCE><FR><NU><IT>4&pgr;t<SUB>N</SUB></IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></C></R></AR></FENCE>  (17)

If we set Z1(tau ) = (ZHR:ZFR), Z2(tau ) = (ZHR:ZFD), xi 1* = (xi HR:xi FR), and xi 2* = (xi HR:xi FD), then for the harmonic regression models SN is
S<SUB>N</SUB>=(y−Z<SUB>i</SUB>(&tgr;)&xgr;<SUP>*</SUP><SUB>i</SUB>)<SUP>t</SUP>&Ggr;(&agr;)<SUP>−1</SUP>(y−Z<SUB>i</SUB>(&tgr;)&xgr;<SUP>*</SUP><SUB>i</SUB>) (18)
where i = 1 and 2 denote, respectively, the free-run and forced desynchrony protocols. Given estimates of tau  and alpha , the maximum likelihood estimate of xi i* is
<A><AC>&xgr;</AC><AC>ˆ</AC></A><SUP>*</SUP><SUB>i</SUB>(&tgr;,&agr;)=[Z<SUP>t</SUP><SUB>i</SUB>(&tgr;)&Ggr;(&agr;)<SUP>−1</SUP>Z<SUB>i</SUB>(&tgr;)]<SUP>−1</SUP>Z<SUP>t</SUP><SUB>i</SUB>(&tgr;)&Ggr;(&agr;)<SUP>−1</SUP>y (19)
For the harmonic regression model, the concentrated log-likelihood in Eq. 12 is maximized by using Newton's procedure to find <A><AC>&tgr;</AC><AC>ˆ</AC></A> and <A><AC>&agr;</AC><AC>ˆ</AC></A>, Gamma (alpha )-1 and  Gamma (alpha )  are computed as above by the Kalman filter, and <A><AC>&xgr;</AC><AC>ˆ</AC></A>i* is computed from Eq. 19. Although the dimension of the full forced desynchrony parameter vector is 18, the dimension of Newton's procedure parameter vector is only two.

We compute the standard errors of the model parameters from the inverse of the observed Fisher information matrices (5, 6, 12, 13). Computing the partial derivatives of the log likelihood with respect to the van der Pol model parameters entails solving an auxiliary differential equation system for each parameter (2). These computations are outlined in APPENDIX B. Explicit analytic formulas for the harmonic regression model parameter information and covariance matrices are given, respectively, in Propositions 1 and 3 of Ref. 12.

Approximate formulas for the variance of the van der Pol parameters can be derived by combining the asymptotic series approximation to the van der Pol model in Eq. 10 with the formula for the asymptotic covariance matrix of the harmonic regression parameter estimates (6, 12). For a harmonic regression model based on Eq. 10, approximate maximum likelihood estimates of gamma  to O(epsilon ) and of epsilon  to O(epsilon 2) are
<A><AC>&ggr;</AC><AC>ˆ</AC></A>=(<A><AC>A</AC><AC>ˆ</AC></A><SUP>2</SUP><SUB>1</SUB>+<A><AC>B</AC><AC>ˆ</AC></A><SUP>2</SUP><SUB>1</SUB>)<SUP>1/2</SUP> (20)

<A><AC>&egr;</AC><AC>ˆ</AC></A>=8[(<A><AC>A</AC><AC>ˆ</AC></A><SUP>2</SUP><SUB>3</SUB>+<A><AC>B</AC><AC>ˆ</AC></A><SUP>2</SUP><SUB>3</SUB>)/(<A><AC>A</AC><AC>ˆ</AC></A><SUP>2</SUP><SUB>1</SUB>+<A><AC>B</AC><AC>ˆ</AC></A><SUP>2</SUP><SUB>1</SUB>)]<SUP>1/2</SUP>
By the Theorem in Ref. 6 and Proposition 3 of Ref. 12, the approximate variances of <A><AC>&tgr;</AC><AC>ˆ</AC></A>, <A><AC>&ggr;</AC><AC>ˆ</AC></A>, and <A><AC>&egr;</AC><AC>ˆ</AC></A> are
Var(<IT><A><AC>&tgr;</AC><AC>ˆ</AC></A></IT>)<IT>=</IT><IT>12&tgr;<SUP>4</SUP> </IT><FENCE><AR><R><C></C></R><R><C></C></R><R><C></C></R><R><C></C></R><R><C></C></R></AR><FR><NU><IT>&ggr;<SUP>2</SUP></IT></NU><DE><IT>f</IT><FENCE><FR><NU><IT>2&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></DE></FR><IT>+</IT><FR><NU><IT>&egr;<SUP>2</SUP>&ggr;<SUP>2</SUP></IT></NU><DE><IT>64f</IT><FENCE><FR><NU><IT>6&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></DE></FR></FENCE><SUP><IT> −1</IT></SUP><IT>&Tgr;<SUP>−3</SUP></IT> (21)

Var(<IT><A><AC>&ggr;</AC><AC>ˆ</AC></A></IT>)<IT>=</IT><IT>4&pgr;f</IT><FENCE><FR><NU><IT>2&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE><IT>&Tgr;<SUP>−1</SUP></IT>

Var(<IT><A><AC>&egr;</AC><AC>ˆ</AC></A></IT>)<IT>=</IT><FR><NU><IT>256&pgr;</IT></NU><DE><IT>&ggr;<SUP>2</SUP></IT></DE></FR> <FENCE><FR><NU><IT>&egr;<SUP>2</SUP></IT></NU><DE><IT>64</IT></DE></FR><IT> f</IT><FENCE><FR><NU><IT>2&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE><IT>+f</IT><FENCE><FR><NU><IT>6&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE></FENCE><IT>&Tgr;<SUP>−1</SUP></IT>
where
f<FENCE><FR><NU>2&pgr;</NU><DE>&tgr;</DE></FR></FENCE>=<FR><NU>&sfgr;<SUP>2</SUP><SUB>&eegr;</SUB>&Dgr;t</NU><DE>2&pgr;<FENCE>1−2e<SUP>−&agr;&Dgr;t</SUP> cos <FENCE><FR><NU><IT>2&pgr;&Dgr;t</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE><IT>+e<SUP>−2&agr;&Dgr;t</SUP></IT></FENCE></DE></FR> (22)
is the spectrum of the discrete AR(1) process.

To evaluate model goodness-of-fit, we report, in addition to graphic summaries, the estimate of the residual variances, Akaike's Information Criterion (AIC), Bayesian Information Criterion (BIC) (4), and analyses of the spectra of the residuals from the model fits. The AIC and BIC are used to help identify the model that is closer to the true model for each subject. With respect to either AIC or BIC, the model with the smallest value of the criterion, i.e., the model with the largest negative value of the criterion, is the one closer to the true model.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Forced desynchrony protocol. We analyze data from five healthy male subjects, ages 21-25, monitored on a 28-h forced desynchrony protocol for 22-30 days, with core-temperature measurements made at 1-min intervals. A subsample taken at 20-min intervals is analyzed for each subject. The harmonic regression model analysis of these data was reported in Ref. 12. We compare those harmonic regression analysis findings with an analysis based on the van der Pol model (Eq. 3) with I(t) = 0, because light levels were 20 lux or less.

It is clear from Fig. 1 that the core-temperature series of each subject consists of a circadian, a forced desynchrony, and a correlated noise component. The strong interaction between the two periodic processes is evidenced by the destructive interference at days 4, 11, and 18 in each subject's data (Fig. 1). This approximate 7-day modulation period is expected, because the two periods are 28 and ~24 h. With both the harmonic regression and van der Pol models, each subject's core-temperature series was readily separated into its circadian, forced desynchrony, and thermoregulatory components, as the data from subject 3 illustrate (Fig. 2). Large values of the residuals occur at the points of maximum excursion of the temperature series (Fig. 3A). The log spectra, computed by periodogram smoothing with a span 100-modified Daniell filter with 20% tapering, were used along with their associated 95% confidence intervals computed by chi 2 approximation (3) to assess goodness-of-fit (Fig. 3B). No subject had any statistically significant frequencies in the spectrum of his residual series for either the harmonic regression or van der Pol models. Both AIC and BIC suggest that the van der Pol model gives a statistically better fit for subjects 1 and 2, whereas the harmonic regression model gives a better fit on the basis of these criteria to the data from subjects 3, 4, and 5 (Tables 1 and 2).


View larger version (38K):
[in this window]
[in a new window]
 
Fig. 1.   Forced desynchrony core-temperature data of subjects 1-5 in order, panels A-E.



View larger version (36K):
[in this window]
[in a new window]
 
Fig. 2.   Van der Pol model analysis of the forced desynchrony core-temperature series of subject 3. A: estimated circadian signal. B: forced desynchrony component. C: the sum of the circadian and forced desynchrony components. Because the circadian signal is not on its limit cycle, the exact pattern of constructive and destructive interference does not repeat every 7 days. D: estimate of the thermoregulatory component.



View larger version (28K):
[in this window]
[in a new window]
 
Fig. 3.   Residual analysis of the fit of the van der Pol model to the core-temperature data of subject 3. A: the residuals; B: the log spectrum. There are large residuals at the points where the core-temperature series has large excursions. None of the estimated spectral ordinates has any statistically significant power, suggesting that the model has captured most of the important systematic structure in the temperature series.


                              
View this table:
[in this window]
[in a new window]
 
Table 1.   Summary of van der Pol equation analysis of forced desynchrony core-temperature series


                              
View this table:
[in this window]
[in a new window]
 
Table 3.   Free-run core-temperature series analysis: subject 3 


                              
View this table:
[in this window]
[in a new window]
 
Table 2.   Summary of harmonic regression analysis of forced desynchrony core-temperature series

With the exception of <A><AC>&egr;</AC><AC>ˆ</AC></A> for subjects 1 and 2, all of the model parameter estimates for each subject are statistically significant, because all estimates are more than two times greater than their associated standard errors. The estimated van der Pol signal (Fig. 2A) is very sinusoidal, consistent with the estimate of this subject's epsilon  being <0.10. Subjects 1, 2, and 3 have estimates of epsilon  between 0.01 and 0.015, whereas for subjects 4 and 5 the estimates are, respectively, 0.15 and 0.18. With the exception of subject 1, the forced desynchrony amplitude estimates for both models are as large or larger than the circadian amplitude estimates. The square-wave shape of the forced desynchrony component is shown in Fig. 2B. The sum of the circadian and forced desynchrony components nicely reproduces the constructive and destructive interference seen in the original data (Fig. 2C). The estimates of alpha  confirm the strong serial dependence in the core-temperature series. The estimates correspond to serial correlation coefficients ranging from 0.866 to 0.920 for the harmonic regression models and from 0.88 to 0.93 for the van der Pol models. The range of alpha -1, the time constant for the thermoregulatory process, is 2.31-3.98 h for the harmonic regression model and 2.74-4.18 h for the van der Pol model.

All of the period estimates are within 15 min of 24 h. The lengths of the approximate 99% Gaussian confidence intervals computed as 5.15 × se(<A><AC>&tgr;</AC><AC>ˆ</AC></A>) suggest that all of the period estimates are different from 25 h. These confidence intervals are 7.1-10.7 min in length for the harmonic regression model and 7.4-12.4 min in length for the van der Pol model. Together with the parameter estimates for <A><AC>&tgr;</AC><AC>ˆ</AC></A>, this finding offers strong evidence that the intrinsic period of the human biological clock is closer to 24 than to 25 h, as recently suggested by Czeisler et al. (18).

Klerman et al. (31) suggested that the light driving effects may bias estimates of intrinsic circadian period unless forced desynchrony studies were >= 20 days long and conducted with light levels in the range of 10-15 lux. Because our experiments were conducted at 20 lux, we investigated the consequences of neglecting the effect of light driving on the estimation of tau .

To investigate the effect of 20-lux light, we simulated the solution to the van der Pol differential equation for 25 days with and without a 20-lux driving term by use of the parameter estimates from subject 4. The parameters of the light driving term were selected to agree with those reported in Ref. 31. Initial conditions were chosen to be identical. Simulations were carried out using the variable step size Cash-Karp Runge-Kutta method (43) of order five with a maximum step size of 1/60 h and an error tolerance of 10-7 to control truncation errors. By day 4, the undriven van der Pol model had decayed to its limit cycle (Fig. 4). The driven solution is also close to an undriven solution over the entire 25 days. The period of the circadian pacemaker, computed by comparing the times of successive minima of the solutions on days 23 and 24, were 24.25 h for the driven solution and 24.18 h for the undriven solution. These differences could not be distinguished by our statistical estimation procedures. By day 25, there is a 6% difference in the amplitudes and a half-hour phase difference in the two solutions. We conclude, therefore, that omitting the light driving term under low light conditions is reasonable.


View larger version (29K):
[in this window]
[in a new window]
 
Fig. 4.   Simulation study of the van der Pol model with and without the light driving terms. A and B together show a 25-day simulation of the van der Pol model equations for parameters and initial conditions chosen to be equal to those of subject 4 in Table 1. Solid line, a simulation in which light driving effects were omitted. Dotted line, light driving effects corresponding to the 20-lux light-dark cycle of the forced desynchrony protocol. The light driving parameters of Eq. 3 were chosen as the circadian light modulation index (m) = 0.333 and a constant of proportionality (C) = 0.018, respectively, to coincide with the estimates given in Ref. 31. Periods of the oscillation differ only slightly, resulting in just a small phase difference near the end of the simulation.

Free-run protocol. To compare the analysis of period estimates obtained from forced desynchrony studies with those from a free-run study, subject 3 underwent both the forced desynchrony protocol and an 18-day free-run analysis. Light intensity during the self-selected lights-on epoch was 150 lux (Fig. 5A). An earlier analysis of this subject's free-run data reported by Czeisler et al. (18) showed that his estimated intrinsic period from the free-run study was 25.13 h. The analysis by Czeisler et al. used the harmonic regression model with correlated noise but without consideration of the effect of light on the pacemaker and of activity on the observed temperature rhythm. Therefore, we reanalyzed this subject's free-run core-temperature series with the van der Pol model, in which the light input is defined by Eq. 3. Light acts directly on the oscillator, as indicated in Eq. 3, and activity acts directly on the observed temperature data during the self-selected lights-on episodes as described by Eq. 6.


View larger version (29K):
[in this window]
[in a new window]
 
Fig. 5.   Harmonic regression model analysis of the free-run core-temperature measurements of subject 3. A: core-temperature series; B: estimated circadian signal; C: estimated thermoregulatory component. Circadian signal is much more square-wave in shape than the corresponding estimate obtained from the van der Pol analysis of the forced desynchrony data in Fig. 2A or the free-run data in Fig. 6B. The first-order autoregressive [AR(1)] estimate of the thermoregulatory response shows a large perturbation between days 6 and 7.

Reanalysis of this subject's free-run data with the harmonic regression model, including the correlated noise, yields a period estimate of 25.13 h and an amplitude estimate of 0.986°F (Table 2). The highly periodic structure of the original data is shown in Fig. 6A. The harmonic regression signal estimate shows an asymmetric waveform that is a square-wave at the crests and narrowed peaks at the troughs (Fig. 6B). The estimated thermoregulatory process shows a prominent perturbation between days 6 and 8 (Fig. 6C). The estimates of tau  and the harmonic regression signal amplitude are statistically significant.


View larger version (28K):
[in this window]
[in a new window]
 
Fig. 6.   Van der Pol model analysis of the free-run core-temperature data of subject 3. A: free-run core-temperature measurements; B: estimated circadian signal; C: beta  (Eq. 6) time-dependent evoked effects of activity during self-selected lights-on intervals; D: estimated thermoregulatory component. The estimated circadian signal is more sinusoidal, with a smaller amplitude than the harmonic regression estimate in Fig. 5B. The pseudoperiodic character of the evoked effects is evident in C. For the van der Pol model, the evoked effect of the perturbation on days 6 and 7 affects both the oscillator and the thermoregulatory process.

The van der Pol model gives estimates of period and amplitude that are much more consistent with those obtained from the analysis of the subject's forced desynchrony data (Tables 1 and 2). The pseudoperiodic effect of activity on the subject's core temperature gives an increase in amplitude during the time in which the lights are on of 0.357°F (Fig. 6C). The sum of the evoked effect of activity on amplitude (0.357°F) and the amplitude of the estimated circadian signal (0.405°F) is approximately equal to the circadian amplitude estimate (0.986°F) from the harmonic regression analysis. The estimated circadian signal from the van der Pol model shows a slightly asymmetrical sinusoidal waveform, which undergoes perturbation between days 6 and 8 (Fig. 6B). Both AIC and BIC suggest that the statistical fit of the van der Pol model is better. A comparison of the van der Pol and harmonic regression analyses from this subject's forced desynchrony study (Table 1) suggests that the van der Pol model provides a more physiologically reasonable fit. The estimate of 0.08 for this subject's internal stiffness parameter shows that, under the free-run conditions, the van der Pol oscillator also behaves like a weakly nonlinear oscillator.

The analytic approximations for the parameter variances presented in Eq. 21 illustrate why the van der Pol and harmonic regression analyses give similar precision for the estimated periods for both the forced desynchrony and free-run analyses. In particular, from Eq. 21 we see that if epsilon  is small, and the circadian signal is close to its limit cycle, the variance of the period estimate is of order T-3.

Constant-routine study. The forced desynchrony study of each subject described in Forced desynchrony protocol was followed immediately by a constant-routine study to estimate the subject's circadian phase without the masking effects of activity and the 28-h-day sleep-wake cycle. Because we established the validity of the unforced van der Pol model at low light levels, and because this model gave a very good description of the forced desynchrony data, we predicted the circadian phase in the constant-routine study from the estimated solution obtained from forced desynchrony data analysis. We did this by evaluating the solution of the estimated van der Pol equation from the beginning of the forced desynchrony protocol until the end of the constant routine. We compared forced desynchrony phase prediction for the constant-routine study with the phase estimates computed from the fits of two harmonic regression models to the constant-routine data. The first model was the two-harmonic model (HR2) for constant-routine core-temperature analyses proposed in Ref. 8, and the second was the limit cycle approximation to the van der Pol model HRvdp in Eq. 10. The HR2 model contains a fundamental and its second harmonic, whereas the HRvdp model contains a fundamental and a third harmonic. We determined the best fit of the harmonic regression models with the period constrained to 24.00-24.30 h as in Ref. 8.

Comparisons of these models are shown in Fig. 7. All three models yield similar amplitude estimates of 0.42, 0.32, and 0.43°F, respectively. None of the three circadian signal estimates lies on the observed core-temperature data, because the thermoregulatory component (not shown) in these data is quite strong. The van der Pol model predicts that the phase of the minimum should occur 36.45 h after the initiation of constant routine, the HR2 model yields an estimate of 35.83 h, and the HRvdp model yields 36.30 h. Although different, all are within the approximate 1.5-h width of the 95% confidence interval for phase of the temperature minimum reported in Ref. 8. The van der Pol model phase estimate agrees more closely with the HRvdp model than with the HR2 model. The HRvdp waveform estimate is less sinusoidal then the van der Pol model prediction, because the HRvdp coefficients are not constrained to satisfy the conditions of the coefficients in Eq. 10. This analysis suggests that the HR2 model used in constant-routine analyses may give a circadian signal estimate different from that estimated by the van der Pol model.


View larger version (19K):
[in this window]
[in a new window]
 
Fig. 7.   Comparison of constant-routine harmonic regression analyses with a prediction based on a van der Pol model fit. Core-temperature data from the post-forced desynchrony constant-routine study of subject 1 are the points in the graph. Superimposed on these data are 3 model fits: the HR2 model (dotted lines), the HRvdp model (dashed lines), and the van der Pol model (undivided line) derived from its fit to the data in the preceding forced desynchrony protocol. Neither of the three fits falls directly on the data because of the strong thermoregulatory component (not shown). The van der Pol model phase prediction agrees most closely with the HRvdp model estimate.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

We have presented a single-equation statistical representation of the human core-temperature circadian rhythm by combining models from previous theoretical and empirical circadian studies of human core temperature. The central elements of our new model are the representation of the circadian pacemaker and its light inputs as a modified van der Pol oscillator, the body's thermoregulatory process as an AR(1) process, and a protocol-specific function to describe the effect of activity on core temperature. We presented both theoretical and computational results for the model.

Our theoretical work links the harmonic regression methods used in analyses of experimental circadian data and the van der Pol differential equation model used in simulation studies of the pacemaker. We showed analytically and with simulation studies that, when there is little or no light input to the pacemaker and the internal stiffness parameter epsilon  is small, the van der Pol behaves like a simple harmonic regression model in the odd harmonics of the period. We combined the perturbation series expansions with previously derived asymptotic formulas for the harmonic regression parameter covariance matrix and showed that, just like the variance of the period estimate of the harmonic regression model (6, 12, 18), the variance of the period estimate for the weakly nonlinear van der Pol model is of the order T-3.

The maximum likelihood model-fitting algorithm (Eqs. 11-16) and the sensitivity algorithm (Eq. B1) for estimating the standard error make it possible to apply the new model in actual core-temperature analyses. The forced desynchrony model can have as many as 18 parameters, whereas the free-run model can have as many as 11. By separating the model parameters into those that are nonlinear and those that are linear, our maximum likelihood algorithm leads to a reliable technique for fitting these high-dimensional models. Because the numbers of data points in the forced desynchrony and free-run studies are large (ranging from 1,548 to 2,160), we used asymptotic theory to compute the standard errors for the maximum likelihood parameter estimates from the estimated Fisher information matrix (5, 6, 12, 13). The sensitivity algorithm gives an efficient way to compute the Fisher information matrix for the differential equation model. The standard errors for the period and amplitude for the van der Pol model agree closely with those derived from harmonic regression model fits. Our results in Eq. 21 illustrate why this should be the case for the forced desynchrony protocol (12). The factors other than sample size that contribute to the higher precision of the period estimates were studied for the harmonic regression forced desynchrony model by Brown et al. (12). Regular use of these standard error algorithms for the van der Pol model in forced desynchrony and free-run analyses will require a systematic study of their accuracy by use of synthetic data to determine the smallest sample size at which the asymptotic theory holds. This analysis will be the topic of a future report.

Application of the new model to experimental data from the three principal circadian protocols offers new insights into the use of human core-temperature data to analyze the circadian pacemaker. First, our work makes it possible to study the dynamic features of the human circadian pacemaker directly from core-temperature series. With a few exceptions (5, 7, 9, 10, 27), the reported values of epsilon , the parameter in the van der Pol model that governs the dynamical properties of the human circadian pacemaker, have been determined almost exclusively through simulation studies. Kronauer and colleagues (29-32, 34, 35) have cited this parameter value as ~0.13.

In a time-zone shift study, Gundel and Spencer (27) used least squares to fit an unscaled, sinusoidally forced version of Eq. 1 to core-temperature series. Eleven of the twelve subjects they studied had values of epsilon  > 0.28 (median: 0.345; range: 0.10-0.74). Gundel and Spencer's approach parallels ours. An important difference between their model and ours is that, in their analyses, the scaling of the van der Pol is not clearly defined. Consequently, the estimated values of epsilon , gamma , and tau  are confounded. This confounding is suggested by the fact that six of the twelve subjects studied had values of epsilon  >= 0.34, a range that makes the van der Pol model less consistent with a weakly nonlinear oscillator. The large values of epsilon  may also reflect confounding of this parameter with the thermoregulatory process, which was not included in their model. Their analysis also does not provide parameter standard error estimates.

Our results showed that, for the forced desynchrony and free-run protocols, the values of epsilon  are ~0.1 and range from 0.015 to 0.19. Fitting the model directly to core-temperature data gives an assessment of the biological variability in the dynamical properties of the human circadian pacemaker. This assessment could not be made from simulation studies. The van der Pol analysis also provides statistical estimates of the other model parameters, such as m and C, which may be used in simulation studies of the human circadian pacemaker. Van der Pol parameters estimated from experimental data can be fed back into simulation studies to yield more accurate investigations of the pacemaker.

Second, our model represents both the dynamical properties of the circadian pacemaker and its response to light. Therefore, it can provide a direct estimate from experimental data of the impact of light on the pacemaker. This point was illustrated in the free-run analysis. By including in the free-run model the 150-lux ambient light input to the oscillator and the evoked effect of activity on core temperature, we extracted a period estimate from the free-run analysis that was completely consistent with the one determined from the forced desynchrony studies. What the harmonic regression analysis of the free-run data estimated as a circadian signal with a period of slightly greater than 25 h, the van der Pol analysis showed to be a combination of circadian signal with a period of 24.25 h, stimulated by a pseudoperiodic process defined by the times at which the subject self-administered light. The van der Pol analysis of the free-run data showed that the harmonic regression amplitude estimate consisted of one component due to the circadian amplitude and a second component attributable to the evoked effect of activity on the observed temperature rhythm. We are currently using our new model to estimate from experimental data the impact of different lighting regimens on the dynamics of the circadian pacemaker.

Third, we reanalyzed 5 of the 10 young forced desynchrony subjects studied by Czeisler et al. (18) with our more comprehensive model and confirmed their recent report that the intrinsic period of the human biological clock is closer to 24 than to 25 h (18). Like the harmonic regression model, the van der Pol model with and without light input yields highly precise period estimates due to the inverse cubic dependence [O(T-3)] of the period variance on the study length (Eq. 21a). This finding suggests that, like the harmonic regression model, the weakly nonlinear van der Pol model may be used to obtain precise period estimates of the circadian pacemaker from forced desynchrony and free-run studies. Furthermore, the covariance matrix of these model parameters makes it possible not only to use single parameter estimates in the simulation studies but also to suggest reasonable conjoined regions of parameters for these studies that are consistent with experimental data.

Finally, our constant-routine analysis identified an important inconsistency between simulation studies of the circadian pacemaker using the van der Pol equation and empirical analyses of core-temperature data based on the harmonic regression model. The harmonic expansion of the weakly nonlinear van der Pol model consists of a first and a third harmonic (Eq. 10), whereas constant-routine core-temperature data are well described by a harmonic regression model consisting of a first and a second harmonic (8). This inconsistency between the two modeling frameworks would not be appreciated without the current analysis. Although the human circadian pacemaker behaves like a weakly nonlinear oscillator, its modeling as a van der Pol oscillator fails to predict the robust second harmonic characteristic of the circadian signal estimated from constant-routine core-temperature data. Core temperature is one of the most reliable markers of the human circadian pacemaker, and the constant-routine study is one of the best protocols for observing the pacemaker output from core temperature with minimal exogenous perturbations. Hence there is compelling need to develop a model that both captures the weakly nonlinear dynamics of the circadian system and accounts for the strong second harmonic present in constant-routine circadian signal estimates. Such a model will almost certainly further our understanding of this oscillator and its dynamics.

The van der Pol model is an analytically tractable, parsimonious equation that describes well many salient features of the circadian signal and its dynamical properties. Aside from the parameter definitions, none of the model components has a specific physiological interpretation. More physiologically consistent core-temperature models must consider the dynamical properties of the circadian signal in a unified framework along with core temperature's strong thermoregulatory component, its protocol-dependent evoked components, and the interactions among these factors. Dynamical models of this category are the topic of our current investigations.


    APPENDIX A
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Asymptotic Series Approximation to the van der Pol Equation

We derive the asymptotic series approximation to the van der Pol equation in Eq. 1 by applying the method of Refs. 36 and 41. The order epsilon K+1 asymptotic series expansion of the solution to Eq. 1 is assumed to have the form
s(t)=&ggr;(t) cos<IT> &psgr;</IT>(<IT>t</IT>)<IT>+</IT><LIM><OP>∑</OP><LL><IT>k=1</IT></LL><UL><IT>K</IT></UL></LIM><IT> &egr;<SUP>k</SUP>u<SUB>k</SUB></IT>[<IT>&ggr;</IT>(<IT>t</IT>)<IT>, &psgr;</IT>(<IT>t</IT>)]<IT>+O</IT>(<IT>&egr;<SUP>K+1</SUP></IT>) (A1)
K = 1,2,..., where gamma (t) and Psi (t) are defined by the differential equations
<FR><NU>d<IT>&psgr;</IT></NU><DE>d<IT>t</IT></DE></FR><IT>=&ohgr;+</IT><LIM><OP>∑</OP><LL><IT>k=1</IT></LL><UL><IT>K</IT></UL></LIM><IT> &egr;<SUP>k</SUP>&psgr;<SUB>k</SUB></IT>(<IT>&ggr;</IT>(<IT>t</IT>))<IT>+O</IT>(<IT>&egr;<SUP>K+1</SUP></IT>) (A2)

<FR><NU>d<IT>&ggr;</IT></NU><DE>d<IT>t</IT></DE></FR><IT>=</IT><LIM><OP>∑</OP><LL><IT>k=1</IT></LL><UL><IT>K</IT></UL></LIM><IT> &egr;<SUP>k</SUP>A<SUB>k</SUB></IT>(<IT>&ggr;</IT>(<IT>t</IT>))<IT>+O</IT>(<IT>&egr;<SUP>K+1</SUP></IT>)
where omega  = 2pi /tau , and for each k, uk(gamma ,psi ) satisfies
<FR><NU>1</NU><DE>2&pgr;</DE></FR> <LIM><OP>∫</OP><LL>0</LL><UL>2&pgr;</UL></LIM> u<SUB>k</SUB>(&ggr;,&psgr;) cos<IT> &psgr;d&psgr;=</IT><FR><NU><IT>1</IT></NU><DE><IT>2&pgr;</IT></DE></FR> <LIM><OP>∫</OP><LL><IT>0</IT></LL><UL><IT>2&pgr;</IT></UL></LIM><IT> u<SUB>k</SUB></IT>(<IT>&ggr;,&psgr;</IT>) sin<IT> &psgr;d&psgr;=0</IT> (A3)
The conditions in Eq. A3 ensure the periodicity of the resulting approximate solution by preventing the appearance of secular terms. The first and second derivatives of Eqs. A2a and A2b up to O(epsilon 3) are
<FR><NU>d<IT>&ggr;</IT></NU><DE>d<IT>t</IT></DE></FR><IT>=</IT><IT>&egr;A<SUB>1</SUB>+&egr;<SUP>2</SUP>A<SUB>2</SUB>+O</IT>(<IT>&egr;<SUP>3</SUP></IT>) (A4)

<FR><NU>d<IT>&psgr;</IT></NU><DE>d<IT>t</IT></DE></FR> <FR><NU>d<IT>&ggr;</IT></NU><DE>d<IT>t</IT></DE></FR><IT>=</IT><IT>&egr;A<SUB>1</SUB>&ohgr;+&egr;<SUP>2</SUP></IT>[<IT>A<SUB>2</SUB>&ohgr;+A<SUB>1</SUB>&psgr;<SUB>1</SUB></IT>]<IT>+O</IT>(<IT>&egr;<SUP>3</SUP></IT>)

<FR><NU>d<IT>&psgr;</IT></NU><DE>d<IT>t</IT></DE></FR><IT>=</IT><IT>&ohgr;+&egr;&psgr;<SUB>1</SUB>+&egr;<SUP>2</SUP>&psgr;<SUB>2</SUB>+O</IT>(<IT>&egr;<SUP>3</SUP></IT>)

<FENCE><FR><NU>d<IT>&ggr;</IT></NU><DE>d<IT>t</IT></DE></FR></FENCE><SUP><IT>2</IT></SUP><IT>=</IT><IT>&egr;<SUP>2</SUP>A</IT><SUP><IT>2</IT></SUP><SUB><IT>1</IT></SUB><IT>+O</IT>(<IT>&egr;<SUP>3</SUP></IT>)

<FR><NU>d<SUP><IT>2</IT></SUP><IT>&ggr;</IT></NU><DE>d<IT>t<SUP>2</SUP></IT></DE></FR><IT>=</IT><IT>&egr;<SUP>2</SUP>A<SUB>1</SUB> </IT><FR><NU>d<IT>A<SUB>1</SUB></IT></NU><DE>d<IT>&ggr;</IT></DE></FR><IT>+O</IT>(<IT>&egr;<SUP>3</SUP></IT>)

<FENCE><FR><NU>d<IT>&psgr;</IT></NU><DE>d<IT>t</IT></DE></FR></FENCE><SUP><IT>2</IT></SUP><IT>=</IT><IT>&ohgr;<SUP>2</SUP>+2&egr;&ohgr;&psgr;<SUB>1</SUB>+&egr;<SUP>2</SUP></IT>(<IT>&psgr;<SUB>1</SUB>+2&ohgr;&psgr;<SUB>2</SUB></IT>)<IT>+O</IT>(<IT>&egr;<SUP>3</SUP></IT>)
Substituting from Eq. A.4 into Eq. 1 gives
<A><AC>s</AC><AC>¨</AC></A>(t)+&ohgr;<SUP>2</SUP>s(t) (A5)

=&egr;<FENCE>−<IT>2&ohgr;A<SUB>1</SUB> </IT>sin<IT> &psgr;−2&ohgr;&psgr;<SUB>1</SUB> </IT>cos<IT> &psgr;+&ohgr;<SUP>2</SUP> </IT><FR><NU><IT>∂<SUP>2</SUP>u<SUB>1</SUB></IT></NU><DE><IT>∂&psgr;<SUP>2</SUP></IT></DE></FR><IT>+&ohgr;<SUP>2</SUP>u<SUB>1</SUB></IT></FENCE><IT>+O</IT>(<IT>&egr;<SUP>2</SUP></IT>)
and
&egr;&ohgr;<FENCE>1−<FR><NU>4</NU><DE>&ggr;<SUP>2</SUP></DE></FR> s<SUP>2</SUP>(t)</FENCE><A><AC>s</AC><AC>˙</AC></A>(t) (A6)

=&egr;<FENCE>−<IT>&ggr;</IT>(<IT>t</IT>)<IT>&ohgr;<SUP>2</SUP></IT><FENCE><IT>1−</IT><FR><NU><IT>1</IT></NU><DE><IT>&ggr;<SUP>2</SUP></IT></DE></FR></FENCE> sin<IT> &psgr;+</IT><FR><NU><IT>&ohgr;<SUP>2</SUP></IT></NU><DE><IT>&ggr;<SUP>2</SUP></IT></DE></FR><IT> &ggr;</IT>(<IT>t</IT>)<SUP><IT>3</IT></SUP> sin<IT> 3&psgr;</IT></FENCE><IT>+O</IT>(<IT>&egr;<SUP>2</SUP></IT>)
The solution to O(epsilon 2) is obtained by equating the coefficients of the epsilon  terms in Eqs. A5 and A6 subject to the conditions in Eq. A3. This yields
<FR><NU>∂<SUP>2</SUP>u<SUB>1</SUB></NU><DE>∂&psgr;<SUP>2</SUP></DE></FR>+u<SUB>1</SUB>=<FR><NU>&ggr;(t)<SUP>3</SUP></NU><DE>&ggr;<SUP>2</SUP></DE></FR> sin<IT> 3&psgr;+</IT><FR><NU><IT>2</IT></NU><DE><IT>&ohgr;</IT></DE></FR><IT> &psgr;<SUB>1</SUB> </IT>cos<IT> &psgr;</IT> (A7)

A<SUB>1</SUB>=<FR><NU>&ggr;(t)&ohgr;</NU><DE>2</DE></FR> <FENCE>1−<FR><NU>&ggr;(t)<SUP>2</SUP></NU><DE>&ggr;<SUP>2</SUP></DE></FR></FENCE>
To satisfy Eq. A3, we require psi 1 = 0, A1 as given in Eq. A7b, and u1 to be
u<SUB>1</SUB>(&ggr;, &psgr;)=−<FR><NU><IT>&ggr;</IT>(<IT>t</IT>)<SUP><IT>3</IT></SUP></NU><DE><IT>8&ggr;<SUP>2</SUP></IT></DE></FR> sin<IT> 3&psgr;</IT> (A8)
When Eq. A7a is used, the solutions for psi (t) and gamma (t) to O(epsilon 2) in Eq. A4 are
&psgr;(t)=&ohgr;t+&psgr;<SUB>0</SUB> (A9)

&ggr;(t)=<FR><NU>2C<SUP>1/2</SUP>e<SUP>&egr;&ohgr;t/2</SUP></NU><DE><FENCE>1+C <FR><NU>4</NU><DE>&ggr;<SUP>2</SUP></DE></FR> e<SUP>&egr;&ohgr;t</SUP></FENCE><SUP>1/2</SUP></DE></FR>
where C is a constant of integration. Substituting Eqs. A8 and A9 into Eq. A1 yields the O(epsilon 2) solution given in Eq. 8.


    APPENDIX B
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Computing the Information Matrix for the van der Pol Parameters

For the van der Pol model, the elements of the information matrix are defined as
E<FENCE><FR><NU>∂<SUP>2</SUP>L</NU><DE>∂&thgr;<SUB>i</SUB>∂&thgr;<SUB>j</SUB></DE></FR></FENCE>=<FR><NU>&sfgr;<SUP>2</SUP><SUB>&eegr;</SUB></NU><DE>(1−e<SUP>−2&agr;&Dgr;t</SUP>)</DE></FR> <FR><NU>∂v<SUB>t<SUB>1</SUB></SUB></NU><DE>∂&thgr;<SUB>i</SUB></DE></FR> <FR><NU>∂v<SUB>t<SUB>1</SUB></SUB></NU><DE>∂&thgr;<SUB>j</SUB></DE></FR>

+<FR><NU>1</NU><DE>&sfgr;<SUP>2</SUP><SUB>&eegr;</SUB></DE></FR> <LIM><OP>∑</OP><LL>n=2</LL><UL>N</UL></LIM> <FENCE><FR><NU>∂v<SUB>t<SUB>n</SUB></SUB></NU><DE>∂&thgr;<SUB>i</SUB></DE></FR>−e<SUP>−&agr;&Dgr;t</SUP> <FR><NU>∂v<SUB>t<SUB>n−1</SUB></SUB></NU><DE>∂&thgr;<SUB>i</SUB></DE></FR></FENCE><FENCE><FR><NU>∂v<SUB>t<SUB>n</SUB></SUB></NU><DE>∂&thgr;<SUB>j</SUB></DE></FR>−e<SUP>−&agr;&Dgr;t</SUP> <FR><NU>∂v<SUB>t<SUB>n−1</SUB></SUB></NU><DE>∂&thgr;<SUB>j</SUB></DE></FR></FENCE>

E<FENCE><FR><NU>∂<SUP>2</SUP>L</NU><DE>∂&thgr;<SUB>j</SUB>∂&xgr;<SUB>k</SUB></DE></FR></FENCE>=<FR><NU>&sfgr;<SUP>2</SUP><SUB>&eegr;</SUB></NU><DE>(1−e<SUP>−2&agr;&Dgr;t</SUP>)</DE></FR> <FR><NU>∂v<SUB>t<SUB>1</SUB></SUB></NU><DE>∂&thgr;<SUB>j</SUB></DE></FR> <FR><NU>∂v<SUB>t<SUB>1</SUB></SUB></NU><DE>∂&xgr;<SUB>k</SUB></DE></FR> (B1)

+<FR><NU>1</NU><DE>&sfgr;<SUP>2</SUP><SUB>&eegr;</SUB></DE></FR> <LIM><OP>∑</OP><LL>n=2</LL><UL>N</UL></LIM> <FENCE><FR><NU>∂v<SUB>t<SUB>n</SUB></SUB></NU><DE>∂&thgr;<SUB>j</SUB></DE></FR>−e<SUP>−&agr;&Dgr;t</SUP> <FR><NU>∂v<SUB>t<SUB>n−1</SUB></SUB></NU><DE>∂&thgr;<SUB>j</SUB></DE></FR></FENCE><FENCE><FR><NU>∂v<SUB>t<SUB>n</SUB></SUB></NU><DE>∂&xgr;<SUB>k</SUB></DE></FR>−e<SUP>−&agr;&Dgr;t</SUP> <FR><NU>∂v<SUB>t<SUB>n−1</SUB></SUB></NU><DE>∂&xgr;<SUB>k</SUB></DE></FR></FENCE>
for i,j = 1,...,7, and where xi k is any non-van der Pol parameter. It follows from Eq. 2 that partial vtn/partial theta i = -partial stn/partial theta i. To compute partial stn/partial theta i, let
F(s)=&egr;<FENCE><FR><NU>2&pgr;</NU><DE>&tgr;</DE></FR></FENCE><FENCE>s(t)−<FR><NU>4s(t)<SUP>3</SUP></NU><DE>3&ggr;<SUP>2</SUP></DE></FR></FENCE>+<FENCE><FR><NU>2&pgr;</NU><DE>&tgr;</DE></FR></FENCE>&ggr;(1−ms(t))CI(t)<SUP>1/3</SUP> (B2)

g(s, z)=−<FENCE><FR><NU><IT>2&pgr;</IT></NU><DE><IT>&tgr;</IT></DE></FR></FENCE><SUP><IT>2</IT></SUP><IT>s</IT>(<IT>t</IT>)<IT>+z</IT>(<IT>t</IT>) <FR><NU>(<IT>1−ms</IT>(<IT>t</IT>))<IT>CI</IT>(<IT>t</IT>)<SUP><IT>1/3</IT></SUP></NU><DE><IT>3</IT></DE></FR>
Substituting Eq. B2 into Eq. 3, and interchanging the order of differentiation with respect to time, and differentiation with respect to theta i, give for each theta i the auxiliary differential equation system (2)
<FR><NU>d</NU><DE>d<IT>t</IT></DE></FR> <FENCE><FR><NU><IT>∂s</IT>(<IT>t</IT>)</NU><DE><IT>∂&thgr;<SUB>i</SUB></IT></DE></FR></FENCE><IT>=</IT><FR><NU><IT>∂z</IT>(<IT>t</IT>)</NU><DE><IT>∂&thgr;<SUB>i</SUB></IT></DE></FR><IT>+</IT><FR><NU><IT>∂F</IT>(<IT>s</IT>(<IT>t</IT>))</NU><DE><IT>∂&thgr;<SUB>i</SUB></IT></DE></FR> (B3)

<FR><NU>d</NU><DE>d<IT>t</IT></DE></FR> <FENCE><FR><NU><IT>∂z</IT>(<IT>t</IT>)</NU><DE><IT>∂&thgr;<SUB>i</SUB></IT></DE></FR></FENCE><IT>=</IT><FR><NU><IT>∂g</IT>(<IT>s</IT>(<IT>t</IT>)<IT>, z</IT>(<IT>t</IT>))</NU><DE><IT>∂&thgr;<SUB>i</SUB></IT></DE></FR>
whose initial conditions are
<FR><NU>∂s<SUB>t<SUB>1</SUB></SUB></NU><DE>∂&thgr;<SUB>i</SUB></DE></FR>=<FENCE><AR><R><C>1</C><C>i=1</C></R><R><C>0</C><C>otherwise</C></R></AR></FENCE> (B4)

<FR><NU>∂z<SUB>t<SUB>1</SUB></SUB></NU><DE>∂&thgr;<SUB>i</SUB></DE></FR>=<FENCE><AR><R><C>1</C><C>i=2</C></R><R><C>0</C><C>otherwise</C></R></AR></FENCE>
because theta 1 = st1 and theta 2 ·t1. Each auxiliary system is integrated using the Runge-Kutta method.


    ACKNOWLEDGEMENTS

The authors are grateful to Brenda Marshall for technical assistance in preparing the manuscript.


    FOOTNOTES

This work was supported by Robert Wood Johnson Foundation Grant 23397; National Institute of Health Grants 1-R01-GM-53559, 1-P01-AG-09975, 1-R01-AG-06072, and 1-R01-MH-45130; National Center for Research Resources General Clinical Research Center Grant M01-RR-02635; and the National Aeronautics and Space Administration through the NASA Cooperative Agreement NCC 9-58 with the National Space Biomedical Research Institute.

Address for reprint requests and other correspondence: E. N. Brown, Statistics Research Laboratory, Dept. of Anesthesia and Critical Care, Massachusetts General Hospital, 55 Fruit St., Clinics 3, Boston, MA 02114-2696 (E-mail: brown{at}srlb.mgh.harvard.edu).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. §1734 solely to indicate this fact.

Received 6 December 1999; accepted in final form 24 March 2000.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
PROTOCOLS, MODEL FORMULATION,...
RESULTS
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

1.   Atkinson, KE. An Introduction to Numerical Analysis (2nd ed.). New York: Wiley, 1989.

2.   Bard, Y. Nonlinear Parameter Estimation. New York: Academic, 1974.

3.   Bloomfield, P. Fourier Analysis of Time Series: An Introduction. New York: Wiley, 1976.

4.   Box, GE, Jenkins GM, and Reinsel GC. Times In: Series Analysis: Forecasting and Control (3rd ed.). Englewood Cliffs, NJ: Prentice Hall, 1994.

5.   Brown, EN. Identification and Estimation of Differential Equation Models for Circadian Data (PhD dissertation). Cambridge, MA: Harvard University, 1987.

6.   Brown, EN. A note on the asymptotic distribution of the parameter estimates for the harmonic regression model. Biometrika 77: 653-655, 1990[ISI].

7.   Brown, EN, Choe Y, and Czeisler CA. Statistical analysis human core-temperature rhythms with differential equation methods (Abstract). J Biol Rhythm Res 26: 38, 1995.

8.   Brown, EN, and Czeisler CA. The statistical analysis of circadian phase and amplitude in constant-routine core-temperature data. J Biol Rhythms 7: 177-202, 1992[ISI][Medline].

9.   Brown, EN, and Luithardt H. Statistical model building and model criticism for human circadian data. J Biol Rhythms 14: 609-616, 1999[Abstract/Free Full Text].

10.   Brown, EN, and Moore-Ede MC. A method for assessing the compatibility of differential equation models with observed circadian data (Abstract). Sleep Res 14: 263, 1985.

11.   Brown, EN, and Schmid CH. Application of the Kalman filter to computational problems in statistics. In: Methods in Enzymology, Numerical Computer Methods (Part B). Orlando, FL: Academic, 1994, p. 171-181.

12.   Brown, EN, Solo V, Choe Y, and Zhang Z. Measuring the period of the human biological clock: an infill asymptotic analysis of harmonic regression estimates. Technical Report 95-03. Boston, MA: Statistics Research Laboratory, Department of Anesthesia and Critical Care, Massachusetts General Hospital, 1997.

13.   Casella, G, and Berger RL. Statistical Inference. Belmont, MA: Duxbury, 1990.

14.   Corradi, C. A note on the computation of maximum likelihood estimates of linear regression models with autocorrelated errors. J Econometr 11: 303-317, 1979.

15.   Czeisler, CA. Human Circadian Physiology: Internal Organization of Temperature Sleep-Wake and Neuroendocrine Rhythms Monitored in an Environment Free of Time Cues (PhD dissertation). Palo Alto, CA: Stanford University, 1978.

16.   Czeisler, CA, Allan JS, Strogatz SH, Ronda JM, Sanchez R, Rios CD, Freitag WF, Richardson GS, and Kronauer RE. Bright light resets the human circadian pacemaker independent of the timing of the sleep-wake cycle. Science 233: 667-671, 1986[ISI][Medline].

17.   Czeisler, CA, Brown EN, Ronda JM, Richardson GS, Kronauer RE, and Freitag WF. A clinical method to assess the endogenous circadian phase (ecp) of the deep circadian oscillator in man (Abstract). Sleep Res 14: 295, 1985.

18.   Czeisler, CA, Duffy JD, Shanahan TA, Brown EN, Mitchell JF, Dijk DJ, Rimmer DW, Ronda JM, Allan JS, Emens JS, and Kronauer RE. Stability, precision, and near-24 hour period of the human circadian pacemaker. Science 284: 2177-2181, 1999[Abstract/Free Full Text].

19.   Czeisler, CA, Johnson MP, Duffy JF, Brown EN, Ronda JM, and Kronauer RE. Exposure to bright light and darkness to treat physiologic maladaption to night work. New Engl J Med 322: 1253-1259, 1990[Abstract].

20.   Czeisler, CA, Kronauer RE, Allan JS, Duffy JF, Jewett ME, Brown EN, and Ronda JM. Bright light induction of strong (Type 0) resetting of the human circadian pacemaker. Science 244: 1328-1332, 1989[ISI][Medline].

21.   Enright, JT. Data analysis. In: Handbook of Behavioral Biology, edited by Aschoff J. New York: Plenum, 1981, vol. 4, p. 21-40.

22.   Enright, JT. Methodology. In: Handbook of Behavioral Biology, edited by Aschoff J. New York: Plenum, 1981, vol. 4, p. 11-20.

23.   Forger, DB, Jewett ME, and Kronauer RE. A simpler model of the human circadian pacemaker. J Biol Rhythms 14: 532-537, 1999[Abstract/Free Full Text].

24.   Gander, PH, Kronauer RE, and Graeber RC. Phase shifting two coupled circadian pacemakers: implications for jet lag. Am J Physiol Regulatory Integrative Comp Physiol 249: R704-R719, 1985[ISI][Medline].

25.   Glotzbach, SF, and Heller HC. Sleep and Thermoregulation. In: Principles and Practice of Sleep Medicine, edited by Kryger MH. Philadelphia, PA: Saunders, 1989, p. 300-309.

26.   Greenhouse, JB, Kass R, and Tsay R. Fitting nonlinear models with ARMA errors to biological rhythm data. Stat Med 6: 167-183, 1987[ISI][Medline].

27.   Gundel, A, and Spencer MB. A circadian oscillator model based on empirical data. J Biol Rhythms 14: 516-523, 1999[Abstract/Free Full Text].

28.   Guyton, AC. Textbook of Medical Physiology (8th ed.). Philadelphia, PA: Saunders, 1991.

29.   Jewett, ME, and Kronauer RE. Refinement of a limit cycle oscillator model of the effects of light on the human circadian pacemaker. J Theor Biol 192: 455-465, 1998[ISI][Medline].

30.   Jewett, ME, and Kronauer RE. Interactive mathematical models of subjective alertness and cognitive throughput in humans. J Biol Rhythms 14: 588-597, 1999[Abstract/Free Full Text].

31.   Klerman, EB, Dijk DJ, Kronauer RE, and Czeisler CA. Simulations of light effects on the human circadian pacemaker: implications for assessment of intrinsic period. Am J Physiol Regulatory Integrative Comp Physiol 270: R271-R282, 1996[Abstract/Free Full Text].

32.   Kronauer, RE. Modeling principles for human circadian rhythms. In: Mathematical Models of the Circadian Sleep-Wake Cycle, edited by Moore-Ede MC, and Czeisler CA. New York: Raven, 1984, p. 105-128.

33.   Kronauer, RE. A quantitative model for the effects of light on the amplitude and phase of the deep circadian pacemaker, based on human data. In: Sleep '90, Proceedings of the Tenth European Congress on Sleep Research, edited by Horne J. Dusseldorf, Germany: Pontenagel, 1990, p. 306-309.

34.   Kronauer, RE, Czeisler CA, Pilato SF, Moore-Ede MC, and Weitzman ED. Mathematical model of the human circadian system with two interacting oscillators. Am J Physiol Regulatory Integrative Comp Physiol 242: R3-R17, 1982[Abstract/Free Full Text].

35.   Kronauer, RE, Forger DB, and Jewett ME. Quantifying human circadian pacemaker response to brief, extended, and repeated light stimuli over the phototopic range. J Biol Rhythms 14: 500-515, 1999[Abstract/Free Full Text].

36.   Krylov, K, and Bogoliubov B. Introduction to Nonlinear Mechanics. Princeton, NJ: Princeton University Press, 1947.

37.   Mills, JN, Minors DS, and Waterhouse JM. Adaptation to abrupt time shifts of the oscillator controlling human circadian rhythms. J Physiol (Lond) 285: 455-470, 1978[Abstract].

38.   Minors, DS, and Waterhouse JM. The use of the constant-routines in unmasking the endogenous component of the human circadian rhythms. Chronobiol Int I: 205-216, 1984.

39.   Moore, RY, and Eichler VB. Loss of a circadian adrenal corticosterone rhythm following suprachiasmatic lesions in the rat. Brain Res 42: 201-206, 1972[ISI][Medline].

40.   Moore-Ede, MC, Sulzman FM, and Fuller CA. The Clocks That Time Us. Cambridge, MA: Harvard Univ. Press, 1982.

41.   Nayfeh, AH. Perturbation Methods. New York: Wiley, 1973.

42.   Nayfeh, AH, and Mook DT. Nonlinear Oscillations. New York: Wiley, 1979.

43.   Press, WH, Teukolsky SA, Vetterling WT, and Flannery BP. Numerical Recipes: The Art of Scientific Computing. Cambridge, MA: Cambridge University Press, 1992.

44.   Priestley, MB. Spectral Analysis and Time Series. London: Academic, 1981.

45.   Rosen, R. Dynamical System Theory in Biology. New York: Wiley, 1970.

46.   Smith, RE. Circadian variations in human thermoregulatory responses. J Appl Physiol 26: 554-560, 1969[Free Full Text].

47.   Stephan, FK, and Zucker I. Circadian rhythms in drinking behavior and locomotor activity of rats are eliminated by hypothalamic lesions. Proc Natl Acad Sci USA 69: 1583-1586, 1972[Abstract].

48.   Strogatz, S. The Mathematical Structure of Human Sleep-Wake Cycle (PhD dissertation). Cambridge, MA: Harvard University, 1986.

49.   Tong, YL. Parameter estimation in studying circadian rhythms. Biometrics 32: 85-94, 1976[ISI][Medline].

50.   Wever, AM. The Circadian System of Man. Berlin: Springer-Verlag, 1979.


Am J Physiol Endocrinol Metab 279(3):E669-E683
0193-1849/00 $5.00 Copyright © 2000 the American Physiological Society