## Abstract

The human sleep process shows dynamic alterations during the night. Methods are needed to examine whether and to what extent such alterations are affected by internal, possibly time-dependent, factors, such as endocrine activity. In an observational study, we examined simultaneously sleep EEG and nocturnal levels of renin, growth hormone (GH), and cortisol (between 2300 and 0700) in 47 healthy volunteers comprising 24 women (41.67 ± 2.93 yr of age) and 23 men (37.26 ± 2.85 yr of age). Hormone concentrations were measured every 20 min. Conventional sleep stage scoring at 30-s intervals was applied. Semiparametric multinomial logit models are used to study and quantify possible time-dependent hormone effects on sleep stage transition courses. Results show that increased cortisol levels decrease the probability of transition from rapid-eye-movement (REM) sleep to wakefulness (WAKE) and increase the probability of transition from REM to non-REM (NREM) sleep, irrespective of the time in the night. Via the model selection criterion Akaike's information criterion, it was found that all considered hormone effects on transition probabilities with the initial state WAKE change with time. Similarly, transition from slow-wave sleep (SWS) to light sleep (LS) is affected by a “hormone-time” interaction for cortisol and renin, but not GH. For example, there is a considerable increase in the probability of SWS-LS transition toward the end of the night, when cortisol concentrations are very high. In summary, alterations in human sleep possess dynamic forms and are partially influenced by the endocrine activity of certain hormones. Statistical methods, such as semiparametric multinomial and time-dependent logit regression, can offer ambitious ways to investigate and estimate the association intensities between the nonstationary sleep changes and the time-dependent endocrine activities.

- multinomial logit model
- transition probabilities
- cortisol
- renin
- growth hormone

eeg data registered and obtained from humans during sleep form the basis of most studies in sleep research. From a statistical point of view, these data, produced by electrophysical measurements (potentials) at various scalp positions, represent continuous time series with continuous state space. Because of the high registration density (generally >240 measurements per second), the volume of sleep EEG data is immense, and, consequently, many difficulties are encountered in statistical evaluation and interpretation of the data. Hence, sleep scientists usually do not work with the raw sleep EEG data but, rather, with their transformations, leading to a state space and to a time scale condensation as well. The most common method for EEG data condensation is based on the criteria of Rechtschaffen and Kales (19). According to these criteria, the original time series of the sleep EEG signals are transformed into new time series with a discrete state space containing about six to eight stages and with a time scale calibration based on units of 30-s intervals (epochs).

Each stage is associated with distinct electrophysical patterns visible in overnight EEG recordings. For example, slow-wave sleep (SWS) is characterized by the presence of slow waves (δ-waves) in the range 1–4 Hz. Thereby, after dividing the EEG time series into epochs of, for instance, 30-s duration, the rater compares the observed 30-s EEG trajectory with the characteristic patterns associated with the different sleep types and assigns the value of the best-matching sleep stage to the epoch. This classification forms the basis of most statistical analyses to explain human sleep.

Two main types can be distinguished: rapid-eye-movement (REM) and non-rapid-eye movement (NREM) sleep. NREM sleep is further divided into four stages: S1 and S2, which correspond to light sleep (LS), and S3 and S4, which correspond to SWS.

To analyze the phenomenon of human sleep as a whole, not only static characteristics of sleep architecture summarizing the overall sleep behavior (such as total duration spent in each sleep stage), but also approaches that are able to take into account the dynamics of the sleep process, should be considered. This also allows exploration of time-dependent factors and their association with the transitions between the sleep stages. In this study, we investigate how alterations in endocrine activity relate to the probability of specific sleep stage transitions. Sleep is a time of considerable activity in various endocrine systems. For example, the major portion of REM sleep occurs during the second half of the night, when plasma cortisol levels are high (23, 25). Furthermore, many studies support the view that the neuropeptides growth hormone (GH)-releasing hormone (GHRH) and corticotropin-releasing hormone (CRH) act as common regulators of sleep architecture and sleep-associated hormone release (7, 11, 17, 22). Renin shows oscillations that are distinctly associated with the REM-NREM sleep cycle. Renin levels increase during NREM periods and decrease during REM periods (4).

A bidirectional interaction exists between sleep EEG and nocturnal hormone secretion (21). Administration of hormones affects sleep EEG. On the other hand, manipulations of sleep (e.g., sleep deprivation) modulate hormone secretion. For example, cortisol levels rise after awakenings (25). The exact causal relationships have not been revealed. However, in this study, we restrict the analysis to the investigation of hormone influence on sleep EEG and propose the use of semiparametric multinomial logit regression models developed by Kneib and Fahrmeir (15). The multinomial logit model explains the probability of selected sleep stage transitions by several influential factors, such as time and hormone secretion, which also varies with time. Extension of purely parametric, linear models developed in the context of generalized linear models (10) to semiparametric versions allows for considerably more flexibility in analyzing temporal developments and hormone effects.

For similar data structures, the association between intensities of sleep stage changes and plasma cortisol concentration has been investigated (27). Event history analysis (EHA), a statistical approach to gain insight into the dynamics of human sleep based on time-continuous data, was applied. Since the 30-s epochs of EEG sleep recordings are very small compared with the overall measurement time in a night (8 h), a time-continuous model seems plausible for analyzing sleep dynamics. As an extension of Cox regression, the aim of EHA is to model and estimate transition intensities as functions of time, taking into account individual sleep history and assessing the influence of factors of interest, such as hormone secretion. A variable containing hormone measurements can be included as a time-varying covariate, provided that it is observed continuously over time. However, the hormone concentration is measured, for technical and ethical reasons, only every 20 min: therefore, interpolation is required to obtain the missing hormone measurements. Although time-continuous models may be the appropriate modeling choice when the dynamics of sleep are of particular interest, the estimated impact of hormone concentration may well be influenced by the choice of a specific interpolation strategy. In particular, the causal interpretation of estimation results can be severely affected by considering interpolation of the hormone values. Since the choice of a particular interpolation strategy will always remain arbitrary, at least to a certain extent, it is difficult to separate true effects from those induced by interpolation. The comparison of results obtained with different interpolation strategies might act as a sort of sensitivity study, but uncertainty about the conclusions and increases in the computational burden remain.

In contrast to EHA, the multinomial logit model works in discrete time and, therefore, offers a natural grouping of the data, which makes interpolation unnecessary. Frequencies of different types of sleep stage transitions can be calculated for each 20-min interval linked to a corresponding hormone measurement. We therefore consider discrete time semiparametric multinomial logit models to be a valuable alternative to EHA when the impact of hormone measurements is the focus of a study. Although we might lose details in the time domain, we obtain a more realistic model for the hormone-sleep relationship, while the reverse implication holds for EHA.

## METHODS

#### Study population.

We analyzed the baseline night of a study reported elsewhere (20). The participants were 47 healthy volunteers (19–69 yr of age): 24 were women [41.67 ± 2.93 (range 19–64) yr] and 23 were men [37.26 ± 2.85 (range 22–69) yr]. A total of 55 subjects entered the study. For technical reasons, eight dropped out and were not included in the analysis. The study was approved by the Ethics Committee for Human Experiments of the University of Munich. After the purpose of the study had been explained to the subjects, they gave their informed consent according to the tenets of the Declaration of Helsinki. All subjects underwent an extensive physical and psychiatric examination; only those with normal results were included in the study. Individuals with sleep disturbances or a personal or family history of psychiatric disorders, as well as those who performed shift work and a time shift in the last 3 mo before the beginning of the study, were excluded. The subjects were free of any medication for ≥8 wk. The subjects abstained from alcohol throughout the study period starting 1 wk before the first (adaptation) night. Caffeine was restricted to 200 ml of coffee in the morning. All premenopausal women (*n* = 17) had regular menstrual cycles and were not taking hormone contraceptives. They were examined 2–6 days after menstruation. Nine of the women were postmenopausal.

#### Experimental procedure.

The study consisted of two consecutive nights: on the first (adaptation) night, the subjects adapted to the laboratory; on the following (examination) night, an indwelling intravenous catheter was inserted at 1930 and connected to plastic tubing that ran through a sound-proof lock into the adjacent room. This arrangement allowed us to collect blood without disturbing the study subject. All subjects were in the supine position from 2000. According to our experience from various similar studies (21) and the literature (12), cannulation does not affect sleep EEG distinctly. Blood samples were collected between 2300 (lights out) and 0700 at 20-min intervals for measurement of cortisol and GH and at 10-min intervals for measurement of renin. The subjects were observed on a television screen in an adjacent room. Without loss of generality, we use a common time scale with 20-min calibration for all hormone measurements between 2300 and 0700.

#### Endocrine recordings.

Free renin, cortisol, and GH concentrations were measured by radioimmunoassay (Advantage, Nichols Institute, San Juan Capistrano, CA; intra- and interassay coefficients of variation <10%). Random samples for each hormone were analyzed in duplicate. According to standard procedures for time series, the remaining specimens were analyzed only once.

#### Sleep EEG recordings.

Electrodes for polysomnographic recordings (Comlab 32 Digital Sleep Lab, Brainlab version 3.3 software, Schwarzer, Munich, Germany) were fixed during the examination night between 2100 and 2200. The subjects were not allowed to sleep until the lights were turned off at 2300. Polysomnographic recordings were obtained from 2300 to 0700 and consisted of two EEGs (C3-A2 and C4-A1; time constant 0.3 s, low-pass filtering at 70 Hz), vertical and horizontal electrooculograms, and electromyogram. The EEG signals were registered at a sampling rate of 256 Hz, then filtered with a digital low-pass filter, and finally stored with a resolution of 128 Hz on a personal computer. The sleep stages [wakefulness (WAKE), stages 1–4 (NREM) sleep, and REM sleep] were visually scored for 30-s epochs according to the criteria of Rechtschaffen and Kales (19). After the scoring, the hypnograms, i.e., the time series of the consecutive sleep stages, serve as samples for our analyses.

#### Statistical analysis.

The data generated in our experimental setup can be characterized as follows. In each 30-s epoch *t*, the current sleep state *Y*_{t} is recorded, taking one of the four values “WAKE,” “REM,” “SWS,” or “LS.” The quantities of interest in event history analysis are the transition intensities (1) where the epoch index *t* is interpreted as a continuous time measurement and the indexes *j* and *k* relate to the sleep state immediately before and after the transition, respectively. Transition intensities are basically a scaled version of the intrinsic transition probabilities *P*(*Y*_{t+Δt} = *k*|*Y*_{t} = *j*), which represent the tendency for moving to state *k* within a very short time interval given that one is in state *j* at time *t*. In a regression model, the transition intensities are related to covariate effects as where η_{jk}(*t*) is an additive predictor comprising time-varying effects, parametric covariate effects, and, possibly, a number of extensions that have been applied in the context of sleep research (16, 27).

Although a time-continuous model is well justified when exploring the time-dependent dynamics of human sleep, it may be questionable when investigating the impact of endocrine activity, which is measured on a much coarser time axis. Therefore, we now consider a discrete time approximation to the transition intensity (*Eq. 1*) based on multinomial logit models (10, 15). For small time increments Δ*t*, *Eq. 1* can be approximated by (2) leading to a regression model for the transition probabilities in terms of the same predictor in the EHA model. In particular, taking Δ*t* = 1 yields a one-step-ahead discrete time approximation.

The regression model corresponding to *Eq. 2* is exactly the multinomial logit model, where the probability π_{jk}(*t*) for the transition from *j* to *k* at time *t* is given by and η_{jk}(*t*) are suitably chosen predictors formed by covariates and *l* = 1,…, *L* denotes all sleep states that are reachable from sleep state *j*, including sleep state *j* itself for the case of no transition to another stage. In the discrete-time multinomial logit model, it is necessary to impose some identifiability constraints on the predictors, since the probabilities π_{jl}(*t*), *l* = 1,…, *L*, have to sum to 1, inducing redundancy in the model formulation. We therefore chose the last category as the reference category and set η_{jl}(*t*) ≡ 0, leading to the final multinomial logit model formulation (3) with (4) Particularly, we chose the reference category to be *L* ≡ *j* for no transition.

In fact, the multinomial logit model defines separate models for each current sleep state *j*. For our analysis, we therefore split the data according to the current state value, obtaining data subsets of observations. For each of these subsets, a number of possible transitions is defined, leading to multinomial logit models with categorical responses of different dimensionality. Table 1 provides an overview of transitions of interest per initial state.

The remaining challenge in setting up a model for investigating the impact of endocrine activity on sleep transitions is to define suitable predictors for *Eq. 3*. Clearly, purely parametric versions of the multinomial logit model requiring all effects to be linear are inappropriate, since the dynamics of human sleep are too variable to be captured in, e.g., a linear time trend. Similarly, if we are interested in temporal changes of hormone effects, i.e., interactions between time trends and hormone secretion effects, flexible surface estimation, instead of parametric interactions, is required. We therefore utilize recently developed semiparametric extensions of multinomial logit models (14, 15) that rely on penalized spline smoothing (see appendix). To be more specific, we consider the following sequence of models for the hormones renin, cortisol, and GH. In *model 1* (*Eq. 5*), for each hormone variable, one model is fit separately (5) In *model 2* (*Eq. 6*), for each hormone variable separately, the interaction with time is assessed by estimating an interaction surface. By use of an interaction term, the possibility that the effect of hormone secretion (its functional form and magnitude) changes with time can be taken into account. Main effects in time and hormone secretion are included in addition, so that the interaction surface can be interpreted as the deviation from the main effects (6) In *model 3* (*Eq. 7*), all three hormones, renin, cortisol, and GH, are analyzed in one model (7)

All models contain a nonparametric age effect and sex as a factor variable as additional covariates. The time dependence of sleep transition intensities is assessed on the basis of the smooth functions *f*_{0,jk}(*t*) and the additional interaction *f*_{3,jk}(*t*, hormone) in *model 2*. Although *models 1* and *3* assume a clear separation between time variation and hormone effects that are time constant, *model 2* allows for interactions between the two. However, *model 2* also induces considerably more complexity and, therefore, requires that more information be estimated. This is also why we did not consider interaction models analyzing the effects of all hormones simultaneously (for more technical details on fitting semiparametric regression models, see appendix).

To account for the fact that the hormone measurements are collected on a time scale that is different from the sleep states, we impose a special structure on the time-varying effects in our models. To be more specific, we assume that the time-dependent functions are defined on the coarser time level provided by the hormone measurements; i.e., η_{jk}(*t*) remains constant over all epochs collected in a 20-min time span of a specific hormone measurement.

The goodness of fit of each model is assessed with Akaike's information criterion (AIC). It is used to determine for each model whether the interaction term between time and hormone variable leads to a reasonable improvement in the (corrected) model fit, so that the inclusion of the term is necessary. In contrast to the log-likelihood criterion, AIC takes into account the increased model complexity that arises when additional terms are included in the predictor. More specifically, AIC combines the log-likelihood criterion as a measure for the goodness of fit with an effective parameter count that measures the model complexity. As a consequence, the combined measure trades off the two mutually contradictory goals to achieve a simple interpretable model that generalizes well, on the one hand, and a model that describes all features of the true model, on the other hand. (For mathematical definition of the AIC, see appendix.)

## RESULTS

#### Model selection.

Table 2 lists the AIC values of models with and without hormone-time interaction. In principle, when two regression models are compared by AIC, the model with the smaller AIC value should be selected. Since, however, some of the AIC values in Table 2 are actually quite close to each other, we have, in some cases, deviated from this rule. In these cases, we combined evidence for the inclusion of an interaction provided by AIC with the assessment of uncertainty reflected in the confidence bands and the interpretability of the estimated surface. The rationale is that we are not in the situation of testing the presence of an interaction but, rather, want to set up a sensible model that describes our data.

For the models with initial state NREM, REM, or SLEEP (i.e., REM or NREM), there is no reduction or only a minor reduction in AIC when the interaction effect is included. Hence, it seems reasonable to assume that the effect of the hormone level does not change with time (and vice versa) in these models. When NREM sleep is divided into LS and SWS, an interaction is detected. It also seems reasonable to include the interaction in the models with the initial state WAKE.

The influence of time and hormone value on the specified transition probabilities (Table 1) is visualized with surface plots in Figs. 1–5. If an interaction is included in the model, the entire interaction surface plus main effects is plotted. Otherwise (models without interaction), a surface obtained from combining the two main effects is plotted. *Equation 3* is used to transform the estimated effects from the scale of the linear predictor to transition probabilities. In particular, the plotted probability surfaces represent the estimated probability of a specific transition within the next 20-min sampling interval (among all observations of the given initial state). For both kinds of plots, the age effect is fixed at a median level of 36 and the sex effect is set to 0 for the male reference. For the female study population and other age levels, the plot is (roughly) proportional to those presented here.

In Fig. 1, the influence of time and cortisol level on the probability of LS-SWS transition is visualized. Comparing the surfaces with and without interaction, one realizes that information is gained when the interaction is taken into account. On the one hand, areas can be recognized where no observations of the specified transition occurred, so that the interaction effect cannot be estimated (high cortisol levels early in the night, low cortisol values late in the night). On the other hand, a refined surface structure can be observed. For low cortisol levels earlier in the night, peaks for the transition probability, which do not exactly coincide with peaks of the curves from the model without interaction, can be found. Hence, structures that cannot be detected if the analysis is restricted to additive main effects emerge. This leads to the conclusion that, in the given example, the interaction surface is a better fit to the data.

#### Models without interaction.

As noted above, for initial state NREM, REM, or SLEEP, AIC implies that it suffices to consider models without interaction (*Eq. 5*). Thus, in these cases, fitting interaction surfaces does not improve the model fit, and the data structure can be explained by the main effects alone. Particularly, this means that there is evidence that the hormone effect does not change with time, and vice versa.

In Figure 2, the results of the model for initial state REM and hormone variable cortisol are visualized. In this model, the set of possible transitions comprises transitions to WAKE, NREM, and no transition. A separate effect is estimated for each transition (category) within the same regression model. The category “no transition” serves as a reference category, and the corresponding probability can be obtained from the estimated effects based on *Eq. 4*. Point-wise subtraction of surfaces in Fig. 2 yields a surface for the probability of no transition shown in Fig. 3. This example is shown to demonstrate how to calculate the probability of the reference category no transition. Subject matter results are, in this case, not easily derived. We observe that the probability of no transition is always >90%. It is leveled at ∼95%, implying that it is very probable to remain in REM sleep under the condition to be currently in REM sleep. This observation is unfortunately no substantial gain in information but, rather, depends on the relation between the typical duration of an REM period and the sampling interval.

Although the interpretation of the probability plots is straightforward, one can use the plots in Fig. 2 (*middle* and *right*) to assess the nonparametric effect of time and cortisol separately. Note that examining the effects separately is only meaningful when no interaction is present. As opposed to the probability surface plot, 95% confidence intervals can easily be attached to the curves to visualize the uncertainty associated with the estimation. Standard error estimates for parameters on the logit scale are a direct result of the optimization procedure used for estimation, whereas standard error estimates on the probability scale require the application of the multivariate Δ rule (see appendix). This is due to the fact that the multinomial probability transform (*Eq. 3*) does not only involve a single predictor but, rather, all predictors of the model.

For the time and cortisol effects, we obtain very smooth, monotonic functions. With increasing cortisol level, the effect of cortisol decreases for the REM-WAKE transition and increases for the REM-NREM transition. The approximate linear form of the effect curves implies that, for explaining REM transitions, it suffices to include cortisol as a linear effect. However, large confidence intervals correspond to large estimation uncertainty, rendering the cortisol effect rather negligible in both models. As one might expect, the time effect for the REM-WAKE transition increases during the night, whereas the time effect for the REM-NREM transition decreases. Hence, with increasing sleep duration, it is more probable to awaken from REM sleep than to pass into NREM sleep.

The increased tendency for the REM-WAKE transition at the very end of the night, however, might be influenced by the lights-on time of 0700. Hence, awakenings from REM sleep induced by this could have led to an artificial increase in the probability for the REM-WAKE transition. Nevertheless, post hoc descriptive statistics show that only a minor proportion of subjects were awakened from REM sleep and that the trend to an increase in REM-WAKE transition toward the end of the night can still be observed when data are restricted to observations made before lights on.

In two further models, the influence of either GH and time or renin and time on transitions from REM sleep is investigated. The implications of the results (not shown) are similar to those of the cortisol model. Changes in the probability for REM-WAKE and REM-NREM transitions are mostly due to time shifts. For increasing time, we observe an increased tendency for the REM-WAKE transition and a decreased tendency for the REM-NREM transition. Taking the uncertainty expressed in point-wise 95% confidence intervals into account, we were not able to detect reasonable effects of hormone variables. Only for renin concentrations <3.4 nmol/l, a slightly negative effect on the REM-WAKE transition can be detected, implying that low renin concentrations inhibit awakening from REM sleep.

In Fig. 4, the influence of time and cortisol concentration on the probability of NREM-REM transition is visualized. In this case, cortisol concentrations ∼45 nmol/l enhance NREM-REM transitions. Also, the course of the time effect exhibits peaks in a periodic manner. The first peak appears at ∼2 h after recording onset, which corresponds to the REM latency (Table 3). Subsequent peaks show a periodicity that is similar to NREM-REM cycles. Although the peaks become less pronounced in the middle of the night, there is a distinct peak at the end of the night, making the NREM-REM transitions very probable. Contrary to possible assumptions, our findings do not support the view that the effect of cortisol on the NREM-REM transition changes with time and, hence, is associated with an enhanced pressure for the NREM-REM transition connected with NREM-REM cycles. This follows from the fact that the interaction term does not improve the model fit.

Analogous conclusions also hold for the models with either renin and time or GH and time (not shown). The nonparametric time effect in either of these models exhibits the same functional form described for cortisol. The interaction term does not add substantial information to explain transitions from NREM during sleep according to AIC. Again, this implies that we have not found evidence that either GH or renin is associated with an enhanced pressure for the NREM-REM transition connected with NREM-REM cycles.

#### Models with interaction.

For initial state LS, SWS, or WAKE, the AIC comparison implies that the data structure cannot be explained by main effects alone, since fitting interaction surfaces does improve the model fit. Specifically, this means that the hormone effect does change with time, and vice versa. Hence, the model formulation of *Eq. 6* is used.

The influence of time and hormone level on the SWS-LS transition is visualized via three-dimensional and image plots in Fig. 5 (*left* and *middle*), along with confidence interval information (Fig. 5, *right*). Significant deviations from the null level are indicated by color gradients toward black (decrease) and white (increase). The null level is defined as the relative frequency of the SWS-LS transition relative to all transitions from SWS observed in the whole night and is estimated to be 17.4%. Additionally, we plotted the location of available observations on the time-hormone level grid (Fig. 6). We observed that, particularly for GH and renin, the scatterplots reveal skewed distributions: the conditional distribution of hormone values given time appears highly skewed to the right. To account for possible fitting instabilities introduced by this data structure, we refitted the models with log-transformed hormone variables. Since the results compare very well with the untransformed model, we chose to present the original results. As noted above, the interaction surface cannot be estimated for regions where no observations are available.

For the cortisol and renin data, relatively rough surfaces are observed. Refitting the models with fewer smoothing spline knots enforcing greater smoothness reveals comparably rough patterns. Hence, we conclude that the roughness is inherent to the data.

For the cortisol model, as one might expect, the scatterplot exhibits a dependency between time and cortisol level: the cortisol concentrations increase with time. Very high concentrations of cortisol are only observed later in the night. Furthermore, the transition probability plots suggest that the probability for the SWS-LS transition is connected with this phenomenon. It is lowest early in the night, when cortisol levels are low. If cortisol concentrations and time values increase, peaks are found, denoting an increase in transition probability. At cortisol concentration ∼90 nmol/l and 3 h after recording onset (onset at 2300), a prominent peak can be observed. Although only a few observations underlay this structure, this peak is significantly different from the null level. Generally, it can be said that there is a considerable increase in probability for SWS-LS transition toward the end of the night, when higher cortisol concentrations can be observed. The reverse is true for the LS-SWS transition (Fig. 1): high probabilities for LS-SWS transition can be observed shortly after sleep onset, when cortisol levels are low. In the middle and toward the end of the night, the transition probability strongly decreases. Hence, the influence of cortisol and time on LS-SWS and SWS-LS transitions shows a complementary trend. Later in the night, fewer LS-SWS transitions occur, and if SWS occurs, the pressure to leave SWS is increased. These findings are compatible with the fact that the amount of SWS is strongly reduced toward the end of the night.

Remarkably, in the SWS-LS transition model, besides the peaks of increased transition probability, dips of decreased transition probability also exist in the immediate vicinity. An intuitive explanation approach is not available for such a complicated structure.

Although the model fit for the hormone variable GH actually does not improve on inclusion of the interaction term (Table 2), for comparison, it is also considered in this discussion. Early in the night, within the first 200 min of analyzed time, some measurements with relatively high concentrations of GH can be observed. The transition probability for this time span remains relatively low, regardless of the GH concentration. For this time span, also larger regions can be found with estimated transition probabilities below the lower limit of corresponding point-wise confidence intervals, implying a significant decrease in the SWS-LS transition for GH concentrations ∼7 and 19 nmol/l. Later in the night, for low GH levels, the transition probability first peaks at ∼4:40 h after recording onset and again at about the end of the night, implying a significant increase in the number of SWS-LS transitions.

For the renin model, the transition probability is highest on the following axis: for renin concentrations ∼15 nmol/l at ∼4 h after recording onset (onset at 2300), with a shift to lower concentrations at ∼7 nmol/l later in the night. These peaks are significantly different from the null level of 17.4%. Moreover, a significant increase in the SWS-LS transition can be found for very high renin concentrations (≥25 nmol/l). Again, the transition probability is lowest early in the night, when the renin level is low.

#### Models with all hormone variables.

When all three hormones, renin, cortisol, and GH, are analyzed in one model, the model formulation of *Eq. 7* is used. It was found that the functional form of the hormone effects is comparable to those of the models for each hormone alone without interaction; therefore, the results are not presented in detail. The fact that there is no substantial difference between the results of the model including all hormone variables simultaneously and models with one hormone variable alone is, nevertheless, interesting. It seems to imply that there is hardly any confounding between the hormones, at least when their impact on the transition probabilities is modeled. Particularly, there is a lack of (linear) dependency between the independent hormone variables, which is supported by only moderate correlations between the hormones in our data. By inspection of both a correlation table for all observations during the night and separate correlation tables per initial state, all correlations with renin are only about ±0.1, whereas correlations between GH and cortisol are about −0.3. The correlations have been computer averaged over all measurements in the data, i.e., using all measurements of all individuals, which may impact the results.

## DISCUSSION

The human sleep process shows many dynamic alterations during the night. To analyze these dynamic structures, it is insufficient to investigate static sleep characteristics summarizing the sleep behavior over larger portions of the night. Methods are needed to examine, on the one hand, how the sleep structure changes (continuously) over time and, on the other hand, whether and to what extent such alterations are affected by internal, possibly time-dependent, factors, such as endocrine activity.

A problem arises because of the different time scales used for measurements of hormones and sleep stages. For technical and ethical reasons, the time between successive hormone measurements is by far greater (20 min) than the time unit (30 s) in which sleep stages are recorded. Depending on the research question at hand, a decision must be made between a statistical model preserving the quasi-continuous time structure of the sleep data and a more realistic model for the hormone-sleep relationship, which works in discrete time. The former model (EHA) was applied by Yassouridis et al. (27) in an investigation of the change in transition intensity of selected sleep stages over time. Hormone measurements were included by fitting a separate time course for high cortisol values (because only high concentrations were expected to influence intensities). By this method, cortisol entered the analysis only as a binary variable, circumventing the interpolation of exact hormone values for each 30-s interval.

In this study, we applied a multinomial logit model that works in discrete time. With this method, a natural grouping of the data recorded within each 20-min interval made interpolation unnecessary. The hormone measurements now entered the analysis as metric variables, enabling us to assess the effect of different cortisol levels. In addition, semiparametric modeling of the time and hormone effect lead to considerably more flexibility in analyzing temporal developments and hormone effects; i.e., with increasing time (hormone level), the linear predictor is not assumed to increase/decrease linearly but, rather, can take a flexible form. By including a hormone-time interaction, we allowed the hormone effect (curve) to vary over time. For models with initial state NREM, REM, or SLEEP (i.e., REM or NREM), the effect does not change with time, rendering an interaction term unnecessary. When NREM is divided into LS and SWS, an interaction is detected. It also seemed reasonable to include the interaction in the models with initial state WAKE.

Results were presented with surface plots showing the effect of endocrine activity and time on selected transition probabilities. As many plots indicate, complicated structures were captured and visualized.

Most of the probability courses of the transitions investigated in the present study show nonstationary shapes. In particular, models that show a significant hormone-time interaction indicate a dynamic structure of the sleep stage changes over the night. However, the aim of the actual study was not verification of nonstationarity of the sleep stage changes but, rather, investigation of the interplay between sleep activity and neuroendocrine activity assigned to the hormones cortisol, renin, and GH.

Some results of the present study are compatible with, and some differ distinctly from, the results obtained from the study by Yassouridis et al. (27) in which a sample of 30 healthy male volunteers (20–30 yr of age) was enrolled and the endocrine activity was focused on cortisol only. The intensity of the SWS-LS transition in the previous study showed an extreme peak at 3 and 7 h after sleep onset. This trend was intensified by high cortisol concentrations, which steadily increased the relative intensity for SWS-LS transition from 2 h after sleep onset to the end of the night. The cortisol effect on the SWS-LS transitions is also present in our study population, which consists of men and women with a wide age range, but we did not observe peaks at 3 and 7 h after sleep onset. In addition to the cortisol effect, there is also a renin effect (Fig. 5) caused by medium and low renin concentrations at ∼4 h (5 and 7 h, respectively) after sleep onset.

Whereas sex (1) and aging (24) affect sleep-endocrine activity, the present study did not consider interaction effects between these variables and hormone measurements (or, alternatively, time) for increasing the computational demand tremendously. The results are adjusted for sex and aging, which have an effect on transition probabilities also in our models. However, interpretation of curves and surfaces in our model does not change for differences in age and sex.

Distinct effects of cortisol on sleep changes were found. Cortisol is the peripheral end point of the hypothalamic-pituitary-adrenocortical (HPA) system. Its central key hormone is CRH. It is thought that a reciprocal interaction of GHRH and CRH plays a key role in sleep-endocrine regulation (6, 21). Changes in the GHRH-to-CRH ratio in favor of CRH appear to result in more shallow sleep. It appears that, during the first half of the night, the influence of GHRH predominates, which results in major portions of SWS and GH during this interval. During the second half of the night, CRH becomes more active, leading to more shallow sleep and low GH and high cortisol levels. Furthermore, the major portion of REM sleep occurs in the second half of the night. Ehlers and Kupfer (6) suggested that GHRH and CRH participate in the two-process model of sleep regulation proposed by Borbély (3). According to this hypothesis, GHRH would represent the process S, a sleep-dependent process related to SWS with an exponential decay during sleep. CRH is thought to represent the process C, the circadian oscillator, which is reflected by a rhythmic variation of sleep propensity in which REM sleep is oscillated to approximate the level of process C during sleep.

Our finding that cortisol concentrations ∼45 nmol/l enhance NREM-REM transitions supports the view that the HPA system promotes REM sleep. During the second half of the night, the longest duration of REM sleep and the highest level of cortisol occur simultaneously. The duration of REM sleep is longer in CRH-overexpressing than wild-type mice (13). After treatment of multiple sclerosis patients with a glucocorticoid, REM sleep was disinhibited (2). Alternatively, it appears possible that an increased probability for NREM-REM transitions is accompanied by increases in cortisol secretion.

Also the decrease in probability for SWS-LS transition early in the night, when GH levels are intermediate or high and cortisol levels are low, fits the hypothesis of a major contribution of GHRH-CRH interaction to sleep regulation. These findings are consistent with results of some other experimental studies (18). Low cortisol and intermediate-to-high GH concentrations early during the night probably reflect a low activity of CRH and a high activity of GHRH related to a long duration of SWS and a low probability for a change in LS. The peak of the transition probability was found at cortisol concentration ∼85 nmol/l and by 0300. At this time, the cortisol morning peak frequently occurs and, probably, a transition toward the preponderance of CRH takes place. The considerable increase in the probability for SWS-LS transition toward the end of the night corresponds with the maximum HPA activity during the early morning hours. On the other hand, the probability of LS-SWS transition is most likely after sleep onset, when HPA activity is low, and decreases in the middle and toward the end of the night. Similarly, high GH concentrations within the first 300 min of the sleep period coincide with a low probability for SWS-LS transition, which declines later in the night.

Renin is the hormone that shows the most distinct ties to the NREM-REM cycle (4). Renin levels show oscillations throughout the night, with maximum values during NREM sleep and decreases during REM sleep. Therefore, it may appear surprising that we did not find an association between renin and the NREM-REM transition. Probably renin is linked to the NREM-REM cycle by an unknown mechanism, whereas it has no effect on the REM-NREM transition. We observed only a slight negative effect for low renin concentrations on the REM-WAKE transition. Furthermore, we found that the probability of SWS-LS transition is lowest early in the night, when the renin levels are low, whereas it increases at ∼4 and 7 h after sleep onset, when renin concentrations increase to ∼14 nmol/l, and then shift to lower concentrations again.

In summary, we can assert that the idea of grouping sleep stage observations around the measurements of hormone concentrations, calculating into them the frequencies of interesting transitions, and using semiparametric extensions of multinomial logit models with flexible surface estimation for evaluating the effects of some presumably influential factors on these transitions provides more flexibility and more promised and more realistic results than other analyses and models. It should be kept in mind that, in addition to the hormones reported here, various neurochemical systems (e.g., monoamines, cytokines, adenosine) participate in sleep regulation.

### Perspectives and Significance

Despite some complexity, the semiparametric multinomial logit models that have been recommended here for studying and quantifying time-dependent hormone effects on the transition probability between certain sleep stages present a powerful way to gain insights into brain function during sleep. Taking into account the innumerable external stimuli to which every conscious human is exposed and noticing that sleep is free of such stimulus-reaction processes, we consider the examination of sleep to be a basic tool for gathering knowledge about the intrinsic functionality of the brain. The EEG is an indispensable device for analyzing sleep, in that it is the major discriminating marker between different sleep stages, e.g., WAKE, REM, and NREM (26). Transitions between these sleep stages are one aspect of a nonstationary sleep behavior. We now hypothesize that these dynamics are caused by themselves or by some internal influential factors. Since almost all observed phenomena in the world have some etiology in the background, we tend to believe that sleep dynamics are caused mostly by certain internal factors, e.g., hormones, and less by themselves. From a mathematical and statistical point of view, the search for associations between dynamic processes such as sleep and the hormone concentrations is naturally a challenge. Any new model, model extension, or discovery in this complex field will provide additional insight into the complicated structure, physiology, mechanism, and function of the most relevant human organ, namely, the brain. We like neither to overvalue our models presented in this study nor overestimate the results obtained with them. However, we see in our models some courageous steps toward a better understanding of human sleep and believe that any extensions of them and further investigations of the objectives in common with them will contribute significantly to a better explanation of brain physiology and functions.

## APPENDIX

#### Implementation.

Semiparametric multinomial logit models as described in *Statistical analysis* are available in the software package BayesX (http://www.stat.uni-muenchen.de/∼bayesx). Although the BayesX documentation employs a Bayesian perspective on semiparametric regression, BayesX is a software tool not only for Bayesian inference but, in general, for semiparametric regression models. In particular, mixed-model-based inference for penalized splines makes perfect sense, and we focused on a frequentist interpretation in our description.

#### Model comparison via AIC.

The general definition of AIC for a regression model with parameters β is given as where l(β) denotes the log likelihood of the model and df denotes the effective degrees of freedom of the model. Although df simply equals the number of regression coefficients in β for usual linear regression models, df in flexible models must take into account the effect of penalizing the coefficients. Generalizing results for linear models, the degrees of freedom are defined as the trace of the hat matrix projecting responses on predicted values in the working model that is employed in model estimation.

#### Penalized splines.

To model the univariate smooth functions *f*_{p,jk}(*x*) within the linear predictors listed in *Statistical analysis*, we consider penalized splines (8), since they provide a flexible and parsimonious representation of nonlinear covariate effects. To simplify notation, we will suppress the category index in the following discussion; i.e., we consider functions *f*(*x*). The basic idea of penalized splines is to represent *f*(*x*) as a polynomial spline of degree *l* based on a moderately large number of B-spline basis functions, i.e. (8) This, in principle, recasts the semiparametric model into a large linear model, since *Eq. 8* is linear in the parameters, but instead of estimating the B-spline coefficients β_{m} unrestricted via maximum likelihood, a penalty term is added to regularize the estimation problem. To be more specific, we want to ensure that the function estimate *f̂*(*x*), despite being flexible, is not too wiggly. From B-spline theory (5), we know that *k*th-order derivatives of B splines depend essentially on *k*th-order differences of the sequence of parameters β_{m}, *m* = 1,…, *M*. Therefore, if we want to ensure smooth function estimates in terms of the *k*th derivative, the log-likelihood has to be augmented by a penalty term constructed on the basis of squared *k*th-order differences, e.g. for first-order differences or for second-order differences. The smoothing parameter τ^{2} controls the trade-off between fidelity to the data (large τ^{2}) and smoothness (small τ^{2}), whereas the additional factor 2 is included merely for numerical convenience.

In vector notation, function evaluations *f*(*x*) can be represented as *f*(*x*) = *z*′β, where *z* = [B_{1}(*x*),…, B_{M}(*x*)] contains the evaluated B-spline basis functions and β = (β_{1},…, β_{M}) collects all regression coefficients. Similarly, the penalty term can be rewritten as a quadratic form (9) where the penalty matrix K = D′D is constructed from difference matrices D of appropriate order. To ensure identifiability within a multinomial logit model, the B-spline coefficients of the reference category *L* within the smooth function *f* ^{(L)}(*x*) = *z*^{(L)′}β^{(L)} are set to 0, i.e., β^{(L)} ≡ 0.

A smooth interaction surface *f*(*x*,*y*) for two continuous covariates *x* and *y* can be constructed via tensor product splines, which result from the pairwise multiplication of univariate splines for *x* and *y*, i.e. This yields the following presentation of *f*(*x*,*y*) where the design vector *z* contains the evaluated basis functions and β the corresponding regression coefficients. Penalty terms can be derived from the application of Kronecker products on corresponding univariate penalty matrices. A detailed discussion of bivariate smoothing splines can be found elsewhere (9).

#### Mixed-model-based inference.

The estimation approach used to fit semiparametric multinomial logit models is described in detail elsewhere (14, 15). Here we present a brief overview of the modeling steps needed to fit the model with algorithms known from mixed-model inference.

The semiparametric multinomial logit model can be fit via a penalized likelihood approach, as used for regularized linear models. This is possible since the semiparametric predictor can be transformed into a linear predictor (see *Penalized splines*). On the basis of generalized linear model theory, the log-likelihood function can be derived from the multinomial distribution induced by the multinomial logit model. To regularize estimation and enforce smooth functions, corresponding penalty terms, such as that presented in *Eq. 9*, are augmented to the log-likelihood function. Iterative optimization schemes, such as a penalized iteratively weighted least-squares algorithm (10), can then be used to obtain the regression coefficient estimates. However, this requires that the smoothness parameters τ^{2} are set in advance. To integrate the determination of appropriate smoothness parameters into an estimation approach, a different perspective on semiparametric regression models that allows determination of the smoothness parameters based on mixed-model methodology is needed.

The general idea is to interpret the penalty term (*Eq. 9*) as a random-effects distribution, where the smoothing parameter turns into a random-effects variance. For mixed models, estimation schemes for the determination of variances are readily available, e.g., maximum likelihood or marginal/restricted maximum likelihood. A truncated power series basis allows for a direct interpretation of semiparametric regression models as mixed models (the regression coefficients of truncated basis functions are defined as random effects), whereas the B-spline basis employed in our approach requires a further reparameterization step. This basically includes finding a proper density function for (parts of) the regression coefficient vector β of one smooth function *f*(*x*) = *z*′β. For simplicity, the category and covariate index are again suppressed.

The first guess for a density function for β reveals itself as partially improper, since *rk*(K) < dim(β), where *rk*(.) denotes the rank of a matrix. To obtain a valid mixed-model representation, it is necessary to reparameterize β to obtain a mixed model with proper random-effects distribution. It has been shown that this can be achieved by decomposing the vector β as (10) where the dimension of *b* is given by the rank of the penalty matrix, i.e., dim(*b*) = *rk*(K), while the dimension of γ is given by dim(γ) = dim(β) − *rk*(K). If appropriate design matrices U_{1} and U_{2} are chosen, we can achieve i.e., γ represents a vector of fixed effects and *b* represents a vector of independent-and-identically distributed random effects with common variance τ^{2}. Details on how to construct the design matrices U_{1} and U_{2} can be found elsewhere (15).

Applying decomposition (*Eq. 10*) to all nonparametric components in our semiparametric multinomial logit yields a mixed-model representation where all random effects have proper Gaussian distributions. Now we can apply mixed-model methodology to estimate not only the regression coefficients but also to automatically determine the smoothness parameters.

#### Confidence intervals for the probability surface.

As outlined (and adopted for the problem at hand) in *Statistical analysis*, the multinomial logit regression model models the dependency of the probability of occurrence of an event category *Y*_{i} = *r*, *r* = 1,…, *q* + 1 at observation *i* on a covariate vector *x*_{i}. For category *r* (*r* ≠ *q* + 1), the probability of occurrence π_{ir} is modeled by where β is the coefficient vector comprising all category-specific vectors β_{s}. An estimate β̂ can be obtained via maximum likelihood and, thus, an estimate for the probability of interest As a standard result of maximum-likelihood theory, β̂ is asymptotically normally distributed with mean β and inverse Fisher information matrix as variance-covariance matrix I(β), i.e. Via the multivariate Δ rule, it holds where ∇*h*_{r}(.) is the vector of partial derivatives of function *h*_{r}(.) with respect to β.

Hence, for the calculation of the asymptotic variance on the probability scale, the partial derivatives of *h*_{r}(β) must be calculated with respect to the category-specific vectors β_{k}, *k* = 1,…, *q*. The vector-valued derivative of *h*_{r}(β) with respect to β_{r} can be calculated as follows On the contrary, the vector-valued derivative of *h*_{r}(β) with respect to β_{k}, *k* ≠ *r* can be calculated as With these results, ∇*h*_{r}(β̂), the corresponding variance-covariance matrix for π̂_{ir} and point-wise confidence intervals based on the corresponding asymptotic normal distribution can be calculated for a given covariate vector *x*_{i} comprising a specific time and hormone value.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2009 the American Physiological Society