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.
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 [
An animal model suitable for assessing the effectiveness of drug candidates in hot flush
treatment is the ovariectomized rat model [
The ovariectomized rat model has been widely applied for the evaluation of substance
effects for treating menopausal syndromes [
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 [
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 [
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
The effects of
In contrast to the experiments described in literature [
As depicted in Figure
1. In the estrogen phase (P_{1}), animals were treated with
2. In order to ensure complete elimination of circulating
3. For the treatment phase (P_{3}), animals were randomized to a test group
(treatment by
Circadian rhythms in body temperature can be studied by dynamic models based on the van
der Pol equation [
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
with a sampling interval Δ
The circadian component
The van der Pol equation
is wellsuited for taking these properties into account [
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
[
The linear trend function
The error terms
The process
where
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
Here,
The covariance matrix
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
The period of oscillation was restricted to
Optimization was then performed only for the remaining 4 parameters contained in
Initial values for start amplitude
In contrast to previous studies where limit cycle amplitudes had been determined prior
to the start of the actual experiment for each individual [
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
Figures
Examples of estimated van der Pol parameters.
Examples of estimated van der Pol parameters  

















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 van der Pol model parameters estimated for animals of the
Examples of estimated trend and variance parameters
Examples of estimated trend and variance parameters  

















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 
Examples of trend and variance parameters estimated for animals of the
In both animals, oscillations vanish after surgery, with the limit cycle being
approached after day two (
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
Estimated van der Pol parameters
Estimated van der Pol parameters  











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  



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 
Mean values and standard deviations of the van der Pol model parameters estimated
for the
Estimated trend and variance parameters
Estimated estimated trend and variance parameters  











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  



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 
Mean values and standard deviations of the trend and variance parameters estimated
for the
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
Boxplots of the limit cycle amplitudes obtained from the time series of each treatment
group (10 animals per group) are shown in Figure
In order to illustrate amplitude reconstruction during treatment in P_{3},
Figure
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
Mean differences between vehicle and treatment phase
Mean differences between vehicle and treatment phase  







1.1265  ( ∞,  0.8247) 

1.6807  ( ∞,  1.3224) 
Mean estimates
In order to compare the performances of both treatments, we computed the ratios
We compared the ratios obtained for each group by a onesided bootstrap test of the
median for independent samples (1000 bootstrap samples, significance level
where
Differences between the amounts of amplitude reconstruction in the E2 and the Tibolone group
Differences between the amounts of amplitude
reconstruction in the 








0.2294  ( ∞, 0.0078)  0.0031 
Median estimates
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 [
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 [
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 [
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 authors declare that they have no competing interests.
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.
We gratefully acknowledge the technical assistance of Tanja Lehmann, Andrea Trieller und Anita Griephan of Female Fertility Control & MM, Bayer HealthCare, Berlin.