HFOs in PFC during REM sleep are associated with theta modulated activity.

(A) CA1 and PFC theta-delta (TD) ratio time locked to either the start or end of a REM sleep bout. Start and end times were determined using CA1 TD ratio. Note the similar increases and decreases in CA1 and PFC TD ratio surrounding REM bout starts and ends, respectively. (B) Raster plot of CA1 and PFC cells during REM sleep with PFC theta and HFO amplitude plotted above. Gray shaded areas indicate PFC HFOs, and arrows denote periods expanded on the right to show single unit spiking during HFO events. (C) (Left) Example PFC HFOs during REM sleep. Note the theta frequency fluctuations in the local field potential. (Right) Example grand average PFC HFO waveform (top) and amplitude (bottom) across all events in an example animal. (D) Example PFC event triggered spectrogram for REM and NREM HFOs/ripples. Note the elevated gamma (40-100 Hz) and theta (6-12 Hz) band activity associated with REM HFOs, in contrast to the elevated slow oscillation (0.1-4 Hz) and spindle (12-16 Hz) band power in during NREM ripples. (E) PFC event rates in REM and NREM sleep (REM = 0.66 ± 0.06, NREM = 0.54 ± 0.05, *p = 0.013, WSR). (F) PFC HFO autocorrelation in REM sleep. Note the peaks in the autocorrelation that repeat approximately every 125 ms. (G) Distribution of inter-event intervals (IEIs) for NREM and REM events. Note the prevalence of IEIs less than 200 ms for REM HFOs (Mode = 0.134 s). Only IEIs up to 1 s are shown for visualization purposes. (H) (Left) REM HFO aligned multiunit activity (MUA) in PFC, and (Right) the power spectral density (PSD) calculated for two example animals. The vertical dashed gray lines indicate the peak theta frequency of the population activity fluctuations surrounding HFO events. (I) Same as in (H) but for NREM ripples. Top and bottom PSDs are from the animals in (H).

Theta, gamma, HFO coupling underlies distinct REM modulation states in PFC.

(A) Averaged modulation across phase-amplitude pair occurrences during high theta power REM periods. The phase of the theta oscillation at a peak frequency of 7-8 Hz strongly modulates the amplitude of higher frequency oscillations. (B) Theta nesting for gamma and HFO frequency oscillations for an example tetrode illustrating coupling of theta, gamma, and HFO oscillations in PFC during REM sleep (see Methods). Briefly, averaged nesting plots were generated by aligning the LFP to the peak power in the preferred theta frequency phase bin for either gamma or HFO oscillations separately. (C) Distribution of preferred phases for PFC cells locked to HFO or gamma events (U2 = 4.74, ***p = 4.56×10−41, Watson’s U2 test). (D) CA1-PFC theta (6-12 Hz) coherence centered on PFC REM HFOs is elevated relative to baseline coherence. (E) Rate of events designated as participating in a chain (<200 ms separation) during NREM and REM sleep (NREM = 0.15 ± 0.02 Hz, REM = 0.25 ± 0.04 Hz, **p = 0.0054, WSR). (F) Example REM HFO chain (doublet) plotted across 4 separate tetrodes with the average HFO power plotted below. (G) PFC theta pairwise phase consistency (PPC) of HFO chains in REM vs NREM (REM = 0.22 ± 0.06, NREM = 0.06 ± 0.02, **p = 0.0097, WSR). (H) HFO chains showed higher PFC (Left) and CA1 (Right) theta PPC compared to isolated HFOs (PFC: Chain = 0.22 ± 0.06, Isolated = 0.09 ± 0.04, *p = 0.013, WSR; CA1: Chain = 0.30 ± 0.08, Isolated = 0.14 ± 0.05, **p = 0.0016, WSR). (I) (Left) Chain and isolated HFO triggered spectrograms. (Right) Power in the gamma (top, 40-100 Hz) and theta (bottom, 6-12 Hz) frequency bands for isolated and chained HFOs. The average power in each band in the ±500 ms window surrounding HFO onset was used for comparison (Gamma: Chain = 0.72 ± 0.06, Isolated = 0.45 ± 0.03, ***p = 8.12×10−7, WSR; Theta: Chain = 0.60 ± 0.20, Isolated = 0.15 ± 0.11, **p = 0.0021, WSR) (J) (Left) CA1-PFC theta (6-12 Hz) coherence during chained vs isolated events and (Right) quantification of the peak coherence in a ±1 second window from HFO onset (Chain = 0.93 ± 0.08, Isolated = 0.65 ± 0.05, **p = 0.0069, WRS). (K) (Left) Phase slope index (PSI) peaked at theta frequency, indicating a CA1 to PFC directed flow of information during HFO chains. (Right) Quantification of PSI in the theta band (6-12 Hz). Each data point corresponds to the averaged PSI across each electrode pair in an epoch. (PSI = −0.068 ± 0.021, **p = 0.0025, t-test against zero).

Differential modulation of PFC neuronal activity during NREM- and REM-sleep prefrontal high-frequency events.

(A) (Left) PFC neurons that were positively modulated (EXC) during REM and NREM PFC HFO/ripple events (REM EXC modulation = 142 out of 922 candidate neurons; NREM EXC modulation = 783 out of 1079 candidate neurons; See Methods section Ripple/HFO aligned modulation). (Right) Average z-scored peri-event time histograms (PETHs) aligned to events. (B) Timing of peak excitation for significantly modulated PFC cells during NREM and REM. Note the broader distribution of peaks during REM sleep (Comparison of peak timing distributions, ***p = 8.29×10−5, WRS). (C) PFC event type (NREM vs REM) can be predicted by population spiking activity (PFC data = 0.68 ± 0.01, shuffle = 0.50 ± 8.00×10−5, ***p = 2.56×10−34, WRS). (D) Modulation of putative CA1 pyramidal neurons aligned to PFC REM HFO events (REM EXC modulation = 66 out of 680 candidate neurons). (E) PFC cofiring and assembly reactivation strength were overall stronger during NREM ripples (Cofiring: NREM = 0.96 ± 0.03, REM = 0.48 ± 0.03, ***p = 3.58×10−27; Reactivation strength: NREM = 2.01 ± 0.05, REM = 1.56 ± 0.07, ***p = 6.91×10−16, WRS). (F) (Left) PFC multiunit activity aligned to chained vs isolated REM HFO events. (Right) HFO aligned multiunit PSD illustrating stronger theta modulated spiking activity in PFC during chained events (mean power in the 6-10 Hz band, *p = 0.029, WRS). (G) The probability of isolated and chained REM PFC HFOs centered on population suppression events detected in REM sleep. Note the strong peak around 0 in the HFO chain distribution, indicating that REM PFC HFO chains are associated with overall decreases in population activity as compared to isolated HFOs. The black dotted lines indicate the 99% confidence intervals based on suppression events extracted from shuffled multiunit activity.

Structured reactivation in PFC during REM HFOs.

(A) Example rasters during HFO chains. Spikes are color-coded based on the HFO number within the chain during which the cell was active. Note the different cell assemblies active during each HFO within a chain. (B) Correlation of PFC population spiking between adjacent REM HFOs and its relationship with inter-event-intervals (r = −0.13, ***p = 6.00×10−8, Pearson correlation). Note that the average IEI for quartile 1 is 0.114 s, which corresponds to 8-9 Hz during HFO chains. (Inset) The Pearson correlation coefficient compared to a distribution of r values generated from shuffling (n = 1000) HFO identity. (C) Rank-order correlation analysis. (Left) Example rank-order template from an example sleep epoch. Templates were generated by concatenating all HFOs in chained events (inter-event periods excluded) and averaging the normalized rank across chains. (Right) The mean rank-order correlation compared to a distribution of means generated by jittering (n = 1000) spike times. (D) (Left) PFC assembly reactivation aligned to NREM ripples or (Right) the first HFO of each REM HFO chain. The assemblies are sorted based on the average peak reactivation bin within a ±500 ms window surrounding the first HFO in a chain. Note the tiling of assembly peaks aligned to HFO events in REM sleep, which manifests in the relatively flat average response in REM vs. a singular peak of synchronous assembly reactivation in NREM as shown in (E). (E) (Left) Averaged event triggered reactivation strength during NREM and REM. (Right) Mean reactivation strength of three example assemblies aligned to high-frequency events in NREM and REM, illustrating the tiling of assembly reactivation relative to events. The example assemblies shown come from the same animal and epoch. (F) Cross-validation of the REM reactivation sequences shown in (D), suggesting a consistent structure of reactivation surrounding REM HFO chains. (Bottom) Distributions of r values calculated from the Pearson correlation between peak reactivation bins across all assemblies for two randomly chosen halves (n = 1000 random splits) of the HFO aligned data (Data = 0.232 ± 0.003, Shuffle = −0.023 ± 0.004, ***p = 1.13×10−264, WRS). (G) Differences in slope of the linear fit between the time of peak reactivation (in milliseconds) and assembly identity for two halves of the HFO aligned data (Data = 0.048 ± 2.21×10−4, Shuffle = 0.070 ± 2.67×10−4, ***p = 2.31×10−302, WRS). (H) Probability of assembly peak displacement between randomly chosen halves of the HFO aligned data. For each assembly, the absolute difference in peak location was calculated between the two halves. To calculate 95% confidence intervals, assemblies were circularly shuffled, and the peak displacement probabilities were calculated over 1000 random splits of the data. Note the peak at small displacement values, indicating consistency of assembly peak activity across splits.

Differential engagement of distinct populations of CA1 neurons during PFC REM HFO chains.

(A) Example raster plot of CA1 and PFC cells during a HFO chain event. Only cells that were active during the plotted time window are shown. Note the coherence of theta band activity in CA1 and PFC with elevated PFC gamma power during each HFO event in the chain. (B) (Left) 2D spatial rate map correlation for non-active, low, and high HFO chain cofiring CA1 pairs (Non-active = 0.123 ± 0.001, Low cofiring = 0.127 ± 0.001, High cofiring = 0.176 ± 0.004, Non vs Low, ***p = 1.26×10−10, Low vs High, ***p = 2.91×10−21, Non vs High, ***p = 2.33×10−51, Kruskal Wallis test with Bonferroni correction). Here, ‘non-active’ indicates cell pairs where cofiring was not assessed because at least one cell was inactive across all events. (Middle) Rate maps and spatial correlation of example low and high cofiring CA1 pairs. (Right) Real probability of high cofiring CA1-CA1 pairs during PFC HFOs (probability = 0.101) compared to a null distribution calculated by shuffling HFO times (n shuffles = 1000). (C) Histograms of CA1-PFC HFO cofiring during CA1-independent NREM PFC ripples and chained REM PFC HFOs. Cofiring during REM HFOs showed a bimodal distribution indicating separable populations of CA1 neurons based on PFC cofiring. (REM, **p = 0.001; NREM p = 0.23, Hartigan’s dip test for bimodality). (D) Firing rate of high cofiring and “other” CA1 cells aligned to HFO chains. Here, ‘other’ includes low cofiring cells and cells that did not have a cofiring value because at least one of the cells in the pair emitted zero spikes across all events (similar to ‘non-active’ above). (Comparison of mean firing rates in the ±100 ms window from HFO onset, ***p = 7.77×10−5, WRS). (E) (Left) High cofiring CA1 neurons exhibited a greater degree of suppression (Modulation index: See Methods section Ripple/HFO aligned modulation) during NREM PFC ripples (Low cofiring = −0.08 ± 0.06, High cofiring = −0.25 ± 0.03, *p = 0.022, WRS; High and low cofiring CA1 neurons were split based on cofiring values > 0 and < 0, respectively). (Right) The cofiring magnitude during HFO chains was correlated with the degree of suppression during NREM PFC ripples (r = −0.23, **p = 0.0045, Pearson correlation). (F) (Left) Changes in CA1 firing rates from the first NREM bout to the last NREM bout within a sleep epoch plotted according to degree of cofiring with PFC (high, > 0; low, < 0), either during CA1-independent NREM ripples or chained REM HFOs (Low NREM = −0.06 ± 0.03, High NREM = −0.02 ± 0.03, p = 0.32; Low REM = −0.12 ± 0.063, High REM = 0.085 ± 0.059, **p = 0.0052, WRS). (Right) Firing rate change separated by cofiring quartile. The change in CA1 firing rates over the course of sleep was correlated with PFC cofiring during REM HFOs but not NREM ripples (NREM, r = 0.125, p = 0.46; REM, r = 0.169, *p = 0.012, Robust linear regression). (G) CA1-PFC HFO cofiring compared for CA1 neurons with firing rate changes higher or lower than the mean (−0.04 Hz) of the distribution (Higher = 0.34 ± 0.04, Lower = 0.22 ± 0.06, *p = 0.038, WRS). Similar to (E), CA1 neurons with positive shifts in firing rate over a sleep epoch had higher PFC cofiring than CA1 neurons with a decrease in rate. (H) (Left) Proportion of CA1-PFC pairs with cofiring greater than 0 during REM HFO chains that contained either non-shifting or REM-shifting CA1 neurons (Non-shifting = 2216/5066 pairs, 43.7%, REM-shifting = 1057/2242 pairs, 47.2%, z = 2.70, **p = 0.0097, Z-test for proportions). (Right) Magnitude of CA1-PFC HFO cofiring (greater than 0) compared for non-shifting and REM-shifting CA1 neurons (Non-shifting = 1.35 ± 0.02, REM-shifting = 1.51 ± 0.03, ***p = 8.30×10−5, WRS). (I) The cofiring magnitude of REM-shifting CA1 neurons was correlated with the degree of suppression during independent NREM PFC ripples (Non-shifting, r = −0.29, p = 0.066, REM-shifting, r = −0.53, **p = 0.0093, Pearson correlation). (J) Burst probability of non-shifting and REM-shifting CA1 neurons during REM HFO chains (Non-shifting = 0.028 ± 0.006, REM-shifting = 0.059 ± 0.017, *p = 0.022, WRS). (K) REM-shifting CA1 neurons showed higher PPC, indicating stronger phase locking to the PFC theta oscillation surrounding HFO chains (Non-shifting = 0.032 ± 0.006, REM-shifting = 0.064 ± 0.011, **p = 0.0026, WRS). Only periods surrounding REM HFO chains were used to evaluate phase locking strength of CA1 neurons.

A simple network model incorporating acetylcholine recreates REM results.

(A) Model schematic. Network of Hodgkin-Huxley (HH) pyramidal and interneurons with an added slow potassium leak conductance. This channel is open in NREM, introducing an adaptation current that lowers cell excitability and bestows resonator properties to all neurons, altering their phase response curves such that connected cells are better able to synchronize spiking. In REM the channel is closed to simulate high ACh concentration, eliminating the adaptation current. (B) Example raster plots of network activity in NREM and REM. Spikes from 400 interneurons shown in red, spikes from 1600 pyramidal neurons shown in black. Shaded regions indicate periods of applied input current (see Methods). (C) Experimental data. (Left) Z-scored theta power surrounding REM PFC HFOs shows an average 3 second increase in theta power. (Right) Normalized MUA aligned to REM PFC HFOs using 100 ms bins shows suppression. (D) Modeling MUA results. (Left) NREM average network firing rate aligned to MUA peaks during stimulus periods, with inset depicting the average stimulus profile. (Center) same as (Left), but for REM. (Right) REM 100 ms binned firing rate aligned to MUA peaks within stimulus periods (i.e., a smoothened version of the Center panel) shows suppression, similar to experimental data. (E) (Left) Experimental data showing the negative correlation between peak theta power during REM HFOs and average MUA during HFOs. (Center) Peak amplitude of applied stimulus vs. average MUA quartiles during stimulus period in the model. (Right) Same as (Center), but for z-score of average theta power (measured from network MUA) during stimulus period. (F) Z-scored cofiring within MUA peaks in stimulus periods in the model. REM peak duration 37 ± 9 ms (mean ± std), NREM peak duration 35 ± 9 ms. (G) Average firing rate surrounding MUA peaks within stimulus periods in the model split between isolated and chained MUA peaks. Isolated peaks defined as having no other peak within a 200 ms.

Stimulation in NREM induces widespread coactivity compared to stimulation in REM.

(A) Model raster plots during example NREM (Top) and REM (Bottom) periods with black trace indicating the percent of all pyramidal cells with at least one spike in 10 ms bins. (B) The same percent active measure (not distinguishing by cell type) calculated on experimental data aligned to NREM and REM events (Left), and in the model aligned to MUA peaks in NREM and REM as a fraction of all pyramidal neurons (Center). Similar plot for all interneurons in the model (Right). (C) Model NREM activity split by cells directly receiving stimulus input (blue) and receiving no external stimulus (pink). Raster of a single example population burst due to stimulation in NREM (Left). Percent of cells with at least one spike in 10 ms bins (as a percent of all pyramidal cells in the network) aligned to stimulus midpoints, the time of peak input amplitude (Center). Similar plot for interneurons (Right). Same as (C), but for REM, for a single stimulation event. Note that non-stimulated cell population (pink) comprise a greater percent of total active neurons outside stimulation times (average stimulation time of 3 sec centered at time of event), as there are more neurons in the non-stimulated group (only 30% of neurons are stimulated per stimulus event).

Amount of NREM and REM sleep (in seconds) across all animals and epochs.

Characterization of REM sleep and prefrontal REM sleep HFOs.

(A) Histological verification of recording locations in dorsal CA1 (Top) and PFC (Bottom). Red arrows indicate electrolytic lesions and scale bars represent 1 mm. (B) Example sleep session showing each behavioral state detected using the sleep state algorithm that utilizes theta-to-delta ratio (TD ratio) and animal velocity (cm/s). (C) (Left) Average velocity during wake and REM sleep (Wake = 1.16 ± 0.06 cm/s, REM = 0.05 ± 0.01 cm/s, *p = 5.39×10−7, WSR). (Right) Theta-to-delta (TD) ratio in NREM and REM sleep (NREM = 0.30 ± 0.11, REM = 2.71 ± 0.10, ***p = 5.39×10−7, WSR). (D) Distributions of intracranial EMG values for wake and REM sleep. All 500 ms binned Pearson correlation values are shown. (E) Comparison of the averaged intracranial EMG for bouts of wake and REM during all sleep epochs (Wake = 0.57 ± 0.01, REM = 0.32 ± 0.02, ***p = 5.26×10−37, WRS) (F) Duration of REM sleep in each epoch (Total duration = 6089 s). Only epochs with at least 30 s of REM sleep were included for analysis. (G) (Left) Average REM sleep duration and (Middle) percentage of REM sleep out of total sleep epoch time (Wake, NREM, and REM) across all 9 sleep epochs. No animals had REM sleep in the first normalized bin, so that point isn’t shown. (Right) Percentage of REM sleep during the first and second halves of the 36 sleep epochs included in analysis (First = 3.85 ± 0.93%, Second = 15.29 ± 0.94%, ***p = 1.93×10−6, WRS). (H) Hypnograms for animals showing sleep state across all 9 concatenated sleep epochs (only sleep epochs are shown; behavior epochs are excluded).

Characterization prefrontal REM sleep HFOs and event aligned multiunit activity.

(A) Average peak aligned amplitude plots across all detected REM HFO events for each animal. (B) Comparison of ripple/HFO duration, amplitude (standard deviations above the mean), and frequency for NREM and REM events (Duration: NREM = 0.111 ± 0.017 ms, REM = 0.126 ± 0.030 ms, p = 0.59; Amplitude: NREM = 4.13 ± 0.12, REM = 4.23 ± 0.20, p = 0.31; Frequency: NREM = 184.96 ± 1.68 Hz, REM = 179.02 ± 1.82 Hz, *p = 0.018, WRS). (C) Average power spectral density calculated from the REM and NREM event aligned multiunit activity in PFC. Note the peak around 8 Hz (theta frequency) for REM HFO aligned multiunit activity. (D) PFC multiunit activity aligned to REM HFO events extracted from a wider frequency range (100-250 Hz). The separation into coordinated and non-coordinated events (Coordinated = 48.75 ± 2.6%) was determined by the overlap of these HFOs with HFOs extracted using the standard frequency range (150-250 Hz). This analysis was performed to demonstrate that the HFOs extracted from the 150-250 Hz range primarily drives the theta frequency multiunit activity in PFC. (E) PFC HFO rates during REM sleep extracted from either CA1 or PFC TD ratio. HFO rates were similar during CA1 and PFC extracted REM sleep, validating the use of CA1-extracted REM sleep for further analysis. (F) PFC multiunit activity and average power spectral density aligned to REM HFO events extracted from the overlap of CA1 and PFC REM sleep. (G) (Left) Peak-to-trough of average spike waveforms and mean firing rate of PFC neurons. We separated populations based on peak-to-trough, since mean firing rate did not give us clear separation. Putative interneurons were identified as cells having a peak-to-trough <0.3 ms. (Right) REM HFO aligned PFC multiunit response when spikes from putative interneurons are excluded. Note the similarity to Figure 1H.

Coupling of gamma and HFOs with PFC theta oscillations.

(A) Example showing bouts of low and high theta power within a REM sleep bout. High theta periods were bouts within REM that exceeded 4 standard deviations above the electrode averaged theta power in REM sleep for at least 900 ms. (B) The CA1 theta inter-peak intervals (IPI) during high PFC theta power periods during REM sleep were lower (corresponding to higher theta frequency) than during low power periods, indicating that these high and low theta power periods may correspond to phasic and tonic substages of REM sleep, respectively (High theta IPI = 0.123 ± 0.002, Low theta IPI = 0.128 ± 8.3×10−4, *p = 0.010, WRS). (C) (Left) Rate of REM HFOs during bouts of high or low theta power (High = 1.19 ± 0.15 Hz, Low = 0.58 ± 0.05 Hz, ***p = 1.33×10−4, WSR). (Right) HFO autocorrelation during periods of high or low theta power. Dots indicate bins where high and low theta power HFO autocorrelations differ (p<0.05, WRS). (D) Modulation indices (MIs) for phase amplitude coupling (PAC) computed for HFO and gamma frequency oscillations compared to MIs calculated from circularly shifted theta phase and HFO/gamma amplitude vectors. Modulation was assessed for a phase frequency range from 2-12 Hz. The amplitude modulated frequencies were 150-250 Hz and 40-100 Hz for HFO and gamma, respectively (HFO = 0.005 ± 3.21×10−4, shuffle = 0.001 ± 8.94×10−5, ***p = 1.13×10−31; Gamma = 0.011 ± 6.49×10−4, shuffle = 0.002 ± 1.67×10−4, ***p = 2.97×10−24, WRS). (E) Frequency of the average theta nested oscillation either aligned to peak HFO or gamma power (See Methods). Note that aligning to peak HFO power at the preferred phase bin yields nested oscillations that are higher in frequency (HFO = 146.48 ± 2.82, Gamma = 128.96 ± 3.091, ***p = 2.85×10−7, WRS). (F) Preferred theta phases of elevated HFO and gamma amplitude (HFO = −0.13 ± 0.08 radians, Gamma = −0.38 ± 0.08 radians, *p = 0.016, WRS). (G) Cross-correlation between PFC gamma and HFO events during REM sleep. The off-center peak at a negative lag (inset) indicates that gamma events tend to precede HFOs, corroborating the phase preference difference in (F). (H) Comparison of nested HFO peaks with a theta phase bin shuffled surrogate (See Methods; HFO = 4.09 ± 0.48 μV, Surrogate = 1.92 ± 0.19 μV, **p = 0.0058, WRS). (I) Same as in (G) but for gamma (Gamma = 11.33 ± 0.62 μV, Surrogate = 2.33 ± 0.20 μV, ***p = 9.46×10−94, WRS). (J) Distribution of preferred phases for PFC cells that are significantly phase locked to either HFOs or gamma. Note the similarity in distribution to Figure 2C in the manuscript (HFO locked = 463 neurons, Gamma locked = 161 neurons, Dual locked = 110 neurons, U2 = 3.56, ***p = 6.37×10-31, Watson’s U2 test) (K) (Left) Gamma triggered spectrograms for events that were uncoupled or coupled to HFOs. (Middle) Comparison of gamma power (40-100 Hz power ±500 ms from event onset) for uncoupled and coupled events. Note the oscillating gamma power during coupled events. (Right) Quantification of average HFO power for coupled and uncoupled events (Coupled = 0.60 ± 0.04, Uncoupled = 0.34 ± 0.03, ***p = 1.03×10−7, WRS). Coupled gamma events were defined as gamma periods that were within 500 ms of a HFO event. (L) (Left) Multiunit PFC activity aligned to gamma events that are coupled with (Left) or uncoupled from HFOs.

Chains of sleep oscillations in NREM (spindles) and REM (HFOs) sleep.

(A) Example isolated and train spindles, similar to previous reports.45 (B) Averaged spindle peak aligned LFP for isolated and train spindles across all animals. (C) Histogram displaying the proportions of NREM spindle train events of specific lengths (Compare with Darevsky et al., 2024). (D) The durations of isolated and train spindles (Isolated = 0.579 ± 0.005 s, Train = 0.550 ± 0.003, ***p = 1.77×10−6, WRS). (E) Histogram displaying the proportions of REM HFO chain events of specific lengths. (F) Histogram of all inter-event intervals (IEI) for HFO chains in (Left) REM or (Right) NREM sleep. Note the tight REM HFO IEI distribution centered around 125 ms, which corresponds to 8 Hz (theta frequency). Black and red vertical lines are displayed at 125 ms for both distributions. (G) PFC multiunit activity aligned to chained and isolated NREM ripple events. Note the absence of theta modulated population activity during ripple chains in NREM sleep. (H) Comparison of HFO amplitude and frequency for chained and isolated PFC REM HFOs (Amplitude: Chain = 4.21 ± 0.26, Isolated = 4.17 ± 0.20, p = 0.995; Frequency: Chain = 172.54 ± 2.71 Hz, Isolated = 181.54 ± 1.51 Hz, *p = 0.018, WRS). (I) (Left) PFC multiunit activity aligned to the preferred theta phase of REM HFOs plotted according to theta power during surrogate events. (Right) Quantification of theta power in each quartile compared to the power during HFO events (Q1 = −1.10 ± 0.007, Q2 = −0.62 ± 0.004, Q3 = −0.20 ± 0.005, Q4 = 0.87 ± 0.024, HFO = 0.26 ± 0.018, ***p = 5.74×10−114, WRS).

REM PFC HFOs are associated with decreases in PFC activity and sequential reactivation as compared to NREM PFC ripples.

(A) PFC neurons that were negatively modulated (INH) during REM PFC HFO events (REM EXC modulation = 114). Although overall modulation in this subset of neurons surrounding PFC HFOs is suppressive, there are transient peaks in activity that occur approximately 125 ms apart (8 Hz theta frequency). (B) Extracted PFC population suppression events (Methods) aligned by their troughs. (Left) All example suppression events in an example epoch. (Right) Averaged PFC multiunit activity aligned to these events. (C) The probability of NREM PFC ripples centered on population suppression events detected in NREM sleep. Unlike REM sleep, there is a decreased probability of NREM ripple occurrence surrounding suppression events, consistent with strong overall excitation during these events (Figure 1I). The slight peak to the right is likely reflective of the suppression in activity prior to the onset of NREM ripples (Figure 1I). The red dotted lines indicate the 99% confidence intervals based on suppression events extracted from shuffled multiunit activity. (D) Rank order correlation analysis on isolated REM PFC HFOs (Left) and pseudo-chain events (Right). (E) Relationship between IRI and (Top) Jaccard index and (Bottom) cosine similarity (Jaccard index, r = −0.094, ***p = 1.30×10−4; Cosine similarity, r = −0.097, ***p = 7.16×10−5, Pearson correlation). (F) PFC multiunit activity aligned to the first HFO in each chained event. Note the peaks in the multiunit activity prior to the initial HFO within a chain, suggesting that fluctuations in PFC multiunit activity precede chained events. Thus, we aligned PFC reactivation strength to the first HFO in each chain for sequence cross-validation analysis. (G) Cross-validation of the REM reactivation sequences surrounding gamma events (Data = 0.147 ± 0.003, Shuffle = −0.051 ± 0.004, ***p = 6.45×10−208, WRS). (H) (Left) Same as in (G) but for isolated REM HFOs (Data = 0.050 ± 0.003, Shuffle = −0.005 ± 0.004, ***p = 4.82×10−4, WRS). (Right) Probability of assembly peak displacement between randomly chosen halves of the isolated REM HFO aligned data. Note the absence of peaks at small displacements as seen in Figure 4H. (I) Same as in (G) but for NREM ripple chains (Data = 0.084 ± 0.003, Shuffle = −0.019 ± 0.004, ***p = 1.37×10−75, WRS). (J) Cross-validation of the NREM reactivation sequences in the post ripple window (+500 ms from first chained ripple onset) (Data = −0.007 ± 0.004, Shuffle = −0.018 ± 0.004, p = 0.062, WRS).

CA1 suppression during NREM PFC ripples and CA1 REM HFO cofiring.

(A) Absolute CA1-CA1 cell pair cofiring during chained and isolated HFOs (Chain = 0.53 ± 0.02, Isolated = 0.48 ± 0.02, ***p = 1.42×10−4, WRS). The magnitude of cofiring (positive and negative) was greater for HFO chains as compared to isolated HFOs, suggesting differential engagement of CA1 neurons. (B) Firing rates of CA1 pyramidal neurons during the first 30 s of NREM, REM, and the last 30 s of NREM sleep within an epoch (First 30 s NREM = 0.54 ± 0.04 Hz, REM = 0.51 ± 0.04 Hz, Last 30 s NREM = 0.51 ± 0.03, First vs REM, ***p = 1.10×10−4, REM vs Last, p = 0.61, First vs Last, *p = 0.013, Friedman test with Bonferroni correction). CA1 firing rates decreased during REM and later in NREM. (C) During NREM sleep, CA1 neurons were strongly suppressed during independent PFC ripples. (D) (Left) There was no difference in the degree of suppression during NREM PFC ripples for low and high cofiring (during isolated REM HFOs) CA1 neurons (Low cofiring = −0.16 ± 0.04, High cofiring = −0.21 ± 0.03, p = 0.40, WRS). (Right) The cofiring magnitude during isolated HFOs was correlated with the degree of suppression during NREM PFC ripples, similar to REM HFO chains (r = −0.19, *p = 0.011, Pearson correlation). (E) Correlation between CA1 modulation during independent NREM PFC ripples and coordinated SWRs. The more suppressed CA1 neurons are during PFC ripples, the stronger they are reactivated during coordinated SWRs (r = −0.85, ***p = 6.45×10−208, Pearson correlation). (F) There was no relationship between CA1 cofiring (with PFC) during REM HFO chains and the degree of modulation during coordinated SWRs in NREM sleep (r = 0.04, p = 0.44, Pearson correlation).

Relationship between REM PFC HFO cofiring, changes in CA1 firing rates over the course of sleep, and CA1 theta phase shifting neurons.

(A) (Left) Mean firing rates of high and low cofiring CA1 pyramidal neurons across the entire sleep epoch or (Right) restricted to REM sleep (Low cofiring mean rate = 0.69 ± 0.10, High cofiring mean rate = 0.79 ± 0.08, p = 0.41; Low cofiring REM mean rate = 0.73 ± 0.12, High cofiring REM mean rate = 0.78 ± 0.10, p = 0.49, WRS). (B) (Left) Changes in CA1 firing rates from the first NREM bout to the last NREM bout within a sleep epoch plotted according to degree of cofiring (cofiring as measured by the z-scored cofiring metric; high, z-cofiring > mean of distribution; low, z-cofiring < mean of distribution) with PFC either during CA1-independent NREM ripples or chained REM HFOs in PFC (Low NREM = −0.06 ± 0.03, High NREM = 0.003 ± 0.048, p = 0.40; Low REM = −0.068 ± 0.056, High REM = 0.122 ± 0.071, *p = 0.013, WRS). (C) Changes in CA1 firing rates from the first NREM bout to the last NREM bout within a sleep epoch plotted according to degree of cofiring (cofiring as measured by the Pearson correlation coefficient; high, r > 0; low, r < 0) with PFC either during CA1-independent NREM ripples or chained REM HFOs in PFC (Low NREM = −0.08 ± 0.03, High NREM = −0.01 ± 0.03, p = 0.080; Low REM = −0.10 ± 0.06, High REM = 0.08 ± 0.06, *p = 0.012, WRS). (D) Same as in (E) but using above or below the mean of the cofiring distribution to designate high or low cofiring CA1 units, respectively (Low NREM = −0.06 ± 0.03, High NREM = 0.01 ± 0.05, p = 0.30; Low REM = −0.08 ± 0.04, High REM = 0.13 ± 0.08, **p = 0.002, WRS). (E) Firing rate change separated by cofiring quartile. The change in CA1 firing rates over the course of sleep was correlated with PFC cofiring during REM HFOs but not NREM ripples (NREM, r = 0.13, p = 0.49; REM, r = 0.17, *p = 0.029, Robust linear regression). (F) Changes in PFC firing rates from the first NREM bout to the last NREM bout within a sleep epoch plotted according to degree of cofiring with CA1 during REM HFO chains in PFC. Results using 0 or the mean of the correlation coefficient distribution as the threshold shown on the left and right, respectively (0 threshold: Low = −0.35 ± 0.20, High = −0.34 ± 0.18, p = 0.90; Mean threshold: Low = −0.30 ± 0.21, High = −0.38 ± 0.21, p = 0.98, WRS). (G) Preferred CA1 theta phase of CA1 neurons during running behavior or REM sleep. A subset of CA1 cells shifted their preferred theta phase during REM sleep (U2 = 0.48, ***p = 1.51×10−5, Watson’s U2 test). Note that only cells that were significantly phase locked to the theta oscillation during both run and REM were used for analysis. (H) (Left) Distributions of preferred CA1 theta phases during run and REM sleep. Only CA1 neurons with significant phase locking during both run and REM sleep were included for analysis. (Right) Histogram of theta phase shifts between run and the subsequent REM sleep session. Red shaded areas indicate neurons that had a phase shift >90°. (I) (Left) SWR aligned firing rates of non-shifting and REM-shifting CA1 neurons. (Right) REM-shifting neurons had higher modulation indices, indicating stronger excitation during SWRs (Non-shifting = 0.60 ± 0.07, REM-shifting = 0.80 ± 0.07, ***p = 8.51×10−5,WRS) (J) The cofiring magnitude during isolated REM HFOs of non-shifting or REM-shifting CA1 neurons was not correlated with the degree of suppression during independent NREM PFC ripples (Non-shifting, r = −0.25, p = 0.085, REM-shifting, r = −0.32, p = 0.13, Pearson correlation). (K) Burst probability of non-shifting and REM-shifting CA1 neurons during isolated REM H (Non-shifting = 0.015 ± 0.002, REM-shifting = 0.021 ± 0.003, p = 0.22, WRS).

Model analyses split by cell type and reversed condition controls demonstrating the specific effects of ACh tone on network response.

(A) Firing rate histograms for NREM and REM simulations split by cell type. Bin size 0.25 Hz. NREM mean rate 1.41 Hz for pyramidal, 1.46 Hz for interneurons. REM mean rate 2.44 Hz for pyramidal, 4.14 Hz for interneurons. (B) Average network firing rate surrounding NREM MUA peaks in stimulus periods, split by cell type. (C) Same as (B), but for REM. (D) Peak stimulus amplitude vs average MUA within the stimulus period split by cell type for NREM (Left) and REM (Right). (E) Reversed stimulus condition for NREM. Same as (B), but swapping the stimulus profile by applying the REM theta-modulated stimulus profile during a low ACh (NREM) simulation. All other parameters are identical. (F) Reversed stimulus condition for REM. Same as (C), but for a simulation using the brief NREM stimulus combined with high ACh (REM). (G) Percent of each cell type with at least one spike in 10 ms bins for the swapped stimulus profile / ACh condition. Note the reduced breadth of pyramidal network activity when applying the brief stimulus in REM compared to NREM, compared to Main Figure 7B. (H) Same as Main Figure 6G, but split by cell type. (I) Swapped stimulus profile / ACh condition, applying the REM theta stimulus during a low ACh NREM simulation. Note the widespread recruitment of stimulated cells relative to non-stimulated cells, compared to REM response shown in Main Figure 7D. (J) Same as (I), but applying the NREM brief stimulus during a high ACh REM simulation. Note the lack of peak in non-stimulated cell recruitment by the stimulus compared to NREM response shown in Main Figure 7C.