mPFC spindle cycles organize sparse thalamic activation and recently active CA1 cells during non-REM sleep

Sleep oscillations in the neocortex and hippocampus are critical for the integration of new memories into stable generalized representations in neocortex. However, the role of the thalamus in this process is poorly understood. To determine the thalamic contribution to non-REM oscillations (sharp-wave ripples, SWRs; slow/delta; spindles), we recorded units and local field potentials (LFPs) simultaneously in the limbic thalamus, mPFC, and CA1 in rats. We report that the cycles of neocortical spindles provide a key temporal window that coordinates CA1 SWRs with sparse but consistent activation of thalamic units. Thalamic units were phase-locked to delta and spindles in mPFC, and fired at consistent lags with other thalamic units within spindles, while CA1 units that were active during spatial exploration were engaged in SWR-coupled spindles after behavior. The sparse thalamic firing could promote an incremental integration of recently acquired memory traces into neocortical schemas through the interleaved activation of thalamocortical cells.


Introduction
Two complementary learning systems acquire and process episodic memories in order to learn adaptive models of the world that generalize beyond specific experiences and continue to incorporate information through the lifetime of an animal [1][2][3][4] . The hippocampus can form memory traces rapidly; however, over time, interactions between the hippocampus and neocortex are thought to extract statistical regularities from the organism's experiences during wakefulness, and update neocortical models or schemas to increase their adaptive value 5,6 . During non-REM sleep, the interactions between the hippocampal and neocortical cell ensembles that underlie these processes are reflected in nested oscillations in the local field potentials (LFPs). The up-states of the slow oscillation (1-4 Hz) provide a processing window in which faster dynamics lead to synaptic plasticity in relevant cell populations. In the neocortex, some cells are preferentially active at the start or at the end of upstates 7,8 , and in the hippocampal CA1 region, sharpwave ripple oscillations (SWRs, 100-275 Hz) occur primarily during up-states 9 . Hippocampal SWRs are brief (up to ~100 ms) oscillations that correlate with the reactivation of spatial representations by hippocampal cell ensembles 10,11 which are thought to contribute the content in the consolidation of episodic memory [12][13][14] . Spindle oscillations (6)(7)(8)(9)(10)(11)(12)(13)(14) are coordinated with the slow oscillation and with SWRs [15][16][17] . Spindles are generated in the thalamus when the inhibitory cells of the thalamic reticular nucleus (TRN) evoke inhibitory post-synaptic potentials and rebound bursting in cells of the dorsal thalamus, which then entrain their postsynaptic targets in neocortex 18,19 . However, the role of spindles in gating and coordinating multi-region communication remains to be explored at the cellular and systems level.
The coupling between these non-REM oscillations appears to be causally linked to memory consolidation, as suggested by the disruption or enhancement of memory following closed-loop manipulations 20,21 , and by the association of spindle dysfunction with perturbation of memory function in some neuropsychiatric disorders 22,23 . However, we are far from understanding the specific mechanisms and single-cell dynamics involved in the hippocampal-neocortical interaction. In particular, information regarding the thalamic involvement is limited. Evidence from anesthetized animals suggests a correlation of thalamic firing with both the neocortical slow oscillation 24,25 and with hippocampal SWRs 26 , while fMRI results from studies in monkeys indicate an overall decrease in thalamic activity during SWRs 27 . A recent study in naturally sleeping rats 28 , found a decrease in activity in cells of the thalamic mediodorsal (MD) nucleus at the time of CA1 SWRs, particularly during SWRs uncoupled from spindles recorded in the skull EEG. This emerging evidence emphasizes the need to study the thalamus, but we have been lacking a multi-region approach (thalamus, neocortex, hippocampus) that would compare in the same preparation the correlations of cells in each region to the sleep oscillations, and that would investigate the dynamics of the thalamic contribution down to the finer timescales relevant for plasticity.
We sought to address this gap by characterizing the fine temporal correlations of thalamic cells with sleep oscillations through extracellular recordings of single units and LFPs in the medial prefrontal cortex (mPFC), hippocampal CA1 region, and nuclei of the midline thalamus in freely behaving rats. We targeted these brain regions because of their anatomical and functional connectivity and their role in memory processes. The limbic thalamus has widespread connections with neocortical regions, and projections from the midline are particularly dense to mPFC [29][30][31] , a critical region for the recall of consolidated memories 32 . The reuniens nucleus in the midline sends collaterals to CA1 and mPFC, making it a compelling candidate to coordinate mPFC-CA1 activity 33,34 . Recording from these areas, we find evidence for a dual role of the thalamus in hippocampal-neocortical sleep interactions; first, an increase in thalamic firing at the down to upstate transition could influence the dynamics of the slow oscillation, and second, the sparse but consistent firing of thalamic cells within individual spindles suggests an interleaved activation of thalamocortical synapses that may promote a flexible integration of recently acquired memory traces into neocortical schemas during systems consolidation.

Results
We used tetrodes to record extracellularly (single units, LFPs) from the midline thalamus, and from areas involved in recent (CA1) and remote (medial prefrontal cortex, mPFC) memory recall 32,35 . During the 1-5 h recording sessions (in freely behaving rats), the animals remained in a quiet, square-shaped enclosure where they cycled through bouts of sleep and wakefulness; in 6 sessions from 3 rats, the animals also explored a radial maze after the recording of sleep, and were brought back to the sleep area after behavior for comparison of sleep activity before and after spatial exploration. LFPs and 317 units were recorded from 7 rats; sessions from an additional animal in which no units were recorded were also included for LFP analyses for a total of n = 8 rats (Figure 1a displays a subset of tetrode locations in brain sections from one rat; Supplementary Figure  1, parts 1 and 2, document the trajectories in the brain for the electrodes and sessions used for analyses). We selected recording sessions in which there was at least one unit recorded in the thalamus and where reliable detection of sleep events could be performed (slow oscillation and spindle activity in mPFC, sharp-wave ripples, SWRs, in CA1; examples of simultaneously recorded LFP traces and units shown in Figure 1b). The majority of recordings in the thalamus were from units in the reuniens (RE), and ventromedial (VM) nuclei, and some from the mediodorsal (MD), paratenial (PT), rhomboid (RH) and centromedial (CM; detailed population numbers in Supplementary Table 1). All of these nuclei are bidirectionally connected to mPFC, and are functionally related to mPFC and other higher order cognitive regions, including CA1, which also receives direct input from reuniens. Here we will use the term limbic thalamus to refer to the thalamic nuclei that we recorded from. See the Methods section for target anatomic coordinates, procedures to identify tetrode location, and additional inclusion/exclusion criteria for the analyses reported below.
To study the modulation of unit spikes in the three regions (mPFC, limbic thalamus, CA1) with non-REM sleep oscillations, we first identified times with both high delta (1-4 Hz) power in mPFC and high SWR (100-275 Hz) power in CA1 (as described in Methods); an example of such non-REM sleep state detection for a representative session is shown in Figure 1c. Within this state, we detected reference features that reflect non-REM sleep events that are key organizers of memory consolidation: Kcomplexes (KCs, marked by red bars in Figure 1d, mPFC LFP), the troughs of spindle cycles in mPFC (red dots in 1d), and SWRs in CA1 (red bars in Figure 1d, CA1 LFP). These LFP events are markers of cell population dynamics, namely, KCs mark the down-states of the neocortical slow oscillation, while SWRs correspond to the high-frequency oscillations that reflect memory-related reactivation of ensembles of CA1 place cells 11,13,14,36 . Spindles (6)(7)(8)(9)(10)(11)(12)(13)(14) are generated through oscillatory bursting in the thalamic reticular nucleus (TRN), which is thought to then engage cells in the dorsal thalamus and, through them, neocortical networks 18,[37][38][39] . We considered that the time window of each spindle cycle has a similar duration to that of SWRs (up to ~ 100ms) and sought to test if individual spindle cycles represent a window for multi-region coordination; we thus detected individual spindle cycles from the mPFC 6-14 Hz filtered LFP. The detection of individual cycles was also motivated by the occasional observation of few spindle cycles in succession (e.g., only 3-4 cycles following the second detected KC in Figure 1d), and by the broader variability in the morphology of the LFP trace in the spindle band compared to the more clearly discrete and consistent signatures of KCs and SWRs (examples in Figure 1d). KCs, spindle cycles, and SWRs were detected from the mPFC and CA1 LFPs after applying an amplitude threshold to the squared, filtered trace (1-4 Hz for KCs; 6-14 Hz for spindles; 100-275 Hz for SWRs), at a level of three (spindles) and four (KCs, SWRs) standard deviations above the average LFP voltage across sleep periods in each recording session (see Methods for details). Although the spindle band overlaps with the theta band, we limited the detection of spindle cycles (and KCs and SWRs) to resting states; peri-spindle spectrograms confirmed that detected spindle cycles are functionally different from theta (Supplementary Figure 2). Peri-event histograms of the LFP and local spiking activity, triggered on the timestamps at the time of largest amplitude of the detected events (trough for KCs and spindles, peak for SWRs), confirmed that the LFP detection identified events that reflect the modulated activity of local cell populations (Supplementary Figure  3). Quantifications are reported as mean ± standard deviation (sd), unless otherwise indicated. Figure 2a illustrates the temporal and phase relations between non-REM sleep events, with examples from the raw and filtered mPFC and CA1 LFPs. The overlaid filtered LFP traces (which have been scaled for display purposes), show that increases in SWR amplitude occur preferentially in the positive half-cycle of the delta-filtered LFP, and in the negative halves of the spindle-filtered trace (inset). The average distribution of delta and spindle phases (estimated with the Hilbert transform, see Methods) at the time of SWRs is shown in Figure 2b for all sessions. Figure 2c shows the peri-event histograms for non-REM detected population events (SWRs, KCs, spindle cycles) across sessions and rats.
Confirming previous results 9,15-17 , we found that KCs, spindles, and SWRs are temporally and phase-coupled; SWRs preferentially occurred following and preceding KCs, phase-locked to the peaks of the delta-filtered LFP, and to the troughs and ascending phase of spindle cycles (circular mean and standard deviation of delta phase at the time of SWRs = -0.04 ± 0.43 rad; spindle phase = -1.72 ± 0.93 rad; n = 23 sessions, 6 rats; Supplementary Figure 4 shows the delta and spindle phase distributions at the time of SWRs for all sessions included in this analysis).

Thalamic cells decrease their overall firing and increase bursting with CA1 SWRs
We next investigated the changes in the firing patterns of single units in the limbic thalamus with CA1 SWRs.
Most thalamic units presented a reduced average firing rate during rest compared to wakefulness. In 28 of 32 sessions the animals (n = 7) alternated between sleep and active wakefulness (in the other 4 sessions rats were only in wakefulness for less than 5 mins and were not included in this comparison), and the firing rate of thalamic units changed from an average of 10.10 ± 4.26 spikes/s during wakefulness to 7.73 ± 3.09 spikes/s during sleep (n = 40 units; p = 1.31 × 10 -6 , Wilcoxon signed rank; Figure 3a). A decrease in firing rate may be due to the hyperpolarization of thalamic cells. Thalamic cells fire in two modes, burst and tonic, depending on the activation state of a transient conductance that relies on T-type calcium channels 40,41 . Bursts correlate with thalamic hyperpolarization, which is a pre-requisite for Tchannel activation, and they can enhance transmission in thalamocortical synapses and induce synaptic plasticity 42,43 . Therefore, analyzing burst firing can provide clues on the mechanism behind the decrease in firing rate in sleep. Bursting can be detected from extracellular recordings because of its characteristic combination of a relatively long interspike interval followed in rapid succession by a few spikes with brief inter-spike intervals 44 , but the correlation between thalamic bursts and SWRs has not been explored.
We detected bursts in the thalamus as groups of action potentials with less than 5 ms of inter-spike interval preceded by at least 70 ms with no spikes, a criterion based on the dynamics of the T-channels underlying burst firing in the dorsal thalamus 41,44,45 . With this criterion, we found that 25 % of the thalamic cells had at least 5 % of their spikes in bursts during wakefulness, while 77.5 % did during sleep (average percentage of spikes in bursts during wakefulness 4.60 ± 7.97 %, and during sleep 18.65 ± 15.3 %; p = 1.8 × 10 -7 ; Wilcoxon signed rank; n = 40 units, 7 rats; Figure 3a). The increase in bursting was more pronounced around the time of SWRs occurrence in CA1 (25.52 ± 17.29 % of spikes in bursts in a 2 s window around SWRs; p = 1.8 × 10 -7 Kruskall-Wallis; Figure 3a). Because the Tchannels responsible for burst firing need a sustained hyperpolarization for about 50-100 ms to be primed for activation 41,46 , the increase in bursting suggests a hyperpolarization of the majority of cells in the limbic thalamus during sleep, and particularly in correlation with hippocampal SWRs.
Indeed, peri-SWR triggered histograms showed that SWRs occur in a background of decreased thalamic spiking. An example of a peri-SWR histogram for a unit from the RE nucleus is shown in Figure 3b, and the peri-SWR for the population of recorded cells in the three areas is shown in Figure 3c; 55.32 % of units in the thalamus showed a reduction in firing of more than 2 standard deviations (sd) from the null (average -4. 23  showed that the firing in CA1 varied inversely with the firing of most thalamic units (61.7 %, p < 0.001; only one unit in the MD nucleus had a significant positive correlation with SWRs rate). For the units that negatively correlated with SWRs, the models predict an average increase of 20.74 % in the rate of SWRs for every 1 spike/s drop in thalamic firing, suggesting anti-correlated activity in the limbic thalamus and CA1 (in contrast, as expected, the rate of SWRs was positively correlated with the spiking activity recorded in CA1; p value < 5 × 10 -189 in all electrodes used for SWR detection). The modulation of thalamic units with SWRs was observed both in RE and in other nuclei. We observed some of the strongest negative correlations with SWRs in RE (unit in Figure 3b), and the size of the effect was marginally significant in RE compared to non-RE units (p = 0.05). However, there was no difference between the proportion of SWR-correlated units in RE compared to non-RE units (p = 0.159), suggesting that the decrease in thalamic firing with SWRs relates to general excitability changes in the state of thalamocortical networks in the limbic thalamus.

Correlation of thalamic units with mPFC delta and spindle oscillations
Although we found a decrease in average thalamic firing and an increase in bursting that suggests a hyperpolarized thalamus during non-REM sleep, it is unclear how the remaining thalamic spikes relate to the mPFC delta (1-4 Hz) and spindle (6-14 Hz) oscillation bands, and if there are specific temporal windows with selective increase in thalamic firing. To investigate this, we used two approaches: we calculated firing rates with respect to local events detected in the mPFC LFP (KCs, spindle troughs); and we estimated the phasepreference of unit spikes to the delta and spindle-filtered mPFC LFP.
Thalamic rebound following the downstate of the mPFC slow oscillation Figure 4a shows peri-KCs histograms for units in the three regions, as well as a representative example of the temporal relations between KCs, SWRs and thalamic activity in one recording session ( Figure 4b). Units in the thalamus often fired asymmetrically with respect to KCs, with stronger firing in the few hundreds of miliseconds after the mPFC KC compared to before (average of 2.48 standard deviations at 274. 76  Previous evidence 8 suggests that the thalamus may have a specific role in the initiation of cortical upstates. To assess the specificity of thalamic firing to the down to up-state transition, we compared the firing of thalamic units before and after mPFC KCs (from the peri-KCs histograms in Figure 4a). For the population of thalamic cells, 38.71 % increased their activity above 2 sd (from the null distribution) following KCs, i.e., in the down to up-state transition, compared to 9.68 % of the cells with a significant peak preceding KCs, i.e., in the up to down-state transition (19.35 % units had a significant peak both before and after KCs). Instead, CA1 and mPFC activity was more similar preceding and following KCs. SWR probability increased within the 250 ms around mPFC KCs (4.65 sd at 208.15 ± 90 ms before KCs, and 3.02 sd at 266.85 ± 100 ms after KCs; n = 23 sessions). This indicates a delta modulated increase in SWR probability in CA1, with the minimum rate of ripples following briefly after the trough of neocortical activity (10.33 ± 5.88 ms). The modulation of SWRs was paralleled by a modulation of firing in CA1 units (n = 111), which increased their activity both before and after KCs (2.32 ± 1.49 sd at -233.7 ± 146 ms before KCs, and 2.35 ± 1.69 sd at 273 ± 140 ms after KCs); on average, CA1 units presented a trough of activity after KCs (19.53 ± 132 ms). Activity in mPFC (n = 46 units) was modulated in a similar way, with increases in unit firing preceding and following the down-state of the slow oscillation (2.33 standard deviations at -267.30 ± 110.40 ms, before KCs; 2.31 standard deviations at 295.20 ± 120 ms, after KCs).
In other words, while units from the three areas showed no difference in the peak of the histograms in the time window following KCs (p = 0.89, Kruskal-Wallis between regions), thalamic units had significantly lower firing before KCs compared to mPFC and CA1 units (p = 0.007, Kruskal-Wallis). This suggest that thalamic cells contribute spikes primarily as the new cycle of slow oscillation starts, while mPFC and CA1 units also fire at up to down transitions. We used an additional metric to quantify the differential firing of units in the three areas with respect to the time of mPFC KCs. For all the units, we calculated the difference between the maximum firing rate in the 0.5 s windows before and after KCs, normalized to the maximum firing in the two windows. This asymmetry index ranges between -1 and 1, with negative values indicating higher firing before KCs. The distribution of this index (Figure 4c) was biased to positive values for thalamic units, but was uniformly distributed for mPFC and CA1 units (average 0.24 ± 0.22 in thalamus; -0.004 ± 0.47 in mPFC; -0.06 ± 0.5 in CA1), consistent with a selective rebound of thalamic firing in the time window following mPFC KCs relative to thalamic activity preceding the down-state.

Sparse thalamic firing and consistent activation of cell pairs during spindles
We then analyzed the firing at the timescale of mPFC spindles ( Figure 5), and found the strongest modulations in the CA1 multi-unit activity (MUA, all spikes recorded by a given electrode). Peri-spindletrough histograms of MUA and single units are shown in Figure 5a (n = 30 tetrodes in mPFC; n = 25 tetrodes in thalamus; n = 52 tetrodes in CA1; n = 46 mPFC units; n = 31 thalamus units; n = 111 CA1 units). Figure 5b shows an example of the histograms for one session in which a unit in the mediodorsal nucleus of the thalamus was recorded at the same time as units in mPFC and CA1. Across the sessions, 84.62 The auto-correlograms of thalamic units did not present peaks at the time intervals of spindle frequency, i.e., in the time window between 50-150 ms, suggesting that the cells did not fire rhythmically in that frequency. Figure 6a illustrates the sparsity of thalamic cell firing during mPFC spindles; the panel shows a trace of mPFC LFP with several cycles of spindle oscillation (6-14 Hz filtered trace overlaid in black); the MUA recorded by one tetrode in the thalamus is shown in orange and the spikes of a single unit in magenta. The autocorrelograms next to the raw data were calculated with the spikes fired during sleep, and show the lack of spindle modulation in the single unit (also when calculated with only burst spikes), in spite of the MUA modulation (which can also be seen in the spike raster). Similar results were obtained for all thalamic units, regardless of the nucleus, and including the units that showed the strongest rebound following KCs (Supplementary Figure  5). Because cells could be firing rhythmically only when they are bursting, we also calculated autocorrelograms with only burst spikes, as well as only for the spikes that occurred in the 0.5 s time window immediately following KCs (when spindle troughs are largest). Thus, these analyses were more restricted to conditions in which rhythmicity in the spindle band would be expected to be strongest, and still failed to reveal evidence of rhythmic firing in the thalamic cells that we recorded from (Supplementary Figure 5).
The lack of rhythmicity in individual thalamic units, together with the modulation of the MUA, raise the possibility that different groups of thalamic units contribute to activity in different spindle cycles. A prediction from this hypothesis is that pairs of units should be more likely to be correlated in brief timescales than individual units with themselves. We tested this in sessions in which we recorded more than one isolated unit in the thalamus (n = 12 pairs, 3 rats), and found that half of the pairs showed a sleep crosscorrelation peak above 3 standard deviations from the average crosscorrelogram (average 2.86 ± 0.87 sd). An example is shown in Figure  6b, which displays the crosscorrelation of two thalamic units during sleep; the plots show that the spikes of 'unit 1' follow reliably in a window of ~50 ms after those of 'unit 2' (reference), while neither of the units show a strong auto-correlation in that period ( Figure 6c). The coordinated activity of cell populations drives amplitude and phase dynamics in the local field potentials 47 , and phase coordination has been suggested to underlie routing mechanisms within and across brain regions [48][49][50] . Therefore, to confirm a functional correlation between thalamic spikes and mPFC non-REM oscillations, we assessed the phase locking of thalamic units to the delta and spindle frequency bands of the mPFC LFP.

Phase-locking of thalamic units to delta and spindle oscillations
We studied the phase-locking of thalamic units to the delta (1-4 Hz) and spindle (6)(7)(8)(9)(10)(11)(12)(13)(14) frequency bands by estimating the LFP phase using two complementary methods, the analytical signal and linear interpolation between consecutive troughs of the spindle and delta filtered LFPs (see Methods for details). Figure 7a shows the estimated phase obtained with both methods for the spindle trace in Figure 6a. Although the unit fires sparsely (raster above the trace), the spikes are more likely to occur on the descending phase of the spindle-filtered and ascending phase of delta-filtered LFP, shown by the phase histograms calculated for the whole sleep session (Figure 7b). We found similar preferential firing for the population of thalamic cells. The phase histograms were fitted for all units to a circular Gaussian (von Mises distribution), with mean and concentration parameters that indicate the direction and depth of the modulation or phase-locking. In order to prevent sampling bias due to non-uniform baseline phase distributions (e.g., due to nonsinusoidal LFP morphology), the unit phase distributions were also normalized to the baseline phase distributions during sleep 51 .
We found that the majority of thalamic units (67.74 % of n = 31 units) showed preferential firing at the descending phases of the spindle-filtered mPFC LFP, as suggested by the positive values of the mean phase distribution (circular mean 0.77 ± 1.15 rad; p < 0.05, Rayleigh test; Figure 7c). Thalamic units were phase-locked to the mPFC delta band, closer to the positive peak of the oscillation (64.52 %, p < 0.05, Rayleigh test; circular mean 0.30 ± 1.04 rad). The amplitude of the thalamic rebound after KCs (reported above) had a tendency to correlate with a preference for earlier phases of delta; although the correlation was not significant (p = 0.08, r = -0.32), we did not observe such a trend when we regressed mean phase preference to the cell firing rate in a time window of the same length but before KCs (p = 0.65, r = 0.08).
Interestingly, we found that the phase lags between the spikes of unit pairs that had the largest temporal cross-correlations (Supplementary Figure 7) were more consistent than the phase lags between uncorrelated pairs. In the sessions with thalamic pairs (n = 12), we calculated the smallest difference between the phase of the spikes during non-REM sleep, and between the phase of the spikes in spindle cycles. The phase difference would be expected to be distributed between -2π and 2π if units fire within one cycle of each other (positive or negative depending on their sequential order), and units that are functionally locked are expected to show a more consistent, less variable phase relation. Figure 7d shows the phase lag distribution for the pair described in Figure 6, calculated both from all non-REM spikes (left panel) and from the spikes that occurred within spindle cycles (right panel); the distributions suggest that most of the spikes occur within less than half a cycle of each other; the slight asymmetry with respect to zero suggests the consistent firing of one of the units before the other. We found similar results (Supplementary Figure 7) for the thalamic pairs that showed significant peaks in their cross-correlogram. The standard deviation of the distribution of phase lags was significantly lower for units with high spike cross-correlations during sleep (pair crosscorrelograms > 3 sd) compared to the units that were not temporally correlated. The standard deviation of phase lag distributions in correlated units was 2.7 compared to 3.05 in uncorrelated units when calculated with all non-REM sleep spikes (p = 0.026, Mann-Whitney U-test), and was 2.46 compared to 2.72 for the spikes within spindles (p = 0.009, Mann-Whitney U-test). These results demonstrate the phase-locked firing of thalamic cells within the period of the spindle oscillation.
Overall, the phase-locking results are consistent with the peri-event histogram analyses in suggesting coordinated thalamic firing with the neocortical delta and spindle oscillations, and demonstrate the phase coordination of pairs of thalamic units within individual spindle cycles.

Spindle activation of thalamic units modulated with CA1 SWRs
Because SWRs were also locked to preferential phases of the spindle oscillation (Figure 2b and Supplementary Figure 4), we investigated the correlation between thalamus and SWRs within the timescale of spindle cycles. We calculated the ripple power within each spindle cycle, and classified cycles in two groups based on the median SWR power (detected spindle cycles and feature distributions for all sessions in Supplementary Figure 8). Cycles with ripple power above the median were slightly shorter in duration (102.17 ± 2.01 ms compared to 103.22 ± 2.35 ms when ripple power was low; p < 0.05 in 65 % of the sessions although not significant for the mean duration across sessions, p = 0.11; Mann-Whitney U-test) and larger in amplitude (229.21 ± 102.14 µV versus 219.47 ± 94.32 µV; p < 0.05 in 85 % of sessions, p = 0.54 across sessions; Mann-Whitney U-test).  in the session on the left panel (Figure 8a), and for the two functionally defined types of spindle cycles, with low and high SWR power, in the middle and right panels. The spindle-trough-triggered histogram of CA1 SWRs for the same session is shown in blue, and confirms that the spindle classification identifies cycles with low or high probability of SWRs in CA1.
The majority of thalamic units showed an increase in spindle cycle modulation in cycles with ripple power above the median (67.7 %, n = 31). A modulation index, calculated as the difference between peak to trough firing relative to the mean in a 200 ms window centered on the peri-spindle-trough histogram, was significantly larger in cells with an increase in modulation compared to those with a decrease in modulation in high-SWR-power spindles. The same metric was not different between the two groups of units during low-SWR-power spindles, suggesting an enhanced recruitment of some thalamic cells in high-SWR-power spindle cycles (the index was 0.34 ± 0.18 in cells with increased modulation, compared to 0.22 ± 0.13 in cells with decreased modulation in cycles with high ripple power; p = 0.04, Mann-Whitney). Cells with increased modulation in high-SWR-power spindles also showed a significantly higher concentration in the distribution of spindle phases at the time of their spikes, indicating that they are more tightly phase-locked to high-SWRpower spindles (concentration parameter 0.23 ± 0.14 compared to 0.13 ± 0.08 in cells with a decrease in modulation in cycles with high ripple power; p = 0.03, Mann-Whitney), without changes in the average preferred phase compared to the cells with decreased modulation (0.35 ± 1.64 rad for units with increased modulation, 0.19 ± 2.04 rad for units with decreased modulation; p = 1), or compared to the average preferred phase of the same units in low-ripple-power cycles (0.12 ± 1.67 rad; p = 0.5, Wilcoxon signed rank). In some of the units, like the example in Figure 8, we observed a strong phase offset between SWRs and unit spikes ( Figure 8b); however, this was not the case in all the units, and some had spindle phase preferences that overlapped with those of SWRs recorded in the same session.
The modulation of thalamic spikes with SWR power was stronger for unit pairs; we calculated the modulation index for the cross-correlation between thalamic pairs recorded in the same session and found that it increased by an average of 26.74 % in spindle cycles with high ripple power compared to low ripple power (average index value 2.11 ± 1.70 in lowripple cycles, 2.51 ± 1.60 in high-ripple cycles; p = 9.77 × 10 -4 , Wilcoxon signed rank; n = 12 pairs), which points to an enhanced co-activation of thalamic cells within spindle cycles that engage hippocampal SWRs.
At broader time scales (seconds), we had found an overall decrease in thalamic firing and increase in bursting with SWRs ( Figure 3). However, it is possible that only spindle-nested SWRs are correlated with thalamic units. To test this possibility, we calculated the cross-correlograms of thalamic units with SWRs that occurred within the time window of detected spindle cycles, and compared them to a sample of equal number of SWRs that did not occur within spindle cycles. We compared the mean spikes/s in the 100 ms central window of these crosscorrelograms between thalamic single units and SWRs within or outside of spindles. We found low firing rates but no significant difference in thalamic firing near SWRs that occurred within detected spindle cycles (mean firing 3.83 ± 3.49 spikes/s with spindle-nested SWRs, compared to 3.99 ± 3.39; p = 0.18, Wilcoxon signed rank; n = 31); only when broadening the window of the cross-correlograms to at least 1 s around SWRs, did we find a slight decrease in thalamic firing with SWRs nested in spindles (mean firing 3.73 ± 3.43 spikes/s with spindle-nested SWRs, compared to 3.94 ± 3.35; p < 10 -2 , Wilcoxon signed rank; n = 31).

Multi-region correlations after spatial exploration
To explore the possibility that awake behavior influences the multi-region coordination and thalamic contribution, we analyzed 6 sessions in which the rats (n = 3; 2 sessions from each rat) explored a radial maze (a 4-arm maze in one rat, 8-arm maze in the other two). The animals explored the maze freely for 30 mins, collecting pellet rewards at the end of each arm. In these sessions, we recorded an average of 52 mins of sleep before and 78 mins of sleep after spatial exploration to compare the non-REM sleep events and their correlation with single units before and after spatial experience. All of the sessions showed an increase in the amplitude (peak-to-trough) and rate (per second) of SWRs during sleep after the maze exploration, as well as an increase in the amplitude of KCs (p = 0.03, Wilcoxon signed rank), while spindles did not change significantly in number or amplitude (p > 0.5).
Although the number and amplitude of detected spindle cycles did not change, there was a significant increase in CA1 SWR power within spindle cycles (p = 0.03, Wilcoxon signed rank), consistent with the overall enhancement of SWR activity in CA1. To find out if there is selective activation of behavior-relevant CA1 cells during spindles, we classified CA1 units recorded in these sessions (n = 42) into those that were active (n = 27) or not (n = 15) when the animal was on the maze (based on visual inspection during manual spike sorting, i.e., unit clusters present both on the maze and during sleep, or only in subsequent sleep; examples from one tetrode in Supplementary Figure 9a). These two groups of CA1 units did not differ in firing rate in the sleep post-behavior (calculated as the mean number of spikes over the total sleep time, or as the reciprocal

Figure 9. CA1 units active during 'run' show stronger activation in spindles with high ripple power, and firing at preferred spindle phases. (A) Peri-spindle trough histograms of the spikes of 2 CA1 units active during exploratory behavior (left) and 2 CA1 units recorded by the same tetrode that were only active during the sleep after maze exploration (right); top plots, spindles with high ripple power; bottom plots, spindles with low ripple power (sorted clusters for these units in Supplementary Figure 9a). (B) Distribution of mPFC spindle phases at the time of unit spikes in SWRs nested in spindles (top) or outside (bottom); p values on top of the plots from Rayleigh test. (C) Population distributions of the percent of spindle cycles in which the CA1 units fired spikes (left panels; Sp LR = spindles with low ripple power in CA1, Sp HR = spindles with high ripple power in CA1), and the mean firing for CA1 units active or not in run (right). (D) Left, distribution of peak firing rate in the peri-SWR histograms for CA1 units active or not during run, separated by SWRs nested or not within spindle cycles (SWRs Sp vs. SWRs not Sp respectively); right, distribution of Rayleigh's Z value suggests preferential spindle phase of CA1 units active in run when firing in SWRs nested in spindles.
of the average inter-spike interval; p > 0.1, Mann-Whitney U-test). However, we found evidence that units that were active during maze exploration are more functionally involved in SWR-spindle oscillations, as suggested by their higher firing and spiking at preferential spindle phases when they participate in SWRs nested in spindles; Figure 9a shows examples of the spindle temporal and phase organization of SWR spikes for units active or not in run. Units that were active during maze exploration were more likely to be active within spindle cycles with high ripple power than units not active on the maze (p = 0.009, Mann-Whitney U-test for the average number of spikes in spindle cycles and p = 0.039 for the percent of spindle cycles with unit spikes; Figure 9). Likewise, the peri-SWR histograms of CA1 cells that had been active during run showed significantly higher peaks (average 5.09 ± 4.23 spikes/s), compared to those units that were only active in the sleep post-behavior (average 2.34 ± 2.09 spikes/s; p < 10 -2 , Mann-Whitney U-test). Although units active in run had higher firing rates within SWRs, we did not find differences between their maximum firing in SWRs that were nested within spindle cycles compared to a sample of an equal number of SWRs not coupled to spindles (p > 0.6; Mann-Whitney U-test). However, the spikes fired by these units in SWRs nested in spindles occurred at preferential phases of the spindle oscillation (circular mean = -1.48 ± 1.33; concentration = 0.63 ± 0.32; average p value = 0.027 ± 0.05; n = 27 units; Rayleigh test; Figure 9). The phase modulation in the same spindle-nested SWRs, was not significant for units not active in run (circular mean = -1.20 ± 1.55; concentration = 0.51 ± 0.29; average p value = 0.218 ± 0.28; n = 14 units, one unit was removed from this analysis due to low firing in SWRs).
These results suggest an enhanced contribution of recently active CA1 units to SWRs, including those SWRs associated with thalamocortical spindles, when the spikes of units active in run are organized according to the phase of the spindle oscillation.
We then characterized post-exploration changes in thalamic activity, and did not find significant changes induced by maze exploration in this sample of thalamic units (n = 9). Figure 10a shows an example of a thalamic unit with a post-KC rebound, which had a similar profile before and after

Figure 10. Modulation of thalamic and CA1 activity before and after free exploration on a radial-arm maze. Peri-event histograms display fairly stable modulation of spiking activity in the thalamus with Kcomplexes before and after spatial exploration; KCtriggered histograms of thalamic spikes are shown in (A) and peri-spindle histograms in (B); instead, a CA1 unit recorded in the same session increases its modulation with spindles (C) in the sleep following behavior. (D) Spindle-triggered SWR histograms show an increase in SWR modulation following behavior. (E) Distributions of the percent of spindle cycles with CA1 unit spikes and the average number of CA1 spikes/spindle cycle before (black) and after exploration (yellow); each circle is one CA1 unit; most CA1 units increased their contribution to spindles post-behavior.
behavior. Thalamic units also had similar correlations with mPFC spindle cycles before and after behavior (Figure 10b), with no change in spindle phase distribution mean or concentration values (p > 0.3, Mann-Whitney U-test), although the modulation index estimated from peri-spindle-trough spike correlograms was on average higher after experience (0.21 ± 0.19 before exploration, 0.37 ± 0.31 after exploration; p = 0.02). This is in contrast to an increase in the spindle correlations with SWRs and with CA1 units that were active on the maze. Figure 10c-d, show examples of CA1 spindlecorrelations from the same session in which the thalamic unit in Figure 10b was recorded. Although the peri-spindle-trough histograms for the thalamic unit are similar before and after exploration (or even show a slight decrease in modulation), the CA1 unit (which was active during run) and the SWRs show enhanced modulation with spindles postexploration. On average, 60 % of the CA1 units fired more spikes in spindle cycles after exploration, and fired in a larger proportion of spindle cycles (distributions for all the CA1 units in Figure 10e).
We tested the possibility of selective modulation of thalamic units with spindle cycles that engage CA1 units that were active during maze exploration (compared to spindle cycles with spikes from CA1 units that were not active during behavior). We calculated the correlation of thalamic spikes with spindle cycles with spikes from the CA1 units that had been active on the maze, and compared it with their correlation with spindle cycles with spikes from CA1 units that were not active during exploration. We found no significant differences between the peak or the average firing values of the cross-correlograms between the two types of spindles and thalamic units (p > 0.3, Wilcoxon signed rank), nor any differences between the mean number of thalamic spikes per spindle cycle or the percent of spindle cycles with spikes from thalamic units (p > 0.1).
Another possibility is that thalamic cells specifically change their modulation with spindles when CA1 units are reactivating representations of spatial trajectories that were experienced while the animal was on the maze. During exploratory behavior, CA1 cells with neighboring place fields become active as part of theta sequences, showing strong phase-locked activation in the theta band when the animal is running (Supplementary Figure 9b); these units then tend to be co-active in resting periods that follow behavior, often during SWRs 11,36 . We selected pairs of CA1 units that had a clear comodulation during run (n = 19 CA1 pairs; examples in Supplementary Figure 9b), and calculated the correlation of thalamic units with spindles in which these pairs were active, compared to cycles in which they were not, and found no difference in thalamic activity (p > 0.9). Because the sample of thalamic units in these sessions was small (n = 9), we repeated these analyses with all the background spikes (MUA) recorded by thalamic tetrodes, finding no significant changes (p > 0.1).
In summary, these analyses suggest that spindles help organize the activation of recently active CA1 units during SWRs but failed to demonstrate a selective modulation of thalamic units with behaviorally relevant CA1 unit activation. It is possible that studying a larger sample of thalamic units, or after behavioral tasks with stronger memory demands, will reveal task specific modulations between these three brain regions. Nonetheless, these results indicate that there is a spindle related organization of the spiking of task-relevant CA1 cells in spindle-nested SWRs. This could promote CA1driven synaptic plasticity in mPFC networks during SWR reactivation during spindles.

Discussion
Here we report on the temporal correlations between cells in the limbic thalamus and the hallmark oscillations of non-REM sleep. We recorded single units and LFPs from functionally connected regions in the mPFC, CA1, and midline nuclei of the thalamus, including reuniens, mediodorsal and ventromedial nuclei. We found that thalamic cells in these nuclei often reduce their firing and are more bursty with CA1 SWRs during non-REM sleep. Although firing in the thalamus becomes sparse, spikes are correlated with other thalamic cells and with LFP population events in mPFC and CA1. Specifically, thalamic units were phase-locked to the delta and spindle frequency bands recorded in deep layers of mPFC, and some thalamic cells showed a strong rebound of activity following the down-state of the mPFC slow oscillation (1)(2)(3)(4). In addition, individual thalamic cells rarely contributed spikes in consecutive cycles of the spindle oscillation, but showed consistent correlations with other thalamic cells within the period of spindle cycles. We also found that CA1 SWRs were strongly correlated with individual cycles of mPFC spindles, a correlation that was enhanced following spatial exploration. Furthermore, CA1 units that were active during behavior were more strongly activated during SWRs compared to units not active during behavior, and including during SWRs nested in spindle cycles, when they phase-locked to the underlying spindle oscillation. The evidence raises the hypothesis that sparse synchronicity in thalamic cell ensembles facilitates CA1-mPFC consolidation processes during non-REM sleep. The enhancement of ripple-spindle correlations after exploratory behavior and preferential engagement of recently active CA1 units in SWRs and phase-locking to spindles, further suggest that individual spindle cycles provide time windows in which the interleaved activation of thalamic cells can promote incremental synaptic plasticity and facilitate the robust integration of recently acquired traces into neocortical networks.

Thalamic inhibition during hippocampal SWRs
The hypothesis of a multi-region interaction that stabilizes and integrates recently acquired episodic memory representations into neocortical schemas is supported by a stream of reports that point to the involvement of the thalamus 9,[15][16][17]21,27,52,53 . However, the investigation of the role of thalamocortical cells in these processes has lagged behind; evidence from unit recordings in naturally sleeping animals is still fairly limited, and we have been lacking simultaneous recordings of LFPs and single units from these brain regions to study the dynamics of their correlations in the same animal. Current results point to a negative correlation between thalamic activity and CA1 SWRs, which is suggested by the decrease in fMRI bold signal in the thalamus of monkeys 27 and the reduction of firing in MD and RE thalamic cells with SWRs 26,28 . Our results extend these reports with the observation of an increase in thalamic bursting, which suggests that inhibition of the thalamus (as opposed to a decrease in excitatory drive) explains the negative correlations between hippocampal and thalamic activity. One possibility is that the GABAergic TRN inhibits the thalamus near the time of SWRs. It is also possible that changes in neuromodulatory tone during non-REM sleep could lead to the hyperpolarization of cells in the dorsal thalamus in coordination with a change in CA1 state that facilitates SWR occurrence. The finding that cells in all the thalamic nuclei we recorded from were coordinated with CA1 and mPFC suggests a non-specific process that may be under neuromodulatory state control, without ruling out the involvement of the TRN. Because RE projects directly to CA1 and to the entorhinal cortex 54 , a change in firing mode in RE cells that project to hippocampus could further influence CA1 state transitions and SWRs occurrence.
Thalamic cells with an increase in firing with SWRs were rare in our sample (~15 %). A recent study reported larger ratios, finding about 40 % of MD multi-unit activity increase with SWRs that were coupled to spindles 28 . In this study, spindles were detected from EEG skull screws using criteria that included at least 0.5 s of duration; spindle-coupled SWRs were then classified as those SWRs that occurred between the onset and offset of spindles. Our LFP recordings in mPFC show that time windows longer than 0.5 s are likely to encompass one or more KCs, which are correlated with SWRs and often followed by a rebound in thalamic firing. In rodents, spindles defined by a duration criterion that overlaps with the delta oscillation are likely to reflect KC-SWRs correlations (in addition to thalamo-SWRs correlations), and could increase the fraction of thalamic electrodes positively modulated with SWRs. The direct recording of LFP and KCs together with the analysis of shorter, within-spindle, timescales helps disambiguate the multi-oscillation correlations.
Computationally, higher order thalamic nuclei (which include those in the midline 55,56 ) have been proposed to gate communication between neocortical areas. This 'gate' may be implemented through direct control of cortical synchronization and effective connectivity [57][58][59] , and through the regulation of the signal-to-noise ratio in resting states 28 . The broad decrease in thalamic activity at the time of potential memory reactivation in CA1, could facilitate effective information transfer between hippocampus and neocortex by reducing interference from the main input (thalamic) to neocortical areas. The signal-tonoise ratio may be particularly favorable to the transmission of novel traces between CA1-mPFC during spindles, due to the enhanced participation of recently active CA1 cells in SWRs. The spindle oscillation may provide a protected, phaseorganized 50 , functional window for select memory consolidation computations, for example those required to extract statistical regularities and promote synaptic plasticity. We find that there is an increase in thalamic bursting near SWRs, a firing mode that is more effective at activating post-synaptic cortical cells 42 . The function of thalamic bursting is not fully understood, but this thalamic response mode promotes better detectability of stimuli during wakefulness [60][61][62] , and has been proposed to serve as a detection signal that can facilitate plasticity and learning in thalamocortical networks 63,64 .

Functional windows for multi-region coordination in non-REM sleep
The up-states of the non-REM slow and delta oscillations are thought to provide discrete and relatively broad windows of hundreds of miliseconds in which hippocampal-neocortical interactions can occur 65,66 . Hippocampal SWRs happen mainly during up-states and nested within the cycles of spindles, and spindles are themselves largest in amplitude following the transition from the down-to-up state 17,67 . Thalamic cells often showed a rebound in firing after mPFC KCs (LFP markers of the slow oscillation trough), suggesting that SWRs that occur during the down-to-up state transition overlap in time with an increase in thalamic drive. The thalamic cells with post-KC rebound could influence CA1-mPFC interactions by influencing the dynamics within the slow oscillation up-state. Several reports have shown that subsets of neocortical and TRN cells are preferentially active at the start or at the end of neocortical slow oscillation up-states, although it is unclear what drives this differential activation 7,68-70 . Our data raise the possibility that cells in the dorsal thalamus (which can drive neocortical and TRN cells) could provide excitation early on in the slow oscillation up-state and influence the ensuing dynamics 8 . Cells in the thalamic nuclei that we recorded from project densely to superficial layers of the neocortex, where they can influence deeper layers through their synapses on the dendritic tufts of pyramidal cells or via layer I interneurons 31,71 .
At finer time scales, individual SWRs events have durations that are strikingly similar to the period of the spindle oscillation (close to 100 ms), within the temporal range of spike-timing dependent synaptic plasticity 72,73 . Indeed, SWR-or spindle-like stimulation patterns induce LTP respectively in CA3-CA1 73 and in neocortical synapses 74 . These cell ensemble and circuit level oscillations and the nesting of SWRs in the spindle cycles can therefore induce the long-term synaptic remodeling behind memory consolidation. Both synaptic enhancement and weakening have been suggested to occur in non-REM sleep and contribute to consolidation [75][76][77][78][79][80] . Although our preparation does not address the synaptic level, our results demonstrate that CA1 cells that were active during exploratory behavior have enhanced correlations with SWRs and are phase-locked to the mPFC spindle oscillation. An enhanced activation of task-specific cells during SWR-spindles, in temporal windows that can promote synaptic plasticity, may explain the increase in memory following spindlelike stimulation 21 , and suggests that spindle oscillations provide windows for synaptic reorganization of recently acquired memories. We observed an increase in spindle-CA1 correlations after behavioral experience even in the absence of significant changes of spindle-thalamic correlations.
It is possible that behavioral tasks with stronger memory demands would reveal a modulation of thalamic cells in addition to the observed CA1 modulation. Still, the differential enhancement of CA1-mPFC relative to thalamo-mPFC correlations could contribute to differential synaptic plasticity changes in these two pathways 77,79 . Spindles could provide a thalamocortical phase-organizing mechanism for the modification of CA1-mPFC synapses following behavioral tagging.

Sparse firing of thalamic cells during spindles and its role in systems memory consolidation
Sparse firing of thalamocortical cells during spindles has been reported before 39,39,[81][82][83][84] , but it was unclear how multiple cells in the dorsal thalamus would correlate with each other, and its relevance for the systems consolidation of memory has not been discussed. We find that pairs of thalamic cells are consistently phase-locked to the mPFC LFP within the period of spindle oscillations, suggesting that individual spindle cycles group and organize the sequential activation of thalamocortical ensembles. The repeated activation of thalamic cell ensembles through non-REM sleep could promote spike-timing dependent synaptic plasticity. Importantly, the sparse firing at the single cell level will allow for the interleaved excitation of different sets of thalamocortical synapses in subsequent spindle cycles, and in coordination with the reactivation of recently active CA1 place cells. The activation of varying combinations of thalamocortical synapses with CA1 SWRs that engage recently active cells (Figure 11) may promote a broad range of synaptic associations at the neocortical level, make the neocortical trace more robust to noise, and facilitate generalization. Critically, the sparse activation could prevent plasticity saturation 85 and contribute to gradual synaptic changes that would allow for extended computations 66 and facilitate continuous learning. One of the challenges in the development of artificial neural networks, has been to overcome the catastrophic interference that can occur during sequential training. While real-world settings require continuous learning, artificial network training algorithms are not yet optimized for these conditions. Sparsity and controlled weight changes have been suggested to facilitate robust representations and learning [86][87][88][89] . Future studies of thalamocortical sleep dynamics can provide clues on the biological implementation of continuous learning algorithms and help improve artificial networks and their applications.
Behaviorally, the systems consolidation of episodic memory leads to the stabilization of memories and to the formation of mental models or schemas that have higher adaptive value for the organism than the original details of the experience 2,3 . The oscillatory interactions between hippocampus and neocortex during sleep support the extraction of statistical consistencies and integration of hippocampal memory traces into gist-like neocortical schemas that are applicable beyond the specific conditions in which a memory was acquired 6,66 . Sparse thalamocortical population synchronicity during spindles may therefore promote incremental synaptic changes that optimize the interleaved extraction and incorporation of implicit information from recently acquired hippocampal memory traces into the neocortical representations.
It is worth noting that in certain states, such as anesthesia or epilepsy 90 , cells in the thalamus may fire at the population rhythm, and sparse synchronicity in the thalamus might be a distinguishing feature of normal resting states 91 . The disruption of thalamocortical oscillations has been proposed as a pathophysiological model for several neuropsychiatric disorders 92,93 . Oscillatory stimulation (mimicking the abnormally increased delta rhythm in schizophrenia) of thalamic projections to CA1 disrupts working memory in rats 94 , and spindle irregularities (reduced or abnormally enhanced spindles) have been observed in children with neurodevelopmental disorders 95 , in schizophrenia 96 , and in absence epilepsy 97 . A careful balance in the excitatory and inhibitory thalamic networks is required to generate the appropriate population level synchronicity and function. Studying spindle dynamics at the cellular and population level will thus be critical to understand the pathophysiology behind these disorders, and to design therapeutic interventions that preserve normal synchronicity in the thalamus as well as the multiregion coordination that underlies cognitive function.

Conclusion
Memory consolidation promotes the stabilization and incorporation of memory traces initially stored in the hippocampus into neocortical networks, contributing to the development and training of neocortical models that allow adaptive behavior. Recording from naturally sleeping rats, we present results that suggest that the individual cycles of the spindle oscillation provide a key processing time window that combines and organizes the activation of recently active CA1 cell ensembles, as well as subsets of cells in thalamocortical networks. The rebound activation of thalamic cells following down-states could facilitate up-state dynamics. As the cycle of the slow oscillation progresses, the enhancement of spindle-CA1 activity relative to spindle-thalamic activity following spatial experience could bias transmission in CA1-mPFC compared to thalamo-mPFC pathways. Lastly, the sparse and interleaved activation of groups of thalamic cells during spindles may provide a critical mechanism to facilitate the continuous integration of novel hippocampal memories into neocortical schemas for generalization and lifelong learning.

Materials and Methods
We used 8 male Long-Evans rats (400-500 g) in the experiments reported here. Animals were housed individually, provided with food and water ad libitum, and monitored in a temperature-controlled room with a 12 h light/dark cycle (lights on/off at 7:00 am/7:00 pm). These experiments were approved by the Committee on Animal Care at the Massachusetts Institute of Technology and conform to US National Institutes of Health guidelines for the care and use of laboratory animals.

Electrode Arrays, Surgical Implantation and Electrophysiology Recordings
Multielectrode arrays (microdrives) with up to 21 nichrome wire tetrodes, 18 tetrodes for recording and 3 references, one per brain region (individual tetrode wires were 12.5 µm in diameter). The microdrives were prepared according to standard procedures in the laboratory, and the base was modified from previous designs ( Figure  1a. The tetrodes spanned the following coordinates from Bregma (Paxinos and Watson atlas, 1986): Anterior 3.65 mm to 2.4 mm for mPFC, posterior 1.4 to 3.2 mm for the thalamus, and P3.1 to 4.4 mm for dorsal CA1. The arrays for mPFC and thalamus were positioned at 16 degrees from the vertical axis to prevent damage to the sagittal sinus. Sterile surgery procedures were performed for chronic microdrive implantation; anesthesia was induced by an intraperitoneal injection of a solution of ketamine (25 mg/kg), xylazine (3 mg/kg) and atropine (0.027 mg/kg), and maintained with 1-2 % inhaled isoflurane. Up to 8 bone screws were secured to the skull for support; three craniotomies were drilled over the target coordinates and the dura mater membrane was removed to allow for tetrode penetration. During surgery, the most anterior thalamic electrode was placed above P1.4 mm and lateral 1.8 mm, and the linear thalamic array was aligned parallel to the sagittal suture; this positioned the rest of the electrodes over their corresponding craniotomies. The implant was secured to the skull with dental cement after surrounding the exposed craniotomies and tetrodes with silicon grease to keep them from being fixed by the cement. Animals were preemptively injected with analgesics (Buprenorphine, 0.5-1 mg/kg, subcutaneous), and monitored by the experimenter and veterinarian staff for three days post-surgery.
The tetrodes were gradually advanced to the target depths over the course of 1.5-2 weeks before recordings started; the slow advancement of electrodes facilitated stable unit and local field potential (LFP) recordings over several hours. Once the target depth was reached, recordings were performed if units were present, and electrodes were lowered further in non-recording sessions to search for new units and allow for electrode drift before the next recording session. Recording sessions typically lasted at least 3 hours, during which the animal was in a quiet squared enclosure of about 50 × 50 cms, were they cycled through several bouts of sleep and wakefulness. Spikes were acquired at 32 kHz after the signal from any of the tetrode channels crossed a preset threshold. Individual units were isolated by manual clustering on peak spike waveform amplitudes across all channels using custom software (Xclust; M.A.W.). LFPs were sampled at 600 Hz after anti-aliasing filtering and downsampling from up to 3.125 kHz. Two infrared LEDs (sampled at 30 Hz) attached to the headstages of the acquisition system (Neuralynx) were used to track the rat's position and estimate velocity.
Spatial exploration in the radial maze. In preparation for the 6 behavioral sessions, 3 of the implanted rats were food restricted for at least 5 days before behavior started. In each session, rats explored a radial maze (with 4 arms for one rat and 8 arms for the other two rats) that was surrounded by consistent visual cues on the walls of the room; rats retrieved one 45 mg sucrose pellet every time the animal reached the end of the arm (provided by the experimenter). The behavioral exploration (two sessions per animal) lasted 30 min, and we also recorded at least 30 mins of sleep pre-and postbehavior. All the rats had less than 7 days of experience exploring the maze at the time of these recording sessions. For one of the rats, the two behavioral sessions were days 1 and 2 of exploration (no previous maze experience before day 1), for rat 2, they were days 2 and 4, and for the last rat they were days 3 and 6. Therefore, all the recordings were performed during early experience in the mazes, when memory traces are expected to be dependent on the hippocampus for retrieval (Varela, 2016). We did not observe differences between the first and second session in each rat, and the data were pooled together.

Selection of Electrodes and Sessions for Analyses
Sessions were selected based on having at least one unit in two of the three regions of interest, clear sleep events (KCs, spindles, SWRs) in the neocortical and hippocampal LFPs, and clear histological confirmation of thalamic tetrode location. While CA1 electrode location can be confirmed via electrophysiological LFP features alone, we find that LFP signals alone do not reliably allow specific localization within the nuclei of the dorsal thalamus and histological confirmation is necessary. For arrays with independently moving tetrodes, this requires the identification of electrode lesions for all the tetrodes in the array to prevent ambiguity when matching signals to electrode position. Thus, for these experiments we only include data from rats in which we successfully confirmed the histological location of all implanted thalamic (and mPFC) tetrodes, and CA1 location was confirmed electrophysiologically. The location of electrodes used for analyses is reported in Supplementary figures 1 and 2. To prevent repeated counts of the same units across sessions, unit analyses include only units for which 1) the tetrode had been moved from the previous recording session at least 80 µm or 2) the amplitude projections used for spike sorting had different profiles compared to the preceding recorded session.

LFP State and Event Detection
For each rat, we selected one wire from one tetrode with clear sleep events (determined by visual inspection) for LFP analyses (one tetrode per brain region). To detect reference events for analyses, the three LFPs were filtered (with zero-phase distortion) in the delta (1-4 Hz), spindle (6)(7)(8)(9)(10)(11)(12)(13)(14), and ripple (100-275 Hz) bands using a band-pass finite impulse response filter designed with a blackman window and filter order 6 times the sampling rate/bandwidth ratio.
We detected reference events during states in which the power of both the delta, in the mPFC LFP, and ripple bands, in CA1, were above the mean for the session for at least two minutes; animal velocity and behavior during these time periods confirmed that the animal was immobile during detected offline states (figure 1). These electrophysiological and behavioral traits are defining features of non-REM sleep. The position information to calculate velocity was obtained from the LEDs attached to the headstages; gaps in the position information (due to occlusion of the LEDs) shorter than 300ms were linearly interpolated using nearby values; gaps longer than 300ms were not considered for further analysis. The space derivative was used to calculate velocity and then smoothed (Gaussian kernel, 250 ms standard deviation) and linearly interpolated to up-sample and match the 600 Hz of the LFP.
Non-REM periods detected with these combined criteria of high delta and SWR power, showed large transient deflections in the mPFC LFP, or Kcomplexes (KCs). To detect KCs, we first rectified the filtered (1-4 Hz) LFP (positive values were made equal to 0), the resulting trace was squared and KCs detected as local maxima, when the amplitude was above the mean plus 4 standard deviations (mean and standard deviation calculated from the LFP in quiet states). We detected spindle troughs as local minima, 3 standard deviations below the mean of the spindleband filtered LFP (6-14 Hz). The SWR detection algorithm detected times when the squared, filtered LFP (100-275 Hz) had an amplitude above the mean plus 3 standard deviations for at least 20 ms (mean and standard deviation calculated for the LFP when the animal was quiet). If two SWRs were closer than 20ms they were considered a single ripple event. The SWR timestamp was selected as the time with the largest absolute value in the ripple filtered LFP.

Correlation quantifications: Peri-event histograms and generalized linear model
We used peri-event histograms (PETHs) to quantify the correlations between sleep events (trough of KCs, spindles, peak of ripples) and the unit activity in each brain region. To compare across sessions and animals, we quantified the change (in units of standard deviation) of the event-triggered average with respect to a null distribution calculated through a permutation test. Specifically, we first calculated the average peri-KC, peri-spindle and peri-SWR histogram across detected events in the recording session (bin 15 ms, followed by smoothing with a Gaussian kernel with 15 ms standard deviation). All the event timestamps were then shifted by a number between 1.5-2.5 s (randomly chosen with uniform probability) to offset the event times with respect to the local spikes, while preserving the inter-event time intervals. The average PETH (15 ms bins, Gaussian smoothed) was calculated from the shifted event timestamps, and this procedure was repeated 100 times to generate a distribution of expected PETHs under the null hypothesis of no correlation between sleep events and local activity (Supplementary Figure 3).
To characterize the relation between thalamic firing and SWR rate more precisely, we implemented generalized linear models (GLM) of CA1 spikes and SWRs as a function of thalamic spikes. Multi-unit spikes from CA1 and the thalamus, and CA1 SWRs, were binned through the non-REM sleep periods in each session (bin 250 ms). A GLM was fit in matlab (equation 1) to estimate the log of counts of spikes or SWRs in CA1 (µ) as a linear function of the spikes in the thalamus (X), assuming a Poisson distribution. From the fitted model, we calculated the estimated change in the response variable (firing in CA1) per unit change of the predictor (firing in thalamus) from the equivalent equation (2). The Wald test was used to test the significance of the thalamic parameter (β1). log(µ) = β0 + β1X (1) µ = exp(β0 + β1X) (2)

Phase-locking analyses
We used some of the phase analyses methods as in Siapas et al., 2005, and matlab´s circular statistics toolbox for quantification. We estimated the instantaneous phase using the Hilbert transform of the mPFC LFP filtered in the delta and spindle bands. Because asymmetry in the LFP shape can bias the phase distributions, the distribution of phases at the time of spikes or SWRs was normalized by the distribution of phases for the full delta or spindlefiltered LFP during non-REM sleep. After normalization, we used the Rayleigh test to determine if the angular phase at the spikes or LFP events timestamps were significantly different from a uniform distribution. A fit to the Von Misses distribution was then used to estimate the mean vector direction (mean preferred phase), and the degree of phase-locking from the concentration parameter kappa. The concentration parameter is a reciprocal measure of dispersion and increases as the density of values near the mean of the distribution increases. We also estimated phase based on linear interpolation between LFP peaks and troughs (arbitrarily assigned -π and π values; equivalent to the 'extrema' method in Siapas et al., 2005); we did not find differences in the phase-locking based on either method and the results are reported for phase estimates based on the Hilbert transform.