1. Neuroscience

REM sleep prefrontal high-frequency oscillation chains mediate distinct cortical – hippocampal reactivation patterns compared to NREM sleep

  1. Justin D Shin
  2. Michael Satchell
  3. Paul Miller
  4. Shantanu P Jadhav  Is a corresponding author
  1. Neuroscience Program, Department of Psychology, and Volen National Center for Complex Systems, Brandeis University, United States
https://doi.org/10.7554/eLife.110795.3
Share this article
Cite this article
  1. Justin D Shin
  2. Michael Satchell
  3. Paul Miller
  4. Shantanu P Jadhav
(2026)
REM sleep prefrontal high-frequency oscillation chains mediate distinct cortical – hippocampal reactivation patterns compared to NREM sleep
eLife 15:RP110795.
https://doi.org/10.7554/eLife.110795.3
  2. Download BibTeX
  3. Download .RIS

eLife Assessment

Shin et al present significant new observations highlighting a novel oscillatory window during REM sleep that impacts neuronal dynamics and cross-regional communication and could have implications for our understanding of sleep's role in memory consolidation. This important study identifies and characterizes high-frequency oscillations (HFOs) in the prefrontal cortex (PFC) during REM sleep, reporting their temporal dynamics and coordination with hippocampal area CA1, and an intriguing dissociation between activation of CA1 neurons during REM HFOs and non-REM ripples that may impact sleep-dependent firing rate decreases. The main claims are supported with convincing evidence including an impressive range of analyses.

https://doi.org/10.7554/eLife.110795.3.sa0

Abstract

REM (rapid eye movement) and non-REM (NREM) sleep stages contribute to systems memory consolidation in hippocampal–cortical circuits. However, the physiological mechanisms underlying REM memory processes remain relatively unclear compared to NREM memory reactivation. Here we report, in rodents, the existence of prefrontal cortical (PFC) high-frequency oscillation (HFO) chains in REM sleep during the consolidation of recently acquired spatial memory. High-density tetrode recordings in hippocampal area CA1 and PFC reveal that REM cortical HFOs occur in characteristic chains that are phase-modulated by theta oscillations, corresponding to increased CA1-PFC theta coherence and delineating periods of enhanced hippocampal–cortical communication. REM HFO chains sequentially organize sparse PFC ensemble reactivation of behavioral activity during periods of local suppression, distinct from widespread reactivation bursts during NREM ripple oscillations. REM HFO chains also preferentially engage CA1 neuronal populations that demonstrate a shift in their preferred theta-phase from behavior to REM sleep. CA1 neuronal activation during REM HFO chains was correlated with CA1 activity suppression during NREM PFC ripples, and linked to differential changes in CA1 firing rates in sleep, suggesting REM-driven regulation of hippocampal excitability. A cortical network model incorporating the effects of acetylcholine can reproduce the distinct REM and NREM activity patterns, providing a mechanistic basis for widespread coactivity during NREM cortical ripples, compared to sparse, temporally extended reactivation on a background of local suppression during REM HFO chains. Overall, these findings establish a role for PFC HFOs in regulating distinct dual sleep stage reactivation patterns.

Introduction

Rapid eye movement (REM) sleep, often referred to as paradoxical sleep due to its electrophysiological similarities to wakefulness, is a relatively transient stage of sleep that occurs after NREM sleep stages in cycles and plays a crucial role in a wide range of brain functions (Blumberg et al., 2020; Peever and Fuller, 2017). Although largely appreciated for its role in emotional processing, recent evidence has highlighted its role in the consolidation of spatial, declarative, and procedural memory, indicating a broader function of REM sleep in mnemonic processing (Tempesta et al., 2018; Mukai and Yamanaka, 2023; Lendner et al., 2023). Importantly, REM dysregulation is associated with memory deficits and is thought to contribute to the progression of neuropsychiatric disorders and neurodegenerative diseases such as depression and dementia, highlighting the critical need for a deeper understanding of this sleep state to benefit human health (Dauvilliers et al., 2018).

Primary models of systems memory consolidation, including the standard two-stage theory of consolidation or trace transformation theory (Squire and Alvarez, 1995; Nadel and Moscovitch, 1997; Buzsáki, 1989; Sekeres et al., 2018), rest under the assumption of cooperative hippocampal–cortical reactivation during sleep to transform initially encoded episodic hippocampal representations into broadly distributed hippocampal–cortical representations for long-term storage and semantic memory formation. In support of this view, several rodent and human studies have confirmed that, in non-rapid eye movement (NREM) sleep, coordination of hippocampal sharp-wave ripples (SWRs) with cortical slow oscillations (1–4 Hz), spindles (12–16 Hz), and ripples (150–250 Hz) underlies the neurophysiological mechanism supporting consolidation (Staresina, 2024; Buzsáki, 2015; Girardeau and Zugaro, 2011; Klinzing et al., 2019). While the role of offline memory reactivation in the hippocampus and prefrontal cortex (PFC) during NREM sleep has been extensively studied in relation to memory function, the precise mechanisms through which REM sleep contributes to these processes, particularly in the context of hippocampal–cortical interactions, are still unclear.

REM sleep is characterized by a high theta-to-delta (TD) ratio (Mizuseki et al., 2011; Tang et al., 2017), and while hippocampal SWRs are primarily reported in NREM states with relatively low theta power, hippocampal theta rhythms during REM sleep have been shown to be important for memory consolidation (Boyce et al., 2016). Correspondingly, phasic theta activity bursts during REM are often coupled with elevated gamma and cortical high-frequency oscillation (HFO) power in hippocampal and cortical networks, potentially organizing neural activity in distributed circuits for mnemonic processing (Montgomery et al., 2008; Brankačk et al., 2012; Bandarabadi et al., 2019; Bueno-Junior et al., 2023; Tort et al., 2013; Ghosh et al., 2022; Arndt et al., 2024; Aleman-Zapata et al., 2022). Furthermore, REM sleep is associated with a unique neuromodulatory landscape that is conducive to the induction of plasticity within and across circuits (Yamada and Ueda, 2019). However, the role of REM sleep cortical oscillations, whether they are coordinated with hippocampal activity patterns or with cortical–hippocampal reactivation, and the relationship between REM and NREM activity for dual sleep state regulation of reactivation and excitability are all open questions.

Cortical ripples (high-frequency or fast oscillations, ~150–250 Hz) have been reported in NREM sleep (Khodagholy et al., 2017; Shin and Jadhav, 2024; Vaz et al., 2019), and our recent study showed that independent PFC ripples that are dissociated from hippocampal SWRs are prevalent in NREM sleep and mediate suppression of hippocampal activity and reactivation, potentially facilitating the selection and consolidation of specific hippocampal memory ensembles, while facilitating local cortical reactivation (Shin and Jadhav, 2024). A related study had similar findings, suggesting that these cortical ripples reflect intracortical processing that occurs independently of hippocampal input, allowing for interference-free consolidation of memory representations within cortical circuits (van Schalkwijk et al., 2023).

In REM sleep, previous studies have established cortical HFOs (in a similar frequency range to high-frequency NREM ripples) along with gamma oscillations – REM HFOs have been detected in different frequency ranges and sometimes reported with different nomenclature. (Brankačk et al., 2012; Bueno-Junior et al., 2023; Tort et al., 2013; Scheffzük et al., 2011) Since recent studies note the possibility of spurious spectral detection in high frequency ranges (van Schalkwijk and Helfrich, 2026), monitoring spiking activity simultaneously with LFP/EEG activity can therefore provide key confirmation about the existence and impact of HFOs on neuronal activity via spiking modulation in local circuits.

Previous studies in rodents have reported a role of REM sleep in homeostatic regulation of firing rates, using recordings lasting from less than an hour to across the 24 hr cycle, but without the context of a behavioral task (Grosmark et al., 2012; Miyawaki and Diba, 2016; Watson et al., 2016). Studies have also reported prominent phenomena such as hippocampal REM-theta phase reversing cells (Poe et al., 2000) and hippocampal reactivation (Louie and Wilson, 2001), using behavioral tasks and examining brief sleep periods in post-task rest sessions. Further, a recent study in mice reported PFC coding and activation of recent experiences and inferred knowledge in NREM and REM sleep using imaging methods, but did not examine any potential role of oscillations (Abdou et al., 2024).

Despite these results, an explicit examination of REM cortical HFO events, and how they impact cortical and hippocampal activity during consolidation of recent experiences is not known, despite the similarity of these HFOs to cortical ripples in NREM. Furthermore, how NREM and REM sleep activity patterns may together complementarily shape the development and refinement of cortical–hippocampal ensembles to support learning is unclear. Therefore, we used high-density electrophysiological recordings in PFC and CA1 during learning and subsequent sleep to investigate whether and how coordinated activity during REM sleep may drive circuit changes underlying systems memory consolidation, and whether prefrontal and hippocampal dynamics differ during cortical ripple and HFO activity patterns in NREM and REM sleep, respectively.

Results

High-frequency events in PFC during REM sleep

We simultaneously recorded local field potential (LFP) and neuronal populations in CA1 (n=1468) and PFC (n=1151) from animals (n=10 rats) during REM sleep (Figure 1—figure supplement 1A, Supplementary file 1, REM duration = 164.6 ± 10.8 s/epoch for the 36 epochs included in analysis) throughout the course of learning a spatial memory task that requires hippocampal–prefrontal interactions (Maharjan et al., 2018; Shin et al., 2019; Figure 1—figure supplement 1B). Activity was continuously acquired across multiple run (8 epochs) and sleep (9 epochs) sessions during learning within a single day (Shin et al., 2019). We separated NREM and REM sleep stages based on TD ratio in CA1 (Mizuseki et al., 2011; Tang et al., 2017), which revealed coherent shifts in TD ratio across CA1 and PFC at the onset and offset of REM sleep (Figure 1A–C). Waking periods showed higher movement and intracranial EMG compared to REM sleep, validating the behavioral state designations (Figure 1—figure supplement 1C–I).

Figure 1 with 1 supplement see all
Download asset Open asset
High-frequency oscillations (HFOs) in PFC during REM sleep.

(A) CA1 and PFC theta-to-delta (TD) ratio time locked to either the start or end of a detected 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) An example sleep session for an animal showing each behavioral state detected using the sleep state algorithm that utilizes TD ratio and animal velocity (cm/s). (C) Average spectrograms across all REM bouts centered on the start (left) and end (right) of each bout. Note the sharp increase and decrease in theta power at the start and end of REM sleep, respectively. (D) 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. (E) (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. (F) Example PFC event triggered spectrogram for REM HFOs and NREM 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. (G) PFC high-frequency event rates in REM (HFOs) and NREM (ripples) sleep (REM = 0.66 ± 0.06, NREM = 0.54 ± 0.05, *p=0.013, WSR). (H) PFC HFO autocorrelation in REM sleep. Note the peaks in the autocorrelation that repeat approximately every 125 ms. (I) Distribution of inter-event intervals (IEIs) for NREM ripples and REM HFOs. 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.

To investigate the high-frequency component of phasic burst activity in PFC, cortical ripples during NREM sleep and HFOs during REM sleep were detected as transient events in the 150–250 Hz range, as previously described (Shin and Jadhav, 2024; Figure 1D and E and Figure 2—figure supplement 1A and B). Similar criteria were used to detect cortical high-frequency events in NREM and REM states; however, to avoid ambiguity and to conform to previous nomenclature (Bueno-Junior et al., 2023; Tort et al., 2013), we refer to cortical NREM events as ripples, cortical REM events as HFOs, and hippocampal SWRs in NREM as SWRs throughout. Cortical event rates during REM were higher than NREM, and REM HFOs were associated with elevated theta and gamma power (Bueno-Junior et al., 2023; Scheffzük et al., 2011), in contrast to spindle- and slow oscillation-coupled ripple events in NREM (Shin and Jadhav, 2024; Staresina et al., 2023; Figure 1F and G). Interestingly, examining the temporal structure of event occurrence, using both the event autocorrelation (the likelihood of finding another HFO at a given time lag) and the distribution of inter-event-intervals (IEIs, the time between consecutive HFOs), we find that REM HFOs often recur at intervals of ~130 ms (Figure 1H and I) as HFO ‘chains’, which is consistent with theta frequency modulation and elevated theta power during these events (Figure 1F).

REM HFOs are associated with fast, theta timescale fluctuations in PFC population/multiunit activity (MUA) during transient periods of locally reduced activity, as evidenced by the overall decrease in population activity during these events (Figure 2A). Furthermore, to isolate the spectral component driving this PFC population response, we re-detected HFOs using a broader 100–250 Hz passband and separated them by co-occurrence with the original 150–250 Hz HFOs (coordinated vs. non-coordinated). HFOs that occurred independently of 150–250 Hz detections elicited little to no phasic PFC spiking modulation, whereas those temporally coordinated with 150–250 Hz events reproduced the phasic theta-modulated response, indicating that this phasic PFC response is selectively associated with HFOs containing prominent 150–250 Hz activity (Figure 2B). This phenomenon was also present when using different parameters for HFO detection or REM sleep designation (Figure 2—figure supplement 1D–F). In contrast, NREM cortical ripples were associated with a singular burst of activity spanning several hundred milliseconds (Shin and Jadhav, 2024, Figure 2C and Figure 2—figure supplement 1C). Thus, these findings establish distinct high-frequency event-associated activity profiles in PFC that may underlie differential processing of memory-related information during NREM and REM sleep.

Figure 2 with 1 supplement see all
Download asset Open asset
HFOs in PFC during REM sleep are associated with theta-modulated population activity and transient suppression.

(A) (Left) REM HFO aligned multiunit activity (MUA) in PFC (n=10 animals, 36 epochs), 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. (B) PFC MUA 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 drive the theta frequency MUA and transient suppression in PFC. (C) Same as in (A) but for NREM ripples. Note the synchronous activity burst time-locked to NREM ripples, in contrast to the distinctive activity pattern during REM HFOs in (A). Top and bottom PSDs are from the same animals as in (A).

Coupling of theta, gamma, and prefrontal HFOs in REM sleep

Theta-gamma phase-amplitude coupling (PAC) is one of the most widely studied cross-frequency coupled phenomena in the brain and involves the coupling of distributed, lower frequency theta oscillations and local higher frequency bursts (Canolty and Knight, 2010). Thus, cross-frequency PAC is considered to be a physiological marker that links distributed brain regions during bouts of coherent activity to integrate information across different spatiotemporal scales (Canolty and Knight, 2010). Furthermore, previous work has proposed that phasic REM sleep may support cortical–hippocampal dialogue through this theta frequency-based coordination (de Almeida-Filho et al., 2021).

Phasic REM is characterized by higher theta power and frequency as compared to tonic REM in both the neocortex and hippocampus (Montgomery et al., 2008; Brankačk et al., 2012). Thus, we detected putative bouts of phasic REM as periods of elevated theta power (Figure 3—figure supplement 1). We further quantified the inter-peak-intervals of the filtered CA1 theta signal during low and high theta power bouts to determine whether they may correspond to putative tonic and phasic substages, respectively. Inter-peak intervals were significantly shorter during high theta power bouts than low power bouts, consistent with these periods reflecting putative phasic REM (Figure 3—figure supplement 1A–C). The rate of REM HFOs, as well as chaining of events, was significantly elevated during periods of high theta power (Figure 3—figure supplement 1D–G), indicating that the prevalence of REM HFO chains is elevated during bouts of putative phasic REM sleep. Consistent with this, we found that theta cycles were more likely to contain embedded HFOs during putative phasic REM than during tonic REM (Figure 3A). Thus, we restricted further analyses for quantifying gamma and HFO coupling to theta oscillations to these phasic bouts of high theta power (Figure 3B, C, Figure 3—figure supplement 2A-F). To further validate the coupling of theta, gamma, and HFOs (illustrated in Figure 1F), we quantified PAC, which is a measure of how strongly the amplitude of fast oscillations varies as a function of the phase of a slower one, between theta phase and gamma to HFO frequency (40–250 Hz) amplitudes (Figure 3B). Examining the modulation indices at two discrete frequency bands, gamma (40–100 Hz) and HFO (150–250 Hz), revealed a higher frequency, HFO component that is also modulated by theta phase (Figure 3—figure supplement 2A–F). The modulation index summarizes this coupling as a single number, with higher values indicating tighter phase-locked amplitude modulation.

Figure 3 with 2 supplements see all
Download asset Open asset
Theta, gamma, and HFO coupling underlies distinct REM modulation states in PFC.

(A) Proportion of theta cycles with HFOs across all, low theta power (putative tonic), and high theta power (putative phasic) cycles (all = 0.075 ± 0.007, low = 0.065 ± 0.006, high = 0.140 ± 0.018, ***p=4.90 × 10–5, WSR). (B) Phase-amplitude coupling (PAC). 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. (C) Theta nesting for HFO and gamma 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 unfiltered LFP to the peak power in the preferred theta frequency phase bin for either gamma or HFO oscillations separately. (D) 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). (E) CA1-PFC theta (6–12 Hz) coherence centered on PFC REM HFOs is elevated relative to baseline theta coherence. The overlaid curve on the spectral plot represents the mean z-scored coherence across epochs, with shading indicating ± SEM. (F) Example REM HFO chain (doublet) plotted across four separate tetrodes with the average HFO power plotted below. (G) HFO chains showed higher PFC (left) and CA1 (right) theta PPC (pairwise phase consistency, a measure of phase-locking strength) 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). (H) (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). (I) (Left) CA1-PFC theta (6–12 Hz) coherence during isolated vs chained events and (right) quantification of the peak theta coherence in a±1 s window from HFO onset (chain = 0.93 ± 0.08, isolated = 0.65 ± 0.05, **p=0.0069, WRS). (J) (Left) Phase slope index (PSI), which is a measure of phase lag consistency across different frequencies, 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).

To further verify that gamma and HFO are separable rhythms, we examined theta nesting of both oscillations, which quantifies the oscillatory components of the average raw LFP waveform centered on peaks of gamma or HFO power within the preferred theta phase bin (see Methods). This revealed two distinct fast oscillatory components embedded within the theta cycle that are tightly coupled, with gamma events leading HFOs (Figure 3C and Figure 3—figure supplement 2C, D). Since the coupling metrics used here quantify relationships across theta cycles in aggregate rather than on a cycle-by-cycle basis, analyses of theta coupling for gamma and HFO activity were restricted to periods of high theta power (putative phasic REM), as mentioned above. Including periods of relatively low HFO occurrence could obscure genuine coupling dynamics. Focusing on epochs with robust theta oscillations, where HFOs are more prevalent, therefore provides a stronger basis for distinguishing gamma from HFOs and yields a more faithful representation of theta, gamma, and HFO cross-frequency coupling.

Furthermore, the phase-locking preferences of PFC neurons to these two oscillatory components differ (Figure 3D, Figure 3—figure supplement 2G), confirming gamma and HFO separation, and potentially reflecting distinct modes of local ensemble memory related processing. Importantly, we observed elevated CA1-PFC theta coherence during these REM HFOs as compared to baseline (Figure 3E), suggesting that increased theta coherence is associated with HFO generation in PFC during REM sleep.

Since we observed short-latency HFO recurrence during REM (Figure 1H and I), we reasoned that these HFO ‘chains’ may be functionally different than HFOs that occur in isolation. Indeed, previous studies have characterized oscillatory events in multiple brain areas that occur in temporal clusters (Davidson et al., 2009; Darevsky et al., 2024; Antony et al., 2018; Mallory et al., 2025). Importantly, chains of spindle events (termed, ‘spindle trains’) in NREM sleep have been shown to underlie persistent reactivation of awake patterns and are associated with precise spatiotemporal coupling of sleep oscillations (Darevsky et al., 2024). Similar to this previous study, we find spindle train events in PFC during NREM sleep (Figure 4—figure supplement 1A–D), supporting the claim that PFC can exhibit temporally clustered events during sleep.

We identified chains of HFOs as events that occur within 200 ms of each other (37.9% of all events). Chains were more prevalent in REM than NREM sleep (Figure 4—figure supplement 1E) and in putative phasic REM than tonic REM (Figure 3—figure supplement 1E and F). Similar to spindle trains, REM HFO chains predominantly occurred in doublets (Figure 3F and Figure 4—figure supplement 1C and G). REM HFOs were more consistently locked to a preferred theta phase than NREM ripples, as quantified by pairwise phase consistency (PPC), which is a bias-free measure of phase-locking strength (Figure 4—figure supplement 1F). Furthermore, REM HFO chains were strongly phase modulated by theta in both PFC and CA1 and had similar properties to isolated HFOs (Figure 3G, Figure 4—figure supplement 1J). The strength of gamma was also elevated during HFO chains and exhibited phasic power fluctuations as well as suppression of population activity (as in Figure 2A) when coupled (within 500 ms) with HFOs, further demonstrating strong temporal coupling (Figure 3H and Figure 3—figure supplement 2I and J).

Corroborating the strong phase locking of chain HFOs to theta oscillations in PFC and CA1 (Figure 3G), we find that hippocampal–prefrontal theta coherence is elevated during chain HFOs as compared to isolated events (Figure 3I), with CA1 theta leading PFC theta, as indicated by the phase slope index (PSI; negative values denote a CA1 to PFC direction of information flow at theta frequencies; Figure 3J). Together, these findings suggest that chains of HFOs can underlie coordinated hippocampal–prefrontal activity that supports REM sleep-mediated memory consolidation.

Differential modulation of PFC neuronal activity during NREM prefrontal ripples and REM prefrontal HFOs

Since we observed theta-modulated population activity aligned to REM HFOs that differed from activity aligned to NREM ripples (Figure 2A–C), we next examined PFC responses to determine how activity surrounding these events are organized at the level of single units and ensembles. Similar to population responses, many PFC neurons exhibited oscillatory spiking aligned to REM HFOs (Figure 4A). Unlike NREM ripples, REM events were associated with a broader distribution of activity peaks, suggesting a more sequential structure of neuronal firing surrounding these events (Figure 4A and B). Furthermore, we trained a binary linear classifier on the spike count vector across PFC neurons during each event and tested its ability to label held-out events as NREM ripples versus REM HFOs (10-fold cross-validation, with event counts equalized across types). PFC neuronal activity reliably distinguished NREM and REM events, confirming differences in spiking content (Figure 4C). In line with elevated CA1-PFC theta coherence surrounding REM HFOs (Figure 3E), we found that a subset of CA1 neurons were modulated by these events (Figure 4D), indicating enhanced interregional coactivity, potentially delineating periods of CA1-PFC ensemble reactivation to support memory consolidation.

Figure 4 with 2 supplements see all
Download asset Open asset
Differential modulation of PFC neuronal activity during NREM prefrontal ripples and REM prefrontal HFOs.

(A) (Left) PFC neurons that were positively modulated (EXC) during REM PFC HFO events and NREM PFC 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. Note the peaks in activity adjacent to REM HFO alignment at theta timescale denoted by the red arrowheads. (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 (MUA) 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) PFC MUA aligned to the first HFO in each chained event. At top is a schematic of the multiunit alignment procedure. The initial HFO in a chain is at time 0. Note the peaks in the MUA prior to the initial HFO within a chain, suggesting that fluctuations in PFC MUA precede chained events. (H) The probability of (top) isolated and (bottom) 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 MUA.

Overall cortical cofiring, which is a measure of coincident activity between neuron pairs during discreet events (Cheng and Frank, 2008), and assembly strength were reduced during REM HFOs relative to NREM ripples (Figure 4E), consistent with overall activity suppression during HFO chains (Figure 4F–H), which can create a permissive background for gating selective reactivation of PFC ensembles in REM sleep. Examining population spiking activity, we found that PFC activity was strongly theta-modulated during REM HFO chains on the background of transient local suppression, compared to both isolated and NREM events (Figure 4F-H, Figure 4—figure supplement 1H-J). This phasic population activity pattern was seen during HFO chains in both putative phasic and tonic REM substates (Figure 3—figure supplement 1F).

To address the possibility that this theta-modulated spiking activity and suppression may simply reflect high-theta-power bouts (rather than HFO occurrence), we generated surrogate event times by randomly sampling time points that matched the preferred theta phase of HFOs, but were at least 250 ms away from any detected HFO. PFC activity aligned to these surrogate events did not exhibit the same theta-modulated activity associated with real HFOs, even when the surrogate events were drawn from the highest-theta-power quartile (Figure 4—figure supplement 1K). Furthermore, we examined population activity as a function of distance from detected HFOs across all REM bouts, or subset by putative phasic REM, and did not observe the same magnitude of phasic response (Figure 4—figure supplement 1L–N), confirming that the transient theta-modulated population activity in PFC is specifically associated with the occurrence of HFOs. Interestingly, only HFO chains were uniquely associated with a transient suppression in overall local population activity (Figure 4F-H, Figure 4—figure supplement 2), which is in contrast to the strong excitation seen in PFC and CA1 during NREM PFC ripples and SWRs, respectively (Shin and Jadhav, 2024). Since REM HFOs predominantly occur in doublets (Figure 4—figure supplement 1G), the suppression of overall population activity is consistent with sublinear integration, in which subsequent events elicit less spiking activity than the first, thus likely creating privileged temporal windows for selective PFC ensemble reactivation and strengthening.

PFC ensemble activity is sequentially organized by REM HFO chains

Since we observed a broader temporal distribution of spiking activity (Figure 4B) as well as multiple peaks in population activity during REM HFO chains (Figure 4F and G), we next examined, in more detail, how ensemble activity is structured around these events. Previous studies have shown that the refractoriness of certain oscillatory events, such as spindles and SWRs, create distinct segregated periods for selective population activation (Davidson et al., 2009; Darevsky et al., 2024; Mallory et al., 2025; Judák et al., 2022). Thus, we first examined population dynamics based on the temporal separation of REM HFOs. We observed selective activation of subsets of PFC ensembles during individual events in HFO chains (Figure 5A). To quantify this structure, we represented each HFO as a binary vector of PFC neurons active during the event and computed the Pearson correlation between the vectors of every pair of consecutive HFOs. Spike content similarity as a function of IEI was examined using quartile ranges for the intervals between HFOs. Quartile 1 had the highest correlation, with the average IEI of this quartile as 114 ms (theta frequency, corresponding to HFO chains), indicating that adjacent HFO events within chains show a high degree of similarity as compared to adjacent isolated events (Figure 5A and B and Figure 6—figure supplement 1A). Furthermore, we examined whether the order in which individual PFC neurons fired during HFOs within a chain was preserved across chains. For each chain, we extracted each cell’s first-spike rank order and compared it to a leave-one-out template constructed from the average normalized rank across all other chains. This procedure revealed preserved sequences of PFC activity across multiple chain events (Figure 5C, Figure 6—figure supplement 1B). These findings suggest that HFO chains stereotypically organize PFC activity, likely to support sequential reactivation of distinct assemblies.

Figure 5
Download asset Open asset
Temporally extended sequential activity in PFC.

(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) Similarity of PFC population spiking between adjacent REM HFOs and its relationship with inter-event-intervals (IEI). Here, the PFC spiking activity during each HFO was binarized across all neurons, and the Pearson correlation coefficient was calculated between adjacent events as a measure of pattern similarity. Then the relationship between pattern similarity and IEI was reported (r = –0.13, ***p=6.00 × 10–8, Pearson correlation). Below each quartile is the average IEI of that quartile in milliseconds. Note that the average for quartile 1 is 114 ms, 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 from the leave-one-out cross-validation procedure. Each event’s rank was correlated with the averaged rank across all other events. The average across all events compared to a distribution of means generated by jittering (n=1000) spike times is shown.

Next, to test our hypothesis of sequential reactivation, we assessed PFC assemblies by examining time-locking of PFC reactivation strength to NREM and REM events. We detected behavioral PFC assemblies that are reactivated in sleep using established methods (Peyrache et al., 2009; van de Ven et al., 2016). Similar to single-unit modulation, the timing of PFC assembly reactivation aligned to NREM ripples was largely restricted to within 200–500 ms after ripple onset, depending on the ripple type examined (Figure 6A and Figure 6—figure supplement 1F). For REM events, only the first HFO within a chain was selected for alignment since we found that peaks in population activity occur prior to the onset of HFO chains (Figure 4G), possibly facilitated by strong phasic gamma surrounding HFO chains (Figure 3H, Figure 3—figure supplement 2I). We found that unlike NREM ripples, where assembly reactivation is synchronous and consistently highest at the onset of events, PFC assemblies were sequentially organized on longer timescales surrounding REM HFOs, similar to single-unit modulation (Figure 4A and B and Figure 4—figure supplement 2A). This temporal tiling of assembly reactivation aligned to REM HFOs manifested in a relatively flat average response, highlighting the differential engagement of PFC assemblies during NREM and REM events (Figure 6B).

Figure 6 with 1 supplement see all
Download asset Open asset
Structured reactivation in PFC during REM HFOs.

(A) (Left) PFC assembly reactivation aligned to NREM ripples or (Right) the first HFO of each REM HFO chain. Reactivation peak locations are shown alongside each plot. The assemblies are sorted based on the average peak reactivation bin within a±500ms 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 (B). (B) (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. (C) Cross-validation of the REM reactivation sequences shown in (A), 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). (D) Differences in slope of the linear fit between the time of peak reactivation (in ms) 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). (E) 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. (F) The distribution of PFC assembly reactivation strength differences between pre and post sleep. The same W-Track session was used as the template for each pair of sleep epochs. Thus, assemblies were matched across pre and post sleep epochs (**p=0.006, t-test against zero). (G) (Left) PFC assembly reactivation strength during sleep over the course of learning. (Right) Comparison of reactivation strength during early (sleep sessions 1–3; epochs prior to grouped average of 75% correct performance on W-Track, see Figure 1—figure supplement 1B) and late (sleep sessions 4–9) sleep epochs (early = 0.475 ± 0.010, late = 0.525 ± 0.007, ***p=1.78 × 10–4, WRS). Except for sleep epoch 1, the preceding W-Track session was used as the template for reactivation. Only neurons that were tracked throughout the entire experiment were included for analysis. (H) Example PFC assembly rate maps. (I) (Left) Distribution of z-scored spatial information values for real assembly maps compared to maps generated from circularly shuffled (n=1000) activation times. Each assembly’s spatial information was expressed as a z-score relative to its shuffled distribution. The red-dotted line at a z=1.65 marks the one-tailed 95th percentile (p=0.05). (Right) Observed spatial information vs. the median of each assembly’s shuffle null distribution (data = 0.535 ± 0.013, shuffle = 0.325 ± 0.005, ***p=2.82 × 10–65, WSR).

To confirm that the sequential tiling of assembly reactivation peaks in Figure 6A was a reliable phenomenon associated with HFO chains rather than an artifact, we used a split-half cross-validation. HFO events were first randomly partitioned into two halves, each assembly’s peak reactivation bin was computed independently within each half, and the similarity of the temporal ordering of assembly peaks between the two halves was calculated by taking the Pearson correlation (see Methods). Consistent with the structured PFC activity during chain events (Figure 5A–C), we found the greatest degree of sequence similarity during REM HFO chains compared to shuffled data (Figure 6C), confirming the tendency for PFC assemblies to be reactivated in a stereotyped, sequential manner as compared to isolated REM HFOs and NREM ripples (Figure 6—figure supplement 1C–F).

Finally, consistent with task-related reactivation, we found stronger PFC ensemble expression during sleep epochs following experience (Figure 6F). Importantly, assembly reactivation in sleep strengthened over the course of learning (Figure 6G, Figure 1—figure supplement 1B), and a large proportion of detected assemblies had behaviorally relevant task representations with high spatial information (Figure 6H–I).

Differential engagement of distinct populations of CA1 neurons by PFC REM HFO chains

Theta coherence and PAC are two prominent modes of interregional coordination thought to underlie mnemonic processing (Canolty and Knight, 2010; Fries, 2005). Since we observed an overall increase in CA1-PFC theta coherence and gamma power during HFO chains (Figure 3H and I), as well as sparse modulation of CA1 activity during REM HFOs (Figure 4D), we reasoned that there may be specific populations of CA1 neurons that are engaged by these events. Examining CA1-CA1 cofiring during these cortical HFO events, we observed a higher degree of spatial rate map correlation for high cofiring pairs, suggesting HFO specific reactivation of behaviorally relevant hippocampal assemblies to support memory consolidation (Buzsáki, 2015; Figure 7A and B). Additionally, the absolute CA1 cofiring, comprising both negative and positive magnitudes of cofiring, was higher during HFO chains, indicating a distinct activity profile compared to isolated HFOs (Figure 7—figure supplement 1A).

Figure 7 with 2 supplements see all
Download asset Open asset
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 (p = 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). For comparisons between cofiring and other metrics (e.g. firing rate), a single cofiring value was calculated for each neuron by averaging the cofiring metric across all neuron pairings. Additionally, high and low cofiring CA1 neurons were split based on average cofiring values >0 and<0, respectively. (E) (Left) High REM PFC HFO cofiring CA1 neurons exhibited a greater degree of suppression during NREM PFC ripples. For a description of modulation index, see Methods section Ripple/HFO aligned modulation. Here, since CA1 neurons exhibit a robust decrease in activity in response to NREM PFC ripples (Shin and Jadhav, 2024), we refer to the modulation as suppressive (Low cofiring = –0.08 ± 0.06, High cofiring = –0.25 ± 0.03, *p=0.022, WRS). (Right) The CA1 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).

Next, we computed CA1-PFC cofiring across all cross-regional cell pairs and found that, during REM HFO chains, the distribution of cofiring values was bimodal, revealing a subpopulation of CA1 cells that are highly engaged with PFC neurons (Figure 7C). This distribution was unimodal during NREM PFC ripples, suggesting that the dual-population structure is specific to REM HFO chains. Analysis of low and high cofiring CA1 neurons during REM HFOs showed that high cofiring neurons exhibited elevated activity during HFO chain events as compared to low cofiring neurons (Figure 7D), independent of baseline firing rates (Figure 7—figure supplement 2A).

Interestingly, we also found a relationship between CA1 activation during REM PFC HFO chains and CA1 suppression by NREM PFC ripples. CA1 neurons that exhibited high cofiring with PFC (termed, ‘high HFO cofiring CA1 neurons’) specifically during REM HFO chains were more strongly suppressed during independent NREM PFC ripples, and the degree of suppression was correlated with the strength of REM CA1-PFC cofiring (Figure 7E and Figure 7—figure supplement 1C and D). Given the sequential occurrence of NREM-REM sleep stages, this finding suggests that NREM ripples may tag specific CA1 populations for subsequent REM-dependent reorganization. However, although there was a strong relationship between CA1 suppression during independent PFC NREM ripples and excitation during coordinated SWRs in NREM sleep, as shown in a previous study (Shin and Jadhav, 2024; Figure 7—figure supplement 1C and E), we did not find correlations between CA1 REM HFO cofiring and CA1 SWR reactivation (Figure 7—figure supplement 1F), indicating that this CA1 subpopulation that we identified as high cofiring with PFC neurons during REM HFOs is specifically engaged by high-frequency PFC events in sleep, with increased activity during REM and suppression during NREM.

Previous studies have shown that NREM and REM sleep are involved in modifying hippocampal excitability, with a net decrease in firing rates over the course of sleep (Grosmark et al., 2012; Miyawaki and Diba, 2016; Watson et al., 2016). To test whether REM HFO engagement predicted shifts in hippocampal excitability across sleep, we computed each CA1 cell’s firing rate change from the first to the last NREM bout within a sleep epoch and asked whether this change was related to how strongly the cell co-fired with PFC during HFOs. Interestingly, we observed that high REM HFO cofiring CA1 neurons differentially shifted their firing rates upward over the course of a sleep session, whereas most of the neuronal population, including low REM HFO cofiring cells, decreased their firing rates (Figure 7F, Figure 7—figure supplement 1B, and Figure 7—figure supplement 2B–E), consistent with the role of REM sleep in hippocampal firing rate diversification (Miyawaki et al., 2019). This effect was hippocampal specific; there was no differential regulation of PFC firing rates over the course of sleep (Figure 7—figure supplement 2G). Furthermore, the degree of cofiring with PFC during REM HFOs predicted overall CA1 firing rate changes during sleep, and this relationship was absent for NREM PFC ripples (Figure 7F and Figure 7—figure supplement 2E). Splitting the population in half based on mean firing rate changes further confirmed this relationship (Figure 7—figure supplement 2F), suggesting a REM-specific role of PFC HFOs in modifying hippocampal excitability.

REM theta-phase shifting CA1 neurons are preferentially engaged during PFC REM HFOs

The dorsal CA1 pyramidal cell layer can be separated into subpopulations that are molecularly defined, have differential engagement with hippocampal SWRs and cortical areas (Valero et al., 2015; Gu et al., 2023; de la Prida, 2020; Harvey et al., 2023), and employ complementary spatial codes for effective coding of environments (Sharif et al., 2021). Notably, a subset of CA1 pyramidal cells preferentially shift their preferred theta phase from run to REM sleep (Mizuseki et al., 2011) – a phenomenon that is experience dependent (Poe et al., 2000). We therefore investigated whether there were differences in REM HFO recruitment for non-shifting and REM-shifting cells. We categorized REM-shifting neurons as cells that had a difference in phase preference greater than 90° (Mizuseki et al., 2011). As expected, there was a subset of CA1 neurons that shifted theta phase preference (n=124/355 neurons significantly locked to both run and REM theta, Shift magnitude = 133.14 ± 2.35°) and were strongly engaged by SWRs (Figure 7—figure supplement 2H–J), as previously reported (Gu et al., 2023). Similar to high cofiring CA1 neurons, these REM-shifting cells had higher cofiring with PFC than non-shifting cells (Figure 8A), and there was a significant relationship between HFO chain cofiring and suppression in the REM-shifting neurons, suggesting a specific modulatory effect on these neurons (Figure 8B and Figure 7—figure supplement 2K). REM-shifting neurons also exhibited strong bursting activity during HFO chains as compared to isolated events and were more strongly phase locked to the PFC theta oscillation surrounding HFO chains (Figure 8C and D and Figure 7—figure supplement 2L). Thus, PFC HFOs during REM sleep preferentially engage a specific population of REM-phase-shifting CA1 neurons, which are suppressed during NREM cortical ripples, and overall show an increase in firing rates during sleep, thus potentially increasing or stabilizing excitability to enhance signal-to-noise for efficient expression during future behavior.

Figure 8
Download asset Open asset
REM theta phase shifting CA1 neurons are preferentially active during PFC REM HFO chains.

(A) (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). Note that all cofiring pairs are assessed here. (B) Similar to Figure 7E, the REM PFC HFO 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). (C) 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). (D) 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.

Modeling high-frequency event-associated PFC activity during NREM and REM sleep

To gain a better understanding of how a cortical network can exhibit different activity patterns during NREM ripples and REM HFOs, we built a simplified network model (Figure 9A). 1600 excitatory pyramidal and 400 inhibitory interneurons were pseudo-randomly connected based on cortical layer 2/3 connectivity probabilities (see Methods). Acetylcholine (ACh) concentration is an especially important factor for differentiating NREM and REM sleep states, with strong impacts on cell-spiking properties and network dynamics (Picciotto et al., 2012; Eban-Rothschild and de Lecea, 2017; Colangelo et al., 2019). To incorporate one of the crucial impacts of ACh concentration on cortical neurons, each neuron in the model used a modified Hodgkin-Huxley (HH) formalism with an additional slow potassium (M-type) leak conductance term (Stiefel et al., 2009). During NREM simulations, this slow potassium channel is opened to simulate low ACh concentration, and the channel is closed for REM simulations to simulate high ACh concentration (Roach et al., 2015; Satchell et al., 2025). By incorporating this additional conductance, a primary effect of the most abundant muscarinic ACh receptor in the brain (the M1 receptor) is included in the model. In vivo and in the model, activation of this potassium leak current alone reduces cell excitability and increases network synchronicity (Czarnecki et al., 2021; Stiefel et al., 2008; Gulledge et al., 2009; Dasari and Gulledge, 2011; Mishra et al., 2024). In the model, modulation of this M-type current also changes the excitatory-inhibitory balance of the network, supporting higher average interneuron rates (Figure 9—figure supplement 1A).

Figure 9 with 2 supplements see all
Download asset Open asset
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 s 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 ± SD), 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.

While the model does not possess the complexity to replicate the same HFO-range spiking frequencies, we found that by applying a stimulus input to a randomly chosen subset of cells of the network, we could recreate some essential characteristics of network activity patterns during NREM and REM cortical events. A brief moderately strong input during NREM, representing input during NREM ripples, was sufficient to induce a larger burst of synchronous spiking in the network, causing a large peak in MUA (Figure 9B and D, left). In order to recreate the MUA surrounding REM PFC HFOs, we found that a theta-modulated Gaussian input with a similar profile to the theta power observed around PFC REM HFOs (representing activity during phasic theta bursts in REM) was able to capture oscillatory and inhibitory features of the MUA (Figure 9C and D). We applied stimulus inputs of a variety of peak amplitudes to the network, and found that for temporally brief, discrete stimuli in NREM, stronger stimulation induced greater MUA (Figure 9—figure supplement 1D). Conversely, the strength of theta-oscillating input (and associated theta power as measured in the model’s MUA) during REM correlated negatively with the MUA level during stimulation (Figure 9E). We found a similar negative relation in the experimental data between theta power surrounding REM HFO times and MUA level, suggesting that increased theta power in PFC is associated with stronger external input.

The cofiring analysis done on experimental spiking data during REM and NREM events (Figure 4E) was implemented for the spikes of a matched number of randomly chosen model neurons during REM and NREM stimulation periods, showing greater cofiring during NREM similar to the experimental data (Figure 9F). In addition, we divided peaks in MUA during REM stimulation periods in the same manner used to separate isolated and chained REM HFOs in the data (minimum 200ms between MUA peaks, Figure 3). Similar to the data, we observed greater theta modulation in the MUA surrounding chained peaks compared to isolated peaks, as well as a larger dip in overall MUA (Figure 9G), with the difference especially clear in the pyramidal population of the model (Figure 9—figure supplement 1H). Overall, these results indicate an ability of the model to capture essential features of PFC activity in NREM and REM sleep, and suggest that REM HFOs reflect transient periods of elevated theta input to PFC neurons, leading to characteristic activity patterns in the network due to high cholinergic tone.

Stimulation in NREM induces widespread network activation compared to stimulation in REM

Analyzing stimulation periods in the model’s NREM and REM states revealed intriguing differences in network activity (Figure 10A). While stimulus periods in REM tended to induce only a slight increase in the total fraction of cells spiking at any given time, stimulation in NREM induced a much wider spread of coactivity throughout the network (Figure 10B). We showed that the wider spread of coactivity depended on the decreased ACh concentration characterizing the NREM state in our model, and not the difference between REM- and NREM-inputs, because application of the brief NREM stimulus during REM did not produce the NREM level of peak coactivity (Figure 9—figure supplement 1G). Analysis of the percentage of PFC cells spiking in periods surrounding NREM and REM events in the model also showed this increased spread of coactivity in NREM, in agreement with the experimental results (Figures 4A, E, 6A, B, 10B). However, it is unclear whether the widespread coactivity in NREM is a result of propagation of spiking due to recurrent excitation, or simply greater efficacy by the stimulus in directly inducing spikes in stimulated neurons during NREM. To investigate this, we analyzed model results by separating neurons within each stimulus event into stimulated and non-stimulated fractions. When stimulating in NREM, a widespread barrage of spiking from stimulated cells was followed by a smaller increase in the percentage of non-stimulated cell firing among both pyramidal and interneuron populations (Figure 10C). In contrast, during REM stimuli, the number of active stimulated neurons increased, but the effect on non-stimulated neurons was suppressive (Figure 10D). Again, swapping stimulus profiles showed that the cholinergic tone in the model plays a large role in creating this difference, especially as the brief NREM-like stimulus applied during REM produces no increase in percentage of coactive non-stimulated neurons (Figure 9—figure supplement 1I and J). These findings nicely align with experiments showing that ACh both suppresses the spread of cortical recurrent excitation and drives inhibition (Kimura et al., 1999; Dasgupta et al., 2018; Zerimech et al., 2020; Chen et al., 2015). Our results suggest that the ACh level regulates the degree to which inputs can induce widespread synchronous network activity, leading to the distinct activity patterns in NREM and REM. Low ACh concentration and the associated reduced inhibitory tone in NREM leads to a rapid spread of activity to many non-stimulated cells upon transient excitatory input to a subset of the network. Such a result suggests that the level of ACh regulates the fraction of cells in a circuit able to respond to a stimulus, accounting for the observation of synchronous bursts of activity among a large portion of cells during NREM (low ACh and low inhibition), whereas coactivity is restricted to smaller subsets of cells in the network during REM (high ACh and high inhibition).

Figure 10
Download asset Open asset
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). (D) Same as (C), but for REM. Note that non-stimulated cell population (pink) comprises a greater percent of total active neurons outside stimulation times (average stimulation time of 3 s centered at time of event), as there are more neurons in the non-stimulated group (only 30% of neurons are stimulated per stimulus event).

Discussion

One of the key challenges in understanding the role of REM sleep in memory consolidation has been disentangling the interregional neurophysiological processes involved in memory-related processing. Previous studies have primarily focused on the role of theta oscillations in REM sleep as a brain-wide pattern that coordinates activity across various brain regions (Boyce et al., 2017), due to their similarity with patterns observed during wakeful behavior. While much attention has been focused on investigating theta oscillations due to their prominence in the hippocampus during REM sleep, the precise dynamics that drive the coordination of multiple brain areas to potentially support systems memory consolidation have been elusive. Furthermore, the prominence of oscillatory phenomena in NREM sleep – including slow oscillations, spindles, cortical ripples, and CA1 SWRs – has enabled systematic dissection of their spatiotemporal dynamics and their relationship to hippocampal–cortical reactivation that supports memory consolidation (Staresina, 2024; Buzsáki, 2015; Girardeau and Zugaro, 2011; Klinzing et al., 2019), leaving REM sleep dynamics poorly understood in comparison.

Here, we show that phasic activity bursts in PFC during REM theta coincide with chains of local HFOs and elevated gamma power. These chains are associated with increased measures of oscillatory coupling between PFC and CA1, and show characteristic theta-modulated spiking activity patterns on the background of local suppression. REM PFC HFO chains organize sparse, sequential reactivation in PFC over extended time periods, and recruit a sub-population of CA1 neurons that exhibit differential regulation of activity over the course of sleep. This contrasts with NREM sleep, where independent PFC ripples are associated with relatively brief, synchronous bursts of local reactivation along with suppression of CA1 activity. We also show that CA1 co-activation during REM PFC HFOs is related to the CA1 activity suppression seen during NREM PFC ripples, with CA1 neurons that are highly coactive with PFC during REM HFO chains being strongly suppressed during NREM ripples, suggesting a dual sleep-stage role of PFC ripples and HFOs in regulating excitability in hippocampal circuits. Thus, these results uncover a novel, high-frequency pattern of PFC activity in REM that can potentially bridge the gap between structured PFC and CA1 firing patterns in NREM-REM stages and changes in hippocampal plasticity during sleep.

In contrast to PFC ripples in NREM sleep, which are involved in suppressing CA1 activity as a possible mechanism for tagging specific memory traces for consolidation (Shin and Jadhav, 2024), here we found that REM cortical HFOs engage specific populations of CA1 neurons to potentially regulate activity homeostasis. Of particular interest regarding this phenomenon of activity regulation is a recent study examining the differential roles of NREM and REM sleep in reorganization of drifting hippocampal assemblies over the course of sleep (Bollmann et al., 2025). This study focused solely on hippocampal assemblies, and the results are in line with our finding of the potential roles of PFC ripples and HFOs in NREM and REM sleep, respectively. We have previously shown that independent NREM PFC ripples that are dissociated from hippocampal SWRs broadly suppress CA1 activity, and that this suppression is related to the strength of CA1 reactivation during SWRs and CA1 assembly reinstatement during subsequent behavior (Shin and Jadhav, 2024). This bidirectional modulation of hippocampal representations during PFC ripples and CA1 SWRs may underlie assembly reorganization during NREM sleep. Relatedly, we report here that CA1 activity during REM PFC HFOs is associated with differential changes in excitability in a subset of hippocampal neurons during sleep, which may, in part, contribute to the stabilization of assemblies in REM sleep and thus drive the sleep-stage dependent reorganization of hippocampal assemblies. This finding is in line with REM sleep’s role in recalibration of neural activity in hippocampal–cortical circuits (Lendner et al., 2023). We speculate that PFC REM HFO chains enable selective linkages between the sparsely reactivated PFC ensembles in coordination with a subset of CA1 neurons.

Cortical HFOs (~120–160 Hz) in REM have been reported before, especially in parietal cortex, although they have been variously termed fast gamma or HFOs (Brankačk et al., 2012; Bueno-Junior et al., 2023; Tort et al., 2013; Scheffzük et al., 2011), with previous reports of strong coupling with theta oscillations in REM (Scheffzük et al., 2011). HFOs have also been reported in parietal cortex and hippocampus during active waking behavior and REM, and it has been argued that HFOs in hippocampus are distinct from hippocampal SWRs (Tort et al., 2013). We detected HFOs in REM sleep using a similar frequency band (150–250 Hz) as cortical NREM ripples. These high-frequency events in NREM are denoted as ‘ripples’, a term used by previous studies for NREM high-frequency events in cortex (Ghosh et al., 2022; Aleman-Zapata et al., 2022; Khodagholy et al., 2017; Shin and Jadhav, 2024; Vaz et al., 2019; Helfrich et al., 2019), although this is likely a separate phenomenon from hippocampal SWRs. PFC NREM ripples that occur independently of hippocampal SWRs are associated with temporally synchronized local PFC reactivation and suppression of hippocampal activity and reactivation (Shin and Jadhav, 2024), while in REM, we show here that they present as PFC HFO chains during theta-modulated activity, associated with sparse, temporally extended local reactivation in the backdrop of suppression, along with co-activation of a hippocampal-subpopulation. These structured cortical sequences seen during HFO chains are gated by theta oscillation bursts in REM sleep, with a higher rate of occurrence in phasic REM than tonic REM substages, and are reminiscent of hippocampal theta sequences (Drieu and Zugaro, 2019), but may have different underlying mechanisms. We examined PFC HFO events during behavior on the W-maze, but did not find similar spiking activity patterns as REM HFO events (Figure 9—figure supplement 2), confirming that this is a REM-exclusive phenomenon.

Based on previous studies, external drive can lead to the generation of PFC ripples and synchronized activity in NREM, whereas high ACh and high inhibitory tone in REM can lead to local suppression in response to external drive (Chen et al., 2015; Sanzeni et al., 2020; Niethard et al., 2016). Notably, a recent computational study demonstrated that differential cholinergic modulation, together with the sequential architecture of sleep (NREM to REM transition), can support memory consolidation through the complementary shaping of ensemble activity during NREM and REM sleep (Satchell et al., 2025). We show here that a simple network model that incorporates the effect of differential ACh tone on excitatory and inhibitory neurons in REM vs. NREM states is able to reproduce the observed experimental activity patterns in response to stimulation, akin to phasic theta bursts during REM and isolated, brief inputs in NREM. It is likely that PFC NREM ripples and REM HFOs are generated locally and modulated by external inputs. We speculate that theta-modulated input to PFC during REM is driven by hippocampal inputs, but it is also possible that theta is driven by inputs arising from regions other than the hippocampal network, a possibility which will require further experimental investigation.

We observed that isolated REM HFOs in PFC were more prevalent than chains, especially in tonic REM, potentially reflecting the need for coordinated activation of distributed circuits to elicit HFO-associated phasic events. One possible mechanism driving this phenomenon could involve PGO waves, which are prominent during phasic REM and are known to play a role in modulating hippocampal plasticity during REM sleep (Ramirez-Villegas et al., 2021). These waves are also coordinated with both SWRs in NREM and theta states in REM (Ramirez-Villegas et al., 2021; Tsunematsu et al., 2023), suggesting that PGO-associated amplification of hippocampal theta oscillations could drive activity above a certain threshold for HFO generation in PFC. Additionally, ACh may also play a role in altering the threshold for HFO generation in PFC in conjunction with PGO coordination. It has been reported that there is coordinated release of ACh in PFC and dorsal hippocampus, potentially elevating theta coherence between the two areas (Teles-Grilo Ruivo et al., 2017). Indeed, a number of studies have shown that ACh boosts theta oscillations (Nuñez and Buño, 2021). Furthermore, ACh induces decorrelation in cortical circuits, which improves cortical representations of relevant stimuli (Goard and Dan, 2009), and our finding that REM HFOs are associated with an overall decrease in PFC population activity, reactivation strength, and cofiring, in stark contrast to NREM ripples, suggests a role for REM HFOs in amplifying local cortical signal-to-noise to enable linking of cortically encoded information with specific populations of hippocampal neurons. Lastly, cholinergic input to the hippocampus also suppresses SWR generation (Zhang et al., 2021; Vandecasteele et al., 2014), which could facilitate the sparse CA1 reactivation during REM sleep for selective changes in excitability and stabilization of salient representations in coordination with cortical networks for subsequent expression during behavior. Future experimental and network modeling studies can investigate how neuromodulatory mechanisms and NREM and REM state-specific cortical activity patterns act in concert to support learning and memory consolidation.

Taken together, these results establish a new cortical LFP pattern in REM sleep (PFC REM HFO chains) associated with a distinctive neuronal population activity signature. The findings are aligned with REM sleep’s role in memory consolidation and provide a neural basis for the cortical–hippocampal modifications that occur during sleep to support activity reorganization and mnemonic processing. The complex interplay between oscillations of varying frequencies revealed here further highlights the importance of segregating sleep into multiple states to uncover how systems memory consolidation can be supported through the multiplexing of distinct phases, both during NREM and REM sleep. The interareal dynamics associated with PFC REM HFOs involve both bottom-up and top-down processes operating in precise temporal coordination, with theta-mediated communication leading to PFC HFO-associated reactivation and modifications in hippocampal excitability. These findings, in conjunction with PFC-mediated suppression of CA1 during NREM sleep, suggest a dual sleep stage processing of mnemonic representations that have implications for understanding the circuit mechanisms underlying systems memory consolidation, and can reveal how neural processes in sleep go awry in disease states.

Limitations

While the present results establish a potential role for REM sleep PFC HFOs, several limitations of our study are acknowledged here. First, we were unable to directly link REM HFO-specific reactivation to learning. Spatial memory tasks in rodents similar to our study have revealed causal links of NREM reactivation driven by CA1 SWRs to learning and memory consolidation (Girardeau et al., 2009; Maingret et al., 2016; Chang et al., 2025; Gridchyn et al., 2020), but reports of REM reactivation are rare (Louie and Wilson, 2001). Investigation of the role of REM reactivation in learning and memory consolidation will need utilization of tasks known to require both REM and NREM states, such as schema-based inference or emotional memories (Abdou et al., 2024; Porter et al., 2025; Zeithamova et al., 2012; Cairney et al., 2015). Second, behavioral task learning and interleaved sleep sessions were conducted in a few hours within a single day, which did not allow us to densely sample the circadian cycle. Third, we did not record eye movements or ponto-geniculo-occipital (PGO) waves, both of which would have allowed for more accurate segregation of tonic and phasic REM sleep states (Simor et al., 2020). Although we observed a bias for isolated and chain HFOs to occur in putative tonic and phasic REM substates, respectively, the scarcity of putative phasic REM bouts made the direct comparison based on substage difficult. Finally, although our model predicts that distinct cell-type activity profiles shape REM sleep HFO dynamics, we did not record a sufficient number of interneurons to test these predictions directly. Future studies using appropriate behavioral tasks, longitudinal sleep recordings, and cell-type-specific opto-tagging will be able to resolve these limitations and further clarify the roles of HFOs in REM sleep.

Resources availability

Materials availability

This study did not generate any new unique reagents.

Experimental model and subject details

All experimental procedures were approved by the Institutional Animal Care and Use Committee at Brandeis University (Protocol #24,001 A) and conformed to US National Institutes of Health guidelines. Ten adult male Long-Evans rats (450–600 g, 4–6 months; RRID:RGD_2308852) were used in this study. Animals were individually housed and kept on a 12 hr regular light/dark cycle.

Methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Chemical compound, drugCresyl VioletAcros OrganicsCat#: AC229630050
Chemical compound, drugFormaldehydeThermo FisherCat#: 50-00-0,67561,7732-18-5
Chemical compound, drugIsofluranePatterson VeterinaryCat#: 07-806-3204
Chemical compound, drugKetaminePatterson VeterinaryCat#: 07-803-6637
Chemical compound, drugXylazinePatterson VeterinaryCat#: 07-808-1947
Chemical compound, drugAtropinePatterson VeterinaryCat#: 07-869-6061
Chemical compound, drugBupivacainePatterson VeterinaryCat#: 07-890-4881
Chemical compound, drugEuthasolPatterson VeterinaryCat#: 07-805-9296
Chemical compound, drugSucroseSigma-AldrichCat#: S8501-5KG
Strain, strain background (Rattus norvegicus)RatCharles RiverCat#: Crl:LE 006; RRID:RGD_2308852
Software, algorithmMATLAB 2022aMathworksMARRID:SCR_001622
Software, algorithmTrodesSpikeGadgets; http://www.spikegadgets.com
Software, algorithmMatclustMattias P. Karlsson; https://www.mathworks.com/matlabcentral/fileexchange/39663-matclust,V1.7
Software, algorithmChronuxPartha Mitra; http://chronux.org/RRID:SCR_005547V2.12
Software, algorithmRobustICAVicente Zarazoso, Pierre Comon; https://webusers.i3s.unice.fr/~zarzoso/robustica.html
OtherElectrophysiology data acquisition systemSpikeGadgets; http://www.spikegadgets.comRRID:SCR_021623
Other12.7 μm NiCr tetrode wireSandvikCat#: PX000004

Surgical implant and electrophysiology

Surgical implantation procedures were as previously described (Shin et al., 2023). Animals (n=10) were implanted with a microdrive array containing 30, 32, or 64 independently moveable tetrodes targeting right dorsal hippocampal region CA1 (–3.6 mm AP and 2.2 mm ML) and right PFC (+3.0 mm AP and 0.7 mm ML). For 64 tetrode recordings (n=1 animal), CA1 and PFC were targeted bilaterally. Tetrodes were split equally between PFC and CA1 (15, 16, or 32 tetrodes in each region). On the days following surgery, hippocampal tetrodes were gradually advanced to the desired depths using characteristic EEG patterns (sharp wave polarity, theta modulation) and neural firing patterns as previously described (Shin et al., 2023). A CA1 tetrode in corpus callosum served as hippocampal reference, and another tetrode in overlying cortical regions with no spiking signal served as prefrontal reference. A ground (GND) screw installed in the skull overlying cerebellum also served as a reference. All spiking activity, gamma, and ripple/HFO-filtered LFPs (gamma: 40–100 Hz; ripple/HFO: 150–250 Hz; see below) were recorded relative to a local reference tetrode. For detection of slower frequency components (i.e. theta and delta), LFP was recorded with respect to the GND screw. Electrodes were not moved at least 4 hr before and during the recording day.

Data were collected using a SpikeGadgets data acquisition system (SpikeGadgets LLC). Spike data were sampled at 30 kHz and bandpass filtered between 600 Hz and 6 kHz. LFPs were sampled at 1.5 kHz and bandpass filtered between 0.1 Hz and 400 Hz. The animal’s position was recorded with an overhead color CCD camera (30 fps) and tracked by color LEDs affixed to the headstage. Additionally, the animals’ speed was calculated based on predetermined cm/pixel values and the position displacement between frames captured at 30 fps.

Behavior

Ten animals were trained on a novel W-maze in a single day, as previously described (Shin and Jadhav, 2024; Shin et al., 2019; Shin et al., 2023). This single-day W-maze alternation learning task is SWR-dependent and requires hippocampal–prefrontal interactions (Maharjan et al., 2018; Jadhav et al., 2012; Fernández-Ruiz et al., 2019). Briefly, following recovery from surgical implantation (~7–10 days), animals were food-deprived and again retrained on a linear track for 2–3 days before exposure to the W-maze task. During the recording day, animals were introduced to the novel W-maze (~80 × 80 cm with ~7 cm wide tracks) for the first time and learned the task rules over eight behavioral sessions during the animals’ light phase between the hours of 9 AM and 6 PM. Each behavioral session lasted 15–20 min and was interleaved with 20–40 min rest sessions in the sleep box for a total of 17 separate sessions (8 W-maze and 9 sleep sessions). The total recording duration spanned ~6–7 hr within a single day. On the W-maze, animals were rewarded for performing a continuous spatial alternation. The task rules are as follows: returning to the center well after visits to either side well (left or right well; inbound trajectories) and choosing the opposite side well from the previously visited side well when at the center well (outbound trajectories). Rewards were automatically delivered to the animal through custom 3D printed reward wells upon nose poke.

Histological verification of recording sites

At the completion of each experiment, electrolytic lesions were made by sending current (30 µA) down 3 electrodes per tetrode for 3–4 s each. 24 hr later, animals were sacrificed by intraperitoneal injection of Euthasol and subsequent intracardial perfusion with 0.9% saline followed by 4% formaldehyde. The skull, with the attached drive, was submerged in 4% formaldehyde for 24 hr before tetrode retraction and brain extraction. Brains were stored in a 4% formaldehyde and 30% sucrose solution for at least 1 week before slicing with a freezing microtome (Leica) at 50 µM. Slices were Nissl stained and imaged at ×4 on a brightfield microscope (Keyence) and stitched together to generate composite images for lesion localization.

Spike sorting

Single units were sorted using a custom manual clustering program (MatClust, M. P. Karlsson) and were distinguished based on spike width, peak and trough amplitude, and principal components. Cluster quality was measured using isolation distance (Schmitzer-Torbert et al., 2005), and only well-isolated neurons were included for analysis.

Sleep state identification

To detect bouts of NREM and REM sleep, animals’ head speed and CA1 LFPs were used as previously described (Mizuseki et al., 2011; Tang et al., 2017). During sleep box sessions, times of awake activity and immobility were defined as periods where head speed was greater or less than 4 cm/s, respectively. Candidate NREM sleep bouts were identified as periods where head speed remained <4 cm/s after an immobility period of at least 60 s. NREM and REM sleep were further identified and separated by using a theta to delta ratio metric. Briefly, CA1 LFP (referenced to GND) for each tetrode was bandpass filtered (6–12 Hz for theta; 1–4 Hz for delta) and averaged, and the TD ratio power was calculated. Periods that exceeded a set threshold (mean +1 SD) for >10 s were categorized as REM bouts within candidate NREM bouts (Zhang et al., 2020; Rothschild et al., 2017). Additionally, candidate REM bouts were visually inspected to further confirm the transition into REM sleep. For comparison of velocity and TD ratio, the average velocity and TD ratio in each epoch for NREM and REM sleep were computed and compared. Additionally, only epochs with at least 30 s of REM sleep were included for analysis. Unless otherwise stated, all analyses were restricted to the 36 (of 90) sleep epochs that satisfied the inclusion criteria.

Intracranial electromyogram

To estimate EMG from the neck, jaw, and face muscles from intracranial LFP signals, a previously established method was used (Watson et al., 2016; Schomburg et al., 2014). Briefly, we extracted the 300–600 Hz filtered signal (Butterworth filter at 300–600 Hz; filter shoulders from 275 to 625 Hz) and calculated the pairwise Pearson correlations in each 500ms bin between pairs of electrodes during the entire sleep epoch. Since we used tetrode recordings and cannot determine the exact separation of electrode sites, we elected to calculate the correlations between the reference tetrode in corpus callosum and CA1 ripple detection tetrodes within the dorsal CA1 pyramidal cell layer (~500 µm separation between corpus callosum and the dorsal CA1 pyramidal cell layer), thus maximizing the separation between the electrode sites. Previous studies have used electrode channels on probes at least two shanks apart (~400 µm). The average of all reference/CA1 layer tetrode pairs was calculated and compared for wake and REM sleep. For the distributions of intracranial EMG values in Figure 1—figure supplement 1D, each 500 ms bin correlation value was reported, and for the comparison between wake and REM sleep, the averages during each wake and REM bout were calculated and compared. A strong relationship between intramuscular EMG and intracranial EMG has been reported (Schomburg et al., 2014).

LFP event detection

PFC ripple/HFO detection: High-frequency events were detected during immobility periods (<4 cm/s) as previously described (Tang et al., 2017). Tetrodes with recorded cells in CA1 and PFC were used as candidate electrodes for event detection as well as downstream LFP analysis. Each LFP from candidate electrodes was bandpass filtered in the ripple/HFO band (150–250 Hz), and the envelope power was calculated using a Hilbert transform and subsequently smoothed with a Gaussian kernel (σ=6). Events were detected as contiguous periods where the ripple/HFO power exceeded 3 SD of the mean on at least one tetrode. A minimum event duration of 15 ms was implemented to exclude spurious noise in the LFP signal. Events were further refined as periods around the detected event that remained above the mean (start and end times). The amplitude of each event was reported as the number of SD above the mean (Shin et al., 2019). To estimate the frequency of each event, the power in each frequency bin was calculated around each event (±500 ms) using multi-taper time-frequency analysis (Chronux toolbox; http://chronux.org/; Mitra and Bokil, 2008) and was z-scored relative to a baseline spectrogram. The frequency of each event was defined as the frequency with the highest power in the ripple/HFO band (100–250 Hz) between the start and end time of each event (Sullivan et al., 2011). Furthermore, HFO chains were defined as 2 or more events that occurred <200 ms apart. Average amplitude plots for each animal were generated by aligning the peaks of each event and averaging.

PFC gamma detection: Bouts of high gamma power were detected as above after bandpass filtering in the gamma band (40–100 Hz).

High theta power detection: Each LFP from candidate PFC tetrodes was used to delineate periods of high theta power during REM sleep. LFPs were bandpass filtered with a theta band filter (6–12 Hz), the envelope power was calculated using a Hilbert transform, averaged across tetrodes, and periods that exceeded 4 SD of the mean during REM sleep were designated as bouts of high theta power. Only high-power periods exceeding 900ms were included as high theta power periods.

PFC spindle detection: Spindle events were detected as previously described (Darevsky et al., 2024). The mean LFP across all candidate PFC electrodes was calculated and the resulting signal was bandpass filtered between 10 and 16 Hz with a zero-phase shifted, third-order Butterworth filter. The envelope power was calculated using a Hilbert transform and subsequently smoothed with a Gaussian kernel (α=2.5). Events that exceeded 2.5 SD of the mean were extracted, and events separated by <300 ms were concatenated. Spindle events were restricted to NREM sleep before analysis. Spindle trains were defined as 2 or more spindle events that occurred ≤2.78 s apart as previously established (Darevsky et al., 2024).

Spectral analysis

Wavelet spectral analyses were utilized to calculate the ripple/HFO aligned power spectra in PFC (2–350 Hz, Morlet wavelets). For each epoch, the PFC tetrode with the greatest number of detected PFC ripples/HFOs was used. The power at each level of the wavelet transform was calculated around each event (±500 ms) and was z-scored relative to a baseline spectrogram calculated for the entire epoch. For each animal, spectrograms from all candidate sleep epochs were averaged.

Event autocorrelation

To calculate the autocorrelation of ripples/HFOs to assess the structure of event occurrence, identical event time vectors were lagged (max time lag 1 s, 10 ms time bins), and the correlations were calculated. The results were smoothed with a Gaussian kernel (σ=3) and averaged across all animals and epochs. Similarly, the cross-correlation between gamma and HFO time vectors was calculated by using the gamma times as the reference vector. Thus, positive and negative lag values indicate that the HFOs tended to occur before and after gamma events, respectively.

Event aligned multi-unit activity

To characterize population activity aligned to LFP events, spike trains of each neuron were binned into 5 ms time windows (Figures 2A and 4F). Multi-unit activity was aligned to PFC ripple/HFO and gamma events and was normalized by the mean population firing rate during either NREM or REM sleep, depending on the events assessed. Subsequently, the power spectral density of the ripple/HFO aligned multi-unit activity was assessed (frequency limits, 2–12 Hz; MATLAB, pwelch).

Control for periods of high theta power

To determine whether theta-modulated population activity in PFC can be largely explained by periods of high theta power, instead of association with HFO events, we extracted surrogate timestamps to align PFC activity to (Figure 4—figure supplement 1K–N). Briefly, for each epoch, the preferred theta phase bin of HFO events was calculated from the PFC tetrode with the highest average theta power and used to randomly select surrogate times to align PFC MUA to. Each preferred phase bin at least 250 ms away from a detected HFO was considered a candidate event. Then, times were randomly chosen, and the number of events selected matched the number of HFO events in that epoch. Theta power at each surrogate time point as well as HFO event was extracted and z-scored to control for any inter-animal variability in theta power. Lastly, the event aligned MUA was split into quartiles based on theta power and compared.

Theta inter-peak intervals during bouts of high and low theta power

Since previous studies have segregated putative tonic and phasic substages of REM sleep based on CA1 theta frequency (Mizuseki et al., 2011), we calculated the inter-peak intervals of the theta oscillation to determine whether bouts of low and high theta power may correspond to putative tonic and phasic substages, respectively (Figure 3—figure supplement 1B). We bandpass filtered and averaged the LFPs from candidate CA1 tetrodes between 5 and 12 Hz, and the peaks were extracted by detecting positive-to-negative zero crossings of the derivative of the filtered, averaged signal. Lastly, the inter-peak intervals were calculated within low and high theta power periods for each epoch and compared. Note that we defined low and high theta periods based on PFC electrodes, whereas we calculated inter-peak intervals based on CA1 electrodes.

Cross-frequency phase-amplitude coupling

PAC was assessed for a range of frequency combinations using a previously established approach (Tort et al., 2010). Since the rate of REM HFOs was highest during bouts of high theta power, we focused on these periods for assessing PAC. For each epoch, the tetrode with the highest HFO power was chosen, and an averaged comodulogram was generated across all animals and epochs for the phase and amplitude frequency ranges from 2 to 12 Hz and 40 to 250 Hz, respectively.

Comodulogram: Briefly, the phases were extracted from the Hilbert transform of the raw LFP of the candidate electrode and binned into 18 bins (intervals of 20°), and the envelope power of the amplitude modulated frequency band was calculated by taking the absolute value of the Hilbert transform. Then, the amplitude of each frequency from 40 to 250 Hz was calculated in each phase bin from 2 to 12 Hz and visualized (Figure 3B).

Modulation index: To quantify the depth of amplitude modulation by the phase frequencies, gamma (40–100 Hz) and HFO (150–250 Hz) signals were separately processed, and the mean gamma or HFO amplitude in each phase bin was calculated. The resulting values were used to calculate the modulation index (Tort et al., 2010), which specifies the strength of phase modulation. Additionally, the phase bin with the maximum amplitude was reported as the preferred phase. For comparison, the phase and amplitude vectors were circularly shifted independently, and the modulation index was calculated.

Theta nesting

To further investigate the coupling of theta, gamma, and HFOs, in PFC during REM sleep, we examined theta nesting (Vaz et al., 2017; Daume et al., 2024) in two different frequency bands (40–100 Hz for gamma and 150–250 Hz for HFO) relative to the theta oscillation (Figure 2C, Figure 3—figure supplement 2E, F). Since we are evaluating the gamma and HFO coupling specifically during these transient HFO events, we restricted our analysis to high theta periods since the rate of REM HFOs was highest during these REM bouts (Figure 3—figure supplement 1D). For each PFC electrode, the theta phase bin for which gamma power was maximal was identified. Then, each instance of this phase bin during high theta periods in REM sleep was detected and the index of the highest gamma power within each preferred phase bin was identified. Only electrodes where there were at least 20 instances of the preferred phase were included. Lastly, the averaged waveform was calculated by averaging the unfiltered LFP centered on every instance of maximal gamma power in the preferred bin. The above procedure was performed separately for HFO band activity. For comparison between the frequency of averaged gamma and HFO waveforms and the preferred phase of theta coupling, electrodes that exhibited at least 5 local maxima within a 100ms window surrounding the preferred phase were analyzed.

To further validate coupling of gamma and HFO to theta, surrogate waveforms were generated by randomly selecting a different phase bin for each theta cycle for alignment of gamma or HFO power. Then, peak-to-trough amplitudes were calculated from the average waveform within the 100 ms window centered on peak power, averaged, and compared between the preferred phase bin aligned data and the random phase bin surrogate.

Phase locking

Phase locking during active behavior on the W-Track and during REM sleep was quantified using previously developed methods (Siapas et al., 2005). The CA1 or PFC tetrode with the greatest ripple/HFO power was filtered in the theta (6–12 Hz), gamma (40–100 Hz), or ripple/HFO (150–250 Hz) range and was used to extract phase information. For HFO and gamma phase locking of PFC cells, each spike that occurred across all extracted events was assigned a phase, and the distributions of the peak phases of all cells were compared between gamma and HFO. During run sessions (for ‘REM theta phase shifting’ analysis below, Figure 7—figure supplement 2H and I), spikes were restricted to CA1 theta periods where the animals’ speed was >5 cm/s and were assigned a phase between 0 and 360°. Only neurons with at least 10 spikes during theta periods were included in this analysis. The Rayleigh test was used to determine whether the phase distribution of each cell deviated from circular uniformity (p<0.05 criterion). Additionally, for phase locking to gamma and HFOs, the preferred phases of each cell were calculated separately for gamma or HFOs. Thus, a cell could be significantly phase locked to both gamma and HFOs. To test for phase preference differences, the distributions were compared using a nonparametric Watson’s U2 test. Subsequently, for each theta-modulated cell, phase-locking strength to the theta oscillation was measured using PPC (Vinck et al., 2010). PPC was calculated by taking the dot product of the angular difference between each pair of phase values then averaging. Similar procedures were performed to determine the HFO and gamma phase preference of PFC cells.

Additionally, PPC was used to assess the strength of theta phase locking of HFO/ripple events in REM or NREM sleep. Only epochs with at least 20 events of each type were included for analysis.

Theta coherence

Theta coherence between pairs of CA1 and PFC tetrodes was calculated during REM sleep periods centered on REM HFOs, and each frequency was normalized by the mean and SD of the baseline coherogram computed across the entire epoch (Figure 3E and I). Coherograms were averaged over CA1-PFC tetrode pairs. For quantification and comparison, theta coherence was measured as the peak coherence between 6 and 12 Hz in a±1 s window around HFO events. For visualization of coherograms, the averaged coherograms across all animals and epochs were smoothed with a 2D Gaussian kernel (σ=4). All coherence-based analyses were conducted using the Chronux toolbox (Mitra and Bokil, 2008). Only epochs with at least 20 events were included for quantification of peak theta coherence.

Phase slope index

Since we observed elevated theta coherence between CA1 and PFC during chained events, we assessed the directionality of information flow at different frequency bands (0–50 Hz) centered on these events (±1 s) using the PSI (Figure 3J; Nolte et al., 2008), which was calculated with the FieldTrip analysis toolbox (Oostenveld et al., 2011). In practice, PSI is used to assess the consistency of phase lag relationships between two signals across different frequency bands and is a measure that is weighted by oscillatory coherence. We opted to use PSI to estimate the directional flow of information instead of other methods, such as Granger causality, since it has been demonstrated that PSI is less prone to false positives (Nolte et al., 2008). Briefly, for each epoch, the PSI was calculated and averaged across each candidate CA1-PFC tetrode pair. Here, positive PSI values indicate that PFC theta leads CA1, while negative values indicate that CA1 leads PFC. Then the distribution of mean PSI values within 6–12 Hz was tested against 0 with a t-test. Only epochs with at least 5 chained events were included for analysis.

Ripple/HFO aligned modulation

Single unit ripple/HFO modulation in CA1 or PFC was calculated as previously described (Figure 4A and D; Shin and Jadhav, 2024; Jadhav et al., 2016). This analysis was restricted to neurons that fired >50 spikes cumulatively across all peri-event windows (±2 s). For each candidate neuron, an averaged event-triggered (NREM or REM event) time histogram (event-PSTH) was calculated. For comparison, 1000 shuffled event-PSTHs were generated by circularly shifting the spikes around each ripple/HFO event by a random amount. To determine whether a neuron was significantly modulated, the squared difference between the real event-PSTH and the shuffled data in the –200 to +200 ms window surrounding event onset was calculated. If this real value exceeded 95% of the shuffled values, the neuron was categorized as significantly modulated. The modulation index for each cell was given by calculating the average difference between the activity in the response window and a background window –600 to –200 ms from event onset. Furthermore, excited (EXC) and inhibited (INH) neurons were separated based on whether the modulation in the event response window was positive (EXC) or negative (INH). Modulation timing distributions for EXC PFC neurons during NREM and REM PFC events were compared by binning (2 ms bins) the peak responses within the modulation window.

Event type prediction

High-frequency event type (NREM or REM) in PFC was predicted by fitting a linear model (MATLAB, fitclinear, Figure 4C) to PFC spike count data. Since the NREM events largely outnumbered the REM events, the event counts were equalized prior to training by randomly excluding a subset of NREM events. The binary decoders were built using 10-fold cross-validation, and prediction accuracy was assessed by comparing performance against models constructed with randomly permuted spike matrices (n=1000). To account for the excluded events, this entire process was repeated 10 times for each animal and epoch.

Ripple and HFO cofiring

To quantify ripple/HFO cofiring during events, a pairwise cofiring z-score was calculated as previously described (Figures 4E and 7B and C; Cheng and Frank, 2008). Briefly, this cofiring measure estimates how likely two cells will fire together during events and is normalized by the incidence of spiking independent of the other cell:

Zcoactivitya,b=nab−nanbNnanb(N−na)(N−nb)/(N2(N−1))

where N is the number of events, na and nb are the number of events where cell a and cell b were active, respectively, and nab is the number of events where both cells were active. Here, cofiring assesses coincident activity between neuron pairs within the start and end times of events, and thus, no explicit temporal threshold was implemented. The same procedure was used to calculate cofiring between pairs in the model simulations.

Detection of population suppression

To validate the suppression of PFC MUA surrounding REM PFC HFOs, periods of decreased PFC activity were separately extracted, and the probability of REM HFOs surrounding these events was calculated (Figure 4H, Figure 4—figure supplement 2B, C). First, MUA was binned (5 ms bins), normalized by the average activity during REM sleep, and smoothed with a Gaussian kernel (σ=3). The population activity was then z-scored, and periods that fell below 0 for 300–500 ms were extracted as population suppression events. Lastly, the probability of HFO occurrence surrounding these events was calculated in 100ms bins. To obtain the 99% confidence intervals, suppression events were extracted from shuffled multi-unit data. Briefly, each cells’ spike count was circularly shifted independently, and events were detected. A total of 1000 shuffles were performed, and the probability of HFO occurrence surrounding these events was calculated for each shuffle.

Relationship between inter-event-interval and content of adjacent HFOs

To assess whether there was a relationship between IEI and population spike content between HFO events, the spiking activity of PFC neurons during HFOs was represented as ones and zeros for either active or inactive during all events (Figure 5B, Figure 6—figure supplement 1A). The Pearson correlation between population activity vectors of two adjacent HFOs was calculated, and the IEI was given by calculating the temporal difference between the centers of the HFO events. Then, the relationship between spike content and IEI was assessed using the Pearson correlation coefficient. To validate this relationship, we employed two other methods for quantifying similarity of activity across two events: the Jaccard index and cosine similarity. The Jaccard index, or the Jaccard similarity coefficient, was used to quantify the similarity between two sets. It is defined as the ratio of the size of the intersection of the two sets to the size of their union. The Jaccard index (J) is expressed as

J(A,B)=|A∩B||A∪B|

where A and B represent the two sets being compared (spike activity for HFO A vs HFO B), ∣A∩B∣ is the number of elements in the intersection of the sets, and |A∪B| is the number of elements in the union of the sets. The resulting value ranges from 0 to 1. Similarly, cosine similarity, which quantifies the cosine of the angle between two vectors, was used to measure the similarity between two spike vectors. Cosine similarity is given by

CS=A∙B∥A∥∥B∥

where A and B are the two vectors being compared, A∙B is the dot product, and ∥A∥ and ∥B∥ represent the magnitudes of the vectors. The values range from –1–1 where 1 indicates identical vectors and –1 indicates that the vectors have an angular difference of 180°. Actual spike counts were used for cosine similarity calculation.

Rank-order correlation

To evaluate the consistency with which populations of neurons in PFC fired in sequence during chained events, we calculated the rank-order correlation between each chain and its template as previously described (Figure 5C, Figure 6—figure supplement 1B; Stark et al., 2015; Valero et al., 2021). The population activity during each HFO in a chained event was concatenated (inter-event activity excluded), and for each concatenated event with ≥5 cells active, the spiking activity was reduced to each unit’s first spike time, chronologically ordered, and normalized from 0 to 1. Then, a cross-validated template-based method was used to calculate sequence similarity between each event and its template. For each event, a template was created by averaging the normalized rank of each cell across all events excluding the event in question. The rank correlation was reported as the mean Pearson correlation coefficient between the event and template across all events. For comparison, a null distribution was generated by jittering spike times (1000 jitters) and performing the rank-order analysis above. For the pseudo-chain control, each HFO chain was time-shifted ±5 s such that each HFO within the chain was shifted the same value (1000 time shifts), thus preserving the IEIs (i.e. the entire chain was shifted coherently). Additionally, the entire shifted chain must have stayed within an REM sleep bout. The above rank-order procedure was performed, and the real correlation was compared to the pseudo-chain null distribution.

Assembly detection

To detect cell assemblies in PFC populations, a previously described method was used (Figure 6; Peyrache et al., 2009; van de Ven et al., 2016). Significant coactivity patterns in PFC were detected using a method based on principal and independent component analyses. For each W-Track epoch, the spike trains of each neuron were binned into 20ms time windows and Z-scored to eliminate any biases due to differences in firing rates. Principal component analysis was then applied to the Z-scored spike matrix (Z). The resulting eigenvalue decomposition of the correlation matrix (C) was given by:

∑i=1nλipipiT=C

where pi is the ith principal component of C and λi is the corresponding eigenvalue. Note that C=1nZZT is the correlation matrix of Z. To determine the number of significant coactivity patterns in the data, the Marcenko-Pastur law was used to calculate a threshold eigenvalue λmax, which was given by λmax=(1+n/B)2, where n is the number of recorded units and B is the total number of bins. Note that an eigenvalue exceeding λmax indicates that the corresponding cell assembly explains more correlation than expected if the activity of neurons was independent of each other. The number of eigenvalues above λmax, NA, represents the number of significant patterns found in the data. These principal components were then projected back onto the original z-scored spike data:

ZPROJ=PSIGNTZ

where PSIGN is the nXNA matrix with NA columns. Independent component analysis (ICA, RobustICA; Zarzoso and Comon, 2010) was then used to identify patterns that were maximally independent from each other. ICA was applied to the matrix ZPROJ to find an NAXNA unmixing matrix, W, which was used to obtain the weights of each cell in each assembly.

V=PSIGNW

Since the sign of the output of ICA is arbitrary, the signs of the weight vector were set such that the highest absolute weight was set to positive.

Assembly reactivation strength

To assess the reactivation of these detected assemblies during sleep sessions (Figure 6), the expression strength of each pattern over time was given by:

Rk(t)=z(t)TPkz(t)

Rk(t) is the reactivation strength of assembly k at time t, z(t) is the z-scored firing rate vector for all neurons at time t, and Pk is the projection matrix for assembly k, constructed by taking the outer product of its weight vector vk. Additionally, the diagonal of the projection matrix was set to zero to reduce the influence of highly active neurons, and the reactivation strength was calculated such that only the positive weight cells were assessed. Thus, only patterns of coactivity are detected as periods of high reactivation strength.

Assembly sequence detection surrounding HFOs

The consistency of assembly reactivation surrounding REM HFOs was assessed by implementing a cross-validation process similar to procedures used for place cell sequences (Figure 6A–D; Sosa et al., 2025; Plitt and Giocomo, 2021). We employed a split-event procedure where HFO events were randomly split into two halves. The average reactivation strength for each assembly was then calculated across all HFOs across both halves separately. This was performed across all animals and epochs. The assemblies aligned to the first half were subsequently sorted by peak reactivation strength within a ±500 ms window around HFO times, and the assemblies aligned to the second half were sorted using the same sort indices such that each row in the two matrices correspond to the same assembly. The consistency of temporal assembly reactivation was quantified as the Pearson correlation between the peak reactivation strengths of the first and second halves. This procedure was conducted 1000 times using different subsets of HFOs for the first and second halves to obtain a distribution of correlation coefficients. For comparison, a distribution of shuffled data was obtained by circularly shuffling assembly reactivation strength prior to HFO alignment, and the above procedure was performed. This procedure was also performed for gamma and NREM PFC ripples for comparison.

Furthermore, the slope and peak differences were calculated as additional metrics of sequence and temporal reactivation consistency between the two halves of data, respectively. For the slope difference metric, the slope of the best fit line between assembly ID and peak reactivation index was taken for the two halves of data and compared. A smaller slope difference compared to shuffled data indicates that assembly reactivation is structured in a more similar manner across the two halves of the real data. For peak reactivation difference, the temporal displacement of the peak reactivation index between the two halves of data was calculated and compared to shuffle. A high probability of small peak differences indicates that the timing of peak reactivation of assemblies relative to HFO chain onset is similar across the two halves of data. Shuffling of assembly strength was carried out as above.

Reactivation strength across epochs

For comparison of pre and post reactivation strengths, the same run-derived templates were projected onto the preceding and following sleep epochs (neurons were matched across pre, run, and post), and the pre and post strengths during periods of sleep were compared. Thus, each assembly was matched across pre and post. Furthermore, to assess reactivation strength changes over learning, the preceding run epoch’s template (following run epoch for sleep epoch #1, since no preceding run) was projected onto the following sleep epoch, and the average reactivation strength was calculated and compared for early (epochs 1–3) and late (epochs 4–9) sleep epochs during periods of sleep. Epochs prior to a grouped average of 75% correct performance on the W-Track task (Figure 1—figure supplement 1B) were defined as early.

2D assembly maps

To generate assembly rate maps, a procedure similar to computing place fields was used (van de Ven et al., 2016). For each assembly, all assembly activation timestamps during W-Track behavioral sessions that exceeded 5 at positional bins where occupancy was >20ms were used for map generation. 2-D assembly maps were calculated by dividing the binned activation count by the occupancy in each 1 cm2 positional bin and smoothed with a Gaussian kernel (σ=2).

Spatial information of assembly fields

Spatial information (bits/event activation) was computed using the following formula:

Ispike=∑i=1Npiλiλlog2λiλ

where pi is the occupancy probability of bin i, λi is the mean assembly activation in bin i, and λ is the occupancy-weighted mean activation. Spatial information is reported in bits per activation event. The significance of each assembly’s spatial information was assessed by a per-assembly circular-shift procedure. The full assembly activation time series was circularly shifted by a random offset >20 s, the activation shuffled 2D assembly field was computed, and the spatial information was calculated. This procedure was repeated 1000 times per assembly to generate a null distribution. The observed spatial information was expressed as a z-score relative to the assembly-specific shuffle distribution, and a one-tailed permutation p-value was computed as p = (1 + #[Ishuffled ≥ IObserved]) / (1+1000). Assemblies with p<0.05 (z>1.65) were classified as showing significant spatial structure (Figure 6H and I).

2-D rate maps

2-D spatial rate maps during active behavior on the W-Track were calculated during periods of high mobility (>5 cm/s; all SWR times excluded) at positional bins where occupancy was >20 ms. Then, the 2-D spatial maps were calculated by dividing the binned spike count by the occupancy in each 1 cm2 positional bin and smoothed with a Gaussian kernel (σ=2).

Calculation of a single cofiring metric and separation into populations of high and low cofiring neurons

Since we wanted to relate the above cofiring metric to other measures, we needed to obtain a single value cofiring metric for each neuron. To do this, we averaged the cofiring values across all pairings for a neuron (e.g. 1 CA1 neuron paired with all PFC neurons) and reported it as the cell’s cofiring. Furthermore, since we observed a bimodal distribution of CA1-PFC cofiring values in REM sleep, we split the population based on whether the average (across all cell-cell combinations) cofiring value or correlation coefficient of a cell was above or below 0 (High cofiring >0; low cofiring <0). Also, since the NREM ripple cofiring distribution was unimodal, we additionally split the populations using the mean of the average cofiring or correlation coefficient distributions (High cofiring >mean; low cofiring <mean).

Change in firing rate across sleep

To assess differences in firing rate changes in CA1 and PFC over the course of sleep, populations were separated into low or high cross-regional cofiring units based on the z-scored cofiring metric above as well as the Pearson correlation between spike count vectors across all REM HFO chains or independent NREM ripples. Here, we are only including independent ripples for NREM analysis since we have shown that CA1-PFC cofiring is elevated when PFC ripples are coordinated with CA1 SWRs (Shin and Jadhav, 2024). Exclusion of coordinated ripple events allows us to assess the influence of strictly PFC ripples in NREM on firing rate changes. For each cell, the change in the firing rate from the first bout to the last bout of NREM sleep was calculated and compared for low and high coactive cells (Figure 7F, Figure 7—figure supplement 2B–E). Only epochs that had first and last NREM bouts exceeding 30 s in length were included for analysis.

REM theta phase shifting

To assess shifts in the preferred CA1 theta phase of CA1 neurons between run sessions and REM sleep, the above procedure was performed, and only neurons that exhibited significant phase-locking during both run and REM sleep were included for analysis. Neurons that had a phase shift of >90° from the preceding run session to the subsequent REM sleep epoch were categorized as REM-shifting.

Single unit PFC HFO bursting

To determine HFO burst probability for each cell, the number of HFOs where burst spikes were detected (at least 2 spikes <6ms apart Harris et al., 2001) was divided by the total number of events.

Network model

Computational simulations were done using a simplified model of a small cortical network. The model contains 1600 excitatory pyramidal and 400 inhibitory interneuron modified HH neurons laid out in a topographical grid. Connections are random, with probability of connection decaying with spatial distance on the grid. Connection probabilities were drawn from known layer 2/3 cortical connectivity rates of pyramidal cells and somatostatin (SST) interneurons. Synaptic weights were assigned randomly using log-normal distributions. Inhibitory cells in the model are regular-spiking with a voltage equation identical to pyramidal cells.

The voltage equation for the neuron model is:

CMdVdt=−m∞3hgNa(V−ENa)−ngKdr(V−EK)−sgKs(V−EK)−gleak(V−Eleak)+Ibase+Inoise−Isyn+Istim

where gx are the maximal conductances, and each of the activation variables h, n, and s have a differential equation of the form:

dydt=1τ(V)(y∞−y)

where y∞ are the steady state values for each of the gating variables, each a function of voltage (Stiefel et al., 2009).

For every presynaptic spike, the synaptic conductance between neurons is calculated as the difference between two decaying exponentials. Model pyramidal neurons and interneurons both elicit fast, large-amplitude synaptic inputs in connected postsynaptic cells modeled after glutamate and GABA-A receptors, respectively. Interneurons also elicit a slow, low-amplitude GABA-B-like conductance in connected neurons. In addition, a white noise current Inoise is applied to each neuron separately to introduce random variation.

Both pyramidal neurons and interneurons in the model feature a slow potassium conductance, gKs, that is active at neuronal resting potential and modeled after the slow potassium M-current (Satchell et al., 2025; Delmas and Brown, 2005). The biological effect of ACh on the M-current through the M1 muscarinic ACh receptor is recreated in the model by raising or lowering gKs to simulate decreases or increases, respectively, in ACh concentration. As in biological neurons, lowering gKs shifts model neurons from Type 2 excitability (resonators) to Type 1 excitability (integrators), raising their excitability but lowering their tendency to synchronize spiking when connected (Stiefel et al., 2008). The result for a connected network of neurons is, in general, lower sensitivity to inputs and an increased propensity for synchronous spiking in simulations modeling a low ACh condition, but increased sensitivity and less correlated spiking in simulations modeling a high ACh condition, as is observed experimentally (McCormick et al., 1993; Meir et al., 2018). Model simulations of NREM sleep use gKs=1.5 for all neurons, while simulations of REM sleep use gKs=0. Ibase is used to slightly hyperpolarize interneurons during REM to maintain balanced excitation and inhibition.

Each external stimulus input to the network Istim is modeled as a depolarizing current injection to a random 30% of neurons. For the NREM brief stimuli, a Gaussian profile that peaks at the midpoint of the stimulus was used, with duration drawn from a normal distribution with a mean of 70 ms and standard deviation 20 ms. Stimulus peak amplitude was calculated similarly, with a mean of 0.75 and SD of 0.3. For REM theta stimuli, the input was a sinusoid with frequency drawn from normal distribution with mean 8.5 Hz and 3 SDs within 7–10 Hz range. Stimulus amplitude was convolved with a Gaussian profile peaking at the midpoint of the stimulus duration. Negative amplitude phases of the sinusoid were set to zero. For REM stimulus duration, mean 3000 ms and standard deviation 250 ms duration were used. Peak amplitude parameters were identical to NREM. No negative durations or amplitudes were allowed and set to zero if drawn from the distribution. Triggering of all stimulus events was determined by a Bernoulli process with a success rate of on average 1 Hz, with overlapping stimuli not permitted. For theta stimuli, a minimum delay of 1000 ms was enforced after the end of each stimulus to allow the network time to return to baseline activity.

Quantification and statistical analysis

All statistical analyses were performed with MATLAB (MathWorks, Natick, MA; R2022a) functions or custom scripts. No specific analysis was performed to predetermine sample size, but the number used in this study is similar to or larger than typically used. Nonparametric, two-tailed statistical tests (Wilcoxon rank-sum [WRS] and signed-rank [WSR], for unpaired and paired data, respectively) were used throughout the paper unless noted. p<0.05 was considered the cutoff for statistical significance. Boxplots show the median (center line), interquartile range (box), and the furthest points not considered outliers (whiskers). Unless otherwise noted, values and errors in the text denote means ± SEM. p-values in the text are reported as follows: p>0.05, *p<0.05, **p<0.01, ***<0.001.

Data availability

Data underlying these results are available in the NWB (Neurodata Without Borders) format on DANDI Archive (ID#000978). Code to replicate these results are available on our lab GitHub (https://github.com/JadhavLab/REMHFOs; copy archived at Shin, 2026).

The following previously published data sets were used

References

    1. de la Prida LM
    (2020) Potential factors influencing replay across CA1 during sharp-wave ripples
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190236.
    https://doi.org/10.1098/rstb.2019.0236

Article and author information

Author details

  1. Justin D Shin

    Neuroscience Program, Department of Psychology, and Volen National Center for Complex Systems, Brandeis University, Waltham, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7959-7772

  2. Michael Satchell

    Neuroscience Program, Department of Psychology, and Volen National Center for Complex Systems, Brandeis University, Waltham, United States
    Contribution
    Software, Formal analysis, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7932-5697

  3. Paul Miller

    Neuroscience Program, Department of Psychology, and Volen National Center for Complex Systems, Brandeis University, Waltham, United States
    Contribution
    Supervision, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9280-000X

  4. Shantanu P Jadhav

    Neuroscience Program, Department of Psychology, and Volen National Center for Complex Systems, Brandeis University, Waltham, United States
    Contribution
    Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    shantanu@brandeis.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5821-0551

Funding

National Institute of Mental Health (R01MH112661)

  • Shantanu P Jadhav

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by the National Institutes of Mental Health (R01MH112661 to SPJ).

Ethics

All experimental procedures were approved by the Institutional Animal Care and Use Committee at Brandeis University (Protocol #24001-A) and conformed to US National Institutes of Health guidelines.

Version history

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.110795. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2026, Shin et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 532
    views
  • 28
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Justin D Shin
  2. Michael Satchell
  3. Paul Miller
  4. Shantanu P Jadhav
(2026)
REM sleep prefrontal high-frequency oscillation chains mediate distinct cortical – hippocampal reactivation patterns compared to NREM sleep
eLife 15:RP110795.
https://doi.org/10.7554/eLife.110795.3

Share this article

https://doi.org/10.7554/eLife.110795

Categories and tags

Research organism