Macroscopic Models for Human Circadian Rhythms

Abstract

Mathematical models have a long and influential history in the study of human circadian rhythms. Accurate predictive models for the human circadian light response have been used to study the impact of a host of light exposures on the circadian system. However, generally, these models do not account for the physiological basis of these rhythms. We illustrate a new paradigm for deriving models of the human circadian light response. Beginning from a high-dimensional model of the circadian neural network, we systematically derive low-dimensional models using an approach motivated by experimental measurements of circadian neurons. This systematic reduction allows for the variables and parameters of the derived model to be interpreted in a physiological context. We fit and validate the resulting models to a library of experimental measurements. Finally, we compare model predictions for experimental measurements of light levels and discuss the differences between our model’s predictions and previous models. Our modeling paradigm allows for the integration of experimental measurements across the single-cell, tissue, and behavioral scales, thereby enabling the development of accurate low-dimensional models for human circadian rhythms.

Keywords

mathematical modeling circadian rhythms chronotype human phase response curve

Introduction

Disrupted circadian rhythms have been implicated in a vast array of both mental and physical health maladies including cancer, diabetes, addiction, depression, and sleep disorders (Abarca et al., 2002; Haus and Smolensky, 2006; Gale et al., 2011; Lewy et al., 2006). Moreover, the efficacy of health treatments has been found to vary in a circadian manner, meaning that knowledge of a patient’s circadian phase could allow for more effective treatments with reduced side effects (Lévi, 2006; Hrushesky, 1985). Therefore, it is a matter of vital importance to understand and predict human circadian rhythms.

The maintenance of healthy circadian rhythms requires them to be synchronized to environmental cycles by outside forces known as zeitgebers. In mammals, the most powerful zeitgeber is the daily light cycle (Pittendrigh and Daan, 1976). Daily light cycles are sensed in the retina and passed directly along the retino-hypothalmic tract to the master circadian clock (Meijer and Schwartz, 2003). The mammalian master circadian clock has been localized to the suprachiasmatic nucleus (SCN), a cluster of 20,000 neurons in the hypothalamus (Moore and Eichler, 1972; Stephan and Zucker, 1972). Each of these thousands of clock neurons in the SCN contains an intricate genetic feedback loop, which cycles with a period close to 24 h (Liu et al., 1997). The emergent rhythm produced collectively by these clock neurons drives peripheral circadian cycles found throughout the body (Dibner et al., 2010).

The study of human circadian rhythms has been approached using 3 principal paradigms. First, through the use of model organisms, which allow for invasive examinations of the master clock, coupling agents, and molecular details of the clock. Second, through careful laboratory-based human studies, which track circadian dynamics through the use of one of several reliable markers for the human circadian rhythm. Finally, recent technological advances have allowed for the collection of large data sets of self-reported survey data from individuals outside the laboratory setting (Walch et al., 2016; Roenneberg et al., 2007; Wirz-Justice et al., 2003). Under this paradigm, sleep-wake cycles (or other correlates of the circadian phase) have been used as a proxy measure of circadian phases.

Studies in model organisms have revealed the details of the genetic feedback loop present in each clock neuron as well as the coupling forces between the neurons that help shape the circadian waveform produced by the SCN. The circadian waveform is known to vary under a variety of conditions, including age, seasonal day length, and light history of the animal (Evans and Gorman, 2016). Laboratory studies of human circadian rhythms have produced an increased understanding of how light input is integrated into the master circadian clock as well as careful measurements of key parameters, such as the human circadian period (Czeisler et al., 1999; St Hilaire et al., 2012; Khalsa et al., 2003). Large data sets of self-reported circadian data have been used to study the variation in circadian dynamics across the human population (Walch et al., 2016; Roenneberg et al., 2007; Wirz-Justice et al., 2003). In particular, these studies have begun to uncover the prevalence of different sleep timing phenotypes or chronotypes in the human population.

From a mathematical modeling perspective, the increased knowledge of the details of circadian timekeeping has led to a divergence in the field. Detailed high-dimensional models have been created to explain and predict the molecular data sets generated through study of model organisms (DeWoskin et al., 2015; Kim and Forger, 2012). In contrast, models of human circadian data have remained phenomenological and low dimensional to avoid overfitting the available data and to reduce the computational burden of simulations. However, these 2 modeling approaches have not been integrated to allow for the exchange of knowledge between the molecular and human paradigms. To incorporate molecular data, models of human circadian dynamics need to be derived systematically from more detailed high-dimensional models of the master circadian clock.

A mathematical technique capable of supporting such a derivation was introduced by Edward Ott and Thomas Antonsen in 2008 (Ott and Antonsen, 2008). Their technique can be used to reduce a large system of heterogeneous coupled phase oscillators to a low-dimensional macroscopic model. Recently, this technique was applied to the study of circadian rhythms directly for the first time (Lu et al., 2016). However, recent evidence has shown that the accuracy of the Ott-Antonsen approach can be improved upon when describing mammalian circadian rhythms (Hannay et al., 2018). In that work, we introduced a new ansatz, the $m^{2}$ ansatz, which provides a systematic procedure for the extraction of low-dimensional macroscopic models for biological networks of coupled oscillators (Hannay et al., 2018).

In this work, we demonstrate the use of the $m^{2}$ ansatz in modeling the light response of human circadian rhythms (Hannay et al., 2018). Starting from a high-dimensional phase oscillator description of the SCN, we derive a low-dimensional model for the human circadian clock and fit the parameters to available data. The flexibility and extensible nature of our approach allow for the derivation of both a single-population and 2-population description of the core mammalian circadian clock. We validate the model parameter fits against 3 additional data sets and compare the predictions of our models against a predominant phenomenological model for the human circadian clock.

Materials and Methods

Previous Models

The most prolific models of human circadian dynamics are based on the van der Pol (VDP) limit cycle oscillator model (Kronauer et al., 1982; Wever, 1972; Forger et al., 1999). The VDP oscillator provides a low-dimensional and well-understood basis to model the overt circadian rhythms as measured by markers such as core body temperature and melatonin levels. As our knowledge of the light response of the human circadian rhythm has grown, a series of modifications have been introduced to the original models (Jewett and Kronauer, 1998; Kronauer et al., 1999; Forger et al., 1999). These progressive modifications have allowed the VDP models to continue to make accurate quantitative and qualitative predictions of the light response of the human circadian rhythm.

For our purposes, we chose to study the simplest VDP model currently used in predicting human circadian rhythms (Forger et al., 1999),

\frac{d x}{d t} = \frac{π}{12} (x_{c} + B (t))

(1a)

\frac{d x_{c}}{d t} = \frac{π}{12} {μ (x_{c} - \frac{4 x_{c}^{3}}{3}) - x [{(\frac{24}{0.99669 τ_{x}})}^{2} + k B (t)]} .

(1b)

The parameters $τ_{x}$ and µ determine the period of the oscillator and the stiffness of the oscillator, respectively. In Eq. 2.1, the variable $B (t)$ is a transformed version of the light stimulus $L (t)$ according to the Process L formalism (Kronauer et al., 1999). The dynamics of Process L add 1 dynamical dimension to the model and are given by

\frac{d n}{d t} = 60 [α (L) (1 - n) - β n]

(2a)

\hat{B} (t) = G (1 - n) α (L)

(2b)

α (L) = {(\frac{L (t)}{I_{0}})}^{p}

(2c)

B (t) = \hat{B} (1 - 0.4 x) (1 - 0.4 x_{c}) .

(2d)

We note that Eq. 2d, called the sensitivity modulation, assumes the amplitude of the transformed light input $B (t)$ varies as a function of the phase of the master circadian clock (Forger et al., 1999). To tie the limit cycle to an experimental circadian marker, the minimum value of the dynamic variable $x$ is taken to coincide with the core body temperature minimum. The parameter values we use are as specified in Serkh and Forger (2014). The accuracy and simplicity of this model have led to its application to many open questions in human chronobiology including jet lag, sleep dynamics, and the treatment of circadian disorders (Phillips et al., 2010; Gleit et al., 2013; Serkh and Forger, 2014; Walch et al., 2016).

Despite the tremendous success of the VDP formalism in modeling human circadian rhythms, the lingering phenomenological basis can limit the model’s utility. The variables ( $x_{c}, x$ ) and parameter µ do not have interpretations that can be tied to the known physiology of the master circadian clock. Thus, the circadian phase and amplitude measured experimentally can be only loosely interpreted within the VDP formalism. In addition, the lack of a physiological origin for the parameters makes the incorporation of many molecular data sets problematic. For example, the overt circadian rhythm described by the VDP models is known to be produced through the aggregation of the rhythmic contributions of thousands of coupled noisy heterogeneous biochemical oscillators. This large ensemble of coupled oscillators produces an intricate circadian waveform within the SCN that varies with age, seasonal day length, light history, and a host of other factors (Farajnia et al., 2012; Evans et al., 2012; Evans and Gorman, 2016). The VDP model formalism has only a limited use in the study of these phenomena.

In addition, large data sets have begun to shed more light on the diversity of chronotypes present in the human population (Walch et al., 2016; Roenneberg et al., 2007; Wirz-Justice et al., 2003). In understanding human chronotypes, the phenomenological basis for the VDP model could cripple the ability of researchers to incorporate differences between an individual’s circadian rhythms beyond variations in the intrinsic circadian period (Phillips et al., 2010). Variations in the circadian period can only partially explain the variation observed in human chronotypes (Duffy et al., 2001; Duffy and Czeisler, 2002). In applications, the diagnosis and treatment of circadian disorders will likely require additional knowledge beyond variations in the intrinsic circadian period (Duffy and Czeisler, 2002).

Computing

Model simulations were run using a custom C++ library, employing a variable step-size fourth-order Runge-Kutta explicit solver for the ordinary differential equation integration. Light schedules for the fitting and validation were reproduced in silico from the various experimental protocols (Khalsa et al., 2003; St Hilaire et al., 2012; Czeisler et al., 1989; Zeitzer et al., 2000). Data points for each experiment were digitized from the literature. The model was entrained to a regular light schedule for 50 days prior to initialization of the experimental protocols. To mimic the experimental circadian phase determination, the core body temperature crossing times (defined to be $ψ = π$ , $ψ_{v} = π$ for our models) were used to determine the phase shifts induced by the light stimulus.

To quantify the model’s adherence to the fitting data sets, we defined a least-squares cost function:

C (\vec{θ}) = \frac{1}{2} \sum_{k = 1}^{4} \sum_{j = 1}^{N_{k}} {(\frac{M_{k j} (\vec{θ}) - D_{k j}}{σ_{k}})}^{2},

where the sum over k enumerates the 4 experimental protocols used in the fitting, and the sum over j defines a mesh of phases/intensities for the experimental stimulus. The $D_{j k}$ are the experimental measurements assumed to have normally distributed errors with standard deviation $σ_{k}$ , and $\vec{θ}$ is the vector of parameter values for the model output $M (\vec{θ})$ . Data ( $D_{j}$ ) values in the cost function were generated by functions fit to the experimental data rather than the raw data points. For the type 1 resetting data, we fit a biharmonic function to the data sets (Khalsa et al., 2003; St Hilaire et al., 2012). For the type 0 data, each branch of the discontinuous phase-response curve was fit to a function of the form $a + b x + c x^{2} + d / (θ - x)^{2}$ (Czeisler et al., 1989). Given the ambiguity present in the data on the exact placement of the discontinuity, θ was allowed to vary in the range $θ \in [8.85, 10.05]$ , and the minimum distance was assumed in each comparison with the model simulations. For light-intensity response curve data, we used the 4-parameter logistic function as specified by the authors in the original work (Zeitzer et al., 2000). The standard deviation of each of the measurements ( $σ_{k}$ ) was estimated using the variation of the sample data about the fitted function. Therefore, the data points within each experimental data set were assumed to have equal noise variance, but different experimental conditions were allowed to have differing noise values.

Optimal parameters were estimated using a genetic algorithm for global optimization of the least squares cost function using a population size of 1000 and running for 200 generations. The optimal parameter sets were then selected from the final population produced by the genetic algorithm. To speed the parameter evaluations, the cost function evaluations for the genetic algorithm were conducted in parallel, resulting in a significant increase in the speed.

Finally, we implemented a Markov Chain Monte Carlo (MCMC) algorithm to explore the cost basin around our optimal parameter sets. The MCMC algorithm allows for the estimation of the posterior distribution $P (θ_{j} | D)$ for each parameter $θ_{j}$ given the data $D$ . A uniform prior was applied to each parameter with the sole exception of the period of the clock in darkness $τ_{D D}$ . For this parameter, we applied a Gaussian prior $N (24.18, 0.20)$ , which is a slightly inflated version of experimental results on this parameter (Czeisler et al., 1999). The MCMC algorithm was implemented as a Metropolis-Hasting’s walk across the cost basin. Best-fit parameter sets for both models and MCMC quantiles may be found in the supplementary materials.

Results

Derivation of a Macroscopic Model

Overcoming the limitations of the VDP modeling formalism will require moving beyond a phenomenological paradigm. In this spirit, we begin with a high-dimensional model that describes the phase ( $ϕ_{i}$ ) of each clock cell in the SCN. We assume the clock neurons are weakly coupled: that is, we make the assumption that deviations off the autonomous limit cycle induced by the coupling forces are sufficiently small to be safely ignored. In addition, we assume the coupling between the oscillators is all-to-all and is sufficiently weak that the contributions may be averaged into a coupling function with a single harmonic (Kuramoto, 1984). We allow for heterogeneity in the natural frequency $ω_{i}$ of each clock cell and include a white noise factor in time progression of each oscillator. These assumptions lead to the following model for each clock neuron $i = 1, 2 \dots N$ ,

\frac{d ϕ_{i}}{d t} = ω_{i} + \sqrt{D} η_{i} (t) + \frac{K}{N} \sum_{j = 1}^{N} \sin (ϕ_{j} - ϕ_{i} + β) + B (t) Q (ϕ_{i})

(3a)

Q (ϕ_{i}) = σ - A_{1} \sin (ϕ_{i} + β_{L 1}) - A_{2} \sin (2 ϕ_{i} + β_{L 2})

(3b)

where $η_{i}$ is a white noise process with $〈 η_{i} 〉 = 0$ , and $〈 η_{i} (t) η_{j} (t^{'}) 〉 = 2 δ (t - t^{'}) δ_{i j}$ , and $δ_{i j}$ is the Kronecker delta. In addition, we allow the individual oscillators to be heterogeneous and assume the natural frequencies of the clock neurons follow a Cauchy (Lorentzian) distribution,

g (ω) = \frac{γ}{π [(ω - ω_{0})^{2} + γ^{2})]},

(4)

where $ω_{0}$ is the median frequency and $γ$ sets the spread of the distribution about the median value. The $Q (ϕ)$ function describes the phase-response curve of the individual clock neuron to a brief light stimulus, and the coefficients $(A_{1}, A_{2})$ scale the first and second harmonic components of the single-cell phase-response curve, respectively. Finally, as in the VDP model, the time-dependent input $B (t)$ gives the form of the light stimulus after processing by the visual system. We model the processing of light information along the retino-hypothalamic tract using a slight modification in the $α (L)$ term of the Kronauer-Jewett Process L formalism (Eq. 2.1; Kronauer et al., 1999). That is, we define the transformation from the raw light input $L (t)$ to the processed light input to the circadian clock $B (t)$ as

\frac{d n}{d t} = 60 [α (L) (1 - n) - δ n],

(5a)

α (L) = \frac{α_{0} L {(t)}^{p}}{L {(t)}^{p} + I_{0}},

(5b)

B (t) = G (1 - n) α (L) .

(5c)

We note that the microscopic model specified in (Eq. 3.1) incorporates both cellular heterogeneity in period (frequency; $ω_{j}$ ) as well as the inevitable stochastic components $η_{j} (t)$ of a biochemical oscillator. We note that our single-cell model is a simple-phase oscillator and thus does not explicitly include the physiological details of the genetic feedback loop. However, a microscopic model of the form of Eq. 3.1 can be extracted from a biochemically motivated higher-dimensional limit cycle model through the mathematical technique of phase reduction (Kuramoto, 1984; Winfree, 2001). Recently, experimental advancements have allowed for the biochemical oscillations of circadian clock proteins to be recorded from intact SCN networks for the first time (Abel et al., 2016; Myung et al., 2012). These oscillations may be used to estimate the phases ( $ϕ_{i}$ ) of the genetic feedback loop of individual circadian neurons within our microscopic model (Eq. 3). Moreover, through experimental manipulations that decouple these neurons, the intrinsic frequencies $ω_{i}$ and noise strengh $D$ may be estimated (Myung et al., 2012). Moreover, by perturbing the system and recording the resynchronization, we may estimate the coupling strength $K$ and even estimate the network structure (Abel et al., 2016; Myung et al., 2015).

While the phase model coarse-grains over the biochemical details of each clock neuron, the large number of neurons in the SCN means it is still a very high-dimensional [ $N = O {(10}^{4})$ ] microscopic model. However, the phenomenon of interest occurs on the collective (macroscopic) scale. The Daido order parameters $Z_{n}$ ,

Z_{n} (t) = \frac{1}{N} \sum_{j = 1}^{N} e^{i n ϕ_{j} (t)},

(6)

describe the oscillator distribution on a collective scale. The first Daido order parameter $Z_{1}$ is of particular importance, as it is often used to measure synchrony. The phase coherence $R_{1} = | Z_{1} |$ varies over the range $[0, 1]$ , where a value of zero indicates desynchrony and $R_{1} = 1$ describes a population that is phase locked in perfect phase and frequency synchrony. The argument (polar angle) of $Z_{1}$ , $ψ = A r g (Z_{1})$ , gives the mean collective phase of the population.

We may transform the microscopic model (Eq. 3) into a system of coupled equations for the Daido order parameters (Hannay et al., 2018; Ott and Antonsen, 2008),

\begin{matrix} \frac{{\dot{Z}}_{n}}{n} = (i ω_{0} + χ) Z_{n} + \frac{K}{2} (e^{i β} Z_{1} Z_{n - 1} - e^{- i β} {\bar{Z}}_{1} Z_{n + 1}) + \\ + \frac{1}{2} (H_{1} Z_{n - 1} - {\bar{H}}_{1} Z_{n + 1}) + \frac{1}{2} (H_{2} Z_{n - 2} - {\bar{H}}_{2} Z_{n + 2}) . \end{matrix}

(7a)

χ = σ B (t) - γ - D n

(7b)

H_{1} (t) = A_{1} e^{- i β_{L 1}} B (t)

(7c)

H_{2} (t) = A_{2} e^{- i β_{L 2}} B (t),

(7d)

where barred quantities indicate the complex conjugate. This system goes on indefinitely (n = 1,2,. . .); thus, we have replaced our high-dimensional microscopic model (Eq. 3) with an infinite dimensional system. However, in examining dynamic measurements of the phase distribution of circadian clock neurons, we have recently reported the emergence of a low-dimensional relation between the Daido order parameters (Hannay et al., 2018). We found that for systems in which most of the oscillators are synchronized, all higher-order Daido order parameters are slaved to the first one $Z_{1}$ , such that $Z_{n} = Z_{1}^{n^{2}}$ . Indeed, we found this relationship is obeyed by the phase distribution of mammalian circadian neurons measured experimentally (Abel et al., 2016; Hannay et al., 2018).

Applying this moment closure and separating the real and imaginary parts for $Z_{1}$ gives the model for the human circadian clock shown in Figure 1. The 2 variables of the model have the benefit of having direct physiological interpretations: $R \in [0, 1]$ measures the collective amplitude of the oscillator population, and $ψ$ gives the mean phase of the population. The terms $L_{R}$ and $L_{ψ}$ give the impact of the light input on the amplitude and mean phase of the circadian clock, respectively. Our systematic derivation also allows the parameters of the macroscopic model to be traced back to the properties of the high-dimensional microscopic model for each phase oscillator (Eq. 3). For example, $\bar{ω}$ is the median angular frequency of the biochemical clock in each circadian neuron, $γ$ gives the dispersion in the frequency distribution, $D$ gives the noise strength, and $β, K$ specify the form of the pairwise coupling function. The remaining parameters $σ, A_{1, 2}, β_{L 1, L 2}$ specify the form of the phase-response curve of each clock neuron to a light pulse.

Figure 1.

The single-population model for the human circadian clock. Note that we have dropped the subscripts on the Kuramoto terms $Z_{1} = R e^{i ψ}$ .

Two-Population Model

Another advantage of our modeling approach is that alterations in the microscopic model may be easily incorporated and a new macroscopic model derived. To illustrate this process, we note that physiological investigations of the mammalian SCN have revealed that it may be functionally clustered into 2 principal regions: the ventral (core) and the dorsal (shell) clusters (Foley et al., 2011). Especially pertinent to the study of light on the mammalian circadian clock is the discovery that light information channeled to the SCN from the eyes along the retino-hypothalamic tract projects mainly onto the ventral cluster of oscillators (Meijer and Schwartz, 2003).

Within our formalism, we may easily generalize the derivation used for a single-population model of the SCN to include the division of neurons into the ventral and dorsal clusters. The addition of the dorsal population of oscillators yields a 5-dimensional macroscopic model for the circadian clock, which describes the collective amplitude of the ventral and dorsal populations $(R_{v}, R_{d})$ and the mean phase of each population ( $ψ_{v}, ψ_{d})$ , along with the same Process L light-processing variable $B (t)$ as the single-population model. The equations for the 2-population model are shown in Figure 2.

Figure 2.

The 2-population model for human circadian rhythms including a ventral and dorsal population of oscillators. Light input is assumed to be presented to only the ventral population. Coupling terms within and between the regions are given in the format $K_{f r o m, t o} .$

A distinct advantage of the 2-population model is that it allows for finer adjustments of the coupling forces within the SCN. Thus, the coupling strengths may be adjusted to allow for variations in the circadian waveform found under various experimental conditions (seasonal light changes, aging, light history; VanderLeest et al., 2009; Myung et al., 2012, 2015; Evans et al., 2013; Buijink et al., 2016). For example, studies in rodents have found that in longer day lengths, the ventral and dorsal populations tend to separate in phase (Myung et al., 2015). Within the 2-population model, this may be included by allowing the interregion coupling strengths $K_{v d}, K_{d v}$ to decrease for increasing day lengths.

Parameter Fitting

To build a predictive model of human circadian rhythms, we fit our model to available data on the human circadian light response. We make use of 3 experimental measurements of the human phase-response curve to light to calibrate our model (Khalsa et al., 2003; St Hilaire et al., 2012; Czeisler et al., 1989). These 3 phase-response curve studies use a similar protocol for the assessment of circadian phase shifts while varying the light stimulus applied considerably. St Hilaire et al. (2012) used a single bright-white light stimulus of 1 h in length, while Khalsa et al. (2003) employed a 6.7-h bright-light pulse. Finally, Czeisler et al. (1989) used a 3-pulse stimulus delivered over a period of 72 h. The single light pulse curves show weak or type 1 phase resetting, meaning the phase-response curve is continuous and the phase transition curve (initial phase plotted against final phase) shows an average slope of 1 (Winfree, 2001). Conversely, the 3-pulse phase-response curve shows type 0 resetting, where the phase transition curve has an average slope of zero (Czeisler et al., 1989; Winfree, 2001). Type 0 resetting is associated with the stimulus driving the system to a phase singularity and produces large phase shifts in the circadian phase (Winfree, 2001).

In addition to the phase-resetting data sets, we also make use of experimental measurements of the human light intensity dose-response curve (Zeitzer et al., 2000). These results study the effects of the differing light intensities on the magnitude of phase delays induced by a light stimulus applied during the early subjective night. The intensity response curve was found to be nonlinear and sigmoidal, with the inflection point near the light intensities typical of indoor lighting (Zeitzer et al., 2000).

To compare our models with experimental data, we define the collective phase $ψ = π$ in the single-population model and $ψ_{v} = π$ in the 2-population model, to correspond with the minimum of the core body temperature in humans. In addition, for the 2-population model, we make the assumption that the core body temperature marker is driven by the ventral SCN. We find this assumption is required for the model to provide good fits to the type 0 resetting data.

A comparison between the VDP model (Forger et al., 1999) and our models shows each of them are capable of describing the phase-response curve data well (Figs. 3 and 4). However, the alteration of the light processing we introduce allows for an improved fit to the light intensity dose-response curve in our models (Figs. 3d and 4d). In addition, we note that to achieve fits to the data for the VDP model, the authors introduced an ad hoc sensitivity modulation function (Eq. 2d), which requires the assumption of a significant variation in the light processing as a function of the circadian phase (Forger et al., 1999). In fitting our models, we find this sensitivity function is not required to describe the phase-shifting data.

Figure 3.

(a-c) Single-population model fit to 3 experimental measurements of the human phase-response curve to light and (d) to the light intensity dose-response curve. Thick curves show the results for the optimal parameter set, and the shading shows the density across the Markov Chain Monte Carlo parameter ensemble. Large dots show experimental data points (Khalsa et al., 2003; St Hilaire et al., 2012; Czeisler et al., 1989; Zeitzer et al., 2000). Dashed and dotted thick lines show results using the van der Pol model (Forger et al., 1999) for parameter values as given in Serkh and Forger (2014). The thin dashed curves show the interpolated function used in the fitting.

Figure 4.

(a-c) Two-population model fit to 3 experimental measurements of the human phase-response curve to light and (d) to the light intensity dose-response curve. Thick curves show the results for the optimal parameter set, and the shading shows the density across the Markov Chain Monte Carlo parameter ensemble. Large dots show experimental data points (Khalsa et al., 2003; St Hilaire et al., 2012; Czeisler et al., 1989; Zeitzer et al., 2000). Dashed and dotted thick lines show results using the van der Pol model (Forger et al., 1999) for parameter values as given in Serkh and Forger (2014). The thin dashed curves show the interpolated function used in the fitting.

Model Validation

To validate our parameter fits, we used the obtained parameter sets to simulate model responses to 3 additional experimental protocols. The first 2 of these data sets consider phase resetting in subjects exposed to intermittent light exposures in the phase delay (Gronfier et al., 2004) and phase advance regions (Rimmer et al., 2000) of the phase response. In the Rimmer et al. (2000) experiments, subjects were exposed to intermittent bright light exposures of 5.3- and 46-min lengths, alternating with episodes of darkness over a total of 5 h. The phase-shifting efficacy of these intermittent exposures was compared with a baseline constant bright-light exposure of 5 h. These intermittent light exposures were found to produce nearly the same magnitude of phase advances as the full light exposure, with the 90-min intermittent light schedule yielding approximately 90% of the baseline phase shift and the 25-min schedule producing 70% of the value of the baseline phase shift measured in terms of the median phase shifts of each group. For this validation data set, we find the 2-population model outperforms both the single population and the VDP model (Fig. 5a), although each of the models has the property that intermittent light exposures retain a large percentage of the phase-shifting capacity of the baseline constant light exposures.

Figure 5.

Comparison of model results for the 3 validation data sets (Gronfier et al., 2004; Chang et al., 2012; Rimmer et al., 2000), (Exp) using the single-population model (SP), 2-population model (TP), and the van der Pol model (VDP). Error bars for the single- and 2-population models are taken from the probability density estimated by the Markov Chain Monte Carlo parameter ensemble.

The second validation data set measured the effects of intermittent light exposures in the phase delay region of the phase response (Gronfier et al., 2004). Subjects were exposed to a intermittent light schedule consisting of six 15-min bright-light pulses separated by 60 min in very dim light. The phase-delaying effects of this intermittent light schedule were compared against a baseline light exposure of 6.5 h of constant bright light. For this data set, we find that each of the models captures the experimental data closely (Fig. 5b).

The third validation data set we considered is a duration-response curve, measuring the phase delays induced by bright-light exposures of different lengths (0.2, 1.0, 2.5, and 4.0 h; Chang et al., 2012). Similar to the intermittent light experiments, we find each model reproduces the qualitative results. Both the single- and two-population models were found to reproduce the mean phase shifts accurately (Fig. 5). While the VDP model matches the experimental results for shorter light pulses, its accuracy degrades for the longer light pulses, although we note that other VDP model variants have been found to match the duration-response curve more closely (Klerman and Hilaire, 2007; Fig. 5d).

Differences in Model Predictions

A major difference between the VDP model and the models we propose here lies in the assumed sensitivity modulation function (Eq. 2d) of the VDP based models. The sensitivity modulation function introduces a strong circadian phase dependence into the amplitude of the processed light input $B (t)$ presented to the circadian oscillator. The sensitivity function allows the VDP model to match the type 0 (3-pulse) phase-resetting experimental data. Functionally, the sensitivity modulation function amplifies the phase shifts in the critical region for stimuli occurring near the core body temperature minimum. This enables the VDP model to demonstrate type 0 resetting for relatively high stiffness values µ, which, in the absence of the sensitivity function, would prevent the oscillator from showing type 0 resetting. Our formalism does not require the introduction of an ad hoc sensitivity function to match the type 0 phase-resetting behavior. To compensate for the loss of the sensitivity function, the amplitude dynamics of our models differ significantly from the VDP model.

We find that the amplitude recovery rates of the models differ significantly in the absence of time cues. In both the single- and 2-population models, circadian amplitude recovery from small amplitudes takes significantly longer than is predicted by the VDP model (Fig. 6a). However, when the light entrainment cues are provided, the amplitude recovery rate speeds up considerably, such that it is comparable with the rate predicted by the VDP model (Fig. 6b). This slower amplitude recovery has been observed in laboratory treatments when participants are exposed to circadian amplitude suppressing bright-light pulses (Jewett et al., 1991, 1994). The one participant who was kept in darkness following an amplitude suppression showed little evidence of amplitude recovery after 4 circadian cycles (Jewett et al., 1994). However, in participants who received an additional light pulse following amplitude suppression, the circadian amplitude was observed to recover to typical levels within 3 circadian cycles (Jewett et al., 1994). This overestimate of the amplitude recovery rate from small amplitudes by the VDP model has been noted previously in the literature (Indic et al., 2005). However, attempts to remedy this deficiency for VDP-based models introduced higher-order terms into the dynamics of the amplitude recovery (Indic et al., 2005). By comparison, our models do not require higher-order terms in the amplitude recovery function but rather make the hypothesis of weaker coupling forces in the SCN. For small-amplitude reductions, both the VDP model and our models predict the amplitude recovery will occur quickly, in accordance with experimental results (Jewett et al., 1994).

Figure 6.

Amplitude recovery from small initial amplitudes in the models. (a) recovery of amplitude for the van der Pol (black), single-population (light gray), and 2-population (dark gray) models in darkness. (b) amplitude recovery when subjected to a regular 16:8 ld light schedule of 100-lux light following the amplitude reduction.

The difference in the amplitude stiffness also manifests itself in the entrainment of the models to regular light schedules. To study the entrainment of the models, we compute the days required to entrain to shifted light schedules in each of the models. To assess entrainment times, we compute the number of days required to entrain to within 0.1 radians or ≈22 min of the final stable entrainment angle beginning from a fine mesh of amplitude states and initial phases. The single- and 2-population models make similar predictions for the number of days required to entrain to the shifted schedules; however, the VDP model predicts significantly longer entrainment times for larger phase shifts (Fig. 7). Observing the entrainment dynamics via stroboscopic plots, we see our models entrain more quickly because of increased circadian amplitude suppression in response to large shifts in the light schedule (Fig. 7). Finally, we note that similar to the VDP model, our models demonstrate an asymmetry between the entrainment times for east and west shifts in the light schedule. Of note, we find this asymmetry in our models for experimentally measured values of the circadian clock period of ≈24.2 h (Czeisler et al., 1999; Carskadon et al., 1999). This contrasts with the mechanism suggested by a recent macroscopic model, using the Ott-Antonsen ansatz, which requires an assumption that the human circadian period exceed 24.5 h (Lu et al., 2016).

Figure 7.

Entrainment times to sudden time zone shifts in the single-population, 2-population, and van der Pol (VDP) models. Colors indicate the days required to entrain to a regular light schedule starting from the amplitude and phase indicated in the circular plot. The amplitudes are normalized between the models to allow comparison, and the orientation of the VDP model plot is reversed for the same reason. Arrows show stroboscopic snapshots of the phase and amplitude at 24-h intervals during the entrainment process.

The differences between the model’s predictions for entrainment time following a sudden shift in the light schedule should alter model predictions of light therapy prescriptions used in the treatment of circadian maladies (Serkh and Forger, 2014). In particular, the increased amplitude malleability of our models will likely have important effects on predictions of optimal light schedules for reentrainment, which typically seek to push the system toward the phase singularity to allow for a faster entrainment (Serkh and Forger, 2014).

Model Comparison in the Wild

To further examine the differences between our models and the VDP-based models, we simulated the 3 models for the Hispanic Community Health Survey data set from the National Sleep Research Resource (Patel et al., 2015; Redline et al., 2014; Zhang et al., 2018; Dean et al., 2016). The Sueño sleep ancillary study contains the daily light exposure schedules measured by a wrist-worn device (Actiwatch Spectrum) for more than 2000 participants for a period of 7 days. We used these measured light schedules as inputs [ $L (t)$ ] to the single-population, 2-population, and VDP models. To measure the differences between the model predictions, we examined the predicted time of day for the dim-light melatonin onset (DLMO) for each of the 3 models. The single- and 2-population models showed similar predictions for the vast majority of the light exposures considered: 95% of the differences in the predicted DLMO timing lie within the interval of (–0.57, 0.21) h.

Most of the light schedules also showed similar predictions between the VDP model and our models. However, we identified 516 of the 1704 light exposures (30.28%) in which the single-population and VDP model predicted that DLMO times differed by more than 1 h (Fig. 8a). We note that the large majority of these discrepancies between the models show the VDP model predicting a later average DLMO time than our models. To examine the causes for these discrepancies, we compared the average daily light exposure schedules between these groups (Fig. 8b). For the discrepancies between the models, we find that the light schedules show overall weaker light intensities. In particular, we note weaker light intensities during the evening hours, when light is expected to phase delay the circadian clock (Fig. 8b).

Figure 8.

Differences in model predictions for experimental light exposure schedules. (a) Average predicted dim-light melatonin onset times for the single-population model plotted against the van der Pol model. The schedules in which the models agree are shown as open circles, and the schedules where the model predictions differ substantially are shown as crosses. Equal predictions are shown as a diagonal solid line. (b) Average light exposure in log₁₀(lux) for the schedules in which the models agree (dashed) and the schedules in which the model predictions differ substantially (solid).

Our models implement a sharper light intensity dosage-response curve than the VDP model (Fig. 3d and Fig. 4d) in accordance with experimental data (Zeitzer et al., 2000). This leads to a prediction of smaller phase delays caused by light exposures of less than ≈100 lux in intensity. We confirmed this by entraining the 3 models to the averaged light schedule for the discrepancies. We found that the VDP model predicts a later DLMO time by (1.21, 1.10) hours for the single- and 2-population models, respectively.

Discussion

The $m^{2}$ ansatz is consistent with experimental measurements of the phase distribution of circadian neurons and allows for the systematic reduction of high-dimensional stochastic phase oscillator models to low-dimensional macroscopic models. In this work, we have shown how the $m^{2}$ ansatz may be applied to derive 2 new models for human circadian dynamics. The $m^{2}$ ansatz and the associated dimension reduction procedure are extremely extensible, enabling our models and assumptions to be updated to incorporate new experimental results.

Since our models are derived from high-dimensional phase models describing the phase of each circadian neuron in the SCN, the variables and parameters have inherent physiological interpretations that can be traced back to a high-dimensional single-cell model. This allows for easier incorporation of new experimental results and greater falsifiability than can be achieved with phenomenological models currently in use (Forger et al., 1999). Moreover, as our knowledge of population variability increases, the physiological interpretations of parameters in our models could allow for personalized models to be constructed based on the specific properties of each individual’s circadian dynamics.

To provide predictive models, we fit our macroscopic models to measurements of the human circadian light response. The parameter fits were validated against 3 additional data sets to evaluate the accuracy of the model predictions. Moreover, we highlight some key differences between our model predictions and a previous phenomenological model based on the VDP oscillator. In particular, we find that the elimination of the sensitivity modulation function used in the VDP model alters the model predictions on the amplitude recovery dynamics. Our models predict a slower amplitude recovery at smaller amplitudes in the absence of time cues, with a much faster recovery of amplitude predicted when rhythmic light input is provided. This is consistent with available experimental data on circadian amplitude recovery (Jewett et al., 1991, 1994). In addition, we find that our models show large differences from the VDP model in predictions of entrainment times. Weaker coupling in our models led to predictions of shorter entrainment times to large shifts in the light schedule than predicted by the VDP model. This difference in predictions could critically change predictions of optimal light-based chronotherapies for circadian misalignment (Serkh and Forger, 2014).

The final comparison we make is to simulate the 3 models using light exposure estimates recorded by wearable devices (Patel et al., 2015; Redline et al., 2014; Zhang et al., 2018; Dean et al., 2016). We find little differences in the predictions between the 2-population and single-population models for these light schedules. However, we find that the predictions between our models and the VDP model differ by more than 1 h for approximately 30% of these schedules. Examination of the light schedules that give rise to these discrepancies revealed that lower light intensities in the evening hours can account for the divergence observed in the model predictions.

For predictions of reentrainment times and of the response to real-world light schedules, we have not determined which models are more accurate. Additional data will need to be collected to resolve this. However, these results illustrate the important point that the current phase-response data used to fit model parameters are not sufficient to constrain these fundamental applications of human circadian models.

Interestingly, we found that fitting the 2-population model to type 0 resetting data necessitated the assumption that the temperature rhythm is primarily driven by the ventral SCN. This opposes the previous experimental results in rats based on a forced desynchrony protocol (Lee et al., 2009). The reasons for this discrepancy will be the subject of future work.

The VDP-based modeling paradigm has undoubtedly served the circadian community well over the past 50 years (Kronauer et al., 1982; Wever, 1972; Forger et al., 1999). However, the phenomenological basis for these models excludes the current understanding of SCN physiology. We have presented a modeling paradigm that allows movement from the single-cell circadian network scale to the behavioral scale in a systematic manner. This results in low-dimensional models of similar dynamical complexity as the VDP-based models. Moreover, our systematic derivation endows the variables and parameters with physiological interpretations and enables experimental data to be included at multiple scales. The models we derive and fit in this work will need to be revised and improved as we learn additional physiological details of the circadian clock. As we demonstrate with the derivation and analysis of the 2-population model, these physiological details may be included in the microscopic model and then pushed through the reduction to obtain a new macroscopic model. Future studies may include, for example, additional details on the coupling network and mechanisms between circadian clocks in the SCN. Moreover, of particular importance in future work will be the integration of mathematical models with large data sets collected by wearable devices to study chronotype variation in the human population. In this application, a simple model with physiological parameters will be crucial to tying chronotype variations to their physiological and genetic origins.

Supplemental Material

supp_materials_hannay_submission – Supplemental material for Macroscopic Models for Human Circadian Rhythms

Supplemental material, supp_materials_hannay_submission for Macroscopic Models for Human Circadian Rhythms by Kevin M. Hannay, Victoria Booth and Daniel B. Forger in Journal of Biological Rhythms

Footnotes

Acknowledgements

We thank O. Walch for useful discussions related to this work. This work was partially supported by NSF DMS-1412119 and DMS-1853506 (K.M.H. and V.B.), Human Frontiers of Science Program grant RPG 24/2012 (K.H. and D.B.F.), and NSF DMS-1714094 (D.B.F.)

Conflict of Interest Statement

The authors have no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

ORCID iD

Kevin M. Hannay

Notes

References

Abarca

Albrecht

Spanagel

(2002) Cocaine sensitization and reward are under the influence of circadian genes and rhythm. Proc Natl Acad Sci U S A 99(13):9026-9030.

Abel

Meeker

Granados-Fuentes

St John

Wang

Bales

Doyle

Herzog

Petzold

(2016) Functional network inference of the suprachiasmatic nucleus. Proc Natl Acad Sci U S A 113(16):4512-4517.

Buijink

Almog

Wit

Roethler

Olde Engberink

AHO

Meijer

Garlaschelli

Rohling

JHT

Michel

(2016) Evidence for weakened intercellular coupling in the mammalian circadian clock under long photoperiod. Plos One 11(12):e0168954.

Carskadon

Labyak

Acebo

Seifer

(1999) Intrinsic circadian period of adolescent humans measured in conditions of forced desynchrony. Neurosci Lett 260(2):129-132.

Chang

Santhi

St Hilaire

Gronfier

Bradstreet

Duffy

Lockley

Kronauer

Czeisler

(2012) Human responses to bright light of different durations. J Physiol 590(13):3103-3112.

Czeisler

Duffy

Shanahan

Brown

Jude

Rimmer

Ronda

Silva

Allan

Jonathan

, et al. (1999) Stability, precision, and near-24-hour period of the human circadian pacemaker. Science 284(5423):2177-2181.

Czeisler

Kronauer

Allan

Duffy

Jewett

Brown

Ronda

(1989) Bright light induction of strong (type 0) resetting of the human circadian pacemaker. Science 244(4910):1328-1333.

Dean

Goldberger

Mueller

Kim

Rueschman

Mobley

Sahoo

Jayapandian

Cui

Morrical

, et al. (2016) Scaling up scientific discovery in sleep medicine: the National Sleep Research Resource. Sleep 39(5):1151-1164.

DeWoskin

Myung

Belle

MDC

Piggins

Takumi

Forger

(2015) Distinct roles for GABA across multiple timescales in mammalian circadian timekeeping. Proc Natl Acad Sci U S A 112(29):E3911-E3919.

10.

Dibner

Schibler

Albrecht

(2010) The mammalian circadian timing system: organization and coordination of central and peripheral clocks. Annu Rev Physiol 72:517-549.

11.

Duffy

Czeisler

(2002) Age-related change in the relationship between circadian period, circadian phase, and diurnal preference in humans. Neurosci Lett 318(3):117-120.

12.

Duffy

Rimmer

Czeisler

(2001) Association of intrinsic circadian period with morningness–eveningness, usual wake time, and circadian phase. Behav Neurosci 115(4):895-899.

13.

Evans

Gorman

(2016) In synch but not in step: circadian clock circuits regulating plasticity in daily rhythms. Neuroscience 320:259-280.

14.

Evans

Elliott

Gorman

(2012) Individual differences in circadian waveform of Siberian hamsters under multiple lighting conditions. J Biol Rhythms 27(5):410-419.

15.

Evans

Leise

Castanon-Cervantes

Davidson

(2013) Dynamic interactions mediated by nonredundant signaling mechanisms couple circadian clock neurons. Neuron 80:973-983.

16.

Farajnia

Michel

Deboer

VanderLeest

Houben

Rohling

JHT

Ramkisoensing

Yasenkov

Meijer

(2012) Evidence for neuronal desynchrony in the aged suprachiasmatic nucleus clock. J Neurosci 32(17):5891-5899.

17.

Foley

Tong

Foley

Lesauter

Welsh

Silver

(2011) Characterization of orderly spatiotemporal patterns of clock gene activation in mammalian suprachiasmatic nucleus. Eur J Neurosci 33(10):1851-1865.

18.

Forger

Jewett

Kronauer

(1999) A simpler model of the human circadian pacemaker. J Biol Rhythms 14(6):533-538.

19.

Gale

Cox

Qian

Block

Colwell

Matveyenko

(2011) Disruption of circadian rhythms accelerates development of diabetes through pancreatic beta-cell loss and dysfunction. J Biol Rhythms 26(5) 423-433.

20.

Gleit

Diniz Behn

Booth

(2013) Modeling interindividual differences in spontaneous internal desynchrony patterns. J Biol Rhythms 28(5):339-355.

21.

Gronfier

Wright

Kronauer

Jewett

Czeisler

(2004) Efficacy of a single sequence of intermittent bright light pulses for delaying circadian phase in humans. Am J Physiol Endocrinol Metab 287(1):E174-E181.

22.

Hannay

Forger

Booth

(2018) Macroscopic models for networks of coupled biological oscillators. Sci Adv 4(8):e1701047.

23.

Haus

Smolensky

(2006) Biological clocks and shift work: circadian dysregulation and potential long-term effects. Cancer Causes Control 17(4):489-500.

24.

Hrushesky

(1985) Circadian timing of cancer chemotherapy. Science 228(4695):73-75.

25.

Indic

Forger

St Hilaire

Dean

Brown

Kronauer

Klerman

Jewett

(2005) Comparison of amplitude recovery dynamics of two limit cycle oscillator models of the human circadian pacemaker. Chronobiol Int 22(4):613-629.

26.

Jewett

Kronauer

(1998) Refinement of a limit cycle oscillator model of the effects of light on the human circadian pacemaker. J Theor Biol 192(4):455-465.

27.

Jewett

Kronauer

Czeisler

(1991) Light-induced suppression of endogenous circadian amplitude in humans. Nature 350:59-62.

28.

Jewett

Kronauer

Czeisler

(1994) Phase-amplitude resetting of the human circadian pacemaker via bright light: a further analysis. J Biol Rhythms 9:295-314.

29.

Khalsa

SBS

Jewett

Cajochen

Czeisler

(2003) A phase response curve to single bright light pulses in human subjects. J Physiol 549(3):945-952.

30.

Kim

Forger

(2012) A mechanism for robust circadian timekeeping via stoichiometric balance. Mol Syst Biol 8(1):630.

31.

Klerman

Hilaire

(2007) Review: on mathematical modeling of circadian rhythms, performance, and alertness. J Biol Rhythms 22(2):91-102.

32.

Kronauer

Czeisler

Pilato

Moore-Ede

Weitzman

(1982) Mathematical model of the human circadian system with two interacting oscillators. Am J Physiol 242(1):R3-R17.

33.

Kronauer

Forger

Jewett

(1999) Quantifying human circadian pacemaker response to brief, extended, and repeated light stimuli over the phototopic range. J Biol Rhythms 14(6):500-515.

34.

Kuramoto

(1984) Chemical Oscillations, Waves, and Turbulence. Vol. 19. Mineola, NY: Dover.

35.

Lee

Swanson

de la Iglesia

(2009) Circadian timing of REM sleep is coupled to an oscillator within the dorsomedial suprachiasmatic nucleus. Curr Biol 19(10):848-852.

36.

Lévi

(2006) Chronotherapeutics: the relevance of timing in cancer therapy. Cancer Causes Control 17(4):611-621.

37.

Lewy

Lefler

Emens

Bauer

(2006) The circadian basis of winter depression. Proc Natl Acad Sci U S A 103(19):7414-7419.

38.

Liu

Weaver

Strogatz

Reppert

(1997) Cellular construction of a circadian clock: period determination in the suprachiasmatic nuclei. Cell 91(6):855-860.

39.

Klein-Cardeña

Lee

Antonsen

Girvan

Ott

(2016) Resynchronization of circadian oscillators and the east-west asymmetry of jet-lag. Chaos 26(9):094811.

40.

Meijer

Schwartz

(2003) In search of the pathways for light-induced pacemaker resetting in the suprachiasmatic nucleus. J Biol Rhythms 18(3):235-249.

41.

Moore

Eichler

(1972) Loss of a circadian adrenal corticosterone rhythm following suprachiasmatic lesions in the rat. Brain Res 42:201-206.

42.

Myung

Hong

DeWoskin

De Schutter

Forger

Takumi

(2015) GABA-mediated repulsive coupling between circadian clock neurons in the SCN encodes seasonal time. Proc Natl Acad Sci U S A 112(29):E3920-E3929.

43.

Myung

Hong

Hatanaka

Nakajima

De Schutter

Takumi

(2012) Period coding of Bmal1 oscillators in the suprachiasmatic nucleus. J Neurosci 32(26):8900-8918.

44.

Ott

Antonsen

(2008) Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18(3):037113.

45.

Patel

Weng

Rueschman

Dudley

Loredo

Mossavar-Rahmani

Ramirez

Ramos

Reid

Seiger

, et al. (2015) Reproducibility of a standardized actigraphy scoring algorithm for sleep in a US Hispanic/Latino population. Sleep 38(9):1497-1503.

46.

Phillips

Chen

Robinson

(2010) Probing the mechanisms of chronotype using quantitative modeling. J Biol Rhythms 25(3):217-227.

47.

Pittendrigh

Daan

(1976) A functional analysis of circadian pacemakers in nocturnal rodents I. The stability and lability of spontaneous frequency. J Comp Physiol A 106:223-252.

48.

Redline

Sotres-Alvarez

Loredo

Hall

Patel

Ramos

Shah

Ries

Arens

Barnhart

, et al. (2014) Sleep-disordered breathing in Hispanic/Latino individuals of diverse backgrounds. The Hispanic Community Health Study/Study of Latinos. Am J Respir Crit Care Med 189(3):335-344.

49.

Rimmer

Boivin

Shanahan

Kronauer

Duffy

Czeisler

(2000) Dynamic resetting of the human circadian pacemaker by intermittent bright light. Am J Physiol 279(5):R1574-R1579.

50.

Roenneberg

Kuehnle

Juda

Kantermann

Allebrandt

Gordijn

Merrow

(2007) Epidemiology of the human circadian clock. Sleep Med Rev 11(6):429-438.

51.

Serkh

Forger

(2014) Optimal schedules of light exposure for rapidly correcting circadian misalignment. PLoS Comput Biol 10(4):e1003523.

52.

St Hilaire

Gooley

Khalsa

SBS

Kronauer

Czeisler

Lockley

(2012) Human phase response curve to a 1 h pulse of bright white light. J Physiol 590(13):3035-3045.

53.

Stephan

Zucker

(1972) Circadian rhythms in drinking behavior and locomotor activity of rats are eliminated by hypothalamic lesions. Proc Natl Acad Sci U S A 69(6):1583-1586.

54.

VanderLeest

Rohling

JHT

Michel

Meijer

(2009) Phase shifting capacity of the circadian pacemaker determined by the SCN neuronal network organization. PloS One 4(3):e4976.

55.

Walch

Cochran

Forger

(2016) A global quantification of “normal” sleep schedules using smartphone data. Sci Adv 2(5):e1501705.

56.

Wever

(1972) Virtual synchronization towards the limits of the range of entrainment. J Theor Biol 36(1):119-132.

57.

Winfree

(2001) The Geometry of Biological Time. New York: Springer.

58.

Wirz-Justice

Roenneberg

Merrow

(2003) Life between clocks: daily temporal patterns of human chronotypes. J Biol Rhythms 18(1):80-90.

59.

Zeitzer

Dijk

Kronauer

Brown

Czeisler

(2000) Sensitivity of the human circadian pacemaker to nocturnal light: melatonin phase resetting and suppression. J Physiol 526(pt 3):695-702.

60.

Zhang

Cui

Mueller

Tao

Kim

Rueschman

Mariani

Mobley

Redline

(2018) The National Sleep Research Resource: towards a sleep data commons. J Am Med Inform Assoc 25(10):1351-1358.