 Research
 Open Access
A dynamic model of circadian rhythms in rodent tail skin temperature for comparison of drug effects
 Dorothee Girbig^{1_93},
 Karsten Keller^{2_93},
 Katja Prelle^{3_93},
 Vladimir Patchev^{4_93},
 Richardus Vonk^{5_93}Email author and
 BerndWolfgang Igl^{5_93}Email author
https://doi.org/10.1186/17403391101
© Girbig et al.; licensee BioMed Central Ltd. 2012
 Received: 19 August 2011
 Accepted: 5 January 2012
 Published: 5 January 2012
Abstract
Menopauseassociated 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 stochasticdynamic 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 wellinterpretable model parameters. Results regarding the statistical evaluation of these parameters are presented.
Background
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 [1]. 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 [2].
The ovariectomized rat model has been widely applied for the evaluation of substance effects for treating menopausal syndromes [4–8] and for investigating the underlying pathophysiology [9].
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 [10], 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 [13], it was mainly applied for studying the human circadian clock and its reaction to external light stimuli [12]. It must be pointed out that, in addition, this model exhibits some properties making it wellsuited 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 timedependent 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 [14] and can thus be used as references to evaluate the capability of the proposed methods for detecting treatment effects.
Methods
The Experiments
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 OVXinduced deprivation of the cognate ligand.
1. In the estrogen phase (P_{1}), 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, P_{2}). 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 (P_{3}), 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.
Mathematical Model
Circadian rhythms in body temperature can be studied by dynamic models based on the van der Pol equation [13]. 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/darkrhythms 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 thermoregulatoryinduced 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 welltempered body core.
with a sampling interval Δt of 3 minutes. y _{ n } is the nth 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 wellsuited 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 [15].
[12]. 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 onesided 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 }).
Model Fit
Here, y _{ obs } is the N dimensional vector of observed measurements, and is the vector of values computed from the estimated model parameters , with .
The covariance matrix contains the estimated individual variances for one of the three previously defined time periods on its diagonal, and all of its offdiagonal entries are zero.
In order to improve performance and runtime of the optimization, model fit was performed by pseudomaximum 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.
The period of oscillation was restricted to hours to account for the fixed lightdark cycle to which the animals were subjected.
Optimization was then performed only for the remaining 4 parameters contained in β _{ 0 }. Therefore, instead of searching for a global minimum of the 10dimensional negative loglikelihood 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 quasiNewton method with box constraints [16].
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 [10] (with experimental phases being denoted as P_{1}  P_{3} 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 surgeryinduced 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 P_{1}. 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 [12], 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 P_{3}, 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 van der Pol parameters  

E2  Tibolone  









 
P_{1}  0.5000  0.4969  0.0100  4.6416  4.6416  0.5000  0.3753  0.2452  3.6932  3.6931 
P_{2}  2.7760  0.2000  0.4055  1.2618  1.2620  3.8264  0.2000  0.1275  1.4968  1.4971 
P_{3}  0.5000  0.4552  0.1423  2.9695  2.9695  0.9150  0.2000  0.2483  3.8367  3.8363 
Examples of estimated trend and variance parameters
Examples of estimated trend and variance parameters  

E2  Tibolone  









 
P_{1}  28.3782  0.0001  11.0133  3.9545  5.0698  29.7567  0.0001  8.9852  3.4665  7.7807 
P_{2}  30.0259  0.0002  7.2154  4.4644  8.1160  28.7178  0.0005  12.5643  4.2927  8.5226 
P_{3}  30.4638  0.0002  9.3479  3.9833  6.6188  31.8452  0.0005  17.5731  3.4213  7.9173 
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 P_{1}. (For this and the following, compare (a), (b) and (c) in Figures 4 and 5) Amplitudes quickly decrease to a new constant level in P_{2} and rise again after substance treatment in P_{3}. 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 P_{2} to P_{3}, indicating positive effects of each substance on the extent of circadian oscillation.
The reduction in the decrease of nighttime temperatures in P_{2}, as well as its subsequent restoration in P_{3} are also reflected by positive or negative trends in these different experimental phases, respectively.
Estimated van der Pol parameters
Estimated van der Pol parameters  

Group  Experimental phase 





E2  P_{1}  0.65 ± 0.32  0.52 ± 0.26  0.24 ± 0.19  3.97 ± 1.48  3.92 ± 1.40 
P_{2}  2.48 ± 0.98  0.29 ± 0.25  0.40 ± 0.34  1.58 ± 0.55  1.58 ± 0.55  
P_{3}  1.43 ± 0.61  0.31 ± 0.26  0.16 ± 0.36  2.70 ± 1.00  2.70 ± 1.00  
Tibolone  P_{1}  0.57 ± 0.21  0.56 ± 0.34  0.25 ± 0.28  3.55 ± 0.75  3.52 ± 0.71 
P_{2}  3.10 ± 1.37  0.36 ± 0.34  0.34 ± 0.37  1.61 ± 0.55  1.61 ± 0.55  
P_{3}  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  

Group  Experimental phase 





E2  P_{1}  29.53 ± 1.68  (0.31 ± 1.99) · 10^{4}  9.87 ± 3.69  4.06 ± 1.28  10.39 ± 2.76 
P_{2}  30.35 ± 1.12  (1.37 ± 1.39) · 10^{4}  9.56 ± 3.25  4.91 ± 3.55  11.17 ± 2.40  
P_{3}  30.71 ± 0.74  (1.46 ± 0.62) · 10^{4}  10.22 ± 3.71  3.55 ± 0.90  9.75 ± 2.07  
Tibolone  P_{1}  30.02 ± 1.01  (1.14 ± 2.85) · 10^{4}  9.83 ± 2.78  5.25 ± 1.79  15.07 ± 5.24 
P_{2}  30.23 ± 1.86  (1.53 ± 2.58) · 10^{4}  10.78 ± 2.51  5.85 ± 1.33  11.97 ± 2.74  
P_{3}  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 P_{3} on the amplitude parameters
After successful treatment, we expect amplitudes to increase from P_{2} to P_{3}. 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.
Mean differences between vehicle and treatment phase
Mean differences between vehicle and treatment phase  

Treatment group 
 CI 
E2  1.1265  ( ∞,  0.8247) 
Tibolone  1.6807  ( ∞,  1.3224) 
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 P_{3} and P_{1} 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 P_{1}, that have been reconstructed during treatment in P_{3}.
Differences between the amounts of amplitude reconstruction in the E2 and the Tibolone group
Differences between the amounts of amplitude reconstruction in the E2 and the Tibolone group  

Estimator 
 CI  p _{ 0.2 } 
θ ^{(E2) } θ ^{(TI)}  0.2294  ( ∞, 0.0078)  0.0031 
Conclusions
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 ttest to compare groups. In doing so, mean temperatures are calculated, for example, for successive 30 minute intervals [5], entire nighttime periods [6] 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 wellinterpretable 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 [10]). The fundamental advantage of our methodology is to involve a dynamic amplitude function which captures the timedependent 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 pseudomaximum 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 shortterm thermoregulatory processes, and thus measurements are correlated. Consequently, knowledge of shortterm dependencies could provide more realistic models. To this end, the use of ARMA (autoregressive moving average) processes has been suggested to model the thermoregulatory component in body core temperature [11].
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.
Declarations
Acknowledgements
We gratefully acknowledge the technical assistance of Tanja Lehmann, Andrea Trieller und Anita Griephan of Female Fertility Control & MM, Bayer HealthCare, Berlin.
Authors’ Affiliations
References
 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 oestrogendeficient rats. Eur J Pharmacol 2001, 419:47–54.PubMedView ArticleGoogle Scholar
 Berendsen H, Kloosterboer H: Oestradiol and mirtazapine restore the disturbed tailtemperature of oestrogendeficient 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 tailskin temperature in two rodent models of estrogen deciencyrelated thermoregulatory dysfunction. Brain Res 2004, 1028:191–202.PubMedView ArticleGoogle Scholar
 Bowe J, Li X, KinseyJones J, Heyerick A, Brain S, Milligan S, O'Byrne K: The hop phytoestrogen, 8prenylnaringenin, reverses the ovariectomyinduced 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 constantroutine coretemperature 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 coretemperature circadian rhythm. Am J Physiol Endocrinol Metab 2000, 279:E669E683.PubMedGoogle Scholar
 Indic P, Brown E: Characterizing the amplitude dynamics of the human core temperature circadian rhythm using a stochasticdynamic 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. Copernicus; 1990.Google Scholar
 Team RDC: [http://www.Rproject.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
Copyright
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.