A probabilistic model for the ultradian timing of REM sleep in mice.

A salient feature of mammalian sleep is the alternation between rapid eye movement (REM) and non-REM (NREM) sleep. However, how these two sleep stages influence each other and thereby regulate the timing of REM sleep episodes is still largely unresolved. Here, we developed a statistical model that specifies the relationship between REM and subsequent NREM sleep to quantify how REM sleep affects the following NREM sleep duration and its electrophysiological features in mice. We show that a lognormal mixture model well describes how the preceding REM sleep duration influences the amount of NREM sleep till the next REM sleep episode. The model supports the existence of two different types of sleep cycles: Short cycles form closely interspaced sequences of REM sleep episodes, whereas during long cycles, REM sleep is first followed by an interval of NREM sleep during which transitions to REM sleep are extremely unlikely. This refractory period is characterized by low power in the theta and sigma range of the electroencephalogram (EEG), low spindle rate and frequent microarousals, and its duration proportionally increases with the preceding REM sleep duration. Using our model, we estimated the propensity for REM sleep at the transition from NREM to REM sleep and found that entering REM sleep with higher propensity resulted in longer REM sleep episodes with reduced EEG power. Compared with the light phase, the buildup of REM sleep propensity was slower during the dark phase. Our data-driven modeling approach uncovered basic principles underlying the timing and duration of REM sleep episodes in mice and provides a flexible framework to describe the ultradian regulation of REM sleep in health and disease.

Introduction

During sleep, the mammalian brain alternates between two distinct states—rapid eye movement (REM) sleep and non-REM (NREM) sleep. The cyclic occurrence of REM sleep constitutes the REM-NREM or sleep cycle, an ultradian rhythm on a minute-to-hour time scale shared by mammals [1], birds [2], and reptiles [3,4]. Although we know in great detail about the neural mechanisms underlying oscillations on a millisecond-to-second timescale [5] or about transcriptional/translational oscillators generating circadian rhythms [6], we lack knowledge of how the brain generates ultradian rhythms.

Transitions from NREM to REM sleep are thought to be controlled by a network of interconnected REM sleep-promoting (REM-on) and REM sleep-suppressing (REM-off) neurons [7-10]. Research in the last decade has identified key populations of REM-on and REM-off neurons in the brainstem and hypothalamus that powerfully promote or suppress REM sleep and has mapped their connectivity at unprecedented detail [11-18]. However, the mechanisms that regulate when the brain state transitions from NREM to REM sleep and thereby determine the duration of the sleep cycle are still largely unknown.

A common statistical feature of mammalian sleep, observed in multiple species including humans is that the duration of REM sleep is positively correlated with the time till the next REM sleep period (inter-REM interval) [19-25]. This correlation is thought to be the manifestation of a homeostatic process that operates on the ultradian time scale [19,20]: A propensity for REM sleep builds up during the inter-REM interval and partially discharges during REM sleep. Longer REM periods, therefore, lead to a stronger discharge of the REM propensity and thus require longer intervals for accumulation to re-enter REM sleep. According to this model, the sleep cycle is not generated by an oscillator circuit, as originally proposed [26], but rather is the consequence of an hourglass-type ultradian process [19,20,27,28], the neural or molecular correlates of which are still largely unknown.

Recently, the hourglass process has been proposed to be the result of a refractory period following REM sleep, during which transitions to REM sleep are effectively suppressed [29-31]. The probability of NREM to REM sleep transitions is reduced particularly after long REM sleep periods (>1 min) in rats [32], which may explain the positive correlation between the duration of REM sleep and the inter-REM interval. However, although possibly being a crucial subunit within the sleep cycle, there is still no quantitative definition of the refractory period, nor is there an understanding of how the quality of NREM sleep and its microarchitecture may be changed during and after the refractory period. Furthermore, as the duration of REM sleep influences the duration of the following inter-REM interval, it is crucial for our understanding of the ultradian sleep cycle to identify factors that regulate the duration of REM sleep episodes.

The sleep pattern in mammals is further complicated by sequences of REM sleep episodes, whose timing deviates from the statistics expected from an hourglass process. In rats, such sequences comprise several temporally close (< 3 min) REM sleep periods [33-35] and their presence results in a bimodal distribution of the inter-REM interval durations. In addition to the rat, REM sleep sequences have been reported in multiple mammalian species including humans [21,36,37], suggesting that they are a conserved phenomenon in mammalian sleep. An increased frequency of REM sleep sequences underlies the homeostatic rebound following a loss in REM sleep induced by exposure to cold temperatures [33,34] and is characteristic of sleep in stressed animals [38,39]. However, the mechanisms underlying the induction of REM sleep sequences and their defining electrophysiological properties are still largely unclear.

Here, we developed a conditional Gaussian mixture model (GMM) that specifies the relationship between the preceding REM sleep duration and the subsequent time spent in NREM sleep. We applied this model to systematically separate short cycles, which form sequences of REM sleep, from long cycles, which are characterized by the positive correlation between REM and subsequent NREM sleep. For long cycles, we defined the duration of the refractory period as a function of the preceding REM sleep duration. Next, we analyzed the EEG and other features of NREM sleep to identify defining properties of NREM sleep during the refractory period. We then used the cumulative distribution function (CDF) of the model as a measure for REM sleep propensity and found that entering REM sleep at higher propensity resulted, on average, in longer REM sleep episodes and reduced the EEG power. Finally, we employed the model to uncover changes in the regulation of REM sleep between the light and dark phase. Altogether, our model-based approach uncovered basic principles underlying the timing and duration of REM sleep episodes.

Results

A probabilistic model relating NREM sleep to preceding REM sleep

Consistent with the terminology introduced in earlier studies, we refer to the time interval between two successive REM periods as the inter-REM interval, while a sleep cycle (also called REM-NREM cycle) comprises one REM episode and the following inter-REM interval () [21,40]. We recorded spontaneous sleep in wildtype mice during the light phase and confirmed in our data set the positive correlation between the preceding REM sleep duration (REMpre) and the subsequent inter-REM interval (). For further analysis, we divided the inter-REM interval into its total duration of NREM sleep, |N|, and total duration of wakefulness, |W| (). Consistent with previous studies [19,20,25], we found that in mice, as in rats, REMpre was more strongly correlated with |N| than with either the total inter-REM interval or |W| ().

Correlation between REM sleep duration and inter-REM interval.

(A) Example EEG power spectrogram, EMG amplitude, and hypnogram with definitions of terms. REMpre, duration of the preceding REM sleep episode; inter-REM, duration of subsequent interval till next REM episode; |W|, sum of the durations of all wake episodes during the inter-REM interval; MA, microarousal (wake bouts ≤ 20 s); |N|, sum of the durations of all NREM episodes (including MAs) during the inter-REM interval. (B) Scatter plots with REMpre on the x-axis and subsequent inter-REM duration (left), |N| (middle), and |W| (right) on the y-axis. The results of linear regression fits are shown in red (P<0.00001 for all 3 slopes, n = 5098 sleep cycles from 72 mice recorded during the light phase).

To systematically describe the interaction between REM sleep and the duration and quality of subsequent NREM sleep, we developed a probabilistic model that specifies the relationship between REMpre and |N|. Applying the natural logarithm (ln) to |N| and plotting the distribution of ln(|N|) separately for increasing 30 s bins of REMpre, we found that for values of REMpre in the range from 30 s to 150 s, ln(|N|) forms a bimodal distribution, reflecting the presence of short and long sleep cycles in the hypnogram (). For REMpre ≥ 150 s, the distribution became unimodal. The distribution of ln(|N|) appeared to be a mixture of two bell-shaped distributions, with the relative weights of each distribution depending on REMpre. Based on this observation, we estimated for each 30 s bin of REMpre the distribution of ln(|N|) using a two component GMM (Methods). More precisely, for each REMpre bin, we modeled the distribution of ln(|N|) as the weighted sum of two Gaussians, kshort ᐧ N(μshort, σ2short) + klong ᐧ N(μlong, σ2long), where k, μ, and σ refer to the weighting factor, mean, and standard deviation of each Gaussian, respectively (). To test whether the GMM is a valid model for the data, we performed the Lilliefors-corrected Kolmogorov-Smirnov test for each 30 s bin and did not find sufficient evidence to reject the null hypothesis that the data are indeed drawn from the estimated Gaussian mixture distributions (). This was true regardless of the specific threshold used to score microarousals (MAs) ().

Conditional GMM to describe the relationship between REMpre and subsequent NREM.

(A) Scatter plot of REMpre vs. ln(|N|). Vertical dashed lines indicate consecutive 30 s bins of REMpre. Solid black lines represent the mean and standard deviation of ln(|N|) for each 30 s bin. (B) Histograms and probability density plots of ln(|N|) for consecutive REMpre bins as indicated on top. Probability densities were computed using a GMM. The notation [a, b) refers to the bin a ≤ REMpre < b. (C) Histogram of ln(|N|) for inter-REM intervals preceded by REM episodes in the range 30 s ≤ REMpre < 60 s. A GMM composed of two Gaussian distributions captures well the bimodal distribution of ln(|N|). The mean and standard deviation of the Gaussian for long and short cycles are referred to as μlong, σlong, and μshort, σshort, respectively. (D) Estimates of GMM parameters as a function of REMpre. The mixture parameter, klong, denotes the probability that a sleep cycle belongs to the long Gaussian distribution. For each parameter, we fitted a linear or logarithmic function describing its dependence on REMpre. (E) Heatmap in which each grid cell (x,y) represents the probability of transitioning from NREM to REM in between |N| - 25 s ≤ y ≤ |N| + 25 s following a REM episode of duration REMpre = x s for x in [10, 15, …, 250]. Each column of the heatmap sums up to 1. (F) Cumulative distribution function (CDF) of the GMM for 7 different values of REMpre. Each line represents, for the given REMpre value, the likelihood of entering the next REM period within |N| s of NREM sleep since the preceding REM episode.

We found that the weight of the long Gaussian distribution, klong, steadily increased with REMpre, indicating that the probability for long cycles was larger the longer the preceding REM episode (). Second, the mean of the long Gaussian distribution, μlong, increased with REMpre, while the standard deviation, σlong, decreased. The mean and standard deviation of the short cycles, μshort and σshort, both decreased with larger REMpre values.

Next, to describe the relationship between REMpre and each of the GMM parameters, we fit either linear or logarithmic functions to the parameter estimates (Figs ). By finding a function that captures the relationship between REMpre and each Gaussian mixture parameter, we were able to explain the amount of subsequent NREM sleep conditional on REMpre using a single probability model (see ).

We visualized the complete probability model using a heatmap (). Each grid cell (x,y) on the heatmap color-codes the probability to transition into REM sleep after |N| = y s of NREM sleep since the last REM sleep period of duration REMpre = x s. As expected from the distribution of ln(|N|) (), we observed two modes along |N|. The lower mode comprises only short inter-REM intervals and exists only for REMpre < 150 s. The second mode contains longer inter-REM intervals with larger |N| values. We refer to sleep cycles with |N| in the lower mode as sequential, as they result in a sequence (or cluster) of closely interspaced REM sleep periods. Sleep cycles that are part of the long distribution are termed single cycles. One characteristic of single cycles is that the mean values of |N| continuously increase with REMpre, i.e. longer REM sleep episodes are followed by larger amounts of NREM sleep before re-entering REM sleep. As further analyzed below, another key characteristic of single cycles is that the preceding REM period is followed by an interval of NREM sleep during which it is extremely unlikely to transition to REM sleep; the duration of this refractory period increases with REMpre.

Next, we computed the CDF of the conditional GMM for different values of REMpre (). Each CDF represents the probability that the animal transitions to REM sleep within |N| s of NREM sleep. In general, the longer the preceding REM period, the more NREM sleep is required to reach the same cumulative probability of re-entering REM sleep. But, irrespective of REMpre, within 2000 s of NREM sleep, the animal will most likely (> 99.46%) transition to the next REM period. For short REMpre values, the CDF immediately rises because of the comparably high probability for the occurrence of sequential cycles. In contrast, for long REM sleep episodes (REMpre ≥ 150 s), the CDF initially stays close to zero and only starts rising once it leaves the refractory period. Altogether, our statistical model to describe the relationship between REMpre and |N| suggests the existence of two different types of sleep cycles: Sequential cycles form sequences of REM sleep episodes, whereas single cycles are characterized by the existence of a refractory period, the duration of which increases with the duration of the preceding REM episode.

A previous study accounted for the presence of sequential sleep cycles and the resulting bimodality of the inter-REM distribution by using a mixture of two Poisson distributions [35]. Consistent with our findings, the authors reported that the mean duration of the inter-REM interval and the probability of observing a long inter-REM interval both increased with REMpre. However, a key property of the Poisson distribution is that the variance increases with the mean, which is not the case for our dataset for both the total inter-REM duration and |N| (). Our conditional GMM instead reproduces the reduction in the variance for increasing values of REMpre.

Sequential vs single cycles

Using our model, we defined a data-driven criterion to separate single from sequential cycles (). For each 2.5 s increment of REMpre, we calculated the probability density functions (PDFs) of both the short and long Gaussian distribution. The intersection point of the two distributions optimally separates single from sequential cycles (). In total, 19.3% of all cycles were sequential (). Since klong = 1 for REMpre ≥ 150 s, sequential cycles exist only for REMpre < 150 s. In contrast, single cycles exist for the entire range of REMpre, and consequently, their REMpre values are on average larger than those of sequential cycles (). The total duration of NREM sleep during sequential cycles reached values up to 222.5 s with a mean of 84.96 s ± 40.92 s (mean ± s.d.) (). In more than 50% of cases, REM sleep sequences comprised only two REM sleep episodes, forming one cycle; in rare cases, they included up to 5 consecutive cycles ().

Sequential vs single cycles.

(A) (Left) Scatter plot of REMpre vs. ln(|N|) with color-coded single and sequential cycles. The threshold optimally separating sequential from single cycles is shown in black. (Right) Illustration of how the threshold was determined for REMpre = 30 s. The probability density functions (PDFs) for the two distributions of the GMM (for REMpre = 30 s) are plotted along the y-axis. The red asterisk indicates the value of ln(|N|) at which the two Gaussians intersect. Values of ln(|N|) below the intersection point are more likely to be drawn from the short distribution and are consequently labeled as sequential cycles. Gray points correspond to cycles with REMpre < 7.5 s for which the conditional GMM is not defined (). (B) Pie chart indicating the percentage of single and sequential cycles. (C) Box plot comparing REMpre for single and sequential cycles. For sequential cycles, REMpre was shorter than for single cycles (Welch’s t-test, t = -35.13, p = 2.59e-228, nsequential = 947, nsingle = 3961). (D) Histogram of |N| for sequential cycles. The vertical dashed line indicates the mean (85.46 s ± 40.92 s; mean ± s.d.). (E) Bar plot showing the percentage of the number of cycles within a REM sleep sequence. Over half of REM sleep sequences contain only one cycle (i.e. comprise two REM periods). (F) Spectral density of parietal (nsequential = 947, nsingle, = 3961) and prefrontal (nsequential = 936, nsingle = 3919) EEG during REM and NREM sleep for both sequential and single cycles. Horizontal lines indicate frequencies at which the spectral density of sequential and single cycles are statistically different at various ɑ levels; (Welch’s t-test, * p<0.05; ** p<0.01; *** p<0.001). One recording did not contain a prefrontal EEG channel. Shadings, 99% confidence interval (CI).

To test for further differences between sequential and single cycles, we examined the EEGs for both types of cycles (). For both the parietal and prefrontal EEG, the general shape of the spectral density during REM sleep was similar, although there were differences in the overall power, which were particularly pronounced for the prefrontal EEG. Further analysis, however, showed that these differences were the result of the differences in the REM episode duration, REMpre, between single and sequential cycles (). We determined the spectral densities for different REM sleep durations () and used these to calculate weighted averages based on the distribution of REMpre for single and sequential cycles, respectively (see ). The weighted averages were very similar to the actual spectral densities, suggesting that the power differences in the REM sleep EEG between single and sequential cycles are the result of the differences in REMpre ().

In contrast to REM sleep, the spectral density for NREM sleep exhibited considerable differences between sequential and single cycles, particularly in the parietal EEG (). Compared with NREM sleep during single cycles, the NREM δ power (0.5–4.5 Hz) during sequential cycles was strongly reduced, while the θ power (5–9.5 Hz) was enhanced. We observed very similar changes in the δ and θ power when varying the threshold used to score MAs (). Thus, the EEG displayed two features resembling the EEG during REM sleep (reduced δ and increased θ power), suggesting that NREM sleep during sequential cycles may constitute an intermediate state between NREM and REM sleep [41].

Refractory vs. permissive periods during single cycles

During the inter-REM interval of single cycles, REM sleep is followed by a refractory period during which the probability of NREM to REM transitions is close to zero (). Using the conditional GMM, we formulated a data-driven definition of the refractory period (. For each 2.5 s increment of REMpre, we calculated the 1st percentile of the long Gaussian distribution and set it as the threshold separating the refractory period from the remaining permissive period (). Notably, the duration of the refractory period is approximately twice the duration of REMpre ().

Refractory and permissive periods during single cycles.

(A) (Left) Scatter plot of REMpre vs. ln(|N|) along with boundary (solid line) separating the refractory from the permissive period within single cycles. (Right) Illustration of how the threshold separating the refractory from the permissive period was determined for REMpre = 60 s. The CDF of the long Gaussian distribution is plotted along the y-axis. The value of |N| for which CDF(ln|N|) = 0.01 (indicated by the red asterisk) corresponds to the duration of the refractory period. (B) Scatter plot of REMpre vs. |N| along with the threshold separating the refractory from the permissive period. Of note, the refractory period is only defined for single cycles (red dots). The black line represents the threshold and the shading indicates the 99% confidence interval (CI) from 10,000 bootstrap iterations. (C) Spectral density of the prefrontal EEG for NREM sleep during the refractory and permissive period. The densities for both periods are statistically different for frequencies in the range 0–15 Hz (Welch’s t-test, *** p < 0.001, nrefractory = npermissive = 3892). Shadings, 99% CI. (D) Box plot comparing the rate of sleep spindles during the refractory and permissive period (Welch’s t-test, t = -36.96, p = 0.0, nrefractory = npermissive = 3908). The rate was calculated as the number of spindles per 1 min of NREM sleep. (E) Box plot comparing the rate of MAs during the refractory and permissive period (Welch’s t-test, t = 50.32, p = 0.0, nrefractory = npermissive = 3908). The rate was calculated as the number of MAs per 1 min of NREM sleep. (F) Progression of θ power, σ power, spindle rate, and MA rate throughout the refractory and permissive period for different values of REMpre. The refractory period is defined as outlined in A. The permissive period comprises the time from the end of the refractory period to the onset of the next REM period. The durations of both the refractory and permissive period were normalized to unit length and subdivided into quartiles of equal normalized duration. The average for all REMpre values is shown in black. Shadings, 99% CI. (G) Progression of θ power (Row 1) and σ power (Row 2) on non-normalized time scale during the first 600 s of NREM sleep during the inter-REM interval for different values of REMpre. The two vertical dashed lines indicate the lowest and highest threshold separating the refractory from the permissive period corresponding to the low and high bound of REMpre. Shadings, 99% CI.

Next, we tested whether the EEG and other properties of NREM sleep, such as sleep spindles and microarousals (MAs), differ between the refractory period and the following permissive period that ranges from the end of the refractory period till the onset of the next REM episode. We computed the spectral density of the EEG during NREM sleep and found that the δ, θ, and σ (10–15 Hz) power of both the prefrontal and parietal EEG were lower during the refractory period (Figs ). In addition to these differences in the EEG, the spindle rate was reduced during the refractory period while MAs were more frequent (). Sleep spindles have recently been implicated in promoting REM sleep and their overall increased rate during the permissive period may thus facilitate transitions to REM sleep [42].

We then analyzed the time course of the different prefrontal EEG power bands and the rate of MAs and spindles throughout the sleep cycle by normalizing the durations of both the refractory and permissive period and dividing them into quarters. The θ power, σ power, and rate of sleep spindles were strongly reduced after REM sleep, increased with downward concavity throughout the refractory period, reached a plateau near its end, and then continued increasing with upward concavity throughout the permissive period (). The overall time course of the θ and σ power did not depend on the exact value of the threshold used to score MAs (). The increase of these values in the final quarter of the permissive period reflects the transition stage between NREM and REM sleep, which is characterized by increased spindle activity, increased θ and reduced δ power () [41]. Importantly, the normalized time course of the θ power, σ power, and spindle rate was consistent regardless of REMpre (), suggesting that, as the duration of the refractory period increased with longer REM periods, the rate at which the θ power, σ power, and frequency of spindles increased was proportionally reduced. Plotting these quantities throughout the refractory period along a non-normalized time axis, we indeed found a slower rise in their time courses following longer REM periods (Figs ). The rate of MAs followed the opposite pattern. It was strongly increased after REM sleep and decreased with upward concavity throughout the refractory period before reaching a plateau, and then continued to decay with downward concavity throughout the permissive period (). We found a similar pattern, irrespective of the used MA threshold (). The time course at which the MA rate declined also depended on REMpre: The longer REMpre, the less steep the decay (). The fact that the inflection point of the time courses of the θ power, σ power, spindle and MA rate all occur at the threshold suggests that our data-driven definition of the refractory and permissive period reflects a natural separation within single sleep cycles.

In contrast, the time course of the δ power was not aligned with the threshold between the refractory and permissive period (). Although its general time course, with an increase at the beginning of the inter-REM interval followed by a decrease, was consistent for all ranges of REMpre, the normalized time bin at which the δ power started decaying varied with REMpre and was not systematically related with the threshold between the refractory and permissive period. Thus, although the time course of the δ power and, in particular, its overall power was strongly influenced by REMpre, consistent with a previous study [14], it did not reflect the probability of NREM to REM sleep transitions.

Finally, to test whether the observed changes in the EEG power and sleep microarchitecture during the refractory period were specifically due to REM sleep, we compared NREM sleep after REM sleep with NREM sleep following wake periods of equal duration (Figs ). In general, the power of all tested EEG bands was less strongly modulated after a period of wakefulness. After REM sleep, the θ and σ power showed a stronger reduction and a steeper increase () compared to their time courses following wake periods of equal duration (). The rate of MAs was also more strongly modulated following bouts of REM sleep (). Of all tested variables, the θ power and rate of MAs were the least affected by preceding wakefulness (). Thus, the observed changes in the NREM EEG power and sleep microarchitecture, especially in the θ power and MA rate, were particularly associated with preceding REM sleep. These findings suggest that REM sleep has different consequences on the subsequent quality of NREM sleep than does wakefulness.

Relationship between wakefulness and NREM sleep

In addition to REMpre, |N| may also be modulated by wake periods in the inter-REM interval. In our data set, only 12.9% of sequential cycles contained wake periods (). 43.0% of single cycles were not interrupted by wakefulness () and 42.7% contained only one or two wake episodes (). In general, |N| was larger for single cycles with larger |W| () and also for cycles with a larger number of wake episodes (). We observed similar relationships, irrespective of the threshold used to score MAs (). By comparing |N| for cycles with different |W| while keeping the number of wake episodes constant, and also comparing |N| for cycles with different numbers of wake episodes while keeping |W| constant, we confirmed that both |W| and the number of wake episodes within an inter-REM interval contribute to longer |N| ().

Relationship between wake episodes and NREM sleep during an inter-REM interval.

(A) (Left) Pie chart indicating the percentage of sequential cycles without (|W| = 0) and with wake episodes (|W| > 0). (Right) Box plot comparing |N| for sequential cycles with and without wake episodes (t-test, t = -3.22, p = 0.0012, n|W|>0 = 122, n|W| = 0 = 825 cycles). (B) (Left) Pie chart indicating the percentage of single cycles without and with wake episodes. (Right) Box plot comparing |N| for single cycles with and without wake episodes (Welch’s t-test, t = -28.64, p = 3.97e-163, n|W|>0 = 2259, n|W| = 0 = 1702). (C) Bar plot showing the distribution of the number of wake episodes during the inter-REM interval of single cycles. Note that 99.63% of single cycles had 8 or fewer wake episodes. (D) (Left) Box plot comparing total NREM duration, |N|, for single cycles with increasing values of total wake duration, |W| (Welch’s ANOVA, F(5,1233.21) = 301.99, p = 3.78e-211). The x-tick q0 corresponds to cycles without wake. The remaining cycles with |W| > 0 were subdivided into quintiles, labeled q1—q5, based on the distribution of |W| for single cycles. (Right) Box plot comparing |N| for single cycles based on the number of wake episodes occurring during the inter-REM interval (Welch’s ANOVA, F(5,475.26) = 246.53, p = 1.65e-129). Note that 97.61% of single cycles contained 5 or fewer wake episodes. (E) Progression of θ power, σ power, spindle rate, and MA rate during NREM sleep before and after a wake episode. Only sequences with at least 1 minute of NREM sleep both before and after wake during the inter-REM interval of single cycles were included. The duration of NREM episodes was normalized. ‘Before’ refers to all NREM sleep in between either the previous REM or wake episode and the current wake episode. ‘After’ refers to all NREM sleep in between the current wake episode and either the next wake or REM episode. Shadings, 99% CI. (F) Bar plot showing average drop (or increase) in θ power, σ power, spindle rate, and MA rate over wake episodes with different durations. All wake episodes for single cycles were divided into five quintiles based on the distribution of their durations. A drop or increase in each variable was calculated by subtracting the average value for 1 min of NREM after wake from the average value for 1 min of NREM before wake (θ: Welch’s ANOVA, F(4,899.90) = 10.48, p = 2.65e-08; σ: ANOVA, F(4,1901) = 1.99, p = 0.093; Spindles: ANOVA, F(4,1901) = 6.22, p = 5.70e-05; MAs: Welch’s ANOVA, F(4,895.63) = 9.13, p = 3.08e-07). Error bars, 95% CI from 1,000 bootstrap iterations.

To test how wake episodes affect the EEG power during NREM sleep, we computed the θ and σ power before and after wake episodes during the inter-REM interval of single cycles. The power in both frequency bands increased during the preceding NREM episode, dropped to levels lower than those at wake onset, and then started rising again throughout the following NREM episode (; see for different MA thresholds). We observed the same pattern for the rate of spindles, whereas the rate of MAs was enhanced after wake and then decayed throughout NREM sleep. The drop in the spindle rate and the increase in the MA rate was larger the longer the intervening wake episode (Figs ). The drop in θ power, on the other hand, first increased and then declined with the duration of the intervening wake episodes. Thus, each wake episode leads to a reduction in quantities that reflect the probability of NREM to REM transitions (), possibly resulting in more NREM sleep, which could explain the positive correlation of |N| with the total duration and number of wake episodes.

Correlation between model CDF and REM sleep duration

Next, we aimed to determine factors influencing the duration of REM sleep. First, we tested whether the REM duration (REMpost) is influenced by the preceding REM duration (REMpre). We found that REMpost is negatively correlated with REMpre (slope = -0.085, R2 = 0.0070, p = 7.17e-09). The slope of the negative correlation was larger for sequential cycles (slope = -0.22, R2 = 0.010, p = 0.0020) than for single cycles (slope = -0.09, R2 = 0.0089, p = 5.07e-09) (). As single cycles are interrupted by longer inter-REM intervals, the intervening periods of NREM sleep and wakefulness may themselves influence REMpost, thereby weakening the correlation between the two variables. Furthermore, the negative correlation supports the notion of a short-term hourglass process: As more REM propensity is discharged during a longer REM episode, the subsequent episode tends to be shorter.

Effects of sleep history on REM episode duration.

(A) (Left) Box plot comparing the duration of REM sleep (REMpost) following sequential cycles for different values of REMpre. Red line, linear regression (slope = -0.22, R2 = 0.010, p = 0.0020). (Right) Box plot comparing REMpost following single cycles for different values of REMpre. Red line, linear regression (slope = -0.091, R2 = 0.0089, p = 5.07e-09). (B) Box plot comparing REMpost following single cycles with different values of |N|. Red line, linear regression (slope = 0.0041, R2 = 9.83e-04, p = 0.052). (C) Box plots comparing REMpost dependent on the CDF value at REM onset for single cycles. The left plot includes all REMpre values; the remaining plots show the correlation for increasing ranges of REMpre. Dashed lines, linear regression (All: slope = 22.82, R2 = 0.013, p = 1.01e-12; [7.5,60): slope = 8.58, R2 = 0.0013, p = 0.090; [60,120): slope = 29.42, R2 = 0.024, p = 1.57e-07; [120,180): slope = 29.64, R2 = 0.039, p = 3.12e-06; [180,240): slope = 43.28, R2 = 0.077, p = 0.0060). (D) Spectral density of prefrontal EEG during REM episodes following single cycles as a function of the CDF at REM onset (Welch’s ANOVA, δ: F(4,1078.85) = 16.18, p = 7.04e-13; θ: F(4,1084.95) = 6.06, p = 8.0e-05; σ: F(4,1089.74) = 9.46, p = 1.61e-07).

Previous studies found that the duration of a REM episode also depends on the amount of NREM sleep preceding REM sleep. But this correlation is typically only weak [19,21,43,44] or was reported to be non-existent [20]. In our dataset, we found no significant correlation between the preceding inter-REM interval duration of single cycles (; slope = 0.0012, R2 = 4.92e-04, p = 0.17) and REMpost, or between the amount of NREM, |N|, and REMpost (; slope = 0.0041, R2 = 9.83e-04, p = 0.052).

The CDF of the conditional GMM describes the probability of entering REM sleep within |N| (s) of NREM sleep since the last REM sleep episode, and it can therefore be interpreted as a measure of the ultradian propensity for REM sleep throughout a single sleep cycle. For a given value of REMpre, there is considerable variation in the values of |N| during single cycles and we tested whether the resulting differences in the CDF values at REM onset may influence the subsequent REM sleep duration. We indeed found that REMpost and the value of the CDF at REM onset were positively correlated (; slope = 22.82, R2 = 0.013, p = 1.01e-12), suggesting that a higher propensity at REM onset leads to longer REM sleep episodes. This correlation was particularly pronounced for cycles with REMpre ≥ 60 s. The finding that the REM duration is more closely correlated with the CDF than with the preceding NREM duration (Williams’ correlation test, t = 17.59, p = 0.0) suggests that it is not the amount of NREM sleep, but rather the propensity for REM sleep accumulated throughout the sleep cycle that influences the REM duration.

Next, to test whether the presence of wake periods affects the correlation between the CDF and REM duration, we calculated the correlation between the CDF and REMpost separately for sleep cycles with and without wake episodes (). Importantly, the dependence of REMpost on the CDF value was not affected by the presence of wake episodes. This finding suggests that the impact of the REM propensity, as quantified by the CDF, on REMpost is not influenced by the presence of wake episodes.

Finally, we analyzed whether the CDF at REM onset also affects the EEG power during REM sleep. Changes in the EEG power were particularly pronounced for the prefrontal EEG (; see for parietal EEG). We found that if the animal transitioned to REM sleep at a low CDF value, the power of the prefrontal EEG in the δ, θ, and σ range was higher than when REM sleep was entered at a high CDF value (). To test whether the negative correlation between the CDF and EEG power may be explained by differences in the REM sleep duration, we calculated linear approximations of the spectral densities for each CDF range based on the distribution of REM durations within this range (). For high values of the CDF (CDF ≥ 0.4), the linear approximation closely matched the observed densities. However, for low CDF values (CDF < 0.4) the approximation considerably differed from the observed densities, suggesting that, at least in part, the magnitude of the prefrontal EEG power during REM sleep reflects REM sleep propensity.

Changes in REM sleep regulation between light and dark phase

Next, we applied our model to test for differences in the regulation of REM sleep between the light and dark phase. As expected [45,46], mice spent less time in sleep during the dark phase (). In addition, the ratio of REM sleep to the total amount of sleep was reduced (). To further understand why the relative amount of REM sleep is reduced during the dark phase, we fit the conditional GMM to the recordings from the dark phase (Figs ). As was the case for the light phase, the Lilliefors-corrected KS test did not provide enough evidence to reject the null hypothesis that |N| conditional on REMpre follows a lognormal Gaussian mixture distribution ().

Changes in REM sleep regulation between light and dark phase.

(A) Bar plots comparing the percentage of REM, NREM sleep, and wake during the light and dark phase (REM: Welch’s t-test, t = 16.14, p = 3.37e-30; NREM: t-test, t = 20.26, p = 4.53e-38; Wake: t-test, t = -20.79, p = 5.19e-39; nlight = 72, ndark = 35 mice). Error bars, 95% CIs from 1,000 bootstrap iterations. (B) Bar plot comparing the ratio REM/(REM+NREM) for the light and dark phase (t-test, t = 3.17, p = 0.0019, nlight = 72, ndark = 35 mice). Error bars, 95% CIs from 1,000 bootstrap iterations. (C) Pie chart showing the percentage of sequential and single cycles during the dark phase. (D) Comparison of GMM parameters for the light and dark phase (Welch’s t-test with Bootstrap, klong: t = -77.26, p = 0.0; μlong: t = -372.77, p = 0.0; σlong: t = -57.93, p = 0.0; Methods). (E) Comparison of μlong (solid lines), and threshold (Ref., dashed lines) separating the refractory from the permissive period for the light and dark phase. Shadings, 95% CIs obtained from 10,000 bootstrap iterations. (F) Comparison of the CDFs of the GMMs for the light and dark phase shown for 4 different values of REMpre.

Based on our model, we identified two major factors underlying the reduced REM/total sleep ratio during the dark phase. First, during the dark phase, the average of klong, the weighting factor for the long Gaussian, was significantly increased (), reflecting the increased percentage of single cycles and reduction in sequential cycles ( compared to the light phase (). As single cycles contain less REM sleep relative to NREM sleep, their reduced frequency during the dark phase, mirrored in the increase of klong, contributes to the decrease in REM sleep. Second, the mean of the long Gaussian distribution, μlong, was larger in the dark phase, indicating a delayed increase in the likelihood of NREM to REM sleep transitions, a second factor contributing to the increase in NREM relative to REM sleep (). Interestingly, the duration of the refractory period did not significantly differ between the light and dark phase () because σlong was also increased during the dark phase (), counteracting the effects of μlong on the duration of the refractory period. Consequently, for large values of REMpre (≥ 150 s), the CDFs for both the light and dark phase stayed close to zero for approximately the same amount of time but then increased at different rates (). Hence, our model suggests that the probability of sequential cycles is decreased during the dark phase and that the rate at which REM propensity increases throughout the sleep cycle is reduced during the dark phase.

Discussion

In this study, we characterized the relationship between REMpre and subsequent NREM sleep using a conditional GMM. The distribution of NREM sleep during sleep cycles is bimodal, suggesting that two different types of inter-REM intervals exist (Figs ). Using our model, we separated short from long inter-REM intervals and found that NREM sleep during short intervals displays reduced δ and increased θ power (). Longer inter-REM intervals begin with a refractory period, during which transitions to REM sleep are highly unlikely (). The refractory period proportionally increases in duration with REMpre and is characterized by a low θ power, σ power, and spindle rate as well as an increased frequency of MAs. The total duration of NREM sleep also depends on the number and total duration of wake periods during the inter-REM interval (). REMpre is negatively correlated with the subsequent REM sleep duration. In addition, the CDF of the conditional GMM is positively correlated with the next REM duration, suggesting that a higher propensity for REM sleep results in longer REM episodes (). The build-up of the REM propensity is delayed during the dark phase (). Altogether, our analysis suggests that three major factors shape the ultradian regulation of REM sleep: The presence of two distinct types of sleep cycles, a refractory period suppressing transitions to REM sleep, and a propensity for REM sleep that influences the next episode duration.

Sequential sleep cycles

Sequential sleep cycles have been observed in rats [32,33,35], cats [21], monkeys [37] and humans [36,47]. The presence of both short and long cycles results in a bimodal distribution of inter-REM intervals in these species. The minimum between the two modes for short and long inter-REM intervals served as a threshold to separate sequential from single cycles. Comparative analysis showed that this threshold and the average duration of inter-REM intervals increase with brain size [1,34]. In rats and, as shown here, in mice the durations of REM sleep episodes at the beginning of sequential cycles are on average shorter than those at the start of single cycles, and the frequency of sequential cycles is reduced during the dark period [33]. Thus, sequential cycles in both species share a lot of statistical similarities suggesting a common physiological mechanism in both species. Infusion of a serotonin receptor agonist into the laterodorsal tegmentum reduced the frequency of REM sleep sequences, while an antagonist increased their occurrence in rats [48], suggesting a role of serotonin in their regulation. We speculate that an increased release of serotonin may contribute to the suppression of sequential sleep cycles during the dark phase [49,50].

Similar to our findings in mice, the EEG δ and σ power in rats is also reduced during sequential cycles [51]. A previous study in cats described NREM sleep during sequential cycles as light slow wave sleep as it is also characterized by low δ power [21]. NREM sleep during sequential cycles is thus reminiscent of the so-called intermediate stage preceding a transition to REM sleep, which shares both features of NREM and REM sleep [41,52]. For future studies, it would be interesting to perform simultaneous local field potential recordings from multiple sites to test whether both NREM and REM sleep states may locally co-exist in different brain areas during sequential cycles as observed for the intermediate stage [41,52].

Current dynamical systems models generate the ultradian alternation between NREM and REM sleep either through feedback from the arousal system to the circadian or homeostatic drive for sleep [53] or by implementing mutually inhibitory interactions between REM-on and REM-off neurons [54,55]. In the latter model type, a homeostatic hourglass process dictates the timing of REM sleep, an assumption supported by the positive correlation between REMpre and the following inter-REM interval. Testing under which assumptions these models can reproduce the lognormal distribution of NREM sleep may further constrain the time course at which REM propensity accumulates. For future studies, it would be interesting to investigate whether these dynamical models can also explain the generation of sequential cycles, possibly by introducing physiologically motivated noise terms mimicking fluctuations in firing rates or neurotransmitter concentrations [54,56].

REM sleep is followed by a refractory period

As the conditional GMM allowed us to separate single from sequential sleep cycles, we could disentangle the refractory period, which only exists for single cycles, from sequential cycles and define it for the whole range of REMpre. Since the refractory period proportionally increases with the duration of the preceding REM sleep episode (by about 2 · REMpre), it may mechanistically explain the positive correlation between REMpre and succeeding NREM sleep [29-31]. With its low σ power, spindle rate, and high rate of MAs, the refractory period likely constitutes a fragile state of NREM sleep in mice. A study in cats similarly found that REM sleep is followed by a stage of light sleep [21]. In humans, the sleep phase N1, which is characterized by the absence of spindles, is most likely to occur during the sleep cycle after REM sleep [57] (published in preprint). Hence, the refractory period may be classified as a substage of NREM sleep in mice that resembles stage N1 in humans.

At present, the refractory period in our study is statistically defined as an interval during which REM sleep is unlikely to occur. The physiological mechanisms underlying this period, however, are unknown. To further characterize this substage of NREM sleep, it would be interesting to test whether activation of known REM sleep-promoting neurons specifically during the refractory period is indeed ineffective in inducing REM sleep. Electrophysiological recordings in the ventrolateral periaqueductal gray (vlPAG) indicated a role of GABAergic REM-off neurons in this area in the ultradian regulation of REM sleep. The activity of vlPAG REM-off neurons gradually decreases during the inter-REM interval, and abruptly rises at the end of REM sleep in a duration-dependent manner: The longer the REM episode lasts, the more these neurons become subsequently activated [13]; an effect which may mediate the dependence of the refractory period on the preceding REM duration [29-31]. Interestingly, vlPAG GABAergic neurons express orexin/hypocretin receptors [58] and may thus be excited by the wake-active orexin/hypocretin neurons known to project to the vlPAG [17]. Such an excitatory drive may underlie the delayed increase in REM propensity during the dark phase suggested by our model: A higher baseline activity of REM-off neurons during the dark phase, when orexin/hypocretin levels are high [59], may delay the time till the activity of these neurons is low enough to allow for a transition to REM sleep.

Our analysis suggests a close association of θ, σ power, and sleep spindles with the probability of NREM to REM transitions. The spindle rate is reduced during the refractory period, when transitions to REM sleep are highly unlikely, while the increased rate of spindles during the permissive period may facilitate REM sleep. In support of this view, a recent study showed that optogenetically triggering spindles through stimulation of the thalamic reticular nucleus enhances the chance of transitions to REM sleep [42]. For future studies, it is therefore important to understand how thalamocortical circuits generating spindles interact with hypothalamic or brainstem circuits controlling REM sleep.

Interaction between REM propensity and REM episode duration

As the REM sleep duration is negatively correlated with the occurrence of sequential cycles and positively correlated with the duration of the refractory period, it plays a crucial role in temporally shaping the sleep cycle. Previous studies showed that during the homeostatic rebound following REM sleep deprivation, REM sleep episodes are elongated [28,60,61], suggesting that an increased propensity for REM sleep is associated with longer episodes. We similarly found that during spontaneous sleep a higher REM propensity, as reflected in the CDF of the conditional GMM, results in longer REM episodes. Our finding suggests that the REM propensity builds up fast enough during the sleep cycle to influence the duration of REM sleep episodes. Consistent with this, short-term REM sleep deprivation (as short as 20 min) induces a detectable rebound in REM sleep, further supporting the notion that REM sleep propensity accumulates at a time scale relevant to modulate the ultradian timing of REM sleep [60-62].

In addition to the episode duration, the CDF was also correlated with the EEG power during REM sleep. Interestingly, sleep recordings in humans also showed a reduction in the REM EEG power within a similar range (α power), after enhancing the REM propensity through REM sleep deprivation [63-65]. A recent study in rats found that the θ power declines throughout REM sleep and this decline is positively correlated with the amount of NREM sleep during the preceding inter-REM interval, also supporting a close association between REM propensity and changes in the EEG power [66].

For future studies, it would be interesting to test how the propensity for REM sleep reflected in the CDF of the conditional GMM relates to the homeostatic long-term process mediating the rebound in REM sleep following long-term total sleep or REM sleep deprivation. The long-term process regulates the daily quota of REM sleep and is thought to accumulate in the absence of REM sleep during both NREM sleep and wakefulness [29,67,68], and has been proposed to be separate from the short-term process regulating the ultradian timing of REM sleep [28]. In vivo recordings of vlPAG GABAergic neurons showed that the firing rates of these REM-off neurons during inter-REM are modulated by the preceding REM sleep duration, suggesting that they mirror the short-term REM propensity [13]. Recording the same neurons throughout long-term REM sleep deprivation and recovery sleep may provide important insights into the extent to which the short- and long-term processes are disjunct or overlap at the neuronal level.

Role of Wakefulness

In addition to the preceding REM sleep, the total amount of NREM sleep during the sleep cycle is also modulated by wake periods. Similar to REM sleep, although to a lesser degree, wake episodes lead to a decrease in the θ power, σ power, and spindle rate and an increase in the MA rate, quantities which are indicative of the likelihood of NREM to REM sleep transitions. Wake episodes may result in a reduction in the probability of REM transitions and thus lead to more NREM sleep, which in turn could explain the positive correlation between the total duration or number of wake episodes and the total NREM duration. Previous studies also suggest that an increased δ power induced by long intervals of wakefulness enforced by sleep deprivation suppresses the initiation of REM sleep [29,69]. On the ultradian time scale, however, we observed no systematic relation between the time course of the δ power and the likelihood of NREM to REM transitions.

Our finding that the presence or absence of wakefulness during the sleep cycle does not affect the correlation between the CDF value and subsequent REM duration (), indicates that, at least on the ultradian timescale, REM propensity and its impact on the REM episode duration is more closely associated with the time spent in NREM sleep than with the combined time in NREM sleep and wake. Thus, while wake periods reduce the following θ power, σ power, and rate of sleep spindles and may thereby lower the subsequent opportunity for NREM to REM transitions [32], they do not appear to interfere with the build-up of the REM propensity reflected in the CDF, reinforcing the notion that the ultradian REM propensity primarily accumulates during NREM sleep [20]. Altogether, our model-based approach provides a flexible framework to study the key factors underlying the ultradian timing of REM sleep and will inform future experimental studies to understand the mechanisms regulating the REM sleep duration, refractory period, and induction of REM sleep sequences.

Methods

Experimental setup

Ethics statement

All experimental procedures were approved by the Institutional Animal Care and Use Committee (IACUC) at the University of Pennsylvania and conducted in accordance with the National Institutes of Health Office of Laboratory Animal Welfare Policy.

Animals

Experiments were performed in male or female C57BL/6J mice (Jackson Laboratory stock no. 000664). Animals were housed on a 12-h dark/12-h light cycle (lights on between 7 am and 7 pm) and were aged 6–12 weeks at the time of surgery. All mice were group-housed with ad libitum access to food and water.

Surgery

All surgeries were performed following the IACUC guidelines for rodent survival surgery. Prior to surgery, mice were given meloxicam subcutaneously (5 mg/kg). Mice were anesthetized using isoflurane (1–4%) and positioned in a stereotaxic frame. Animals were placed on a heating pad to maintain the body temperature throughout the procedure. Following asepsis, the skin was incised to gain access to the skull. For EEG recordings, stainless steel wires were attached to two screws, one on top of the parietal and one on top of the prefrontal cortex. The reference screw was inserted on top of the cerebellum. For EMG recordings, two stainless steel wires were inserted into the neck muscles. All electrodes, screws, and the mini-connector holding the EEG, EMG wires were secured to the skull using dental cement. After the injection and implantation were finished, bupivacaine (2 mg/kg) was administered at the incision site.

Sleep recordings

Sleep recordings were performed in the animal’s home cage or in a cage to which the mouse was habituated for 3 days, which was placed within a sound-attenuating box. EEG and EMG signals were recorded using an RHD2000 amplifier (intan, sampling rate 1 kHz). EEG and EMG signals were referenced to the common ground screw on top of the cerebellum. During the recordings, EEG and EMG electrodes were connected to flexible recording cables using a small connector. To determine the brain state of the animal, we first computed the EEG and EMG spectrogram with consecutive fast Fourier transforms (FFTs) calculated for sliding, half-overlapping 5 s windows, resulting in a 2.5 s time resolution of the hypnogram. Next, we computed the time-dependent delta, theta, sigma, and high gamma power by integrating frequencies in the range 0.5 to 4 Hz, 5 to 12 Hz, 12 to 20 Hz, and 100 to 150 Hz, respectively. We also calculated the ratio of the theta and delta power (theta/delta) and the EMG power in the range 50 to 500 Hz. For each power band, we used its temporal mean to separate it into a low and high part (except for the EMG and theta/delta ratio, where we used the mean plus one standard deviation as threshold). REM sleep was defined by high theta/delta ratio, low EMG, and low delta power. A state was set as NREM sleep, if delta power was high, the theta/delta ratio was low and EMG power was low. In addition, states with low EMG power, low delta, but high sigma power were scored as NREM sleep. Wake encompassed states with low delta power and high EMG power and each state with high gamma power (if not otherwise classified as REMs). Of the bins classified as wake periods, those forming sequences of 20 s or less were classified as microarousals. Finally, we manually rescored the automatic classification to correct for errors using a graphical user interface, visualizing the raw EEG, EMG signals, spectrograms, and hypnogram. The software for automatic brain state classification and manual scoring was programmed in Python. The light-phase data contained 5098 sleep cycles from 125 recordings of 72 mice. The dark-phase data contained 1242 sleep cycles from 55 recordings of 35 mice.

Data analysis

Gaussian mixture model parameters

Following the definition in [40], a sleep cycle comprises a REM sleep episode and the following inter-REM interval. All sleep cycles with REM sleep episode duration REM < 240 s were divided into 8 groups based on REM. We chose 30 s non-overlapping bins ([0,30), [30,60), …, [210,240)) to ensure that each group contained enough cycles to reliably estimate the model parameters while being able to capture the change in these parameters conditional on REM. For each group, we used the Expectation-Maximization algorithm on ln(|N|) to find the maximum likelihood estimates for klong, kshort, μlong, μshort, σlong, and σshort. In the case of REM ≥ 150 s, the distribution of ln(|N|) was unimodal; consequently, klong = 1 and we only estimated μlong and σlong. For consistency, we assumed that the apparent unimodality of the distribution for REM < 30 s results from the blending of the distributions for short and long cycles. Computations were performed using the python package scikit-learn [70].

Lilliefors-corrected KS test

The standard Kolmogorov-Smirnov (KS) test compares data to a pre-defined distribution by comparing the CDF to the empirical cumulative distribution function (ECDF). However, it is invalid when the parameters for the CDF are estimated using the data [71]. Therefore, to perform a valid goodness-of-fit test with the estimated parameters, we applied the Lilliefors-corrected KS test as follows:

For each consecutive 30 s bin of REMpre, using the estimated parameters, we draw a random sample with the same size as the original data and calculate the KS-statistic.

We repeat (1) 10,000 times and form a distribution of KS-statistics, KSsim.

We calculate the KS-statistic using the observed data (KSobs) and find how extreme KSobs is when compared to the distribution of the simulated KS-statistics.

If KSobs is too large (KSobs > 95% of KSsim), we reject the null hypothesis that our data comes from the specified distribution.

Performing the Lilliefors-corrected KS test, we found that there is not enough evidence to reject the hypothesis that our data is drawn from a Gaussian mixture distribution (). This was the case for both the light and dark phase and held true for different thresholds to score MAs ().

Conditional Gaussian mixture model

After estimating all 6 parameters of the GMM for each 30 s bin of REMpre, we used either a linear (y = ax+b) or logarithmic (y = a·ln(x+b)+c; b≥0) function to model the relationship between each parameter and REMpre (Figs ). kshort was implicitly defined as kshort = 1-klong. Because klong = 1 for REMpre > 150 s, we only used the first 6 values of klong for fitting the function. We calculated the residual sum of squares (RSS) for both the linear and logarithmic fits and chose the function for which the RSS was lower. The fitting was performed using the Trust Region Reflective algorithm implemented in SciPy [72]. The logarithmic fit was generally better for all parameters except for σshort. In total, the model comprises 14 parameters. The complete probability model can be expressed as follows: where f is the probability density function of a Gaussian distribution, and

Note that k(x) is defined as the minimum of the result of the log function and 1, as k(x) is a probability and cannot exceed 1.

For the dark phase, the data for 150 ≤ REM < 180 contained only one single ln(|N|) value falling into the short Gaussian distribution, resulting in a value for k below 1. To avoid that a single data point changes the entire relationship between REM and k, we fit the logarithmic and linear functions on the first 4 values of k (Figs ).

Of note, the complete conditional GMM for the light phase is defined only for sleep cycles with REM ≥ 7.5 s, because REM = 7.5 s is the lowest value of REM for which the intersection point of the short and long Gaussian distributions, x, satisfies μshort < xintersect < μlong (). Similarly, the complete model for the dark phase is only defined for sleep cycles with REM ≥ 12.5 s ().

Model simulation

To assess the goodness-of-fit of the model to the data, we performed a simulation as follows:

For each REM in the data set with 7.5 s ≤ REM < 240 s:

Calculate klong, kshort, μlong, μshort, σlong, σshort using the conditional GMM.

Based on these parameters, define the Gaussian mixture distribution for the given REM and sample one data point (ln(|N|)) from this distribution.

Perform step (1) 10,000 times.

Compare the simulated distribution of ln(|N|) to the actual distribution of ln(|N|) for sleep cycles with REM in the range 7.5 s ≤ REM < 240 s.

We found that the simulated distribution of ln(|N|) closely overlaps with the distribution of the actual data ().

Single and sequential cycles

For each sleep cycle with REM in the range 7.5 s ≤ REM < 240 s, we label the cycle as either single or sequential using the following procedure:

Using the conditional GMM, determine for REM the probability distributions for the short and long Gaussian distributions.

Find the intersection point, x, along the x-axis (ln(|N|) axis) of the PDFs of the two Gaussian distributions.

If ln(|N|) < x, the cycle is sequential, otherwise it is labeled as single.

We use the criterion in step (3) because if ln(|N|) < x, the PDF of the short Gaussian evaluated at ln(|N|) is larger than the PDF of the long Gaussian evaluated at ln(|N|) and vice versa. i.e. where f and f correspond to the PDFs of the short and long Gaussian distributions.

Refractory and permissive period

For each value of REM in the range 7.5 s ≤ REM < 240 s, we determine the threshold separating the refractory from the permissive period as such:

Using REM and the conditional GMM, calculate μ and σ.

Using these parameters, define the CDF of the long Gaussian distribution, F(x).

The ln(|N|) value for which F(ln(|N|)) = 0.01 is set as the threshold.

Sleep spindle detection

For spindle detection, we first calculated the spectrogram for the prefrontal EEG. The spectrogram was computed for consecutive 600 ms windows with 500 ms overlap, resulting in a 100 ms temporal resolution. The spindle detection algorithm used two criteria to determine for each 100 ms time bin whether it was part of a spindle or not: The first criterion was that the height of the maximum peak in the sigma frequency range (10–16.67 Hz) exceeds a threshold, which corresponded to the 96th percentile of all maximum peaks in the sigma frequency range of the sleep recording. We determined the optimal percentile value by maximizing the performance of the algorithm on a manually annotated control data set. Second, the power value of the peak in the sigma range (10–16.67 Hz) had to be greater than half of the peak value in the range 0–10 Hz. The optimal value for this ratio (sigma peak ratio) was again determined on the control data set. Next, the algorithm merged spindle events that were temporally close to each other. First, spindle events in adjacent bins were considered as part of the same spindle. Second, we fused together sequences of spindle events that were interrupted by gaps of less than 300 ms. The optimal value for the gap was again determined on the control data set. Finally, we discarded spindles with duration ≤ 200 ms. Of all the potential spindles, we only considered those as spindles where for at least half of the time bins the peak frequency lied in the range of 10–16.7 Hz. The parameters of the spindle detection algorithm (sigma percentile threshold, sigma peak ratio, and minimum fusing distance) were optimized using a manually annotated data set. The algorithm correctly identified 88.8% of spindles present in the annotated data set with a false positive rate of 6.9%.

Linear approximation of spectral densities

To test if differences in the EEG spectral density between single and sequential REM sleep episodes were due to differences in their respective durations (), we first discretized the distribution of REM for both types of cycles into eight 30 s bins and calculated for each bin the corresponding REM sleep duration-dependent spectral density, P. Then, for single and sequential cycles, we found the proportion of cycles, w and w, with REM falling into bin i. Note,

Using these proportions as weights, we calculated a weighted average of the duration dependent spectral densities to give us linear approximations, and of the true spectral densities as follows:

Similarly, we tested whether differences in the EEG during REM sleep were due to differences in the CDF at REM onset or the result of differences in REM. For each of the five groups determined by the CDF value at REM onset, we found the proportion of REM periods falling into each bin and used these proportions as weights for the weighted averages ().

Spectral density and power estimation

The EEG and EMG signals were sampled at 1000 Hz. The hypnogram was binned in 2.5 s epochs. Spectral densities of the EEG were computed using Welch’s method with a Hanning window for 3 seconds long, half overlapping intervals, resulting in a frequency resolution of 1/3 Hz. Frequency bands were defined as follows: δ: 0.5–4.5 Hz, θ: 5–9.5 Hz, σ: 10–15 Hz. To calculate the power for each frequency band, we approximated the corresponding area under the spectral density curve using a midpoint Riemann sum.

Bootstrap comparison of light and dark phase GMM parameters

We performed bootstrapping with 10,000 iterations on the entire data set for both the light and dark phase and estimated the GMM parameters using the same procedure as explained above. To compare differences between the two phases, we computed for klong, μlong, and σlong the log-function relating each parameter to REMpre and then determined the average value of that function over the entire range of REMpre for which the conditional GMM is defined. Repeating this procedure for each bootstrap iteration resulted in two distributions for each parameter corresponding to the light and dark phase, respectively. Welch’s t-test revealed significant differences between the light and dark phases (klong: t = -77.26, p = 0.0; μlong: t = -372.77, p = 0.0, σlong: t = -57.93, p = 0.0).

Statistics

Statistical tests were performed using the python modules scipy [72] and pingouin [73] and the R package cocor [74]. Linear regressions were performed with the python module statsmodels [75]. For comparisons of quantities between two groups, we used the Levene test to check the homoscedasticity of the data and performed either t-tests or Welch’s t-tests. To compare quantities between multiple groups, the data sets were compared using either one-way or two-way ANOVA followed by multiple comparisons tests.

Conditional GMM for the light phase and model validation.

(A) Comparison of logarithmic (solid lines) and linear (dashed lines) fits for the functions describing the relationship between REMpre and each GMM parameter. For klong, we only fitted the functions to the first 6 REMpre bins (filled circles). For the remaining REMpre bins, klong was set to 1. (B) Estimated probability density functions (PDFs) of the Gaussian distributions for short and long cycles for each 2.5 s increment of REMpre in the range 2.5 s—12.5 s. (C) Histogram over |N| on the normal (left) and natural log scale (right). The histogram for the actual data is compared with the prediction by the GMM (10,000 model simulations).

(PDF)

Click here for additional data file.

Mean and standard deviation of inter-REM and |N| as a function of REMpre.

(A) Mean and standard deviation of ln(inter-REM) as a function of REMpre (n = 5090). (B) Mean and standard deviation of ln(|N|) as a function of REMpre.

(PDF)

Click here for additional data file.

Prefrontal EEG during REM sleep for sequential and single cycles.

(A) Spectral density of prefrontal EEG during REM sleep for different REM sleep durations (REMpre). (B) Spectral density of prefrontal EEG during REM sleep for sequential and single cycles. Solid lines represent the actual densities with shadings indicating the 99% CIs. The dashed lines represent the weighted averages of the duration-dependent densities in A. The weighted averages were calculated based on the proportion of REMpre values falling into each 30 s bin (Methods).

(PDF)

Click here for additional data file.

Comparison of different MA thresholds (1).

(A) Pie chart showing the percentage of REM, NREM, and Wake for different MA thresholds. All wake episodes with duration ≤ 30 s or ≤ 10 s were scored as MA; for the threshold of 0 s, no MAs were scored. (B) Spectral density of parietal (top) and prefrontal EEG (bottom) during NREM sleep for both sequential and single cycles using different MA thresholds. Horizontal lines indicate frequencies at which the spectral densities for sequential and single cycles are statistically different (Welch’s t-test, *** p<0.001). Shadings, 99% CI. (C) Progression of θ power, σ power, and MA rate throughout the refractory and permissive period for different MA thresholds. The durations of both the refractory and permissive period were normalized to unit length and subdivided into quartiles of equal normalized duration. The average for all REMpre values is shown in black. Shadings, 99% CI.

(PDF)

Click here for additional data file.

EEG power, spindle rate, and MA rate throughout the refractory and permissive period.

(A) Spectral density of the parietal EEG for NREM sleep during the refractory and permissive period. Horizontal lines indicate frequencies at which the spectral densities for sequential and single cycles are statistically different; (Welch’s t-test, *** p<0.001, nrefractory = npermissive = 3934). (B) Progression of δ power throughout the refractory and permissive period for different ranges of REMpre. The average for all REMpre values is shown in black. Shadings, 99% CI. (C) Progression of θ power, σ power, rate of spindles, and rate of MAs on normalized time scale throughout NREM sleep. Similar to but with the last 40 s of NREM sleep before the onset of the next REM episode excluded. (D) Progression of spindle rate and MA rate on the non-normalized time scale throughout the first 600 s of NREM sleep during the inter-REM interval for different values of REMpre. The two vertical dashed lines indicate the lowest and highest threshold separating the refractory from the permissive period derived from the low and high bounds of the corresponding REMpre range. Shadings, 99% CI. (E) Progression of θ power, σ power, spindle rate, and MA rate on the non-normalized time scale throughout the first 600 s of NREM sleep after a wake period. Different ranges of wake episode durations were selected to match the corresponding durations of REMpre in D. Shadings, 99% CI.

(PDF)

Click here for additional data file.

Comparison of different MA thresholds (2).

(A) Box plots comparing total NREM duration, |N|, for single cycles with increasing values of total wake duration, |W|, for different MA thresholds (30 s: Welch’s ANOVA, F(5,832.43) = 320.67, p = 2.87e-191; 10 s: Welch’s ANOVA, F(5,1850.87) = 359.09, p = 2.22e-269; No MA: Welch’s ANOVA, F(4,1951,53) = 1113.21, p = 0.0). The x-tick q0 corresponds to cycles without wake. The remaining cycles with |W| > 0 were subdivided into quintiles, labeled q1—q5, based on the distribution of |W| for single cycles. For the threshold of 0 s, there was no cycle without wake. (B) Box plots comparing |N| for single cycles based on the number of wake episodes occurring during the inter-REM interval (30 s: Welch’s ANOVA, F(5,249.10) = 281.68, p = 2.54e-100; 10 s: Welch’s ANOVA, F(8,511.64) = 202.86, p = 2.55e-153; No MA: Welch’s ANOVA, F(19,681.55) = 448.18, p = 0.0). (C) Progression of θ power, σ power, spindle rate, and MA rate during NREM sleep before and after a wake episode for different MA thresholds. Only sequences with at least 1 minute of NREM sleep both before and after wake during the inter-REM interval of single cycles were included. The duration of NREM episodes was normalized. ‘Before’ refers to all NREM sleep in between either the previous REM or wake episode and the current wake episode. ‘After’ refers to NREM sleep in between the current wake episode and either the next wake or REM episode. Shadings, 99% CI.

(PDF)

Click here for additional data file.

Relationship between wake episodes and NREM sleep.

(A) Scatter plots for |W| vs. |N| during single cycles for an increasing number (1–10) of wake episodes during the inter-REM interval. For each number of wake episodes, |W| and |N| are positively correlated. Red lines, linear regression fits. (B) Scatter plots for the number of wake episodes vs. |N| during single cycles for increasing ranges of |W|. In each case, the number of wake episodes and |N| was positively correlated. Red lines, linear regression fits.

(PDF)

Click here for additional data file.

Variables influencing REM episode duration and EEG.

(A) Box plot comparing the duration of REM sleep (REMpost) following single cycles for different values of inter-REM. Red line, linear regression (slope = 0.0012, R2 = 4.92e-04, p = 0.17). (B) Spectral density of parietal EEG during REM episodes following single cycles for different values of the CDF at REM onset (δ: ANOVA, F(4,3794) = 1.27, p = 0.27; θ: Welch’s ANOVA, F(4,1158.22) = 16.18, p = 6.65e-13; σ: Welch’s ANOVA, F(4,1161.57) = 16.27, p = 5.62e-13). (C) Comparison of the true and estimated spectral densities of REM periods for different ranges of the CDF values at REM onset. Solid lines represent the true spectral densities with shadings representing the 99% CIs. Dashed lines indicate estimated weighted averages (Methods). (D) Comparison of the relationship between the CDF (REM propensity) at REM onset and REMpost for single cycles with wake (|W| > 0) and single cycles without wake (|W| = 0). The first row shows box plots comparing REMpost of single cycles with |W| > 0 with different CDF values at REM onset. The first column shows the correlation for the full range of REMpre; the remaining columns display the results for different bins of REMpre. The second row shows the same comparisons for single cycles with |W| = 0. Dashed lines, linear regression (‘s’ indicates the slope of the regression line). The third row shows the linear regression results of CDF vs. REMpost for single cycles with |W| > 0 and single cycles with |W| = 0. Shadings indicate 95% CIs obtained from 10,000 bootstrap iterations. Columns correspond to different bins of REMpre.

(PDF)

Click here for additional data file.

Conditional GMM for the dark phase.

(A) Comparison of logarithmic and linear fits for the functions describing the relationship between REMpre and each GMM parameter for the dark phase. For klong, the functions were only fitted to the first 4 REMpre splits (filled circles). (B) Estimated PDFs of the short and long Gaussian distributions for each 2.5 s increment of REMpre in the range 5 s—15 s for the dark phase.

(PDF)

Click here for additional data file.

Inter-individual variability.

Mean and standard deviation across animals of key variables in the sleep pattern of mice during the light and dark phase. Note that the R2 values for REMpre vs. |N|, |W|, or inter-REM differ from those in , as linear regression was performed for each animal individually before averaging, instead of computing R2 values for the whole data distribution from all animals.

(PDF)

Click here for additional data file.

GMM parameters for the light phase.

The columns indicate for each range of REMpre the number of sleep cycles, the weight of the long Gaussian distribution, mean and standard deviation of the short and long Gaussian distribution and the p-value obtained for the Lilliefors-corrected Kolmogorov-Smirnov (KS) test (Methods).

(PDF)

Click here for additional data file.

P-values from Lilliefors-corrected KS-test for GMMs estimated for different MA thresholds.

The GMM was fitted on the light phase data set for varying thresholds used to score MAs. The columns show the p-values of the Lilliefors-corrected KS-test for the different MA thresholds.

(PDF)

Click here for additional data file.

Coefficients of the conditional GMM for the light phase.

Each column shows the coefficients (a,b,c) for the logarithmic or linear functions describing each GMM parameter as a function of REMpre (Methods). For all parameters, we used a logarithmic fit except for σshort.

(PDF)

Click here for additional data file.

GMM parameters for the dark phase.

The columns indicate for each range of REMpre the number of sleep cycles, the weight of the long Gaussian distribution, mean and standard deviation of the short and long Gaussian distribution and the p-value obtained for the Lilliefors-corrected KS test for the dark phase (Methods).

(PDF)

Click here for additional data file.

Coefficients of the conditional GMM for the dark phase.

Each column shows the coefficients (a,b,c) for the logarithmic or linear functions describing each GMM parameter as a function of REMpre for the dark phase (Methods). For all parameters, we used a logarithmic fit except for σshort.

(PDF)

Click here for additional data file.

30 Apr 2021

Dear Dr Weber,

Thank you very much for submitting your manuscript "A probabilistic model for the ultradian timing of REM sleep in mice" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Paul Franken

Guest Editor

PLOS Computational Biology

Wolfgang Einhäuser

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note that the review of reviewer #2 is uploaded as attachment.

Reviewer #1: This paper addresses an important unsolved problem in sleep research: the genesis of REM/NREM sleep cycles. The authors find that the dynamics in mice can be adequately explained in terms of a bifurcation between sequential REM bouts (i.e., essentially a continuation of the REM state) and a refractory period that separates distinct REM bouts, with the length of this period being dependent on the length of the preceding REM bout. The paper provides a powerful treatment of this issue. I applaud the authors for being so incredibly thorough with their analysis of the data. This approach allows a number of alternate interpretations to be directly tested. Overall I enjoyed reading the paper. Please see my suggestions below.

1) Regarding the sequential bouts, the authors muse on the idea that the interceding NREM bouts appear to have some features of REM and thus may be a sort of intermediate state. Have the authors considered this pattern in the context of local sleep?

2) The genesis of REM/NREM sleep cycles has been considered previously using dynamical systems models. These models have helped to propose and test potential theoretical mechanisms. Specifically, see:

Diniz Behn CG, Ananthasubramaniam A, Booth V. Contrasting existence and robustness of rem/non-rem cycling in physiologically based models of rem sleep regulatory networks. SIAM Journal on Applied Dynamical Systems. 2013;12(1):279-314.

Phillips AJ, Robinson PA, Klerman EB. Arousal state feedback as a potential physiological generator of the ultradian REM/NREM sleep cycle. Journal of theoretical biology. 2013 Feb 21;319:75-87.

It would be valuable to discuss the present findings in the context of these candidate models. I can imagine that the sequential states could be readily explained within a dynamical model by introducing stochasticity, as this could cause temporary transit between REM and NREM sleep, even before a REM sleep bout has really concluded. The positive association between REMpre and N is clearly in line with models that assume a homeostatic balance of the stages (and probably not in line with models that assume a clock-like mechanism). However, it is less clear to me whether the findings for W are consistent with existing dynamical models, or whether models would need to be amended to explain these results.

3) I would avoid using the word "impact" in cases where you do not establish causality or direction of effects (which is not established simply by temporal precedence), such as on line 137.

4) I think that the treatment of REMpre<30 s requires some more consideration/justification. In the text, the authors say that "we found that for values of REMpre < 150 s,142 ln(|N|) forms a bimodal distribution". But this is clearly not the case for REMpre < 30 s, so this statement is not accurate. The argument for effectively treating this as bimodal on Line 723 is not strong. A distribution can be unimodal without being Gaussian, as appears to be the case here. A decision on the basis of the Shapiro-Wilk test therefore does not make sense. My interpretation looking at Fig. 2E is that the unimodality for short REM bouts is effectively due to blending of the two types of NREM bouts (short and long), which would arguably justify treating it as a bimodal distribution.

5) Line 752: "We calculated the residual sum of squares (RSS) forlong753 both the linear and logarithmic fits and chose the function for which the RSS was lower." - This is a 2-parameter vs. 3-parameter fit comparison. Since the slope of the logarithmic function is a/(x+b)+c, it can also be made to approximate a linear fit for very large |b|. Without bounds on a, b, c, this does not seem to be a meaningful comparison (i.e., the log function is theoretically guaranteed to do better), as the logarithmic function could always be made to fit the data at least as well as the linear fit based on RSS. The comparison would make sense only with adjustment for the number of parameters (e.g., with AIC). It may be simplest to just drop the linear fits.

Reviewer #2: uploaded as an attachment

Reviewer #3: This manuscript is an important and stimulating contribution to the understanding of the ultradian organization of sleep that gives further support to notion that hourglass processes may significantly contribute to shape the architecture of the sleep –wake cycle. It also recognizes the centrality of REM sleep episodes (and its associated hourglass mechanism) in the timing of ultradian sleep cycles. The analytic strategy is elegant, solid and highly reliable, and is an invitation to be replicated in other animal models including humans. This work also confirms in mice the existence of two separate categories of REM sleep episodes (sequential and single) that were described by the Bologna group more than twenty years ago (Amici et al 1994 ). It also confirms that the factor “time of day” (see below) affects the timing of the REMs hourglass process (as described in ref 24: Vivaldi et al. 2005).

In my opinion, the main theoretical contribution the authors are: (i) incorporates the recently proposed REMs dependent “refractory period” as a relevant timing mechanism within the ultradian cycle of mice as has been described in rats; (ii) that the REMs “refractory period” is an exclusive attribute of single REMs episodes.; (iii) that REMs episodes duration may be determined within ultradian cycles by REMs propensity estimated by the Cumulative Density Function of NREM.

I have the following observations:

1. The definition of REM sleep propensity as the CDF is clear but the relationship between N an W in the buildup of REM sleep propensity is in my opinion not clearly stated. It is a complicated issue because it could be discussed form the point of view of the relationship between REMpre vs. inter-REMinterval, and between inter-REMinterval vs REMpost.

The paper discussed mostly the second alternative, as stated in line 653: “REM propensity and its impact on the REM episode duration i s more closely associated with the time spent i n NREM sleep than with the combined time i n NREM sleep and wake.” Authors limited the analysis of REMpre vs. inter-REMinterval to the results presented in Fig 1B. The problem of direct REMs episode duration with following amount of wakefulness (in inter-REM-interval) is the interference (masking) of feeding, drinking, exploring, etc, activities that obviously are outside sleep modulatory proceses, introducing loud noise in the timing of wakefulness with respect to the ultradian temporal frame.

I would suggest to explore the relationship of REMpre and REMpost with following interval duration and W duration by setting CDF constant.

2. Detailed information should be presented regarding the inter-individual variability. On the other hand, no details are presented with respect to the validation of automated sleep scoring (line 709: “we manually verified the automatic classification using a graphical user i nterface, visualizing the raw EEG, EMG signals, spectrograms and hypnogram”).

3. Results regarding night-day comparisons were not included in Discussion section? On the other hand I consider inappropriate the term “circadian” when referring to night-day comparisons. Circadian processes are estimated in circadian time or under specific experimental conditions that exclude local or zeitgeber time as main determinant, and is not the case in this experiment.

4. Presented results (particularly those related to sequential vs. single REMs) should be discussed in the light of evidences obtained in rats. In this sense I would suggest authors to review the following reference (Amici, et al. https://doi.org/10.1142/9781860947186_0008, in Parmeggiani &Velluti Eds. The Physiologic Nature of Sleep) where rat neurophysiological and vegetative aspects of REM sleep physiology are summarized. In this sense, there are interesting similarities between mice and rats in the neurophysiological description of sequential episodes.

Minor:

Place references 31 (Amici 1994) before reference 31 (Zamboni et al. 1999) , as the original description correspond of sequential/single episodes correspond to Amici.

Reviewer #4: Summary: This is a very interesting paper that uses a wide range of statistical approaches to provide a highly detailed analysis of REM sleep in mice. As the authors suggest, their analyses establish a useful framework for characterizing “the ultradian regulation of REM sleep in health and disease.” In addition, these analyses identify three major factors shaping ultradian regulation of sleep that have implications for the physiological basis of REM sleep regulation. My specific questions and concerns are listed below.

Major concerns:

1. p.4, line 94: more detail/explanation is needed to support the claim that REM sleep sequences have been observed in multiple species including humans. The timescales of these sequences (and ultradian cycles in general) vary widely across species, so these terms should be defined carefully.

2. P. 25: the authors describe a semi-automatic scoring approach for determining sleep/wake behavior. Has this scoring approach been published previously? Validated against other rodent scoring?

3. Related to the question about scoring, how was the 20 sec threshold for microarousals determined? This choice significantly affects the analysis through the definition of |N| and |W|.

4. The authors should be careful with the language used in the interpretation of the model, particularly in the last paragraph of the Results on P. 18 and 19. Since the model has been fit to the data, the model parameters describe the data but do not explain features such as an increased percentage of single cycles or contribute to an increase in NREM relative to REM sleep. Similarly, on p. 2, line 24: it seems too strong to say that the model “confirms” the existence of two types of sleep cycles; it would be better to say “supports” here.

Minor concerns:

1. Some of the citations in the Introduction seem a little sparse; e.g.,

a. P.3, line 67: citations to REM networks should include work from Yang Dan’s group such as Chen et al., Neuron, 2018

b. P. 4, line 79: citations to REM homeostat should include Franken et al. (ref [51])

2. P. 6, line 141: more detail is needed for the construction of the probability densities in Figure 2B; I am interpreting the terminology “REM_pre splits” to refer to subsets of REM_pre, but I am unclear on the significance of these subsets. How were the 30 second bins chosen? Were there similar numbers of REM_pre episodes in each subset? Approximately how many?

3. P. 8, line 198-201: the interpretation of the duration of the refractory period as a function of REM_pre needs clarification. Does this correspond to the increased black region below the red band in the heat map in Fig 2E?

4. P. 9, lines 216 – 223: this comparison with previous work may be better in the Discussion.

5. P. 9, line 229: what is the |N| duration associated with optimal separation in Figure 3A (right)? How does this threshold vary with |N|, and using this threshold, what percentage of sequential cycles are misclassified as long and vice versa (for different values of |N|)?

6. P. 9, line 235: the interpretation of a sequence of REM sleep with only one cycle as comprising 2 REM periods is given in the caption of Figure 3 caption but not in the text. It would be useful to define in both places.

7. P. 17, line 498: the authors claim that the REM duration is more closely correlated with the CDF than the preceding NREM duration. Is this conclusion supported by a formal analysis of the correlations (e.g., Williams correlation test)?

8. P 18, line 508: careful language should be used when talking about circadian differences; although circadian modulation of REM sleep has been established in other work, these results focus on light/dark differences.

9. P. 18, line 537-538: it would help to remind the readers that k_long is the weighting parameter for the long REM bouts in the Gaussian mixture model.

10. P. 20, line 583: the sleep phase N1 is increased in narcolepsy, and, although this light phase of sleep is associated with transitions to wake, fragmentation of N1 is not increased in narcolepsy (type 1) compared to controls (Maski et al., SLEEP, 2021).

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

Reviewer #4: No: I may have missed it, but it looked like only the raw data and not the code were available in the repository the authors indicated.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: H. Craig Heller

Reviewer #3: No

Reviewer #4: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

Submitted filename: Review of REM-NREM paper.docx

Click here for additional data file.

9 Jun 2021

Submitted filename: Response_to_reviewers.pdf

Click here for additional data file.

21 Jul 2021

Dear Dr Weber,

Thank you very much for submitting your manuscript "A probabilistic model for the ultradian timing of REM sleep in mice" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Dear Franz:

Sorry for the delay in responding. Reviewer 4 still found a minor issue that needs your attention. Reviewer 2 (Craig Heller) maintains that your use of the term refractory is incorrect. I disagree with Craig as I introduced this term in the 2002 paper which you already cited and also the more recent work by Le Bon (2021) and Vivaldi (2020), both already cited, use this terminology (there are probably other publications). Even the algorithm that Joel used to determine NRTs (Benington, Kodali, Heller SLEEP 1994), and I used in the 2002 paper, has a 40s delay build-in, in which NRTs cannot occur. I leave it up to you whether you want to add a sentence qualifying your use of the term refractory. To conclude, consider your manuscript as provisionally accepted and I hope you can still revise your manuscript by addressing these last two issues.

with kind regards,

Paul

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Paul Franken

Guest Editor

PLOS Computational Biology

Wolfgang Einhäuser

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Thank you for carefully addressing all my comments. This is an excellent paper and a significant contribution to our understanding of REM/NREM sleep cycles.

Reviewer #2: Thanks for the interesting discussion. I re-iterate that I am not scoring your manuscript even though I think your data are great. I just don't have the expertise to evaluate your modeling that is a major portion of the paper. However, I was trying to convince you that you have a perfectly good physiological explanation of your data without invoking a mysterious refractory phase, which the data do not support. That explanation is classical homeostatic regulation, which your data do support. The definition of refractory is to be insensitive to a stimulus that was previously effective. A neuron is refractory to depolarizing input immediately following an action potential. A man is refractory to sexual stimulation following orgasm. If you went out into the cold without a coat and started to shiver, but then put on a coat and stopped shivering, are you refractory to cold? No, you have just eliminated the stimulus. At one point you also seem to accept someone else's description of the REM/NREM process as being like an hour glass --- no it is not. Maybe an old fashioned ice box is like an hour glass as the ice melts, but a modern refrigerator is a regulated system. The NREM REM cycles are regulated homeostatically as shown by both our results. If you are looking for a function of REM, you are more likely to find it by looking for the feedback information relating NREM to REM, and not in a supposed mechanism of refractoriness.

As we showed, the brief REM episodes are likely to occur at any time in the classical sleep cycle, but are increasingly likely as the cycle continues. Also, a brief REM episode is likely to be followed soon by another one, and so on with increasing rapidity. These brief REM episodes are failures to maintain REM, and when they fail they do not dissipate much REM propensity, so they come again sooner and sooner until a REM bout is sustained and dissipates the REM need.

Physiologically, we need to be looking for what process is exhausted or built up during the hyperpolarized NREM state that has to be reversed during REM. That is investigating a homeostatic process and not a refractory mechanism.

You can interpret your data anyway you want in your paper, of course. But, I just want you to try thinking of your beautiful data in a physiological way.

Reviewer #3: I have no further concerns.

Reviewer #4: Summary: This is a very interesting paper that uses a wide range of statistical approaches to provide a highly detailed analysis of REM sleep and the impact of REM sleep on NREM sleep in mice. The authors have addressed my previous comments in the revised manuscript, and the new Discussion subsection on Sequential Sleep Cycles is a great addition. One few minor suggestion:

Lines 612-615: The authors state that “In rats and, as shown here, in mice the average duration of REM sleep episodes preceding a sequential cycle is on average shorter than that before a single cycle, and the frequency of sequential cycles is reduced during the dark period [33].” Is this comment actually restricted to the REM sleep episodes preceding these events or should it describe REM sleep episodes within sequential cycles compared to REM sleep episodes in a single cycle? If the authors are highlighting a feature of the first REM sleep episode in a sequential cycle, it would be helpful to clarify this point.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

Reviewer #4: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Andrew J K Phillips

Reviewer #2: Yes: H. Craig Heller

Reviewer #3: Yes: Adrián Ocampo-Garcés, MD, PhD

Laboratorio de Sueño y Cronobiología

ICBM, Facultad de Medicina

Universidad de Chile, Santiago

Chile

Reviewer #4: Yes: Cecilia Diniz Behn

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

28 Jul 2021

Submitted filename: response_to_reviewers.pdf

Click here for additional data file.

29 Jul 2021

Dear Dr Weber,

We are pleased to inform you that your manuscript 'A probabilistic model for the ultradian timing of REM sleep in mice' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Paul Franken

Guest Editor

PLOS Computational Biology

Wolfgang Einhäuser

Deputy Editor

PLOS Computational Biology

***********************************************************

20 Aug 2021

PCOMPBIOL-D-21-00423R2

A probabilistic model for the ultradian timing of REM sleep in mice

Dear Dr Weber,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Livia Horvath

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol