Abstract
Rapid eye movement (REM) sleep is a recurring state throughout the sleeping period. Based on the examination of 45 sleep records of 3moold male rats during the middle of the light phase, a stochastic model is proposed for the sequenceX _{1},Y _{2},X _{2},Y _{2}, . . . of REM sleep durations X and interREM sleep waiting times Y experienced by a rat during a sleeping period. In our model the probability distribution of any variable in the sequence, given the past, is allowed to depend on only the immediately previous variable. The conditional distributionsf(y_{i} ‖x_{i} ) andg(x _{i+1} ‖y_{i} ) do not depend on the indexi. It is shown that the marginal distributions tend to stationarity. Aggregations of the data on a discrete time scale suggest that the conditional distributions be formulated as twocomponent mixtures. These component distributions are modeled as Poisson and their means are called the means of short and long waiting time and the means of short and long REM sleep duration. Associated with each mean is a probability weight. Parametric forms are given to the means and probability weights. The model estimated by maximum likelihood shows a good fit to data of the 3moold rats. The model fit to a smaller data set obtained from rats aged 15–22 mo shows a significant shortening of the means for both short and long REM sleep bout durations compared with the means of the 3moold rats. Neuronal correlates for the behavior of the model are discussed in the context of the reciprocal interaction model of REM sleep regulation.
INTRODUCTION
The phenomenon of rapid eye movement (REM) sleep was first described by Aserinsky and Kleitman (1953) and was further characterized by Dement and Kleitman (1955). Since then, this state has been found and studied in a variety of mammalian and avian species (see reviews,Siegel 1995; Tobler 1995) but it has not been found unequivocally in nonhomeothermic species (seediscussion) (Hartse 1994). In all species where it has been found, REM sleep is an episodic state whose episode duration and interepisode interval are species dependent (Jones 1991). However, the REM sleep episode durations and the interREM sleep intervals are highly variable from episode to episode, even within a given individual. Previous research has attempted to explain the episodic nature of REM sleep in mammals both at the neural substrate level (Hobson et al. 1975) and in mathematical terms (McCarley and Hobson 1975; McCarley and Massaquoi 1986). These descriptions have tried to model a smooth cyclical function to the recurrence of REM sleep. The problem with such descriptions is that the high level of variation seen from one REM sleep episode to the next is not considered but instead is smoothed out.
Various statistical models have been developed to describe the sleep–wake cycle and statetostate transitions in adult humans (Achermann et al. 1993; Borbely 1982;Yang and Hursch 1973) and in infants (HolditchDavis et al. 1998). Markovian models have also been developed to describe the effects of temazapam on human sleep (Karlsson 2000). However, all these statistical models deal more with nonREM sleep and waking and do not investigate the REM sleep architecture in particular.
Previous studies have focused mainly on REM sleep latency (the period from first sleep to the first REM sleep bout), interREM sleep interval (the period from the beginning of one bout of REM sleep to the beginning of the next), the average duration of REM sleep bouts, and the proportion of total REM sleep. Structurally, these studies assume that the REM sleep cycle is sinusoidal, or at least highly periodic, with a specific frequency, and that variations from the basic cycle are due to small perturbations. In this article we present a descriptive probability model of REM sleep architecture that considers the variations in the previously mentioned quantities as important parts of the underlying mechanisms regulating REM sleep and not simply as perturbations. Furthermore, we compare the REM sleep of young and old rats, demonstrating that there are important changes to the REM sleep architecture of these animals with aging.
METHODS
Animals and the implantation of EEG and EMG electrodes for sleep recordings
The 31 male Sprague–Dawley rats weighing between 300 and 375 g (approximately 3 mo old) and 11 male Sprague–Dawley rats weighing between 600 and 670 g (aged between 15 and 22 mo) used for this analysis were obtained from Harlan (Indianapolis, IN). The animals were housed individually with a 12L:12D lighting schedule (lights on at 06:30) and food and water were available ad libidum. Animals were handled on a daily basis.
To implant the EEG and electromyogram (EMG) electrodes, animals were first anesthetized with an intramuscular injection of 2 ml/kg of a cocktail containing 6.3 mg/ml xylazine, 0.25 mg/ml acepromazine maleate, and 25 mg/ml ketamine HCl. Three screw EEG electrodes (1 frontal and 2 parietal) (Plastics One, Roanoke, VA), one screw reference electrode (occipital), and two EMG electrodes (Plastics One) (threaded into the nuccal muscles) were permanently implanted. The electrode ends were secured to a sixelectrode pedestal (Plastics One) and the whole assembly was secured to the calvarium using cranioplastic cement (Plastics One). Wounds were treated with lidocaine and a topical antibiotic (containing neomycin, polymixin B, and bacitracin) following surgery and animals were monitored until they regained consciousness. Wounds were cleaned and received antibiotic ointment daily until healed.
Sleep recordings and sleep record scoring
Animals were given 2 weeks to recover from surgery. During this time, animals were acclimated to the recording chamber and the presence of the recording cable (Plastics One). Following the recovery period, one or two baseline measures of sleep were recorded for each animal using a Grass Model 8 EEG recorder. Sleep was recorded on paper for 4 h starting at 10:00 ± 00:30. Records were manually scored and quantified for wakefulness, drowsy, slow wave sleep, and REM sleep in 6s epochs. Wakefulness was scored when the frontal/parietal electrodes showed >50% highfrequency alpha and beta waves (>8 Hz) of lowamplitude accompanied by a high muscle tone; drowsy was scored when the frontal/parietal electrodes showed a mixture of waves but contained <20% of delta waves (0.5–2 Hz) accompanied by a lower muscle tone; slow wave sleep was scored when the frontal/parietal electrodes showed >20% delta wave activity accompanied by low muscle tone; and REM sleep was scored when the frontal/parietal electrodes showed theta activity (4–7 Hz) accompanied by muscle atonia or very decreased muscle tone.
MODEL DEVELOPMENT
Conditional probability model
The measures of REM sleep latency, average interREM sleep interval, and average REM sleep duration give important information about what may be happening to REM sleep but they fail to describe the REM sleep architecture accurately. No studies have identified the notion of a waiting time between bouts of REM sleep (the time between the end of one REM sleep bout and the beginning of the next REM sleep bout) or considered the serial characteristics of the sequence of REM sleep bout durations and waiting times. We introduce the concept of the interREM sleep bout waiting time and consider the series of REM sleep bout durations and interREM sleep bout waiting times. Our quantification of a sleep record will begin with the first bout of REM sleep and will be denoted by the series
In our twostate probability model the probability distribution of any variable in the series (I) conditional on all previous variables in the series is a function of model parameters and the value of only the immediately previous variable. This is a Markovtype property and in a generalized sense the model could be called a Markov model. Let g(y‖x) denote the distribution of a waiting time y, given the valuex of the previous REM sleep duration, andf(x‖y) denote the distribution of a REM sleep duration x, given the value of the previous waiting time y. We take the distribution of the initial REM sleep duration to be of the formf(x
_{1}‖y
_{0}) where y
_{0} = ∞. Using a standard probability argument the joint probability distribution of Eq.1
given n is
λ_{1} = mean of short waiting time
λ_{2} = mean of long waiting time
1 − w = probability of short waiting time
w = probability of long waiting time
θ_{1} = mean of short REM sleep duration
θ_{2} = mean of long REM sleep duration
1 − v = probability of short REM sleep duration
v = probability of long REM sleep duration
We have therefore the statistically motivated construct of two states, short REM sleep or long REM sleep durations, having probabilities 1 − v and v, respectively. Similarly, a state of waiting is assumed to be either a state of short waiting or long waiting time with probabilities 1 − wand w.
Baseline behavior of 3moold rats
Thirtyone Sprague–Dawley rats, aged approximately 3 mo, were used in a study to determine baseline characteristics of the pattern of recurrence of REM sleep. Fortyfive sleep records of the form inEq. 1 were obtained. The average series length was 13. To induce more regularity into the data we discounted short REM sleep durations of <12 s, since these may actually represent an intermediate state and perhaps failed attempts to enter into REM sleep. Such durations were ignored, added to the waiting time, and the other REM durations were decreased by 12 s. These 45 sleep records were obtained from 31 different rats, 14 of which had duplicate records. We have treated these 45 records as mutually independent with the same probability distribution. This resulted from 1) a visual inspection of time plots that did not suggest any type of association between duplicate records and 2) a nonparametric analysis of REM durations and waiting times. The nonparametric analysis involved the four variables X _{1},Y _{1},X _{2},Y _{2}. For each variable the 45 values were ranked from 1 to 45 using average ranks when ties occurred. The ordered pairs of ranks for the 14 duplicate records were plotted in a triangular region. Under the hypothesis that all records are mutually independent with the same probability distribution, the points have approximately a uniform probability distribution over the triangle. With positive association between the pairs of values for duplicate records, points would tend to cluster in the acute corners of the triangle. The general impression for each of the four graphs suggests a uniform distribution with no clustering. We do not take this to be proof of the independence of multiple records for the same rat. Indeed, one may safely assume that they are not independent. Our conclusion is merely that the dependence is likely to be rather subtle. This issue bears further scrutiny. In this study, whose primary goal is to ascertain the characteristics of population distributions, we feel that more is to be gained by using as much data as possible than by adhering strictly to the independence of samples and suffering a loss of sample size.
The basic unit of measurement for REM sleep duration (X) and the waiting time (Y), is in units of 6 s. Since the underlying scale is continuous, it is natural to consider continuous time models. However, to induce additional regularity in this first treatment of a complex modeling problem, we impose a discretization of time measurements, with different units forX and Y. Henceforth, X values of 0,1,2 . . . will represent the time intervals (in seconds) (0,24], (24,48], etc., andY values of 0,1,2 . . . will represent the time intervals (in seconds) (0,96], (96,192], (192,288], etc. The widths of these intervals are smoothing parameters and were subjectively chosen to provide smoothing while retaining major features of the distributions. Observed counts for the combined Xs in the data and for the combined Ys in the data are shown in Figs.1 and 2. The major features of the X and Ydistributions are the relative concentrations at X = 0 and 4, and Y = 0 and 4. The presence of two modes for REM sleep duration and waiting time suggest mixture distributions. Using the simplest type of serial association in the sequence (Eq. 1 ), firstorder Markov, the mixture formulation (Eq. 3 ) arises.
To gain insight into the types of functions that should be considered for the weights and the component means in (Eq. 3 ), we examined the distribution of observed counts for the combined (X,Y) pairs and the combined (Y,X) pairs. These are shown in Figs. 3 A and4 A respectively. Inspecting Fig. 3 A sheds light on the nature ofg(y‖x), while inspecting Fig.4 A sheds light on the nature off(x‖y). In Fig. 3 A the second waiting time mode increases with REM duration and also becomes more prominent. In Fig. 4 A the REM duration modes do not change so much with waiting time and the prominence of the second mode may decrease. Following are the suggested characteristics:
the mean of short waiting time λ_{1} is constant,
the mean of long waiting time λ_{1} = λ_{2}(x) an increasing function ofx,
the probability of long waiting time w =w(x) is an increasing function of x,
the mean of short REM duration θ_{1} is constant,
the mean of long REM duration θ_{2} is constant,
the probability of long REM duration v =v(y) is a decreasing function of y.
The following parametric forms were adopted for those parameters that were not constant functions
The parameters of this model have been estimated by maximum likelihood. This procedure selects the parameter values that maximize the likelihood or probability of observing the actual data. The likelihood function L is the product of probabilities (Eq.2
) over all independent series of the form (Eq. 1
). This function is maximized over the parameter space by minimizing −2lnL, which is
It is of interest to test the following hypotheses:
H_{1}: the probability of a long waitingw(x) is a constant function,
H_{2}: the mean of long waiting λ_{2}(x) is a constant function,
H_{3}: the probability of long REM sleep durationv(y) is a constant function.
If H_{1} and H_{2} are true, thenX and the following Y are independent. If H_{3} is true, then Y and the followingX are independent. To test each hypothesis, we minimize −2lnL under the restriction of each hypothesis, to obtain
According to established principles of largesample hypothesis testing (Kendall and Stuart 1961), these may be compared to
We have now illustrated, for a onesample situation, estimation and hypothesis testing in the twostate stochastic model. This model was suggested by observing gross features in the data.
The calculation of exact expected numbers from the model for comparison with the observed numbers from data are not straightforward. The first issue is the nonstationary quality of the model: the marginal distributions of the Xs in an observed sequence (Eq.1 ) are not all the same and neither are the marginal distributions of the Ys. However, it is shown inappendix that the marginal distributions of theXs approach a stationary distribution, as do the marginal distributions of the Ys. These stationary limits can effectively be achieved very quickly. Our model is then a pseudostationary model.
The expected counts shown in Figs. 14 were obtained from stationary model probabilities using parameter values estimated from the data. In each figure an expected count is a model probability times the corresponding total observed count (number of X values for Fig. 1, number of Y values for Fig. 2, number of (X,Y) pairs for Fig. 3, and number of (Y,X) pairs for Fig. 4). We have also obtained expected counts from simulations. Using parameter values from the estimated full model, 5,000 sequences (Eq. 1 ) were simulated with n = 13 (the average value in the data). Plots of counts from the simulation adjusted to have the same totals as the data are indistinguishable from the expected counts in Figs. 14. The level of agreement between the expected counts and the observed counts seems to be generally supportive of our model. We have not calculated goodnessoffit statistics to test the model, principally because we believe that the model needs further refinement. The challenge for future work is to consider model distributions other than Poisson, even continuous ones, and other devices, to obtain a closer fit of expected counts to observed counts. At this juncture, we feel that the fit is close enough to make our model useful for data summary and interpretation, as well as hypothesis testing. Since there is some lack of fit, we have made the likelihood procedure more robust by replacing the function ln(probability) in (Eq. 5 ) by ln(probability + 0.0001).
Table 1 presents serial correlations in the actual 45 sequences and in the simulated 5,000 sequences. Values from the data are weighted averages of serial correlations for individual rats, weighted by the number of data pairs present for the calculation of a given serial correlation. Significance at the 0.001 level for the sign test of zero population correlation, applied to the individual values that went into a weighted average is indicated. This table has a twofold purpose: 1) to establish the existence of serial correlation in a nonparametric setting and 2) to offer further comparison between the data and the model.
Comparison with rats aged 15 to 22 mo
The previous section established the relevancy of the twostate stochastic model in describing the recurrence of REM sleep in 3moold rats. In this section we consider additional data consisting of 11 sleep records for 11 different rats aged 15–22 mo, which were treated identically to the younger rats. We assume that the model, with possibly different parameter values, applies to these data. The purpose of this section is to illustrate the use of our model in making group comparisons by comparing rats in the two age groups.
Maximum likelihood estimation of the model using the data from the older rats produced the results shown in appendix
. For the older rats, the probability of long REM sleep duration depends only slightly on the previous waiting time. Since this is the only quantity linking REM sleep duration to the previous waiting time, it is seen that REM sleep duration is almost independent of the previous waiting time, in the model for the older rats. Figures 5 and6, discussed earlier with respect to the model for the young rats, also present the model for the older rats. The means of long and short REM sleep, as well as the weighted average, is less for the older animals, but it is smoother, as more of it is in the long REM sleep state. To test this apparent shortening of REM sleep duration in the older animals, a parametrization for the combined groups needs to be chosen. For the younger group we set w
_{0} = 1 and retain the other 11 parameters. In the older group, we setw
_{0} = 1,v
_{1} =v
_{2} = 0 and retain the other 9 parameters. The full model for the combined sample has 20 parameters. Since −2lnL for the combined sample is additive over groups
DISCUSSION
We have developed a statistical model for the recurrence of REM sleep in the rat that does not describe sleep as a sinusodal function with a given frequency or as a limited cycle (McCarley and Massaquoi 1986). We have introduced the concept of a waiting time and feel that this concept is more valuable for the study of sleep architecture and is more in line with the reciprocal interaction model of REM sleep regulation (McCarley and Hobson 1975) than the concept of a repeating interval. The model presented here can be understood within the context of the reciprocal interaction model of REM sleep control proposed by McCarley and Hobson (1975). This model explains how the REM sleep cycle, both in terms of REM sleep duration and waiting time, may be generated by several pontine nuclei (see review Datta 1995). In this neuronal model, neurons of the dorsal raphe (DR) and locus coeruleus (LC) nuclei inhibit the firing of the laterodorsal tegmental (LDT) and pedunculopontine (PPT) neurons, which are critical to the generation of REM sleep. The model proposes that the neurons of the DR and LC have an inhibitory collateral autofeedback that eventually stops their own activity and allows the neurons of the LDT/PPT to gain activity and generate REM sleep (Datta 1995; McCarley and Hobson 1975). The strength of these autofeedback connections accounts for the length of the waiting time. The model further proposes that the neurons of the LDT/PPT send excitatory projections to the DR and LC. These excitatory projections eventually cause the neurons of the DR and LC to regain activity, reestablishing the inhibition of the LDT/PPT neurons and thus ending the REM sleep bout. In this neuronal model, the function of these excitatory projections determines the length of the REM sleep bout.
Our model is based on the observation that REM sleep bout durations are clustered around either a short duration time (12 + 24θ_{1} = 18.6 s) or a longer duration time (12 + 24θ_{2} = 101.5 s), suggesting that the excitatory connections from the LDT/PPT to the DR and LC may exist in one of two general states. A strong connectivity would account for the short durations and a weak connectivity would account for the long durations. Further work on the state of these cholinergic synapses should help to understand if two states predominate, which is an important prediction of our model.
In the 3moold rats, the probability of having a long REM sleep bout duration is a decreasing function of the waiting time, having an asymptotic value of 1 − v _{0} = 0.59. The probability of having a short REM sleep bout duration increases to the asymptotic value of 0.41. This indicates that the long REM sleep bout duration state of the excitatory connections is more stable than that of the short duration.
The model suggests a more complex regulation of the autofeedback connections of the DR and LC, which control the length of the waiting time. The waiting times are distributed about a short waiting time mean of 48 + 96λ_{1} s and a long waiting time mean that increases from 48 + 96a _{0}(1 − a _{1}) = 409 s to 48 + 96a _{0} = 728 s as the previous REM sleep bout duration lengthens. The probability of having a short waiting time is highest [1 −w _{0}(1 −w _{1}) = 0.51] when the preceding REM sleep bout duration is the shortest and decreases to 0.10 as the preceding REM sleep bout duration increases. In the case of the long waiting time the opposite is true, with its probability increasing from 0.49 to 0.90 as the preceding REM sleep bout duration increases. Because a short waiting time should result from having a strong collateral inhibition of the DR and LC, it is possible that a more generalized residual inhibition occurs in the cells of the DR and LC and that this residual inhibition dissipates with the length of the REM sleep bout duration time. Thus short REM sleep bout durations would affect this residual inhibition less, making the following waiting time shorter since activation of the cells from the DR and LC would be less efficacious. Furthermore, since the longer a REM sleep bout lasts the long waiting time mean also increases, the model suggests that the strength of the inhibitory feedback weakens during the actual REM sleep bout. This weakening of the inhibitory connection must be somewhat specific to the serotonergic and noradrenergic connections since the two means of the REM sleep bout durations are not affected by timing dynamics. If the weakening in the inhibitory connections were more globally caused, then one would expect to see the length of means of REM bout durations also affected.
The magnitude and statistical significance to changes in the amount of REM sleep seen with aging has been unclear. Several studies have reported changes to the amount, duration, and interREM sleep intervals of REM sleep in aging rats (ArankowskySandoval et al. 1992; Crespi 1999; Markowska et al. 1989; Rosenberg et al. 1979; Stone et al. 1989, 1992; Witting et al. 1993). These studies report significant decreases in the amount of REM sleep that aged rats have. However, other studies have failed to find these REM sleep changes in older rats (Li and Satinoff 1995;Mendelson and Bergman 1999; Ruigt et al. 1989). Our model indicates (see Figs. 5 and 6) that both REM sleep bout duration and waiting time are reduced in the older rats. Moreover, we show that this reduction in REM sleep bout duration is significant (P = 0.00005). Thus our work indicates that, while older rats have shorter bouts of REM sleep than younger rats, they return to REM sleep more quickly, so that the percentage of sleep that is REM sleep may not be much different for the older animals than it is for the younger ones.
The shortening of the short and long REM sleep bout duration means with age would suggest that the strength of cholinergic connections at the DR and LC increase in late life while maintaining the existence of two predominant states. In the aged rats, however, the long REM sleep bout duration state is even more stable than that of young animals, since the probability of having a long REM sleep bout duration is stable at 0.804 [1 −v _{0}(1 −v _{1})], regardless of the length of the preceding waiting time. Although the long REM sleep bout duration predominates, the mean REM sleep bout duration is always shorter in the aged rats than in the young rats (see Fig. 6).
We feel that the proposed twostate stochastic model of REM sleep in the rat has great utility both in terms of suggesting ways that one might study the regulation of REM sleep architecture and in terms of understanding experimental manipulations that alter REM sleep. The site of action of specific pharmacological agents may be gleaned by observing what factors of the model are affected and in which direction. Although we have used the reciprocal interaction model of REM sleep regulation to add light to the probabilistic model, the present model does not depend on the accuracy of the neuronal model proposed by McCarley and Hobson (1975). Many other factors play a role in how these pontine centers behave to regulate REM sleep [Orexin, Somatostatin, vasoactive intestinal peptide (VIP), Anandimide], however, we feel that the use of this probabilistic model will aid in understanding how these other factors and brain centers regulate the REM sleep architecture.
Acknowledgments
We acknowledge the technical assistance of C. Moreno and we thank J. JosephVanderpool and J. Arteaga for the kind use of Model 8 Grass EEG machines. Careful reading by referees has resulted in a much improved paper.
This work was supported by Division of Research Grant G12RR08124 and National Institute of General Medical Services Grant SO6GM–08012.
Footnotes

Address for reprint requests: G. G. Gregory, Department of Mathematical Sciences, The University of Texas at El Paso, El Paso, TX 79968 (Email: gavin{at}math.utep.edu).
 Copyright © 2002 The American Physiological Society
Appendix
Let the distribution of X
_{1} be