JN Journal of Neurophysiology
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


J Neurophysiol 88: 2589-2597, 2002; doi:10.1152/jn.00861.2001
0022-3077/02 $5.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via Web of Science (7)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Gregory, G. G.
Right arrow Articles by Cabeza, R.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Gregory, G. G.
Right arrow Articles by Cabeza, R.

J Neurophysiol (November 1, 2002). 10.1152/jn.00861.2001
Submitted on 19 October 2001
Accepted on 24 July 2002

A Two-State Stochastic Model of REM Sleep Architecture in the Rat

Gavin G. Gregory1 and Rafael Cabeza2

 1Department of Mathematical Sciences and  2Department of Biological Sciences, The University of Texas at El Paso, El Paso, Texas 79968


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Gregory, Gavin G. and Rafael Cabeza. A Two-State Stochastic Model of REM Sleep Architecture in the Rat. J. Neurophysiol. 88: 2589-2597, 2002. Rapid eye movement (REM) sleep is a recurring state throughout the sleeping period. Based on the examination of 45 sleep records of 3-mo-old male rats during the middle of the light phase, a stochastic model is proposed for the sequence X1,Y2, X2,Y2, . . . of REM sleep durations X and inter-REM 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 distributions f(yixi) and g(xi+1 yi) do not depend on the index i. 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 two-component 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 3-mo-old 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 3-mo-old rats. Neuronal correlates for the behavior of the model are discussed in the context of the reciprocal interaction model of REM sleep regulation.


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

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 (see DISCUSSION) (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 inter-REM 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 state-to-state transitions in adult humans (Achermann et al. 1993; Borbely 1982; Yang and Hursch 1973) and in infants (Holditch-Davis 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 non-REM 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), inter-REM 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
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

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 six-electrode 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 6-s epochs. Wakefulness was scored when the frontal/parietal electrodes showed >50% high-frequency alpha and beta waves (>8 Hz) of low-amplitude 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
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Conditional probability model

The measures of REM sleep latency, average inter-REM 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 inter-REM sleep bout waiting time and consider the series of REM sleep bout durations and inter-REM 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
(<IT>X</IT><SUB><IT>1</IT></SUB><IT>, </IT><IT>Y</IT><SUB><IT>1</IT></SUB><IT>, </IT><IT>X</IT><SUB><IT>2</IT></SUB><IT>, </IT><IT>Y</IT><SUB><IT>2</IT></SUB><IT>,…, </IT><IT>X<SUB>n</SUB></IT><IT>, </IT><IT>Y<SUB>n</SUB></IT>) (1)
where the random variable Xi is the REM sleep duration of the ith bout, the random variable Yi is the waiting time until the next REM sleep bout, and n +1 is the random number of REM sleep bouts in a period of sleep. In this paper we introduce a parametric two-state probability model for the stochastic series (Eq. 1) conditional on the value of n. Then maximum likelihood significance testing is utilized for assessing group or treatment effects in experimental data of this kind. This methodology is illustrated later in the text in the comparison of 3-mo-old rats with similar rats aged 15 to 22 mo.

In our two-state 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 Markov-type 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 value x of the previous REM sleep duration, and f(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 form f(x1|y0) where y0 = infinity . Using a standard probability argument the joint probability distribution of Eq. 1 given n is
<IT>f</IT>(<IT>x</IT><SUB><IT>1</IT></SUB><IT>‖</IT><IT>y</IT><SUB><IT>0</IT></SUB>)<IT>g</IT>(<IT>y</IT><SUB><IT>1</IT></SUB><IT>‖</IT><IT>x</IT><SUB><IT>1</IT></SUB>)<IT>f</IT>(<IT>x</IT><SUB><IT>2</IT></SUB><IT>‖</IT><IT>y</IT><SUB><IT>1</IT></SUB>)<IT> … </IT><IT>f</IT>(<IT>x<SUB>n</SUB></IT><IT>‖</IT><IT>y</IT><SUB><IT>n</IT><IT>−1</IT></SUB>)<IT>g</IT>(<IT>y<SUB>n</SUB></IT><IT>‖</IT><IT>x<SUB>n</SUB></IT>) (2)
The bimodality of the data presented in Baseline behavior of 3-mo-old rats led us to consider a mixture model for the conditional distributions
<IT>g</IT>(<IT>y</IT><IT>‖</IT><IT>x</IT>)<IT>=</IT>(<IT>1−</IT><IT>w</IT>)<IT>p</IT><SUB><IT>1</IT></SUB>(<IT>y</IT><IT>‖</IT><IT>x</IT>)<IT>+</IT><IT>wp</IT><SUB><IT>2</IT></SUB>(<IT>y</IT><IT>‖</IT><IT>x</IT>)

f(x‖y)=(1−v)<IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT><IT>‖</IT><IT>y</IT>)<IT>+</IT><IT>vq</IT><SUB><IT>2</IT></SUB>(<IT>x</IT><IT>‖</IT><IT>y</IT>) (3)
where the means of the distributions p1 and p2 will be denoted lambda 1 and lambda 2, which along with weight w may depend on x and the means of the distributions q1 and q2 will be denoted theta 1 and theta 2, which along with weight v may depend on y. We take lambda 1 < lambda 2 and theta 1 < theta 2. The parameters are described as follows:

lambda 1 = mean of short waiting time

lambda 2 = mean of long waiting time

- w = probability of short waiting time

w = probability of long waiting time

theta 1 = mean of short REM sleep duration

theta 2 = mean of long REM sleep duration

- 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 - w and w.

Baseline behavior of 3-mo-old rats

Thirty-one Sprague-Dawley rats, aged approximately 3 mo, were used in a study to determine baseline characteristics of the pattern of recurrence of REM sleep. Forty-five sleep records of the form in Eq. 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 X1, Y1, X2, Y2. 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 for X and Y. Henceforth, X values of 0,1,2 . . . will represent the time intervals (in seconds) (0,24], (24,48], etc., and Y 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 Y distributions 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), first-order Markov, the mixture formulation (Eq. 3) arises.



View larger version (16K):
[in this window]
[in a new window]
 
Fig. 1. Frequency distributions of both observed and model predicted REM sleep bout durations. REM sleep bout durations are plotted in seconds and by category with the solid circles representing the observed distribution and the empty circles representing the model predicted distributions. Each category interval represents 24 s and points are plotted in the middle of each interval.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 2. Frequency distributions of both observed and model predicted inter-REM sleep bout waiting times. Waiting times are plotted in seconds and by category with the solid circles representing the observed distribution and the empty circles representing the model predicted distributions. Each category interval represents 96 s and points are plotted in the middle of each interval.

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. 3A and 4A respectively. Inspecting Fig. 3A sheds light on the nature of g(y|x), while inspecting Fig. 4A sheds light on the nature of f(x|y). In Fig. 3A the second waiting time mode increases with REM duration and also becomes more prominent. In Fig. 4A 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:



View larger version (51K):
[in this window]
[in a new window]
 
Fig. 3. Distribution of observed and model predicted inter-REM sleep waiting times following a given REM sleep bout duration. The surface plot shown in A represents the observed distributions while the plot shown in B represents the model predicted distribution.



View larger version (58K):
[in this window]
[in a new window]
 
Fig. 4. Distribution of observed and model predicted REM sleep bout durations following a given inter-REM sleep waiting time. The surface plot shown in A represents the observed distributions while the plot shown in B represents the model predicted distribution.

the mean of short waiting time lambda 1 is constant,

the mean of long waiting time lambda 1 = lambda 2(x) an increasing function of x,

the probability of long waiting time w = w(x) is an increasing function of x,

the mean of short REM duration theta 1 is constant,

the mean of long REM duration theta 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
&lgr;<SUB>2</SUB>(<IT>x</IT>)<IT>=</IT><IT>a</IT><SUB><IT>0</IT></SUB>(<IT>1−</IT><IT>a</IT><SUB><IT>1</IT></SUB><IT>e</IT><SUP><IT>−</IT><IT>a</IT><SUB><IT>2</IT></SUB><IT>x</IT></SUP>)

<IT>w</IT>(<IT>x</IT>)<IT>=</IT><IT>w</IT><SUB><IT>0</IT></SUB>(<IT>1−</IT><IT>w</IT><SUB><IT>1</IT></SUB><IT>e</IT><SUP><IT>−</IT><IT>w</IT><SUB><IT>2</IT></SUB><IT>x</IT></SUP>)

<IT>v</IT>(<IT>y</IT>)<IT>=1−</IT><IT>v</IT><SUB><IT>0</IT></SUB>(<IT>1−v<SUB>1</SUB>e</IT><SUP><IT>−</IT><IT>v</IT><SUB><IT>2</IT></SUB><IT>y</IT></SUP>) (4)
where all new parameters are nonnegative and v0, v1, w0, w1, and a1 are less than or equal to one. In this paper, the component distributions p1, p2, q1, and q2 were taken to be Poisson. This was done for two reasons: 1) the Poisson distribution is one of the most utilized discrete distributions for general purpose modeling and 2) the Poisson distribution has but a single parameter, its mean, and thus keeps the total number of parameters to a minimum. There are altogether 12 parameters in the parameter space of the model defined by Eqs. 2, 3, and 4.

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>2 </IT><LIM><OP>∑</OP><LL><IT>i</IT><IT>=1</IT></LL><UL><IT>n</IT></UL></LIM> [<IT>ln </IT><IT>f</IT>(<IT>x<SUB>i</SUB></IT><IT>‖</IT><IT>y</IT><SUB><IT>i</IT><IT>−1</IT></SUB>)<IT>+ln </IT><IT>g</IT>(<IT>y<SUB>i</SUB></IT><IT>‖</IT><IT>x<SUB>i</SUB></IT>)] (5)
summed over all independent series. Function minimization was via the FORTRAN subroutine BCONF in the IMSL collection, available with Compaq Visual Fortran (1999). The estimated model parameters, along with the minimum value of -2lnL, are given in APPENDIX B. A FORTRAN program for obtaining such values is available from the first author. In the estimated model, the association between X and the following Y is positive. The association between Y and the following X is negative. The conditional mean waiting time is (1 - w)lambda 1 + wlambda 2 and is shown in Fig. 5 along with the means lambda 1 and lambda 2 of short and long waiting. The conditional mean REM sleep bout duration is (1 - v)theta 1 + vtheta 2 and is shown in Fig. 6 along with the means theta 1 and theta 2 of short and long REM sleep bout duration. Such a pair of graphs provides a useful summary of model characteristics. At any abscissa the probability of the long state is the distance between the short and the average curve, expressed as a proportion of the distance between the short and long state curves. In the next section we compare these graphs to ones obtained by fitting our model to data from older rats.



View larger version (24K):
[in this window]
[in a new window]
 
Fig. 5. Model mean waiting times as a function of previous REM sleep duration. The dark symbols show the model means for the young rats while the empty symbols show the information for the old rats. The circles are the model means for the short waiting times, the squares are the model means for the long waiting times, and the triangles are the model means for the expected inter-REM sleep waiting time following a given length of REM sleep.



View larger version (27K):
[in this window]
[in a new window]
 
Fig. 6. Model mean REM sleep durations as a function of previous inter-REM waiting time. The dark symbols show the model means for the young rats while the empty symbols show the information for the old rats. The circles are the model means for the short REM sleep bouts, the squares are the model means for the long REM sleep bouts, and the triangles are the model means for the expected REM sleep bout following a given length of waiting time.

It is of interest to test the following hypotheses:

H1: the probability of a long waiting w(x) is a constant function,

H2: the mean of long waiting lambda 2(x) is a constant function,

H3: the probability of long REM sleep duration v(y) is a constant function.

If H1 and H2 are true, then X and the following Y are independent. If H3 is true, then Y and the following X are independent. To test each hypothesis, we minimize -2lnL under the restriction of each hypothesis, to obtain
−<IT>2 ln </IT><IT>L</IT>(<IT>w</IT><SUB><IT>1</IT></SUB><IT>=</IT><IT>w</IT><SUB><IT>2</IT></SUB><IT>=0</IT>)<IT>=5184.8</IT>

−<IT>2 ln </IT><IT>L</IT>(<IT>a</IT><SUB><IT>1</IT></SUB><IT>=</IT><IT>a</IT><SUB><IT>2</IT></SUB><IT>=0</IT>)<IT>=5222.5</IT>

−<IT>2 ln </IT><IT>L</IT>(<IT>v</IT><SUB><IT>1</IT></SUB><IT>=</IT><IT>v</IT><SUB><IT>2</IT></SUB><IT>=0</IT>)<IT>=5169.6</IT>

According to established principles of large-sample hypothesis testing (Kendall and Stuart 1961), these may be compared to
−<IT>2 ln </IT><IT>L</IT>(<IT>full model</IT>)<IT>=5162.2</IT>
to obtain test statistics for the respective hypotheses. The bigger a difference in the value of -2lnL the more untenable the hypothesis is. Under a hypothesis, the associated difference in the value of -2lnL will have a distribution which is approximately chi 2 with 2 df the reduction in the dimensionality of the parameter space under the hypothesis. P values obtained from the chi 2 distribution are
<IT>P</IT>(<IT>H<SUB>1</SUB></IT>)<IT><0.00005</IT>

<IT>P</IT>(<IT>H<SUB>2</SUB></IT>)<IT><0.00005</IT>

<IT>P</IT>(<IT>H<SUB>3</SUB></IT>)<IT>=0.025</IT>
The data, therefore possess sufficient evidence to reject H1, H2, and H3.

We have now illustrated, for a one-sample situation, estimation and hypothesis testing in the two-state 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 in APPENDIX A that the marginal distributions of the Xs 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. 1-4 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. 1-4. The level of agreement between the expected counts and the observed counts seems to be generally supportive of our model. We have not calculated goodness-of-fit 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.


                              
View this table:
[in this window]
[in a new window]
 
Table 1. Serial correlations in the sequence of REM sleep durations and waiting times

Comparison with rats aged 15 to 22 mo

The previous section established the relevancy of the two-state stochastic model in describing the recurrence of REM sleep in 3-mo-old 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 B. 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 and 6, 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 w0 = 1 and retain the other 11 parameters. In the older group, we set w0 = 1, v1 = v2 = 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
−<IT>2 ln </IT><IT>L</IT>(<IT>full model</IT>)<IT>=5162.2+1354.7=6516.9</IT>
Consider testing the hypothesis that the population means of long and short REM sleep durations for the older animals are identical to the values for the younger animals. Under this hypothesis, the dimensionality of the parameter space is reduced by two. Minimization of -2lnL under this hypothesis produces
−<IT>2 ln </IT><IT>L</IT>(<IT>same REM means</IT>)<IT>=6537.6</IT>
Therefore (6537.6 - 6516.9) = 20.7 is the value of a variate that has, under the hypothesis, a distribution that is approximately chi 2 with 2 df degrees of freedom. The P value determined from the chi 2 distribution is <0.00005 and the hypothesis is rejected.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

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 + 24theta 1 = 18.6 s) or a longer duration time (12 + 24theta 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 3-mo-old 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 - v0 = 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 auto-feedback 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 + 96lambda 1 s and a long waiting time mean that increases from 48 + 96a0(1 - a1) = 409 s to 48 + 96a0 = 728 s as the previous REM sleep bout duration lengthens. The probability of having a short waiting time is highest [1 - w0(1 - w1) = 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 inter-REM sleep intervals of REM sleep in aging rats (Arankowsky-Sandoval 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 - v0(1 - v1)], 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 two-state 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.


    APPENDIX A
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES

Let the distribution of X1 be
<IT>f</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)<IT>=</IT>(<IT>1−</IT><IT>v</IT><SUP>(<IT>1</IT>)</SUP>)<IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)<IT>+</IT><IT>v</IT><SUP>(<IT>1</IT>)</SUP><IT>q</IT><SUB><IT>2</IT></SUB>(<IT>x</IT>)<IT>=</IT><IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)<IT>+</IT><IT>v</IT><SUP>(<IT>1</IT>)</SUP>[<IT>q</IT><SUB><IT>2</IT></SUB>(<IT>x</IT>)<IT>−</IT><IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)]
The distribution of Y1 is
<IT>g</IT><SUB><IT>1</IT></SUB>(<IT>y</IT>)<IT>=</IT><LIM><OP>∑</OP><LL><IT>x</IT></LL></LIM> <IT>g</IT>(<IT>y</IT><IT>‖</IT><IT>x</IT>)<IT>f</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)

=<LIM><OP>∑</OP><LL><IT>x</IT></LL></LIM> <IT>g</IT>(<IT>y</IT><IT>‖</IT><IT>x</IT>)<IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)<IT>+</IT><IT>v</IT><SUP>(<IT>1</IT>)</SUP> <LIM><OP>∑</OP><LL><IT>x</IT></LL></LIM> <IT>g</IT>(<IT>y</IT><IT>‖</IT><IT>x</IT>)[<IT>q</IT><SUB><IT>2</IT></SUB>(<IT>x</IT>)<IT>−</IT><IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)]
The distribution of X2 is
<IT>f</IT><SUB><IT>2</IT></SUB>(<IT>x</IT>)<IT>=</IT><LIM><OP>∑</OP><LL><IT>y</IT></LL></LIM> <IT>f</IT>(<IT>x</IT><IT>‖</IT><IT>y</IT>)<IT>g</IT><SUB><IT>1</IT></SUB>(<IT>y</IT>)

≡<IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)<IT>+</IT><IT>v</IT><SUP>(<IT>2</IT>)</SUP>[<IT>q</IT><SUB><IT>2</IT></SUB>(<IT>x</IT>)<IT>−</IT><IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)]
where
<IT>v</IT><SUP>(<IT>2</IT>)</SUP><IT>=</IT><LIM><OP>∑</OP><LL><IT>y</IT></LL></LIM> <IT>v</IT>(<IT>y</IT>)<IT>g</IT><SUB><IT>1</IT></SUB>(<IT>y</IT>)

=<LIM><OP>∑</OP><LL><IT>x</IT></LL></LIM> <LIM><OP>∑</OP><LL><IT>y</IT></LL></LIM> <IT>v</IT>(<IT>y</IT>)<IT>g</IT>(<IT>y</IT><IT>‖</IT><IT>x</IT>)<IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)<IT>+</IT><IT>v</IT><SUP>(<IT>1</IT>)</SUP> <LIM><OP>∑</OP><LL><IT>x</IT></LL></LIM> <LIM><OP>∑</OP><LL><IT>y</IT></LL></LIM> <IT>v</IT>(<IT>y</IT>)<IT>g</IT>(<IT>y</IT><IT>‖</IT><IT>x</IT>)[<IT>q</IT><SUB><IT>2</IT></SUB>(<IT>x</IT>)<IT>−</IT><IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)]

≡<IT>a</IT><IT>+</IT><IT>bv</IT><SUP>(<IT>1</IT>)</SUP>
We have 0 < a < 1 and 0 < a + b < 1. Now we can proceed recursively. The sequence of weights v(1), v(2), v(3), . . . converges to a/(1 - b). The mean of X converges to &thgr;<SUB>1</SUB>+<FR><NU><IT>a</IT></NU><DE>1−<IT>b</IT></DE></FR> (&thgr;<SUB>2</SUB>−&thgr;<SUB>1</SUB>). The mean of Y converges to
<LIM><OP>∑</OP><LL><IT>x</IT></LL></LIM> [(<IT>1−</IT><IT>w</IT>(<IT>x</IT>))<IT>&lgr;<SUB>1</SUB>+</IT><IT>w</IT>(<IT>x</IT>)<IT>&lgr;<SUB>2</SUB></IT>(<IT>x</IT>)]<IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)<IT>+</IT><FR><NU><IT>a</IT></NU><DE><IT>1−</IT><IT>b</IT></DE></FR> <LIM><OP>∑</OP><LL><IT>x</IT></LL></LIM> [(<IT>1−</IT><IT>w</IT>(<IT>x</IT>))<IT>&lgr;<SUB>1</SUB>+</IT><IT>w</IT>(<IT>x</IT>)<IT>&lgr;<SUB>2</SUB></IT>(<IT>x</IT>)][<IT>q</IT><SUB><IT>2</IT></SUB>(<IT>x</IT>)<IT>−</IT><IT>q</IT><SUB><IT>1</IT></SUB>(<IT>x</IT>)].


    APPENDIX B
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES


                              
View this table:
[in this window]
[in a new window]
 


    ACKNOWLEDGMENTS

We acknowledge the technical assistance of C. Moreno and we thank J. Joseph-Vanderpool 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 G12RR-08124 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 (E-mail: gavin{at}math.utep.edu).


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
METHODS
MODEL DEVELOPMENT
DISCUSSION
APPENDIX A
APPENDIX B
REFERENCES


0022-3077/02 $5.00 Copyright © 2002 The American Physiological Society




This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via Web of Science (7)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Gregory, G. G.
Right arrow Articles by Cabeza, R.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Gregory, G. G.
Right arrow Articles by Cabeza, R.


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Visit Other APS Journals Online