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 |
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 |
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
|
(1)
|
where s(t) is the signal from the circadian
pacemaker,
is the approximate period of the oscillation,
is the
approximate limiting amplitude, and
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
, the slower the approach of the oscillator to its limit cycle and the more sinusoidal its oscillations become. Weakly nonlinear means 0 <
1. If
= 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
,
, and
. 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 |
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 = n
t, n =
1,2,...,N, and N
t =
T. For each protocol, the core-temperature measurement may
be expressed as
|
(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
|
(3)
|
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
,
, and
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
|
(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
|
(5)
|
(12, 18). For the self-selected
timing of activity during the free-run protocol,
xtn is
|
(6)
|
This term models the effect of activity during free run as
an increase in mean core temperature by an amount
for
> 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
|
(7)
|
where 
1 is the time constant for the
thermoregulatory response, and the
tn
are independent, Gaussian random variables with zero mean and variance

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 <
1 and
the harmonic regression models can be made explicit by expanding Eq. 1 in an asymptotic series in powers of
(41). The nonlimit cycle perturbation solution of
Eq. 1 to O(
2) is
|
(8)
|
where
|
(9)
|
and C0 is a constant of
integration. Equations 8 and 9 are derived in
APPENDIX A. Because
= t
lim
(t), the
asymptotic or limit cycle solution of Eq. 8 to
O(
2) is
|
(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 2
/
. 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
|
(11)
|
where constants not depending on the model parameters or
the data have been omitted,
(
) 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

2 = SN/N is the maximum likelihood
estimate of 
2, the concentrated log likelihoods
are
|
(12)
|
Let y = (yt1,...,ytn),
and for the van der Pol model let
= (st1,
·t1,
,
,
,m,C),
and s(
) = [st1(
),
...,stn(
)].
For the free-run protocol, define the regression
coefficient
FR = (µ,
) and the N × 2 design matrix
|
(13)
|
where X(t) is given in Eq. 6.
For the forced desynchrony protocol, define the regression
coefficient
FD = (µ,C1,D1,...,C4,D4), and the N × 9 design matrix
|
(14)
|
Then SN has the form
|
(15)
|
where
is either
FR or
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
, s(
) can be computed by a
Runge-Kutta algorithm (1), whereas for a given value of
,
(
)
1 and
(
) are computed using the
Kalman filter (11). The Kalman filter algorithm exploits
the Markov structure in the correlated noise model to invert
(
)
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
and
, the maximum likelihood estimate of
is
|
(16)
|
where ^ denotes the estimate. Therefore, the maximum
likelihood estimates of
and
can be computed using Newton's
procedure to maximize Eq. 12 with respect to
and
at
each iteration and computing
(
,
) from Eq. 16.
The dimension of the parameter vector in the Newton's method is 8, i.e., the dimension of
plus
.
For the harmonic regression model, we define for the free-run protocol
the regression coefficient
HR = (A1,B1,A2,B2)t
and the N × 4 design matrix
|
(17)
|
If we set Z1(
) = (ZHR:ZFR),
Z2(
) = (ZHR:ZFD),
1* = (
HR:
FR), and
2* = (
HR:
FD), then for
the harmonic regression models SN is
|
(18)
|
where i = 1 and 2 denote, respectively,
the free-run and forced desynchrony protocols. Given estimates of
and
, the maximum likelihood estimate of
i* is
|
(19)
|
For the harmonic regression model, the concentrated
log-likelihood in Eq. 12 is maximized by using Newton's
procedure to find
and
,
(
)
1 and
(
) are computed as above by the Kalman filter, and
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
to O(
) and of
to
O(
2) are
|
(20)
|
By the Theorem in Ref. 6 and Proposition 3 of Ref. 12, the approximate variances of
,
, and
are
|
(21)
|
where
|
(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 |
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
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 (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.
|
|
With the exception of
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
being <0.10. Subjects 1, 2, and 3 have estimates of
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
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

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(
) 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
, 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
.
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
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:
(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
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 |
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
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
, 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
> 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
,
, and
are confounded.
This confounding is suggested by the fact that six of the twelve
subjects studied had values of
0.34, a range that makes the van
der Pol model less consistent with a weakly nonlinear oscillator. The
large values of
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
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 |
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
K+1
asymptotic series expansion of the solution to Eq. 1 is
assumed to have the form
|
(A1)
|
K = 1,2,..., where
(t) and
(t) are defined by the differential equations
|
(A2)
|
where
= 2
/
, and for each k,
uk(
,
) satisfies
|
(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(
3) are
|
(A4)
|
Substituting from Eq. A.4 into Eq. 1 gives
|
(A5)
|
and
|
(A6)
|
The solution to O(
2) is
obtained by equating the coefficients of the
terms in Eqs.
A5 and A6 subject to the conditions in Eq. A3. This yields
|
(A7)
|
To satisfy Eq. A3, we require
1 = 0, A1 as given in
Eq. A7b, and u1 to be
|
(A8)
|
When Eq. A7a is used, the
solutions for
(t) and
(t) to
O(
2) in Eq. A4 are
|
(A9)
|
where C is a constant of integration.
Substituting Eqs. A8 and A9 into Eq. A1 yields the O(
2) solution given in
Eq. 8.
 |
APPENDIX B |
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
|
(B1)
|
for i,j = 1,...,7, and where
k is any non-van der Pol parameter. It
follows from Eq. 2 that
vtn/
i = 
stn/
i. To
compute
stn/
i,
let
|
(B2)
|
Substituting Eq. B2 into Eq. 3, and
interchanging the order of differentiation with respect to time, and
differentiation with respect to
i,
give for each
i the auxiliary differential
equation system (2)
|
(B3)
|
whose initial conditions are
|
(B4)
|
because
1 = st1 and
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 |
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