A dynamic model of circadian rhythms in rodent tail skin temperature for comparison of drug effects
© Girbig et al.; licensee BioMed Central Ltd. 2012
Received: 19 August 2011
Accepted: 5 January 2012
Published: 5 January 2012
Menopause-associated thermoregulatory dysfunction can lead to symptoms such as hot flushes severely impairing quality of life of affected women. Treatment effects are often assessed by the ovariectomized rat model providing time series of tail skin temperature measurements in which circadian rhythms are a fundamental ingredient. In this work, a new statistical strategy is presented for analyzing such stochastic-dynamic data with the aim of detecting successful drugs in hot flush treatment. The circadian component is represented by a nonlinear dynamical system which is defined by the van der Pol equation and provides well-interpretable model parameters. Results regarding the statistical evaluation of these parameters are presented.
The rapid drop in hormone production during the climacteric implies a variety of changes in the physiology of the female body. This can lead to hot flushes which occur in around 70% of all women in Europe and North America . Hot flushes are characterized by a sudden sensation of intense heat, accompanied by flushing of certain peripheral skin parts as well as intensive sweating. The pathophysiology of hot flushes is not yet completely understood, but their occurrence is assumed to originate in disturbances of the thermoregulatory processes in the hypothalamus, which acts as the body's thermostat .
So far, there have been limited efforts in developing sophisticated techniques based on characteristic properties of the rhythms when evaluating the effects of a tested substance on circadian rhythms. Important characteristics are (1) a dynamic amplitude function monitoring the adaptation of oscillations to new treatments and (2) the capability of the circadian pacemaker to maintain limit cycles after sufficient adaptation to the treatment.
It should be noted that circadian rhythms in the body core have been subject to intensive research in the last few decades, and appropriate mathematical models have been formulated. These approaches comprise, for example, harmonic regression models based on weighted sums of multiple sines and cosines of varying frequencies , or differential equation models [11, 12].
In this work, we present an approach for analyzing circadian rhythms in tail skin temperature, which is based on a nonlinear dynamical model defined by the van der Pol equation. Originally proposed for the analysis of circadian rhythms within body core temperature data , it was mainly applied for studying the human circadian clock and its reaction to external light stimuli . It must be pointed out that, in addition, this model exhibits some properties making it well-suited for the investigation of the tail skin temperature of the rat. In particular, it describes sinusoidal oscillations in body temperature weighted by a continuous amplitude function, which allows monitoring of rhythm adaptations to a new treatment.
Here, the van der Pol model is fundamental for analyzing time-dependent amplitudes. The model parameters intuitively provide interpretable measures of the properties characterizing the oscillation. For example, they contain information about the start amplitude, the final amplitude at the end of a treatment period and about the individual pace of amplitude adaptation in each subject.
We apply the van der Pol equation to investigate the effects of the natural estrogen 17β Estradiol (E2) and of the synthetic hormonally active steroid Tibolone on circadian rhythms in the tail skin temperature. Both substances are known to reduce hot flushes in affected women successfully  and can thus be used as references to evaluate the capability of the proposed methods for detecting treatment effects.
In contrast to the experiments described in literature [3–5, 7, 8], measurements were started directly after ovariectomy in order to prevent desensitization of estrogen receptors due to the OVX-induced deprivation of the cognate ligand.
1. In the estrogen phase (P1), animals were treated with E2 for 6 days in order to simulate the physiological conditions prior to OVX. This allowed circadian rhythms to be maintained after the elimination of estrogen production. Their magnitude serves as a reference of the animals' original rhythms prior to OVX.
2. In order to ensure complete elimination of circulating E2, no medication was administered for 6 days (vehicle phase, P2). Nevertheless, a vehicle was applied each day to keep habituation to the experimental procedure. The aim of this phase was to guarantee that all effects observed in the subsequent test phase were evoked by the test substance.
3. For the treatment phase (P3), animals were randomized to a test group (treatment by Tibolone), and a control group (treatment by E2), each comprising 10 animals. In each group, treatment was administered for a period of 10 days.
Circadian rhythms in body temperature can be studied by dynamic models based on the van der Pol equation . So far, applications have mainly been focused on the analysis of measurements in the human body core [11, 12]. In doing so, reactions of the circadian pacemaker in the hypothalamus to external light stimuli or changes in light-/dark-rhythms are investigated in order to get deeper insights about the underlying physiological processes. The application of the model to circadian rhythms in tail skin temperature requires some adjustments because of the different properties of human body core temperature and rat tail skin temperature and because of the different experimental protocols.
Measurements of body core temperature usually have far less variability than tail skin temperature recordings since the latter are strongly influenced by thermoregulatory-induced variations in vasodilatation. These can result in high frequency perturbations with peaks that often exceed the circadian amplitude making the intrinsic circadian component difficult to detect. Furthermore, tail skin temperature is more sensitive to changes in ambient temperature (evoked for example by rats sitting on their tails) than the well-tempered body core.
with a sampling interval Δt of 3 minutes. y n is the n-th temperature measurement, where we start with measurement 0. Furthermore, c(t) = c 0 + c 1 · t with c 0 > 0 and c 1 ∈ ℝ describes a general baseline and a linear trend, while e n is a random error term for each n.
The circadian component x(t) needs to incorporate the most important properties of circadian rhythms in the living organism, which are (1) the capability to be sustained without external stimuli and (2) a dynamic amplitude function to monitor the adaptation of the extent of circadian oscillation to different treatment modes in different experimental phases. Similar to the circadian pacemaker maintaining steady oscillations in the undisturbed organism, the model should incorporate oscillations approaching a constant amplitude.
is well-suited for taking these properties into account [11, 13]. Here, is the frequency of oscillation (with τ describing the circadian period), and ε is a flexibility parameter which describes the extent to which the oscillation deviates from the harmonic oscillator model describing sinusoidal oscillations with constant phase and amplitude. The parameter γ describes a limit cycle amplitude. The nonlinear dynamical system (2) contains a supercritical Hopf bifurcation at the bifurcation point ε = 0, which enables limit cycles for positive values of ε. The larger the value of ε, the faster the limit cycle is approached but the less it resembles a sinusoidal oscillation .
. Here, k is a constant of integration and ψ 0 describes an angular phase shift. The initial amplitude a 0 := a(0) depends on k and γ according to , showing also dependence of k on a 0 and γ. The parameters of the circadian component can therefore be summarized by the parameter vector β 0 = (a 0, γ, ε, ψ, τ).
The linear trend function c(t) is added to the oscillation in order to enable adaptation of the circadian component to the overall changes in temperature baseline. The van der Pol oscillator describes periodic oscillations with approximately sinusoidal shape that are symmetric around the baseline. Consequently, a decrease in overall amplitude would imply that temperatures decrease about the same amount during the day as they are elevated at night. Under estrogen withdrawal, however, temperature is elevated during nighttime, but remains unmodified at daytime. The linear trend accounts for this one-sided decrease in amplitudes by enabling the actual baseline of the oscillatory component to remain centered between the upper and lower peaks of the sinusoidal waves.
The error terms e n describe all additional influences on the oscillation, resulting from various causes such as movement, electrical artifacts due to the measurement device or rapid thermoregulatory processes in the tail skin. Different levels of activity during day- and nighttime suggest different variance in these time spans. Even if the variance caused by thermoregulatory processes can be assumed to be small during the quiet phase during the day, husbandry and treatment each morning can lead to excitement periods that largely increase variability of measurements. Furthermore, injection can cause reactions which result in distinct changes in tail skin temperature occurring shortly after the treatment period each morning. In order to cope with these varying scenarios, different variance parameters are assumed for the night period between 6 p.m. and 6 a.m., for the morning period between 6 a.m. and 12 p.m., and for the afternoon period between 12 p.m. and 6 p.m. The random error terms are assumed to be independently and normally distributed, where the distributions are identical within each period. This tributes to the fact that an exact modeling of the 'noise' containing observational errors and treatment effects seems to be very complicated. The assumptions leading to usual least square estimates provide reasonable results as shown below.
where T day 1 is the period of points in time of 6 - 12 a.m. (distortions due to husbandry and treatment may happen), T day 2 is the period between 12 a.m. and 6 p.m. (quiet period, least variance), and T night is the nighttime period between 6 p.m. and 6 a.m., characterized by increased activity. Thus, model (1) depends on the vector containing a total of 10 parameters (note that the 5 parameters of the van der Pol model are comprised in β 0 ).
In order to improve performance and runtime of the optimization, model fit was performed by pseudo-maximum likelihood estimation. In this context, a subset of the model parameters was estimated prior to each step of the optimization iteration. In a first step, the linear trend parameters (mean c 0 and slope c 1) were estimated by linear regression for each phase. These values stayed fixed in the subsequent optimization. Variance parameters for the residuals were derived from the empirical variances of the residuals in each step of the iteration.
Optimization was then performed only for the remaining 4 parameters contained in β 0 . Therefore, instead of searching for a global minimum of the 10-dimensional negative log-likelihood function, the minimum of this function with respect to a subset of these parameters was assessed, while the other parameters were fixed. Optimization was performed in R, using the BFGS quasi-Newton method with box constraints .
Initial values for start amplitude a 0 were obtained by fitting a harmonic regression model y(t) = A 1 cos(ωt) + A 2 sin(ωt) to the observations of the first day in each experimental phase, as described in the literature  (with experimental phases being denoted as P1 - P3 as described in the The Experiments in the Methods section). Starting values were then given by . Starting values for the limit cycle amplitude γ were obtained in similar manner from the last day of each experimental phase. Since aftereffects of surgery can prevent regular rhythms at day one of measurements, but surgery-induced hypothermia can lead to unrealistically large estimates of a 0, small start amplitudes were promoted by starting optimization with initial value a 0 = 0.5 in P1. Upper boundaries for parameter a 0 (or γ) were computed as half of the maximum span of values occurring during the first (or last) day in each phase.
In contrast to previous studies where limit cycle amplitudes had been determined prior to the start of the actual experiment for each individual , here, the true limit cycle amplitudes were unknown and had to be estimated together with the other model parameters simultaneously. To ensure realistic estimates for γ, we had to bound ε from below. Otherwise, extremely large estimates for γ could occur, which were then compensated by nearly vanishing values of ε. To this end, we restricted ε to a lower boundary of 0.2. In those cases where no dynamical effects on the amplitudes occurred, i.e. due to unsuccessful treatment in P3, estimates of start and limit cycle amplitude were roughly similar. Consequently, the lower boundary on ε could not artificially induce amplitude growth in those cases. The upper boundary for ε was set to 1 in order to restrict the deviation of the perturbation solution from the true van der Pol oscillator to small values. The phase shift ψ was constrained to the interval [-π,π].
When subjecting the initial values of each parameter to small normally distributed perturbations of variance 0.1, the average absolute differences between the corresponding parameter estimates and the estimates obtained using unperturbed initial values lay between 6.2 · 10-6 (parameter ε) and 2.5 · 10-4 (parameter a 0). This indicates that (1) the large fluctuations in the data apparently prevent the optimization routine from finding the global minimum, but (2) that nevertheless, the model fit should be sufficiently robust against small alterations in the initial values, so that changes in starting value do not drastically influence the result of the fit.
Results and discussion
Exemplary illustration of the model fit
Examples of estimated van der Pol parameters.
Examples of estimated trend and variance parameters
In both animals, oscillations vanish after surgery, with the limit cycle being approached after day two (E2 group animal) or day three (Tibolone group animal) in P1. (For this and the following, compare (a), (b) and (c) in Figures 4 and 5) Amplitudes quickly decrease to a new constant level in P2 and rise again after substance treatment in P3. In each phase, limit cycle amplitudes γ closely coincide with the final values of the estimated amplitude function (a(T)), indicating that the limit cycles amplitude provides a good representation of the amplitude at the end of experiments. For both animals, limit cycle amplitudes approximately double from P2 to P3, indicating positive effects of each substance on the extent of circadian oscillation.
The reduction in the decrease of nighttime temperatures in P2, as well as its subsequent restoration in P3 are also reflected by positive or negative trends in these different experimental phases, respectively.
Estimated van der Pol parameters
Estimated van der Pol parameters
0.65 ± 0.32
0.52 ± 0.26
0.24 ± 0.19
3.97 ± 1.48
3.92 ± 1.40
2.48 ± 0.98
0.29 ± 0.25
0.40 ± 0.34
1.58 ± 0.55
1.58 ± 0.55
1.43 ± 0.61
0.31 ± 0.26
0.16 ± 0.36
2.70 ± 1.00
2.70 ± 1.00
0.57 ± 0.21
0.56 ± 0.34
0.25 ± 0.28
3.55 ± 0.75
3.52 ± 0.71
3.10 ± 1.37
0.36 ± 0.34
0.34 ± 0.37
1.61 ± 0.55
1.61 ± 0.55
1.45 ± 0.67
0.20 ± 0.00
0.05 ± 0.21
3.29 ± 0.99
3.29 ± 0.99
Estimated trend and variance parameters
Estimated estimated trend and variance parameters
29.53 ± 1.68
(0.31 ± 1.99) · 10-4
9.87 ± 3.69
4.06 ± 1.28
10.39 ± 2.76
30.35 ± 1.12
(1.37 ± 1.39) · 10-4
9.56 ± 3.25
4.91 ± 3.55
11.17 ± 2.40
30.71 ± 0.74
(-1.46 ± 0.62) · 10-4
10.22 ± 3.71
3.55 ± 0.90
9.75 ± 2.07
30.02 ± 1.01
(1.14 ± 2.85) · 10-4
9.83 ± 2.78
5.25 ± 1.79
15.07 ± 5.24
30.23 ± 1.86
(1.53 ± 2.58) · 10-4
10.78 ± 2.51
5.85 ± 1.33
11.97 ± 2.74
30.86 ± 1.28
(-3.15 ± 1.15) · 10-4
14.77 ± 3.50
4.57 ± 1.34
11.28 ± 2.68
Effects of the treatment in P3 on the amplitude parameters
After successful treatment, we expect amplitudes to increase from P2 to P3. When searching for treatment effects, we thus focus on the limit cycle amplitudes parameter γ which indicates the state of the system at the end of these experimental phases.
Comparing amplitude reconstruction between the E2 and the Tibolone group
In order to compare the performances of both treatments, we computed the ratios of limit cycle amplitudes between P3 and P1 for all animals i = 1, ..., 10 of each group k (k either refers to E2 orTibolone). These ratios provided a measure for the proportion of the natural amplitude values, as observed at the end of P1, that have been reconstructed during treatment in P3.
Usually, telemetric data analysis is based on simple graphical methods to illustrate time series, with compounds being assessed subjectively afterwards. Sometimes, the information within the huge amount of longitudinal measurements is condensed to a few mean values and then, statistical inference is performed by using simple methods like a t-test to compare groups. In doing so, mean temperatures are calculated, for example, for successive 30 minute intervals , entire nighttime periods  or time spans of several hours after substance injection [4, 8].
Here, we propose a new statistical strategy to analyze stochastic dynamic data in order to detect promising compounds for hot flush treatment using the whole time series. Our methodology is based on a nonlinear dynamical system defined by the van der Pol equation. The resulting differential equation model consists of sensible and well-interpretable model parameters to fit and to analyse tail skin temperature measurements including circadian rhythms. The mechanism presented here is a modified version of the models for human body core temperature series [11–13]. As an improvement to these models, an overall trend in temperature baseline has been included, and, most importantly, separate day- and nighttime variances are used to account for varying levels of activity and excitement at different times of the day. Moreover, common approaches model the circadian component by sinusoidal oscillations with constant amplitudes (i.e. Cosinor and harmonic regression models ). The fundamental advantage of our methodology is to involve a dynamic amplitude function which captures the time-dependent changes in amplitude during the treatment phase.
Two ideas for slight improvements for future work are suggested in the following. The first is related to the fact that the model fit was performed by pseudo-maximum likelihood fitting based on the assumption of independent residuals. This assumption does not represent true characteristics of the residual data, since high frequency oscillations are mainly based on short-term thermoregulatory processes, and thus measurements are correlated. Consequently, knowledge of short-term dependencies could provide more realistic models. To this end, the use of ARMA (auto-regressive moving average) processes has been suggested to model the thermoregulatory component in body core temperature .
Secondly, substance effects have only been investigated with respect to their implications on amplitude parameters. However, other van der Pol model parameters exist which describe different characteristics of circadian rhythms, e.g. the amplitude flexibility ε, the phase shift ψ or the differences between start amplitudes and final amplitudes of a certain experimental phase. The influence of different treatment modes on these parameters have to be investigated. It could, for example, be examined whether significant differences in these types of parameters occur in the E2 group. Observed effects could then be compared to those effects probably occurring in test groups of various substances.
We gratefully acknowledge the technical assistance of Tanja Lehmann, Andrea Trieller und Anita Griephan of Female Fertility Control & MM, Bayer HealthCare, Berlin.
- Andrikoula M, Prelevic G: Menopausal hot flushes revisited. Climacteric 2009, 12:3–15.PubMedView ArticleGoogle Scholar
- Rossmanith WG, Ruebberdt W: What causes hot flushes? The neuroendocrine origin of vasomotor symptoms in the menopause. Gynecol Endocrinol 2009, 25:303–314.PubMedView ArticleGoogle Scholar
- Berendsen H, Weekers A, Kloosterboer H: Tibolone and raloxifene on the tail temperature of oestrogen-deficient rats. Eur J Pharmacol 2001, 419:47–54.PubMedView ArticleGoogle Scholar
- Berendsen H, Kloosterboer H: Oestradiol and mirtazapine restore the disturbed tail-temperature of oestrogen-deficient rats. Eur J Pharmacol 2003, 482:329–333.PubMedView ArticleGoogle Scholar
- Sipe K, Leventhal L, Burroughs K, Cosmi S, Johnston G, Deecher D: Serotonin 2A receptors modulate tail-skin temperature in two rodent models of estrogen deciency-related thermoregulatory dysfunction. Brain Res 2004, 1028:191–202.PubMedView ArticleGoogle Scholar
- Bowe J, Li X, Kinsey-Jones J, Heyerick A, Brain S, Milligan S, O'Byrne K: The hop phytoestrogen, 8-prenylnaringenin, reverses the ovariectomy-induced rise in skin temperature in an animal model of menopausal hot flushes. J Endocrinol 2006, 191:399.PubMedView ArticleGoogle Scholar
- Opas EE, Gentile MA, Kimmel DB, Rodan GA, Schmidt A: Estrogenic control of thermoregulation in ER α KO and R β KO mice. Maturitas 2006, 53:210–216.PubMedView ArticleGoogle Scholar
- Opas EE, Scafonas A, Nantermet PV, Wilkening RR, Birzinc ET, Wilkinson H, Colwell LF, Schaeffer JM, Towler DA, Rodan GA, Schmidt A: Control of rat tail skin temperature regulation by estrogen receptor- β selective ligand. Maturitas 2009, 64:46–51.PubMedView ArticleGoogle Scholar
- Dacks PA, Rance NE: Effects of Estradiol on the Thermoneutral Zone and Core Temperature in Ovariectomized Rats. Endocrinology 2010, 151:1187–1193.PubMedView ArticleGoogle Scholar
- Brown E, Czeisler C: The statistical analysis of circadian phase and amplitude in constant-routine core-temperature data. J Biol Rhythms 1992, 7:177–202.PubMedView ArticleGoogle Scholar
- Brown E, Choe Y, Luithardt H, C C: A statistical model of the human core-temperature circadian rhythm. Am J Physiol Endocrinol Metab 2000, 279:E669-E683.PubMedGoogle Scholar
- Indic P, Brown E: Characterizing the amplitude dynamics of the human core temperature circadian rhythm using a stochastic-dynamic model. J Theor Biol 2006, 239:499–506.PubMedView ArticleGoogle Scholar
- Brown E: Identification and estimation of differential equation models for circadian data. PhD thesis, Harvard University 1987.Google Scholar
- Modelska K, Cummings S: Tibolone for postmenopausal women: systematic review of randomized trials. J Clin Endocr Metab 2002, 87:16–23.PubMedView ArticleGoogle Scholar
- Guckenheimer J, Holmes P: Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Coper-nicus; 1990.Google Scholar
- Team RDC: [http://www.R-project.org] R A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. 2011. [ISBN 3–900051–07–0]Google Scholar
- Efron B, Tibshirani R: An introduction to the bootstrap. New York: Chapman & Hall; 1994.Google Scholar
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.