A dynamic model of circadian rhythms in rodent tail skin temperature for comparison of drug effects

Abstract

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.

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 menopause-like 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 heat-dissipating 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.

Figure 1

Principles of the ovariectomized rat model. The effects of ovariectomy (OVX) on circadian rhythms in the tail skin temperature of the rat. In the intact animal, the temperature baseline is lowered during the active nocturnal phase. In the ovariectomized animal, the temperature in this episode is elevated, which results in a reduction of the distance to the baseline of the quiescent phase. The graphics have been adapted from [3].

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 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 [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-/night-rhythm 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 husbandry-related 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.

Figure 2

Telemetric device. Measurement device and transmitter used in the telemetry experiments (source: Bayer Pharma AG).

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.

As depicted in Figure 3, the experiment was split up into three consecutive phases:

Figure 3

Different phases and treatment groups. Experimental protocol of the OVX experiments.

1. In the estrogen phase (P₁), 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₂). 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₃), 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-/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.

We assume that the observed time series of some length N is given by the temperature values

(1)

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 ₀ + c ₁ · t with c ₀ > 0 and c ₁ ∈ ℝ 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

(2)

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 [15].

Since no exact solution of the van der Pol oscillator is available, we currently use the second-order perturbation solution which is given by

(3)

with the dynamic amplitude function

(4)

being a solution of the differential equation

[12]. Here, k is a constant of integration and ψ ₀ describes an angular phase shift. The initial amplitude a ₀ := a(0) depends on k and γ according to , showing also dependence of k on a ₀ and γ. The parameters of the circadian component can therefore be summarized by the parameter vector β ₀ = (a ₀, γ, ε, ψ, τ).

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.

The process 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 containing a total of 10 parameters (note that the 5 parameters of the van der Pol model are comprised in β ₀ ).

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:

(5)

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 off-diagonal entries are zero.

Equation (5) can be reformulated to the negative log-likelihood, which serves as the objective function for optimization:

(6)

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 ₀ and slope c ₁) 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 light-dark cycle to which the animals were subjected.

Optimization was then performed only for the remaining 4 parameters contained in β ₀ . 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 [16].

Initial values for start amplitude a ₀ were obtained by fitting a harmonic regression model y(t) = A ₁ cos(ωt) + A ₂ 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₁ - P₃ 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 ₀, small start amplitudes were promoted by starting optimization with initial value a ₀ = 0.5 in P₁. Upper boundaries for parameter a ₀ (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₃, 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 ₀). 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 short-time 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.

Figure 4

Example of model fit (E2 group). Examples of estimated model components belonging to an animal of the E2 group, separated by experimental phase (P₁ - P₃). (a) Tail skin temperature time measurements (y _n ); (b) circadian components ; (c) amplitude functions ; (d) trend curves ; (e) combination of estimated circadian components and trend curves , embedded in the time series. Daytime periods are shown in yellow; nighttime periods are shown in blue. The model parameter estimates are given in the Tables 1 and 2.

Figure 5

Example of model fit (Tibolone group). Examples of estimated model components belonging to an animal of the Tibolone group, separated by experimental phase (P₁ to P₃). (a) Tail skin temperature time measurements (y _n ); (b) circadian components ; (c) amplitude functions ; (d) trend curves ; (e) combination of estimated circadian components and trend curves , embedded in the time series. Daytime periods are shown in yellow; nighttime periods are shown in blue. The model parameter estimates are given in the Tables 1 and 2.

Table 1 Examples of estimated van der Pol parameters.

Table 2 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 P₁. (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₂ and rise again after substance treatment in P₃. 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₂ to P₃, indicating positive effects of each substance on the extent of circadian oscillation.

The reduction in the decrease of nighttime temperatures in P₂, as well as its subsequent restoration in P₃ 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.

Table 3 Estimated van der Pol parameters

Table 4 Estimated trend and variance parameters

Effects of the treatment in P₃ on the amplitude parameters

After successful treatment, we expect amplitudes to increase from P₂ to P₃. 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₂ tend to be located at lower levels than those at the end of P₁. In P₃, an overall rise in amplitude parameters can be observed relative to P₂.

Figure 6

Estimated limit cycle amplitudes per phase and treatment group. Limit cycle amplitudes occurring at the end of each experimental phase. (a) E2 group, (b) Tibolone group.

In order to illustrate amplitude reconstruction during treatment in P₃, Figure 7 provides boxplots of the pairwise differences in limit cycle amplitudes of vehicle and treatment phase (P₂ - P₃). As expected, all pairwise difference values are located below zero, indicating increase in amplitudes from P₂ to P₃ and thus successful amplitude restoration in both groups.

Figure 7

Changes in estimated amplitudes between vehicle and treatment phase. Pairwise differences (P₂ - P₃) of limit cycle amplitudes estimated in both treatment 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₂ to P₃ in both treatment groups. In the E2-group, an average amplitude rise of approximately 1.13°C occured. Amplitude rise due to Tibolone amounted to an average of 1.68°C.

Table 5 Mean differences between vehicle and treatment phase

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₃ and P₁ 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₁, that have been reconstructed during treatment in P₃.

We compared the ratios obtained for each group by a one-sided bootstrap test of the median for independent samples (1000 bootstrap samples, significance level α = 0.05). The corresponding hypothesis pair was

(7)

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₃ as the E2 group.

Table 6 Differences between the amounts of amplitude reconstruction in the E2 and the Tibolone group

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 t-test 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 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 [10]). 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 [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.