Statistical model building and model criticism for human circadian data

Mathematical models have played an important role in the analysis of circadian systems. The models include simulation of differential equation systems to assess the dynamic properties of a circadian system and the use of statistical models, primarily harmonic regression methods, to assess the static properties of the system. The dynamical behaviors characterized by the simulation studies are the response of the circadian pacemaker to light, its rate of decay to its limit cycle, and its response to the rest-activity cycle. The static properties are phase, amplitude, and period of the intrinsic oscillator. Formal statistical methods are not routinely employed in simulation studies, and therefore the uncertainty in inferences based on the differential equation models and their sensitivity to model specification and parameter estimation error cannot be evaluated. The harmonic regression models allow formal statistical analysis of static but not dynamical features of the circadian pacemaker. The authors present a paradigm for analyzing circadian data based on the Box iterative scheme for statistical model building. The paradigm unifies the differential equation-based simulations (direct problem) and the model fitting approach using harmonic regression techniques (inverse problem) under a single schema. The framework is illustrated with the analysis of a core-temperature data series collected under a forced desynchrony protocol. The Box iterative paradigm provides a framework for systematically constructing and analyzing models of circadian data.     

A simpler model of the human circadian pacemaker

Numerous studies have used the classic van der Pol oscillator, which contains a cubic nonlinearity, to model the effect of light on the human circadian pacemaker. Jewett and Kronauer demonstrated that Aschoff's rule could be incorporated into van der Pol type models and used a van der Pol type oscillator with higher order nonlinearities. Kronauer, Forger, and Jewett have proposed a model for light preprocessing, Process L, representing a biochemical process that converts a light signal into an effective drive on the circadian pacemaker. In the paper presented here, the authors use the classic van der Pol oscillator with Process L and Jewett and Kronauer's model of Aschoff's rule to model the human circadian pacemaker. This simpler cubic model predicts the results of a three-pulse human phase response curve experiment and a two-pulse amplitude reduction study with as much, or more, accuracy as the models of Jewett and Kronauer and Kronauer, Forger, and Jewett, which both employ a nonlinearity of degree 7. This suggests that this simpler cubic model should be considered as a potential alternative to other models of the human circadian system currently available.     

The statistical analysis of circadian phase and amplitude in constant-routine core-temperature data.

Accurate estimation of the phases and amplitude of the endogenous circadian pacemaker from constant-routine core-temperature series is crucial for making inferences about the properties of the human biological clock from data collected under this protocol. This paper presents a set of statistical methods based on a harmonic-regression-plus-correlated-noise model for estimating the phases and the amplitude of the endogenous circadian pacemaker from constant-routine core-temperature data. The methods include a Bayesian Monte Carlo procedure for computing the uncertainty in these circadian functions. We illustrate the techniques with a detailed study of a single subject's core-temperature series and describe their relationship to other statistical methods for circadian data analysis. In our laboratory, these methods have been successfully used to analyze more than 300 constant routines and provide a highly reliable means of extracting phase and amplitude information from core-temperature data.     

Revised limit cycle oscillator model of human circadian pacemaker

In 1990, Kronauer proposed a mathematical model of the effects of light on the human circadian pacemaker. This study presents several refinements to Kronauer's original model of the pacemaker that enable it to predict more accurately the experimental results from a number of different studies of the effects of the intensity, timing, and duration of light stimuli on the human circadian pacemaker. These refinements include the following: The van der Pol oscillator from Kronauer's model has been replaced with a higher order limit cycle oscillator so that the system's amplitude recovery is slower near the singularity and faster near the limit cycle; the phase and amplitude of the circadian rhythm in sensitivity to light from Kronauer's model has been refined so that the peak sensitivity to light on the limit cycle now occurs approximately 4 h before the core body temperature minimum (CBTmin) and is three times as great as the minimum sensitivity on the limit cycle; the critical phase (at which type 1 phase response curves [PRCs] can be distinguished from type 0 PRCs) that occurs at CBT,n now corresponds to 0.8 h after the minimum of x (x(min) in this refined model rather than to the exact timing of x(min) as in Kronauer's model; a direct effect of light on circadian period was incorporated into the model such that as light intensity increases, the period decreases, which is in accordance with Aschoff's rule.     

Characterizing the amplitude dynamics of the human core-temperature circadian rhythm using a stochastic-dynamic model

Two measures, amplitude and phase, have been used to describe the characteristics of the endogenous human circadian pacemaker, a biological clock located in the hypothalamus. Although many studies of change in circadian phase with respect to different stimuli have been conducted, the physiologic implications of the amplitude changes (dynamics) of the pacemaker are unknown. It is known that phase changes of the human circadian pacemaker have a significant impact on sleep timing and content, hormone secretion, subjective alertness and neurobehavioral performance. However, the changes in circadian amplitude with respect to different stimuli are less well documented. Although amplitude dynamics of the human circadian pacemaker are observed in physiological rhythms such as plasma cortisol, plasma melatonin and core temperature data, currently methods are not available to accurately characterize the amplitude dynamics from these rhythms. Of the three rhythms core temperature is the only reliable variable that can be monitored continuously in real time with a high sampling rate. To characterize the amplitude dynamics of the circadian pacemaker we propose a stochastic-dynamic model of core temperature data that contains both stochastic and dynamic characteristics. In this model the circadian component that has a dynamic characteristic is represented as a perturbation solution of the van der Pol equation and the thermoregulatory response in the data that has a stochastic characteristic is represented as a first-order autoregressive process. The model parameters are estimated using data with a maximum likelihood procedure and the goodness-of-fit measures along with the associated standard error of the estimated parameters provided inference about the amplitude dynamics of the pacemaker. Using this model we analysed core temperature data from an experiment designed to exhibit amplitude dynamics. We found that the circadian pacemaker recovers slowly to an equilibrium level following amplitude suppression. In humans this reaction to perturbation from equilibrium value has potential physiological implications.     

A phase dynamics model of human circadian rhythms

Nonphotic entrainment of an overt sleep-wake rhythm and a circadian pacemaker-driving temperature/melatonin rhythm suggests existence of feedback mechanisms in the human circadian system. In this study, the authors constructed a phase dynamics model that consisted of two oscillators driving temperature/melatonin and sleep-wake rhythms, and an additional oscillator generating an overt sleep-wake rhythm. The feedback mechanism was implemented by modifying couplings between the constituent oscillators according to the history of correlations between them. The model successfully simulated the behavior of human circadian rhythms in response to forced rest-activity schedules under free-run situations: the sleep-wake rhythm is reentrained with the circadian pacemaker after release from the schedule, there is a critical period for the schedule to fully entrain the sleep-wake rhythm, and the forced rest-activity schedule can entrain the circadian pacemaker with the aid of exercise. The behavior of human circadian rhythms was reproduced with variations in only a few model parameters. Because conventional models are unable to reproduce the experimental results concerned here, it was suggested that the feedback mechanisms included in this model underlie nonphotic entrainment of human circadian rhythms.     

A novel approach to the van der pol oscillator: Natural frequency and entrainment

Abstract

The physics of the van der Pol oscillator as realized by the Meissner circuit is discussed by analogy to the beat phenomenon and by a consequent analysis of current balance. The current balance method leads to a new, very accurate equation for the dependence of the oscillator frequency on the feedback parameter. Several aspects of entrainment (existence, limited frequency range, dependence on parameters, phase shift) can be explained, too. Numerical results are presented which have been obtained by solving the homogeneous and inhomogeneous van der Pol equation with a Runge‐Kutta method.     

Control systems models for the circadian clock of the New Zealand weta, Hemideina thoracica (Orthoptera: Stenopelmatidae)

The New Zealand weta, Hemideina thoracica, is a nocturnal orthopteran insect which emerges from holes in trees or from under bark soon after sunset to forage for several hours on plant and animal material before returning to its refuge before dawn. In tests of the internal clock hypothesis it exhibits clear circadian locomotor rhythms in which the period is initially somewhat less than 24 h, but frequently spontaneously increases to over 25 h. The rhythms are entrainable by light and temperature cycles, obey Aschoff's Law and are temperature compensated. A single oscillator feedback model accounts for these basic properties of the weta clock, but does not explain a variety of examples of rhythm lability, such as day skipping, spontaneous change in period, scalloping and desynchrony typically found in the real data. To account for these characteristics the model is expanded into two linked populations of oscillators, which retain the basic properties of the simple model and in addition interact through their coupling to show the various types of free-run lability. To make these control systems models compatible with the molecular interpretation of circadian biology, each of the components in the feedback loop is matched with molecular function and structure.     

Global bifurcation structure of a Bonhoeffer-van der Pol oscillator driven by periodic pulse trains

The Bonhoeffer-van der Pol (BVP) oscillator is a valuable dynamical system model of pacemaker neurons. Isochrons, phase transition curves (PTC), and two dimensional bifurcation diagrams served to analyze the neuron's response to periodic pulse stimuli. Responses are described and explained in terms of the nonlinear dynamical system theory. An important issue in the generation of spikes by pacemaker neurons is the existence of both slow and fast dynamics in the state point's trajectory in the phase plane. It is this feature in particular that makes the BVP oscillator a faithful model of living pacemaker neurons. Comparison of the model's responses with those of a living pacemaker was based also on return maps of interspike intervals. Analyzed in detail were the complex discharges called stammering which involve interspike intervals that arise unpredictably and exhibit histograms with several modes separated by the equal intervals.Supported by Trent H. Wells Jr. Inc.     

The method of harmonic balance applied to coupled asymmetrical van der Pol oscillators for intestinal modelling

The spontaneous electrical rhythms recorded from the gastro-intestinal tract of humans and animals have been successfully modelled by an array of interconnected van der Pol oscillators. To account for asymmetry in the recorded waveforms (with particular reference to the human small intestine) an additional term in the van der Pol dynamics has been included. It is shown that the method of harmonic balance can be used to give analytical results for this asymmetrical condition. The non-linear algebraic equations are solved by hill-climbing to give values of d.c., fundamental and second harmonic amplitudes together with the entrained frequency. The results correlate well with actual measurements made on an analogue simulation by three different methods for waveshape factors of 0.1 and 1.0     

A Mathematical Model of the Human Circadian System and Its Application to Jet Lag

A mathematical model of the circadian system is described that is appropriate for application to jet lag. The core of the model is a van der Pol equation with an external force. Approximate solutions of this equation in which the external force is composed of a constant and an oscillating term are investigated. They lead to analytical expressions for the amplitude and period of free-running rhythms and for the frequency limits of the entrainment region. The free-running period increases quadratically with stiffness. Both period and amplitude depend on the value of the constant external force. The width of the range of entrainment is mostly determined by the external force, whereas the relative position of this range follows the intrinsic period of the oscillator. Experiments with forced and spontaneous internal desynchronization were evaluated using these analytical expressions, and estimates were obtained for the intrinsic period of the oscillator, its stiffness, and the external force. A knowledge of these model parameters is essential for predictions about circadian dynamics, and there are practical implications for the assessment of the adaptation after rapid transmeridian travel.     

A Mathematical Model of the Human Circadian System and Its Application to Jet Lag

A mathematical model of the circadian system is described that is appropriate for application to jet lag. The core of the model is a van der Pol equation with an external force. Approximate solutions of this equation in which the external force is composed of a constant and an oscillating term are investigated. They lead to analytical expressions for the amplitude and period of free-running rhythms and for the frequency limits of the entrainment region. The free-running period increases quadratically with stiffness. Both period and amplitude depend on the value of the constant external force. The width of the range of entrainment is mostly determined by the external force, whereas the relative position of this range follows the intrinsic period of the oscillator. Experiments with forced and spontaneous internal desynchronization were evaluated using these analytical expressions, and estimates were obtained for the intrinsic period of the oscillator, its stiffness, and the external force. A knowledge of these model parameters is essential for predictions about circadian dynamics, and there are practical implications for the assessment of the adaptation after rapid transmeridian travel.     

Revisiting spontaneous internal desynchrony using a quantitative model of sleep physiology

Early attempts to characterize free-running human circadian rhythms generated three notable results: 1) observed circadian periods of 25 hours (considerably longer than the now established 24.1- to 24.2-hour average intrinsic circadian period) with sleep delayed to later circadian phases than during entrainment; 2) spontaneous internal desynchrony of circadian rhythms and sleep/wake cycles--the former with an approximately 24.9-hour period, and the latter with a longer (28-68 hour) or shorter (12-20 hour) period; and 3) bicircadian (48-50 hour) sleep/wake cycles. All three results are reproduced by Kronauer et al.'s (1982) coupled oscillator model, but the physiological basis for that phenomenological model is unclear. We use a physiologically based model of hypothalamic and brain stem nuclei to investigate alternative physiological mechanisms that could underlie internal desynchrony. We demonstrate that experimental observations can be reproduced by changes in two pathways: promotion of orexinergic (Orx) wake signals, and attenuation of the circadian signal reaching hypothalamic nuclei. We reason that delayed sleep is indicative of an additional wake-promoting drive, which may be of behavioral origin, associated with removal of daily schedules and instructions given to participants. We model this by increasing Orx tone during wake, which reproduces the observed period lengthening and delayed sleep. Weakening circadian input to the ventrolateral preoptic nucleus (possibly mediated by the dorsomedial hypothalamus) causes desynchrony, with observed sleep/wake cycle period determined by degree of Orx up-regulation. During desynchrony, sleep/wake cycles are driven by sleep homeostasis, yet sleep bout length maintains circadian phase dependence. The model predicts sleep episodes are shortest when started near the temperature minimum, consistent with experimental findings. The model also correctly predicts that it is possible to transition to bicircadian rhythms from either a synchronized or desynchronized state. Our findings suggest that feedback from behavioral choices to physiology could play an important role in spontaneous internal desynchrony.     

Analysis of limit-cycle conditions in intercoupled relay oscillators with reference to gastrointestinal modelling

This paper presents an exact analysis of a mutually coupled relay oscillator based on a method orignated by Tsypkin. Limit-cycle frequencies and phases can be determined exactly using this method, unlike other approximate methods based on describing functions and harmonic balance techniques. A new method of exact determination of limit-cycle stability is also shown to give excellent agreement with simulation studies. Different types of intercoupling are shown to give different stability conditions, and these are discussed in relation to gastrointestinal (GI) smooth muscle modelling. GI tract electrical activity has previously been modelled using bidirectionally coupled nonlinear oscillators, and the results of the present analysis of relay oscillators is compared with other studies using van der Pol dynamics.     

On the complex dynamics of intracellular ganglion cell light responses in the cat retina

We recorded intracellular responses from cat retinal ganglion cells to sinusoidal flickering lights, and compared the response dynamics with a theoretical model based on coupled nonlinear oscillators. Flicker responses for several different spot sizes were separated in a smooth generator (G) potential and corresponding spike trains. We have previously shown that the G-potential reveals complex, stimulus-dependent, oscillatory behavior in response to sinusoidally flickering lights. Such behavior could be simulated by a modified van der Pol oscillator. In this paper, we extend the model to account for spike generation as well, by including extended Hodgkin-Huxley equations describing local membrane properties. We quantified spike responses by several parameters describing the mean and standard deviation of spike burst duration, timing (phase shift) of bursts, and the number of spikes in a burst. The dependence of these response parameters on stimulus frequency and spot size could be reproduced in great detail by coupling the van der Pol oscillator and Hodgkin-Huxley equations. The model mimics many experimentally observed response patterns, including non-phase-locked irregular oscillations. Our findings suggest that the information in the ganglion cell spike train reflects both intraretinal processing, simulated by the van der Pol oscillator, and local membrane properties described by Hodgkin-Huxley equations. The interplay between these complex processes can be simulated by changing the coupling coefficients between the two oscillators. Our simulations therefore show that irregularities in spike trains, which normally are considered to be noise, may be interpreted as complex oscillations that might carry information.To the memory of Prof. Otto-Joachim Grusser     

Coupled oscillators utilised as gait rhythm generators of a two-legged walking machine

The gait of current two-legged walking machines differs from that of humans, although the kinematic structures of these machines' legs frequently imitate human limbs. This paper presents a method of generating the trajectories of hip and knee joint angles resulting in a gait pattern similar to that of a human. For this purpose the solutions of coupled van der Pol oscillator equations are utilised. There is much evidence that these equations can be treated as a good model of the central pattern generator generating functional (also locomotional) rhythms in living creatures. The oscillator equations are solved by numerical integration. The method of changing the type of gait by changing appropriate parameter values in the oscillator equations is presented (change of velocity and trajectory of leg-ends). The results obtained enable enhanced control of twolegged walking systems by including gait pattern generators which will assume a similar role to that of biological generators.     

A model for the speed of memory retrieval

The memory retrieval process of number problems with external noise is studied with the use of the Bonhoeffer–van der Pol oscillator model. Three cell assembly responses are simulated, coding one true number and two neighboring erroneous. The time of a correct response, T _c, was averaged over statistical assemblies of numerous trials. It is demonstrated that T _c takes a minimum value for a certain noise intensity. This result correlates well with experimental data by Usher and Feingold (2000). The location of the minimum as a function of the time delay between two consecutive simulation trials is investigated.     

The circadian locomotor rhythm of hemideina thoracica (Orthoptera; Stenopelmatidae): A population of weakly coupled feedback oscillators as a model of the underlying pacemaker

A control systems model consisting of a population of weakly-coupled feedback oscillators has been developed to simulate the circadian locomotor rhythm of the insect, Hemideina thoracica (Orthoptera; Stenopelmatidae). The model is an extension of a previously published single oscillator feedback model (Gander and Lewis, 1979) which successfully simulates entrainment, phase response curves, temperature compensation and Aschoff's Rule for Hemideina activity rhythms. The population model described here has the additional properties of predicting some of the free-run period lability (Pavlidis, 1978a, b) observed in the Hemideina rhythm (Christensen and Lewis, 1982) which is unexplained by single oscillator systems. Model behaviour is compared with the experimental data derived from the insect activity rhythms.     

Comparisons of the variability of three markers of the human circadian pacemaker

A circadian pacemaker within the central nervous system regulates the approximately 24-h physiologic rhythms in sleep cycles, hormone secretion, and other physiologic functions. Because the pacemaker cannot be examined directly in humans, markers of pacemaker function must be used to study the pacemaker and its response to environmental stimuli. Core body temperature (CBT), plasma cortisol, and plasma melatonin are three marker variables frequently used to estimate the phase of the human pacemaker. Measurements of circadian phase using markers can contain variability due to the circadian pacemaker itself, the intrinsic variability of the marker relative to the pacemaker, the method of analysis of the marker, and the marker assay. For this report, we compared the mathematical variability of a number of methods of identifying circadian phase from CBT, plasma cortisol, and plasma melatonin data collected in a protocol in which pacemaker variability was minimized using low light levels and regular timing of both the light pattern and the rest/activity schedule. We hoped to assess the relative variabilities of the different physiological markers and the analysis methods. Methods were based on the crossing of an absolute threshold, on the crossing of a relative threshold, or on fitting a curve to all data points. All methods of calculating circadian phase from plasma melatonin data were less variable than those calculated using CBT or cortisol data. The standard deviation for the phase estimates using CBT data was 0.78 h, using cortisol data was 0.65 h, and for the eight analysis methods using melatonin data was 0.23 to 0.35 h. While the variability for these markers might be different for other subject populations and/or less stringent study conditions, assessment of the intrinsic variability of the different calculations of circadian phase can be applied to allow inference of the statistical significance of phase and phase shift calculations, as well as estimation of sample size or statistical power for the number of subjects within an experimental protocol.