## 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

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
Equation 1where *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*y*
_{t1},…,*y*
_{tn
}are collected in the interval [0,*T*], where *t _{n} = nΔt*,

*n*= 1,2,…,

*N*, and

*N*Δ

*t*=

*T*. For each protocol, the core-temperature measurement may be expressed as Equation 2where μ is the core-temperature mean,

*s*

_{tn }is the circadian oscillation,

*x*

_{tn }is the evoked effect on core-temperature of the subject's activity level, and

*v*

_{tn }is the fluctuation in core-temperature measurements due to the body's thermoregulatory response. The variable

*s*

_{tn }is represented as the solution to the modified van der Pol equation, defined in Ref. 31 as Equation 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,

*s*

_{tn }is represented as Equation 4where

*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
depends on the circadian protocol. For the forced desynchrony protocol with a 28-h light-dark cycle, the regular, square-wave shape of
is well described by the harmonic regression model
Equation 5(12, 18). For the self-selected timing of activity during the free-run protocol,*x*
_{tn
} is
Equation 6This 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*x*
_{tn
} = 0 for all*t _{n}
*.

The random variable *v*
_{tn
}is a discrete sample from a continuous first-order autoregressive [AR(1)] process and is defined as
Equation 7where α^{−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
Equation 8where
Equation 9and *C*
_{0} is a constant of integration. *Equations 8
* and *
9
* are derived inappendix
. Because γ =_{t→∞}
^{lim}γ(*t*), the asymptotic or limit cycle solution of *Eq. 8
* to*O*(ε^{2}) is
Equation 10The 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*v*
_{tn
}, the log likelihood for each model has the form
Equation 11where 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*S _{N}
* is a quadratic form that depends on the data, the model, and the circadian protocol. Becauseς̂

_{η}

^{2}=

*S*/

_{N}*N*is the maximum likelihood estimate of ς

_{η}

^{2}, the concentrated log likelihoods are Equation 12Let

*y = (y*

_{t1},…,

*y*

_{tn }), and for the van der Pol model let θ = (

*s*

_{t1},

*s˙*

_{t1},τ,γ,ε,

*m,C*), and

*s(θ*) = [

*s*

_{t1}(θ), …,

*s*

_{tn }(θ)]. For the free-run protocol, define the regression coefficient ξ

_{FR}= (μ,β) and the N × 2 design matrix Equation 13where

*X(t)*is given in

*Eq. 6*. For the forced desynchrony protocol, define the regression coefficient ξ

_{FD}= (μ,

*C*

_{1},

*D*

_{1},…,

*C*

_{4},

*D*

_{4}), and the

*N*× 9 design matrix

Equation 14

Then *S _{N}
* has the form
Equation 15where ξ is either ξ

_{FR}or ξ

_{FD}, and

*Z*is either

*Z*or

_{FR}*Z*, depending on whether the protocol is free run or forced desynchrony, respectively. For a given value of θ,

_{FD}*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 Equation 16where ^ 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} = (*A*
_{1},*B*
_{1},*A*
_{2},*B*
_{2})^{t}and the *N *× 4 design matrix
Equation 17

If we set *Z*
_{1}(τ) = (*Z _{HR}
*:

*Z*),

_{FR}*Z*

_{2}(τ) = (

*Z*:

_{HR}*Z*), ξ

_{FD}_{1}

^{*}= (ξ

_{HR}:ξ

_{FR}), and ξ

_{2}

^{*}= (ξ

_{HR}:ξ

_{FD}), then for the harmonic regression models

*S*is Equation 18where

_{N}*i*= 1 and 2 denote, respectively, the free-run and forced desynchrony protocols. Given estimates of τ and α, the maximum likelihood estimate of ξ

_{i}

^{*}is Equation 19For 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
. 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
Equation 20
By the *Theorem* in Ref. 6 and *Proposition 3*of Ref. 12, the approximate variances of τ̂, γ̂, andε̂ are
Equation 21
where
Equation 22is 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.3
*A*). 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. 3
*B*). 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 and2).

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. 2
*A*) 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.2
*B*. The sum of the circadian and forced desynchrony components nicely reproduces the constructive and destructive interference seen in the original data (Fig. 2
*C*). 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.

#### 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. 5
*A*). 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
*.

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.6
*A*. The harmonic regression signal estimate shows an asymmetric waveform that is a square-wave at the crests and narrowed peaks at the troughs (Fig. 6
*B*). The estimated thermoregulatory process shows a prominent perturbation between *days 6* and *8* (Fig. 6
*C*). The estimates of τ and the harmonic regression signal amplitude are statistically significant.

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. 6
*C*). 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. 6
*B*). 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 (*HR _{2}
*) for constant-routine core-temperature analyses proposed in Ref. 8, and the second was the limit cycle approximation to the van der Pol model

*HR*in

_{vdp}*Eq. 10*. The

*HR*

_{2}model contains a fundamental and its second harmonic, whereas the

*HR*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.

_{vdp}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 *HR*
_{2} model yields an estimate of 35.83 h, and the *HR _{vdp}
*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

*HR*model than with the

_{vdp}*HR*

_{2}model. The

*HR*waveform estimate is less sinusoidal then the van der Pol model prediction, because the

_{vdp}*HR*coefficients are not constrained to satisfy the conditions of the coefficients in

_{vdp}*Eq. 10.*This analysis suggests that the

*HR*

_{2}model used in constant-routine analyses may give a circadian signal estimate different from that estimated by the van der Pol model.

## 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.EB1
*) 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.

## Acknowledgments

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

## Appendix

### 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
Equation A1
*K *= 1,2,…, where γ(*t*) and Ψ(*t*) are defined by the differential equations
Equation A2
where ω = 2π/τ, and for each *k*,*u _{k}
*(γ,ψ) satisfies
Equation A3The conditions in

*Eq. EA3*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 Equation A4 Substituting from

*Eq. A.4*into

*Eq.1*gives Equation A5 and Equation A6 The solution to

*O*(ε

^{2}) is obtained by equating the coefficients of the ε terms in

*Eqs.EA5*and

*EA6*subject to the conditions in

*Eq.EA3*. This yields Equation A7 To satisfy

*Eq. EA3,*we require ψ

_{1}= 0,

*A*as given in

_{1}*Eq. A7b,*and

*u*to be Equation A8When

_{1}*Eq. A7a*is used

*,*the solutions for ψ(

*t*) and γ(

*t*) to

*O*(ε

^{2}) in

*Eq. EA4*are Equation A9 where

*C*is a constant of integration. Substituting

*Eqs. EA8*and

*EA9*into

*Eq.EA1*yields the

*O*(ε

^{2}) solution given in

*Eq. 8*.

## Appendix

### 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
Equation B1
for *i,j* = 1,…,7, and where ξ_{k} is any non-van der Pol parameter. It follows from *Eq. 2
* that ∂v_{tn}/∂θ_{i} = − ∂s_{tn}/∂θ_{i}. To compute ∂s_{tn}/∂θ_{i}, let
Equation B2
Substituting *Eq. EB2
* 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)
Equation B3
whose initial conditions are
Equation B4
because θ_{1} =*s*
_{t1
} and*θ*
_{2} s˙_{t1}. Each auxiliary system is integrated using the Runge-Kutta method.

## 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.

- Copyright © 2000 the American Physiological Society