Research
A dynamic model of circadian rhythms in rodent tail skin temperature for comparison of drug effects
 Bioinformatics Group, Max Planck Institute of Molecular Plant Physiology, Potsdam, Germany
 Institute of Mathematics, University of Lübeck, Germany
 Toxicology, Global Early Development, Bayer HealthCare, Wuppertal, Germany
 TRG Women's Healthcare, Global Drug Discovery, Bayer HealthCare, Berlin, Germany
 Global Drug Discovery Statistics and Experimental Medicine Statistics, Bayer HealthCare, Berlin, Germany
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].
An animal model suitable for assessing the effectiveness of drug candidates in hot flush treatment is the ovariectomized rat model [3]. Here, a menopauselike state is simulated by ovary removal (ovariectomy, OVX), resulting in estrogen deprivation. It is assumed that the lack of estrogens lowers the set point for the activation of heatdissipating mechanisms. As shown in Figure 1, this results in a lack of decrease in tail skin temperature during the rats' active phase at night and, consequently, flattening of the circadian oscillations in this parameter. Then, the effectiveness of a substance in compensating the effects of estrogen deprivation can be assessed by its capability to restore the amplitude of circadian rhythms in the tail skin temperature.
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
The effects of E2 and Tibolone were tested using the ovariectomized rat model. During ovariectomy, a temperature sensor was tunneled under the tail skin at a distance of 4  5 cm from the base of the tail, and a transmitter was implanted (see Figure 2 for details). Temperature values were averaged over 10 seconds for each record; records being taken every 3 minutes. A day/nightrhythm was simulated by a 'daytime period' of 6 a.m.  6 p.m. in which room light was switched on, and by a 'nighttime period' of the remaining 12 hours where the light was switched off. Each morning at 6 a.m., animal husbandryrelated activities began. Treatment was administered between 8:30  9:30 a.m. each day. In the remaining time, animals were undisturbed in order to minimize the amount of perturbing influences on the measurements.
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.
As depicted in Figure 3, the experiment was split up into three consecutive phases:
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.
We assume that the observed time series ${\left({y}_{n}\right)}_{n=0}^{N1}$ of some length N is given by the temperature values
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.
The van der Pol equation
is wellsuited for taking these properties into account [11,13]. Here, $\omega =\frac{2\pi}{\tau}$ 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 $\ddot{x}\left(t\right)+{\omega}^{2}x\left(t\right)=0$ 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].
Since no exact solution of the van der Pol oscillator is available, we currently use the secondorder perturbation solution which is given by
with the dynamic amplitude function
being a solution of the differential equation
[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 ${a}_{0}={\left(\frac{1}{4k}+\frac{1}{{\gamma}^{2}}\right)}^{\frac{1}{2}}$, 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.
The process ${\left({e}_{n}\right)}_{n=0}^{N1}$ is then defined by
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 $\beta =\left({\beta}_{0},{c}_{0},{c}_{1},{\sigma}_{day1}^{2},{\sigma}_{day2}^{2},{\sigma}_{night}^{2}\right)$ containing a total of 10 parameters (note that the 5 parameters of the van der Pol model are comprised in β_{0}).
Model Fit
A separate model was fitted to the time series obtained from each animal in each experimental phase, aiming at a maximization of the observations' likelihood. Assuming normally distributed and independent error terms e_{n}, this likelihood is provided by the density function of the multivariate normal distribution:
Here, y_{obs }is the N dimensional vector of observed measurements, and $\mathrm{y}\left(\widehat{\beta}\right)$ is the vector of values computed from the estimated model parameters $\widehat{\beta}=\left({\widehat{\beta}}_{0},{\u0109}_{0},{\u0109}_{1},{\widehat{\sigma}}_{day1}^{2},{\widehat{\sigma}}_{day2}^{2},{\widehat{\sigma}}_{night}^{2}\right)$, with ${\widehat{\beta}}_{0}=\left({\widehat{a}}_{0},\widehat{\gamma},\widehat{\epsilon},\widehat{\psi},\widehat{\tau}\right)$.
The covariance matrix $\mathrm{\Sigma}\left(\widehat{\beta}\right)$ 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.
Equation (5) can be reformulated to the negative loglikelihood, which serves as the objective function for optimization:
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 $\widehat{\tau}=24$ 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 ${a}_{0}=\sqrt{{A}_{1}^{2}+{A}_{2}^{2}}$. 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
Figures 4 and 5 show examples of the circadian components x(t), the corresponding amplitude functions a(t) and the linear trend components detected in the tail skin temperature time series of two exemplary animals. The parameter estimates obtained for each animal are provided in Table 1 (van der Pol parameters) and in Table 2 (trend and variance parameters). Note that the residuals obtained for the considered time periods are unimodally distributed with a slight difference to a normal distribution. However, analysis of the autocorrelation function revealed only very shorttime dependencies in the residuals. This is not in complete accordance with our assumptions, but does not contradict the point that our model serves as a rough description of the considered processes.
Examples of estimated van der Pol parameters  

E2  Tibolone  
$\widehat{a}\left(0\right)$  $\widehat{\epsilon}$  $\widehat{\psi}$  $\widehat{\gamma}$  $\widehat{a}\left(T\right)$  $\widehat{a}\left(0\right)$  $\widehat{\epsilon}$  $\widehat{\psi}$  $\widehat{\gamma}$  $\widehat{a}\left(T\right)$  
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  

E2  Tibolone  
${\u0109}_{0}$  ${\u0109}_{1}$  ${\widehat{\sigma}}_{day1}^{2}$  ${\widehat{\sigma}}_{day2}^{2}$  ${\widehat{\sigma}}_{night}^{2}$  ${\u0109}_{0}$  ${\u0109}_{1}$  ${\widehat{\sigma}}_{day1}^{2}$  ${\widehat{\sigma}}_{day2}^{2}$  ${\widehat{\sigma}}_{night}^{2}$  
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.
In total we fitted the three phases for 20 animals, with AIC values between 12000 and 23000, each being slightly lower than the value obtained for the corresponding fit of the harmonic regression model. Average parameter estimates obtained for both treatment groups are provided in Table 3 and 4.
Estimated van der Pol parameters  

Group  Experimental phase  $\widehat{a}\left(0\right)$  $\widehat{\epsilon}$  $\widehat{\psi}$  $\widehat{\gamma}$  $\widehat{a}\left(T\right)$ 
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 estimated trend and variance parameters  

Group  Experimental phase  ${\u0109}_{0}$  ${\u0109}_{1}$  ${\widehat{\sigma}}_{day1}^{2}$  ${\widehat{\sigma}}_{day2}^{2}$  ${\widehat{\sigma}}_{night}^{2}$ 
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.
Boxplots of the limit cycle amplitudes obtained from the time series of each treatment group (10 animals per group) are shown in Figure 6. As expected, amplitudes at the end of P_{2 }tend to be located at lower levels than those at the end of P_{1}. In P_{3}, an overall rise in amplitude parameters can be observed relative to P_{2}.
In order to illustrate amplitude reconstruction during treatment in P_{3}, Figure 7 provides boxplots of the pairwise differences in limit cycle amplitudes of vehicle and treatment phase (P_{2 } P_{3}). As expected, all pairwise difference values are located below zero, indicating increase in amplitudes from P_{2 }to P_{3 }and thus successful amplitude restoration in both groups.
In order to assess whether this increase is statistically significant and to obtain quantitative estimates of the amount of amplitude reconstruction in each treatment group, Table 5 provides mean estimates and confidence intervals of the pairwise differences. To account for the uncertainty of the parameters' distributions, the estimates and their confidence intervals have been computed by bootstrapping of the mean with 1000 bootstrap samples being drawn in each simulation [17]. All confidence intervals in Table 5 are located below zero, indicating a significant increase in amplitudes from P_{2 }to P_{3 }in both treatment groups. In the E2group, an average amplitude rise of approximately 1.13°C occured. Amplitude rise due to Tibolone amounted to an average of 1.68°C.
Mean differences between vehicle and treatment phase  

Treatment group  $\left(\widehat{\mu}\right)$  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 ${Q}^{\left(k\right)}:=\frac{{\gamma}_{3}^{\left(k\right)}}{{\gamma}_{1}^{\left(k\right)}}$ 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}.
We compared the ratios obtained for each group by a onesided bootstrap test of the median for independent samples (1000 bootstrap samples, significance level α = 0.05). The corresponding hypothesis pair was
where θ^{(E2) }and θ ^{(TI) }are the location parameters of the ratio distributions in the E2 or Tibolone group, respectively. As shown in Table 6, the upper confidence interval boundary is slightly above zero, but far below the threshold of 0.2. This indicates that the Tibolone group exhibits a similar recovery rate in P_{3 }as the E2 group.
Differences between the amounts of amplitude reconstruction in the E2 and the Tibolone group  

Estimator  $\left(\widehat{\mu}\right)$  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.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
The work presented here is based on a collaboration between all authors working extensively on the experiments described in this paper. KP and VP designed and supervised the experiments. DG, BWI, KK and RV provided the mathematical and statistical framework and developed analytical tools. DG did the programming and the implementation. DG and BWI analyzed the data and DG wrote the main parts of the paper. All authors discussed the results and implications and commented on the manuscript at all stages. All authors read and approved the final manuscript.
Acknowledgements
We gratefully acknowledge the technical assistance of Tanja Lehmann, Andrea Trieller und Anita Griephan of Female Fertility Control & MM, Bayer HealthCare, Berlin.
References
 Andrikoula, M Prelevic, G Menopausal hot flushes revisited. Climacteric 12: 3–15, DOI: http://dx.doi.org/10.1080/13697130802556296 [PubMed]
 Rossmanith, WG Ruebberdt, W What causes hot flushes? The neuroendocrine origin of vasomotor symptoms in the menopause. Gynecol Endocrinol 25: 303–314, DOI: http://dx.doi.org/10.1080/09513590802632514 [PubMed]
 Berendsen, H Weekers, A Kloosterboer, H Tibolone and raloxifene on the tail temperature of oestrogendeficient rats. Eur J Pharmacol 419: 47–54, DOI: http://dx.doi.org/10.1016/S00142999(01)009669 [PubMed]
 Berendsen, H Kloosterboer, H Oestradiol and mirtazapine restore the disturbed tailtemperature of oestrogendeficient rats. Eur J Pharmacol 482: 329–333, DOI: http://dx.doi.org/10.1016/j.ejphar.2003.09.061 [PubMed]
 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 1028: 191–202, DOI: http://dx.doi.org/10.1016/j.brainres.2004.09.012 [PubMed]
 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 191: 399.DOI: http://dx.doi.org/10.1677/joe.1.06919 [PubMed]
 Opas, EE Gentile, MA Kimmel, DB Rodan, GA Schmidt, A Estrogenic control of thermoregulation in ERαKO and RβKO mice. Maturitas 53: 210–216, DOI: http://dx.doi.org/10.1016/j.maturitas.2005.04.006 [PubMed]
 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 64: 46–51, DOI: http://dx.doi.org/10.1016/j.maturitas.2009.07.007 [PubMed]
 Dacks, PA Rance, NE Effects of Estradiol on the Thermoneutral Zone and Core Temperature in Ovariectomized Rats. Endocrinology 151: 1187–1193, DOI: http://dx.doi.org/10.1210/en.20091112 [PubMed]
 Brown, E Czeisler, C The statistical analysis of circadian phase and amplitude in constantroutine coretemperature data. J Biol Rhythms 7: 177–202, DOI: http://dx.doi.org/10.1177/074873049200700301 [PubMed]
 Brown, E Choe, Y Luithardt, H C, C A statistical model of the human coretemperature circadian rhythm. Am J Physiol Endocrinol Metab 279: E669–E683. [PubMed]
 Indic, P Brown, E Characterizing the amplitude dynamics of the human core temperature circadian rhythm using a stochasticdynamic model. J Theor Biol 239: 499–506, DOI: http://dx.doi.org/10.1016/j.jtbi.2005.08.015 [PubMed]
 Brown, E Identification and estimation of differential equation models for circadian data. PhD thesis, Harvard University,
 Modelska, K Cummings, S Tibolone for postmenopausal women: systematic review of randomized trials. J Clin Endocr Metab 87: 16–23, DOI: http://dx.doi.org/10.1210/jc.87.1.16 [PubMed]
 Guckenheimer, J Holmes, P Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Copernicus.
 Team RDC. R A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, http://www.Rproject.org [ISBN 3900051070].
 Efron, B Tibshirani, R An introduction to the bootstrap. New York: Chapman & Hall.
How to cite: Girbig, D, Keller, K, Prelle, K, Patchev, V, Vonk, R and Igl, B 2012. A dynamic model of circadian rhythms in rodent tail skin temperature for comparison of drug effects. Journal of Circadian Rhythms 10:1, DOI: http://dx.doi.org/10.1186/17403391101 

This is an Open Access article distributed under the terms of the Creative Commons Attribution License


This article has been peer reviewed (journal peer review policy). 

Published on 5 January 2012. 
ISSN: 17403391  Published by Ubiquity Press  This work is licensed under a Creative Commons License.