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.^1,2 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.^3–5 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.⁶

Primary models of systems memory consolidation, including the standard two-stage theory of consolidation or trace transformation theory,^7–10 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.^11–14 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, is still unclear.

REM sleep is characterized by high theta-delta ratio,^15,16 and 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,¹⁷ and corresponding phasic theta activity bursts during REM are often coupled with elevated gamma and cortical high-frequency oscillation (HFO) power, potentially organizing neural activity in distributed circuits for mnemonic processing.^18–25 Furthermore, REM sleep is associated with a unique neuromodulatory landscape that is conducive for the induction of plasticity within and across circuits.²⁶ 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,^27–29 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.²⁸ A related study had similar findings, suggesting that these ripples reflect intracortical processing that occurs independently of hippocampal input, allowing for interference-free consolidation of memory representations within cortical circuits.³⁰

Previous studies have established cortical HFOs (in a similar frequency range to high-frequency NREM ripples) along with gamma oscillations in REM sleep – REM HFOs have been detected in different frequency ranges and sometimes reported with different nomenclature.^19,21,22,31 Recent studies also note the possibility of spurious spectral detection in high frequency ranges.³² Monitoring spiking activity simultaneously with LFP/EEG activity can therefore provide key confirmation about the existence and impact of high-frequency oscillations 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 a 24-hour cycle, but without the context of a behavioral task.^33–35 Studies have also reported prominent phenomena such as hippocampal REM-theta phase reversing cells³⁶ and hippocampal reactivation,³⁷ using behavioral tasks and examining sleep periods in post-task rest sessions. 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.³⁸

An explicit examination of REM HFO events, and how they impact cortical and hippocampal activity during consolidation of recent experiences is not known, despite the similarity of these high-frequency oscillations 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/HFO activity patterns in NREM and REM sleep.

Results

High-frequency events in PFC during REM sleep

We simultaneously recorded local field potential (LFP) and neuronal populations in CA1 (n = 1,468) and PFC (n = 1,151) from animals (n = 10 rats) during REM sleep (Figure S1A, Table S1, 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^39,40 (Figures S1A,B). Activity was continuously acquired across multiple run (8 epochs) and sleep (9 epochs) sessions during learning.⁴⁰ We separated NREM and REM sleep based on theta-to-delta (TD) ratio,^15,16 and find coherent shifts in TD ratio across CA1 and PFC at the onset and offset of REM sleep (Figures 1A and S1B,C), as well as higher movement and intracranial EMG during wake as compared to REM sleep, validating the behavioral state designations (Figures S1C-H). To investigate the high-frequency component of phasic burst activity in PFC, cortical ripples and HFOs during NREM and REM sleep respectively were detected as transient events in the 150-250 Hz range, as previously described²⁸ (Figures 1B,C and S2A,B), with similar criteria used to detect cortical events in NREM and REM states. Cortical event rates during REM were higher than NREM (Figure 1E), and REM HFOs were associated with elevated theta and gamma power,^21,31 in contrast to spindle- and slow oscillation-coupled ripple events in NREM^28,41 (Figure 1D). Interestingly, examining the difference in autocorrelation peaks and inter-event-intervals (IEI), we find that REM HFOs often recur at intervals of ∼130 ms (Figures 1F,G) as HFO “chains”, which is consistent with theta frequency modulation and elevated theta power during these events (Figure 1D). REM HFOs are associated with fast, theta timescale fluctuations in PFC population/multiunit activity (MUA) during periods of locally reduced activity, as evidenced by the overall decrease in population activity during these events (Figure 1H). This phenomenon was also present when using different parameters for HFO detection or REM sleep designation (Figures S1D-G). Conversely, NREM cortical ripples were associated with a singular burst of activity spanning several hundred milliseconds²⁸ (Figures 1I and S1C). 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 1.

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.⁴² 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.⁴² Furthermore, previous work has proposed that phasic REM sleep may support cortical-hippocampal dialogue through this theta frequency-based coordination.⁴³ Phasic REM is characterized by higher theta power and frequency as compared to tonic REM in both the neocortex and hippocampus.^18,19 Thus, we detected putative bouts of phasic REM as periods of elevated theta power (Figure S3). 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 tonic and phasic substages, respectively, and find that high theta power is also associated with higher theta frequency, based on lower inter-peak intervals (Figures S3A,B). The rate of REM HFOs, as well as chaining of events, was elevated during periods of high theta power (Figure S3C), indicating that REM HFOs preferentially occur during bouts of putative phasic REM sleep. We restricted further analyses for quantifying gamma and HFO coupling to theta oscillations to these phasic bouts of high theta power. To further validate the coupling of theta, gamma, and HFOs (Figure 1D), we quantified PAC between theta phase and gamma to HFO frequency (40-250 Hz) amplitudes (Figure 2A). 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 (Figures S3D-I). Additionally, we quantified nesting of gamma and HFO to theta oscillations (Methods) and confirmed the presence of separate gamma and HFO frequency components that are tightly coupled, with gamma events leading HFOs (Figures 2B and S3F,G). Furthermore, the phase locking preferences of PFC neurons to these two oscillatory components differ (Figures 2C and S3J), 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 2D), suggesting that increased coherence is associated with HFO generation in PFC during REM sleep.

Figure 2.

Since we observed short-latency HFO recurrence during REM (Figures 1F,G), 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.^44–47 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.⁴⁵ Similar to this study, we find spindle train events in PFC during NREM sleep (Figures S4A-D), supporting the claim that PFC exhibits 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 2E) and similar to spindle trains, predominantly occurred in doublets (Figures 2F and S4C,E). REM HFOs exhibited higher consistency of phase locking to the theta oscillation as compared to NREM ripples (Figure 2G). Furthermore, REM HFO chains were strongly phase modulated by theta in both PFC and CA1 and had similar properties to isolated HFOs (Figures 2H and S4H). 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 1H) when coupled with HFOs, further demonstrating strong temporal coupling (Figures 2I and S3K,L). Corroborating the strong phase locking of chain HFOs to theta oscillations in PFC and CA1 (Figure 2H), we find that hippocampal-prefrontal coherence is elevated during chain HFOs as compared to isolated events (Figure 2J), with CA1 leading PFC (Figure 2K), highlighting a mechanism through which chains of HFOs reflect hippocampal-prefrontal coordination to potentially support REM sleep-mediated memory consolidation.

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

Since we observed theta modulated population activity aligned to REM HFOs that differed from activity aligned to NREM ripples (Figures 1H,I), 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 3A). 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 (Figures 3A,B). Furthermore, PFC neuronal activity was able to distinguish NREM and REM events, indicating differences in content (Figure 3C). In line with elevated CA1-PFC theta coherence surrounding REM HFOs (Figure 2C), we found that a subset of CA1 neurons were modulated by these events (Figure 3D), indicating enhanced interregional coactivity, potentially delineating periods of CA1-PFC ensemble reactivation to support memory consolidation. Overall cortical coactivity and assembly strength were reduced during REM HFOs relative to NREM ripples (Figure 3E), consistent with activity suppression during HFO chains (Figures 3F-G), which can gate 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, compared to both isolated and NREM events (Figures 3F and S4F-H), and this phasic population response was largely due to the association with HFOs and not necessarily periods of high theta power (Figure S4I). Interestingly, only HFO chains were uniquely associated with a decrease in overall local population activity (Figures 3F,G and S4), which is in contrast to the strong excitation seen in PFC and CA1 during NREM PFC ripples and SWRs, respectively.²⁸ Since REM HFOs predominantly occur in doublets (Figure S4E), sublinear integration may underlie the population suppression during HFO chains and can reflect privileged temporal windows for selective PFC ensemble reactivation and strengthening.

Figure 3.

PFC ensemble activity is sequentially organized by REM HFO chains

Since we observed a broader distribution of activity (Figure 3B) as well as multiple peaks in population activity during REM HFO chains (Figure 3F), we next examined, in more detail, how activity is structured around these events. Previous studies have shown that the refractoriness of certain oscillatory events, such as spindles and SWRs, elicits differential population responses.^44,45,47,48 Thus, we first examined population dynamics based on the temporal separation of REM HFOs and found that HFO adjacency determines the similarity of the reactivation content in PFC across HFOs (Figures 4A,B and S5D). Reactivation similarity as a function of Inter-Event-Interval (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 in chains show a high degree of similarity as compared to adjacent isolated events. Furthermore, rank order analysis on the content of concatenated HFO chains demonstrated a preserved sequential firing of PFC neurons across events (Figures 4C and S5D), suggesting a process through which HFO chains organize PFC activity in a stereotypic manner that allows for sequential reactivation of distinct assemblies.

Figure 4.

To test this sequential reactivation, we next assessed PFC assemblies by time-locking PFC reactivation strength to NREM and REM events. Similar to single unit modulation, the timing of PFC assembly reactivation aligned to NREM ripples was largely restricted to 200-500 ms after ripple onset, depending on the ripple type examined (Figures 4D and S5J). 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 S5F), possibly facilitated by strong phasic gamma surrounding HFOs (Figures 2I and S3K). We found that unlike NREM ripples, where assembly reactivation is synchronous and consistently highest at the onset of events, PFC assemblies were sequentially organized surrounding REM HFOs, similar to single unit modulation (Figures 3A,B and 4D). 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 4E). To confirm the sequential activity during REM HFOs, we quantified the consistency of the assembly reactivation using a cross-validated approach (see Methods). Consistent with the high temporal coordination between chain events and gamma (Figure 2H), we found the greatest degree of sequence similarity during REM HFO chains and gamma events compared to shuffled data, confirming the tendency for PFC assemblies to be reactivated in a stereotyped, sequential manner as compared to isolated REM HFOs and NREM ripples (Figures 4F-H and S5G-J).

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

Theta coherence and phase amplitude coupling are two prominent modes of interregional coordination thought to underlie mnemonic processing.^42,49 Since we observed an overall increase in CA1-PFC theta coherence and gamma power during HFO chains (Figures 2I,J), as well as sparse modulation of CA1 activity during REM HFOs (Figure 3D), we reasoned that there may be specific populations of CA1 neurons that are engaged by these events. Examining CA1-CA1 cofiring during these events, we observed a higher degree of spatial correlation for high cofiring pairs, suggesting HFO specific reactivation of behaviorally relevant hippocampal assemblies to support memory consolidation¹² (Figures 5A,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 S6A). Next, we evaluated interregional activity and observed a bimodal distribution of CA1-PFC cofiring during REM HFOs, indicating a subpopulation of CA1 cells that are highly engaged with PFC neurons (Figure 5C). Examining low and high cofiring CA1 neurons during REM HFOs showed that indeed high cofiring CA1 neurons exhibited elevated activity during chained events as compared to low cofiring CA1 neurons (Figure 5D).

Figure 5.

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 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 (Figures 5E and S6C,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 ripples and excitation during coordinated SWRs in NREM sleep, as shown in a previous study²⁸ (Figures S6C,E), we did not find correlations between CA1 REM HFO cofiring and CA1 SWR reactivation (Figure S6F), indicating that these CA1 sub-populations are specifically engaged by high-frequency PFC events, 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.^33–35 Interestingly, we observed that these 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 (Figures 5F, S6B, and S7A-E), consistent with the role of REM sleep in hippocampal firing rate diversification.⁵⁰ This effect was hippocampal specific; there was no differential regulation of PFC firing rates over the course of sleep (Figure S7F). Furthermore, the degree of cofiring with PFC during REM HFOs predicted CA1 firing rate changes during sleep, and this relationship was absent for NREM PFC ripples (Figures 5F and S7E). Splitting the population in half based on mean firing rate changes further confirmed this relationship (Figure 5G), suggesting a REM specific role of PFC HFOs in modifying hippocampal excitability.

The dorsal CA1 pyramidal cell layer can be separated into subpopulations that are molecularly defined, have differential engagement with hippocampal SWRs and cortical areas,^51–54 and employ complementary spatial codes for effective coding of environments.⁵⁵ Notably, a subset of CA1 pyramidal cells preferentially shift their preferred theta phase from run to REM sleep¹⁵ – a phenomenon that is experience dependent.³⁶ 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°.¹⁵ 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 (Figures S7G-I), as previously reported.⁵² Similar to high cofiring CA1 neurons, these REM-shifting cells had higher coactivity with PFC than non-shifting cells (Figure 5H), 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 (Figures 5I and S7I). 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 (Figures 5J,K and S7K). 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.

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 6A). 1,600 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.^56–58 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.⁵⁹ 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.^60,61 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.^62–66 In the model, modulation of this M-type current also changes the excitatory-inhibitory balance of the network, supporting higher average interneuron rates (Figure S8A).

Figure 6.

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 6B,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 6C,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 S8D). 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 6E). 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 3E) 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 6F). 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 200 ms between MUA peaks, Figure 2). 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 6G), with the difference especially clear in the pyramidal population of the model (Figure S8H). 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 7A). 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 7B). 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 S8G). 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 3A,E, 4D,E, and 7B). 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 7C). In contrast, during REM stimuli, the number of active stimulated neurons increased, but the effect on non-stimulated neurons was suppressive (Figure 7D). 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 S8I,J). These findings nicely align with experiments showing that ACh both suppresses the spread of cortical recurrent excitation and drives inhibition.^67–70 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 7.

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,⁷¹ 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,^11–14 leaving REM sleep poorly understood in comparison.

Here, we show that phasic theta bursts in PFC during REM coincide with chains of local HFOs and elevated gamma power. These chains are associated with increased measures of functional connectivity between PFC and CA1, organize sparse, sequential reactivation in PFC over extended time periods, and recruit a sub-population of CA1 neurons for reactivation that exhibit differential regulation of activity over the course of sleep. This is in contrast to NREM sleep, where 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/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,²⁸ 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.⁷² This study focused solely on hippocampal assemblies, and the results are in line with our finding of the potential roles of PFC ripples/HFOs in NREM and REM sleep. 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.²⁸ This bidirectional modulation of hippocampal representations during PFC and CA1 ripples 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.⁵ We speculate that PFC REM HFO chains enable selective linkages between the sparsely reactivated PFC ensembles in coordination with a subset of CA1 neurons.

A limitation of our study is the inability to investigate the relationship of the REM reactivation phenomena 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,^73–76 but reports of REM reactivation are rare.³⁷ Investigation of the role of REM reactivation in learning and memory consolidation will need utilization of tasks known to require REM and NREM states, such as inference or emotional memories.^38,77–79 Another limitation is that behavioral task learning and interleaved sleep sessions were conducted in a few hours within a single day, and not across the circadian cycle.

Cortical high-frequency oscillations (∼120-160 Hz) in REM have been reported before, especially in parietal cortex, although they have been variously termed fast gamma or HFOs (high-frequency oscillations),^19,21,22,31 with previous reports of strong coupling with theta oscillations in REM.³¹ 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.²² 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,^{23,25,27–29,80} although this is likely a separate phenomenon from hippocampal sharp-wave ripples (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,²⁸ while in REM, we show 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 oscillations in phasic REM, and are reminiscent of hippocampal theta sequences,⁸¹ but may have different underlying mechanisms.

External drive can lead to the generation of PFC ripples and synchronized activity in NREM, whereas high acetylcholine (ACh) and high inhibitory tone in REM can lead to local suppression in response to external drive.^70,82,83 Notably, a recent study demonstrated that differential cholinergic modulation, together with the sequential architecture of sleep (NREM to REM transition), supports memory consolidation through the complementary shaping of ensemble activity during NREM and REM sleep.⁶¹ 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 ripples/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, potentially reflecting the need for coordinated activation of distributed circuits to elicit HFO-associated phasic events. One possible mechanism driving this phenomenon could involve ponto-geniculo-occipital (PGO) waves, which are prominent during phasic REM and are known to play a role in modulating hippocampal plasticity during REM sleep.⁸⁴ These waves are also coordinated with both SWRs in NREM and theta states in REM,^84,85 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 coherence between the two areas.⁸⁶ Indeed, a number of studies have shown that ACh boosts theta oscillations.⁸⁷ Furthermore, ACh induces decorrelation in cortical circuits, which improves cortical representations of relevant stimuli,⁸⁸ and our finding that REM HFOs are associated with an overall decrease in PFC population activity, as well as reactivation strength and cofiring, in contrast to NREM ripples, indicates a role for REM HFOs in amplifying local signal to noise to link cortically encoded information with specific populations of hippocampal neurons. Lastly, cholinergic input to the hippocampus also suppresses SWR generation,^89,90 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 support learning and memory consolidation.

Taken together, these results 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 highlight 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, suggests 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.

Method

Behavior

Ten animals were trained on a novel W-maze in a single day, as previously described.^28,40,91. This single-day W-maze alternation learning task is SWR-dependent, and requires hippocampal-prefrontal interactions.^39,92,93 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 between the hours of 9 AM and 6PM. 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 hours 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.

Surgical implant and electrophysiology

Surgical implantation procedures were as previously described.⁹¹ Animals (n = 10) were implanted with a microdrive array containing 30-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. 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.⁹¹ 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 hours 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 and running speed were recorded with an overhead color CCD camera (30 fps) and tracked by color LEDs affixed to the headstage.

Histological verification of recording sites

At the completion of each experiment, electrolytic lesions were made by sending current (30 uA) down 3 electrodes per tetrode for 3-4 seconds each. 24 hours 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 hours 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 uM. Slices were Nissl stained and imaged at 4x 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,⁹⁴ 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.^15,16 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 ratio of theta to delta (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.^95,96 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 (EMG)

To estimate EMG from the neck, jaw, and face muscles from intracranial LFP signals, a previously established method was used.^35,97 Briefly, we extracted the 300-600 Hz filtered signal (Butterworth filter at 300 – 600 Hz; filter shoulders from 275 – 625 Hz) and calculated the pairwise Pearson correlations in each 500 ms 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 S1D, 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.⁹⁷

LFP event detection

PFC ripple/HFO detection: High-frequency events were detected during immobility periods (<4 cm/s) as previously described.¹⁶ 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 were 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.⁴⁰ 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/)⁹⁸ 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.⁹⁹ 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 were 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 900 ms were included as high theta power periods.

PFC spindle detection: Spindle events were detected as previously described.⁴⁵ The mean LFP across all candidate PFC electrodes was calculated and the resulting signal was bandpass filtered between 10-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.⁴⁵

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 second, 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. 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. 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 multiunit activity 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 multiunit activity was split into quartiles based on theta power and compared.

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

Since it has been previously shown that tonic and phasic substages of REM sleep can be segregated based on CA1 theta frequency,¹⁵ we calculated the inter-peak intervals of the theta oscillation to determine whether bouts of low and high theta power may correspond to tonic and phasic substages, respectively. We bandpass filtered and averaged the LFPs from candidate CA1 tetrodes between 5-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

Phase-amplitude coupling was assessed for a range of frequency combinations using a previously established approach.¹⁰⁰ Since the rate of REM HFOs was highest during bouts of high theta power, we focused on these periods for assessing phase-amplitude coupling. 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-12 Hz and 40-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-250 Hz was calculated in each phase bin from 2-12 Hz and visualized.

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,¹⁰⁰ 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^101,102 in two different frequency bands (40-100 Hz for gamma and 150-250 Hz for HFO) relative to the theta oscillation. 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 S3C). 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 100 ms 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.

Ripple/HFO aligned modulation

Single unit ripple/HFO modulation in CA1 or PFC was calculated as previously described.^28,103 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) 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.

Theta coherence

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. Coherograms were averaged over CA1-PFC tetrode pairs. For quantification and comparison, coherence was measured as the peak coherence between 6-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.⁹⁸ Only epochs with at least 20 events were included for quantification of peak 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 phase slope index (PSI),¹⁰⁴ which was calculated with the FieldTrip analysis toolbox.¹⁰⁵ 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.

Phase-locking

Phase-locking during active behavior on the W-Track and during REM sleep was quantified using previously developed methods.¹⁰⁶ 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 were 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), spikes were restricted to CA1 theta periods where the animals’ speed was >5 cm/s and were assigned a phase between 0-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 was 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 U² test. Subsequently, for each theta modulated cell, phase locking strength to the theta oscillation was measured using pairwise phase consistency (PPC).¹⁰⁷ 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.

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 <6 ms apart¹⁰⁸) was divided by the total number of events.

Detection of population suppression

To validate the suppression of PFC multiunit activity surrounding REM PFC HFOs, periods of decreased PFC activity were separately extracted, and the probability of REM HFOs surrounding these events was calculated. First, multiunit activity 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 100 ms bins. To obtain the 99 percent 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 (IEI) 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. 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

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

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 to 1 where 1 indicates identical vectors and −1 indicates that the vectors have an angular difference of 180 degrees. Actual spike counts were used for cosine similarity calculation.

Assembly detection

To detect cell assemblies in PFC populations, a previously described method was used.^109,110 Significant cofiring 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 20 ms 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:

Where p_i is the i^th principal component of C and λ_i is the corresponding eigenvalue. Note that 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 , 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, N_A, represents the number of significant patterns found in the data. These principal components were then projected back onto the original z-scored spike data:

where P_SIGN is the nXN_A matrix with N_A columns. Independent component analysis (ICA, RobustICA)¹¹¹ was then used to identify patterns that were maximally independent from each other. ICA was applied to the matrix Z_PROJ to find an N_AXN_A unmixing matrix, W, which was used to obtain the weights of each cell in each assembly.

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, the expression strength of each pattern over time was given by:

R_k(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 P_k is the projection matrix for assembly k, constructed by taking the outer product of its weight vector v_k. 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. Only sleep sessions with a preceding run session to use as a template were included for reactivation analysis (sleep sessions 2-9).

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.^112,113 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 or 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.

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.^114,115 The population activity during each HFO in a chained event was concatenated (inter-event activity excluded), and for each concatenated event with 25 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 seconds such that each HFO within the chain was shifted the same value (1000 time shifts), thus preserving the inter-event-intervals (ie. The entire chain was shifted coherently). Additionally, the entire shifted chain must have stayed within a REM sleep bout. The above rank-order procedure was performed, and the real correlation was compared to the pseudo-chain null distribution.

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 cm² positional bin and smoothed with a Gaussian kernel (σ = 2).

Ripple and HFO cofiring

To quantify ripple/HFO cofiring during events, a pairwise coactivity z-score was calculated as previously described.¹¹⁶ Briefly, this coactivity measure estimates how likely two cells will fire together during events and is normalized by the incidence of spiking independent of the other cell:

where N is the number of events, n_a and n_b are the number of events where cell a and cell b were active, respectively, and n_ab is the number of events where both cells were active. The same procedure was used to calculate cofiring between pairs in the model simulations.

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 coactivity 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 coactivity is elevated when PFC ripples are coordinated with CA1 SWRs.²⁸ Exclusion of coordinated ripple events allows us to assess the influence of strictly PFC ripples in NREM on firing rate changes. Since we observed a bimodal distribution of CA1-PFC cofiring values, 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. Since the NREM ripple cofiring distribution was unimodal, we additionally split the populations using the mean of the average cofiring or correlation coefficient distributions. 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. Only epochs that had first and last NREM bouts exceeding 30 s in length were included for analysis.

Network model

Computational simulations were done using a simplified model of a small cortical network. The model contains 1,600 excitatory pyramidal and 400 inhibitory interneuron modified Hodgkin-Huxley 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:

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

Where y_∞ are the steady state values for each of the gating variables, each a function of voltage.⁵⁹

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 I_noise is applied to each neuron separately to introduce random variation.

Both pyramidal neurons and interneurons in the model feature a slow potassium conductance, g_Ks, that is active at neuronal resting potential and modeled after the slow potassium M-current.^61,117 The biological effect of acetylcholine (ACh) on the M-current through the M1 muscarinic ACh receptor is recreated in the model by raising or lowering g_Ks to simulate decreases or increases, respectively, in ACh concentration. As in biological neurons, lowering g_Ks 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.⁶³ 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.^118,119 Model simulations of NREM sleep use g_Ks = 1.5 for all neurons, while simulations of REM sleep use g_Ks = 0. I_base is used to slightly hyperpolarize interneurons during REM to maintain balanced excitation and inhibition.

Each external stimulus input to the network I_stim 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 standard deviation 0.3. For REM theta stimuli, the input was a sinusoid with frequency drawn from normal distribution with mean 8.5 Hz and 3 standard deviations 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.

Experimental model and subject details

All experimental procedures were approved by the Institutional Animal Care and Use Committee at Brandeis University 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.

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, ***p < 0.001.

Supplementary material

Table S1.

Figure S1, related to Figure 1.

Figure S2, Related to Figure 1.

Figure S3, Related to Figure 2.

Figure S4, Related to Figure 2.

Figure S5, Related to Figures 3 and 4.

Figure S6, Related to Figure 5.

Figure S7, Related to Figure 5.

Figure S8, Related to Figures 6 and 7.

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 will be available on our lab GitHub (TBD) upon acceptance. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Acknowledgements

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

Additional information

Author contributions

J.D.S. and S.P.J. conceived the study. J.D.S. performed the experiments. J.D.S. processed and analyzed the data. M.S. conceptualized the model and P.M. advised on model implementation. J.D.S., M.S., and S.P.J. wrote the manuscript. S.P.J. provided supervision and funding for all aspects of the study.

Funding

HHS | NIH | National Institute of Mental Health (NIMH) (R01MH112661)

• Shantanu P Jadhav