Differential Spike Timing and Phase Dynamics of Reticular Thalamic and Prefrontal Cortical Neuronal Populations during Sleep Spindles

The 8–15 Hz thalamocortical oscillations known as sleep spindles are a universal feature of mammalian non-REM sleep, during which they are presumed to shape activity-dependent plasticity in neocortical networks. The cortex is hypothesized to contribute to initiation and termination of spindles, but the mechanisms by which it implements these roles are unknown. We used dual-site local field potential and multiple single-unit recordings in the thalamic reticular nucleus (TRN) and medial prefrontal cortex (mPFC) of freely behaving rats at rest to investigate thalamocortical network dynamics during natural sleep spindles. During each spindle epoch, oscillatory activity in mPFC and TRN increased in frequency from onset to offset, accompanied by a consistent phase precession of TRN spike times relative to the cortical oscillation. In mPFC, the firing probability of putative pyramidal cells was highest at spindle initiation and termination times. We thus identified “early” and “late” cell subpopulations and found that they had distinct properties: early cells generally fired in synchrony with TRN spikes, whereas late cells fired in antiphase to TRN activity and also had higher firing rates than early cells. The accelerating and highly structured temporal pattern of thalamocortical network activity over the course of spindles therefore reflects the engagement of distinct subnetworks at specific times across spindle epochs. We propose that early cortical cells serve a synchronizing role in the initiation and propagation of spindle activity, whereas the subsequent recruitment of late cells actively antagonizes the thalamic spindle generator by providing asynchronous feedback.


Introduction
Sleep spindles are a highly conserved signature of mammalian non-REM (NREM) sleep, appearing in neocortical electroencephalographic (EEG) and local field potential (LFP) recordings as discrete 0.5-3 s oscillatory events at 8 -15 Hz. During NREM sleep, lowered neuromodulatory tone in the thalamus causes thalamic reticular nucleus (TRN) and thalamocortical (TC) cells to become hyperpolarized and fire low-threshold spike bursts due to their expression of T-type Ca 2ϩ channels (Crunelli et al., 1989;Huguenard and Prince, 1992). The reciprocal connectivity between GABAergic TRN cells and glutamatergic TC cells under these conditions forms a resonant circuit that is the basis for generating spindle oscillations (Steriade et al., 1993;von Krosigk et al., 1993).
The basic mechanisms of spindle generation between TRN and TC cells have been characterized in detail in cats, ferrets, and rats, under anesthesia and in slice preparations (Steriade, 2005). Evidence that spindle-like activity is seen in the decorticated thalamus but not in isolated cortex confirms that the essential pacemaker for the spindle oscillation is thalamic. However, the loss of long-range coordination between thalamic spindles in the decorticated ventrolateral thalamus in anesthetized cats suggests that the cortex is essential for coordinating long-range spindle synchrony (Contreras et al., 1996). Furthermore, another study also conducted in anesthetized cats revealed that electrical stimulation of motor or somatosensory cortex could evoke spindles, suggesting that synchronous corticothalamic feedback is able to trigger spindle initiation (Contreras and Steriade, 1996). However, a more recent hypothesis posits that the cortex may also contribute to spindle termination by depolarizing TC and TRN cells and thus desynchronizing rhythmic activity in the thalamus (Timofeev et al., 2001;Bonjean et al., 2011). Although the precise mechanistic bases of the contrasting cortical roles in spindle initiation and termination are unclear, single unit recordings from the medial prefrontal cortex (mPFC) in rats have shown that neuronal subtype and depth affect firing rate modulation and spindle-firing phase relative to the cortical LFP, both under anesthesia (Puig et al., 2008;Hartwich et al., 2009) and during natural sleep . This raises the possibility that different classes of cortical cells could fulfill distinct functional roles over the time course of individual spindles.
Given the likely importance of cortical influences on the temporal dynamics of spindles, we used multisite recording techniques to probe concurrent activity in the TRN and mPFC, a cortical region central to cognitive function. To date, such approaches have typically been restricted to small numbers of simultaneously recorded cells in primary motor or sensory regions and have been performed during anesthesia rather than natural sleep (Contreras et al., 1996;Bonjean et al., 2011;Ushimaru et al., 2012). Here, through the use of simultaneous tetrode recordings in mPFC and TRN in naturally sleeping rats, we reveal novel dynamic processes that occur during individual sleep spindles, thus gaining insight into the network mechanisms for the initiation and termination of the oscillation.

Materials and Methods
All procedures were performed in accordance with the UK Animals Scientific Procedures Act of 1986 and were approved by the University of Bristol Ethical Review Group.
Surgery and electrophysiology. Six adult male Sprague Dawley rats (300 -500 g, Charles River Laboratories) were anesthetized with either sodium pentobarbital or isoflurane and implanted with 16-tetrode custom-built microdrives targeting the medial prefrontal cortex (ϩ3.2 mm from bregma, 0.6 mm lateral, 1.5-3.0 mm ventral to brain surface; see Fig. 2A) and the rostral sector of the TRN (Ϫ1.5 mm from bregma, 1.8 mm lateral, 4.5-7.0 mm ventral to brain surface). The TRN rostral sector was targeted because it receives direct projections from mPFC and from thalamic nuclei that project to mPFC (Cornwall et al., 1990). Tetrodes were made by twisting bundles of 13 m polyimide-insulated nichrome wire (Kanthal).
Electrophysiological data were acquired using Digital Lynx hardware and Cheetah software (Neuralynx) while rats were at rest or asleep in a dimly lit, sound-attenuating chamber. Behavior was continuously monitored via four video cameras. Prefrontal tetrode signals were referenced to a local reference tetrode located superficially in the prefrontal cortex. Thalamic tetrode signals were referenced to a prefrontal or thalamic tetrode that registered no spike activity. The polarity of all signals was inverted, such that depolarization of neurons resulted in upward deflection of the recorded extracellular voltage. Threshold-triggered extracellular action potentials were sampled at 32 kHz and band-pass filtered at 600 -6000 Hz; simultaneous local field potential recordings were sampled at 2000 Hz and band-pass filtered at 0.2-475 Hz.
For targeting the thalamus, tetrodes were sometimes advanced by ϳ4 mm immediately after surgery. After a 24 h postsurgical recovery period, tetrodes were gradually advanced toward their target structures. TRN units were recognized by their characteristic "accelerando-decelerando" burst-firing activity during NREM sleep (see Fig. 2C; Domich et al., 1986) and narrow spike widths (see Fig. 2G). Occasionally, units without canonical burst-firing activity were classified as TRN units if they were recorded on the same tetrode as confirmed TRN units and also had narrow spike widths.
Perfusions, lesioning, and histology. At the end of each experiment, rats were deeply anesthetized with sodium pentobarbital (60 mg/kg) and electrolytic lesions were made by passing an anodal current (30 A for 10 s) through each tetrode. Rats were then perfused transcardially with 0.9% saline followed by 4% paraformaldehyde. Brains were removed and refrigerated for ϳ48 h in 4% paraformaldehyde and then for several days in 30% sucrose. 50 m sections were cut on a freezing microtome and Nissl-stained with cresyl violet or thionin.
Electrophysiological data analysis. All analyses were performed in MATLAB (MathWorks). Multitapered spectral analyses (Mitra and Pesaran, 1999) were performed using the Chronux toolbox (www. chronux.org), and circular statistics were computed using the CircStat toolbox (Berens and Velasco, University of Tübingen, Germany). Other analyses used custom-written code.
Spikes were sorted semiautomatically on the basis of waveform characteristics (energy and first principal component) using KlustaKwik (K.D. Harris, http://klustakwik.sourceforge.net/), and clusters were then checkedmanuallywithMClustinMATLAB(A.D.Redish,http://redishlab. neuroscience.umn.edu). Occasionally, when data files were very large due to high spike counts, clustering was performed entirely manually due to the infeasibility of automatic sorting. After clustering, only units with a mean peak amplitude Ͼ50 microvolts and with cluster isolation distances Ͼ15 (Schmitzer-Torbert et al., 2005) were selected for further analysis. Units with Ͼ1% of interspike intervals (ISIs) Ͻ2 ms were also excluded.
PFC units were classified in a multistep process: first, putative fastspiking interneurons (here referred to as narrow-spiking [NS] cells) were identified using an unsupervised hierarchical clustering algorithm (Ward's method) on the basis of peak-to-valley amplitude ratio, spike width, and mean firing rate (see Fig. 2 B, D,E), similarly to a previously described method (Quirk et al., 2009). The NS group was likely to contain fast rhythmic bursting or "chattering" cells (Gray and McCormick, 1996), which we eliminated by discarding any units with non-unimodal ISI distributions. To classify putative pyramidal (Py) cells, we applied a spike width threshold of 255 microseconds, which split the bimodal spike width distribution of the whole mPFC cell population approximately halfway between its two modal values (cf. Benchenane et al., 2010;Phillips et al., 2012). All cells with spike widths above this threshold (apart from a small number that had been classified as NS cells) were classed as Py cells.
mPFC principal pyramidal cells have been shown to possess different intrinsic burst-firing properties in vitro (Yang et al., 1996). Although our extracellular recordings did not permit us to measure intrinsic firing properties, we separated Py cells into subtypes according to their burstfiring characteristics (see Fig. 2F ) across the whole recording session. Cells with ISI distributions that contained at least two peaks and had Ͼ10% of ISIs Ͻ15 ms were classified as "bursting" cells. Cells with Ͻ1% of ISIs Ͻ15 ms were classified as "regular spiking" cells (Jones and Wilson, 2005). All remaining Py cells were classified as "intermediate." Spindle detection and other analyses. Spindles were detected using a previously described algorithm (Phillips et al., 2012). Briefly, mPFC local field potential signals were band-pass filtered at 6 -18 Hz (to avoid the filter roll-off compromising 8 -15 Hz spindle power) and the amplitude envelope of the rectified filtered signal was calculated. If the envelope exceeded a detection threshold of 3.5 SDs and remained above a start-toend threshold of 2.0 SDs for longer than 0.4 s, the event was classed as a spindle, with onset and offset times defined by the start-to-end threshold crossings. If two spindle events were separated by Ͻ0.5 s, they were amalgamated into one. Spindle events lasting Ͼ3 s were discarded.
Spindle phase was calculated by linear interpolation between troughs (phase 0) and peaks (phase ) of the 6 -18 Hz band-pass-filtered LFP. To eliminate any bias in the detection of spindle onset and offset toward particular spindle phases, the onset and offset times were pseudorandomly jittered by Ϯ0.5 phase cycles in spike-field phase analyses.
Because the precise timing of spindle activity may vary across mPFC subregions and between different layers of cortex, we compared spindles from all pairs of mPFC tetrodes in each rat (total n ϭ 29 pairs). Across tetrode pairs, we found that spindle troughs were temporally aligned, as indicated by a strong central peak in the mean cross-correlogram of spindle trough times (Fig. 1C). Furthermore, when we compared concurrent spindles on pairs of mPFC tetrodes, there was a strong correlation in spindle duration between tetrode pairs (mean r ϭ 0.65, Spearman rank). We therefore used LFP from one designated mPFC tetrode that showed robust spindle activity in each rat for subsequent analyses of spindle phase. The spindle phase relationships of unit populations were not biased by this approach (see Fig. 4 F, G).
Instantaneous frequency was calculated by measuring the time interval between the nearest pair of peaks or troughs to a particular perispindle time in the 6 -18 Hz band-pass-filtered LFP signal. If the evaluation time was closer to a peak than a trough, the nearest pair of troughs was used; if the time was nearer to a trough, the nearest pair of peaks was used.
To account for the variable duration of spindles, we used a measure of normalized time, here termed "relative spindle time" (RST; Fig. 1C). According to this measure, RST "0" was defined as the onset time of each detected spindle, and RST "1" as the termination time, with all other values of RST linearly interpolated around these 2 time points. RSTnormalized spectrograms were computed by resampling the timefrequency spectral matrices for individual spindle events such that they were aligned at RST 0 and RST 1, after which the aligned array of spectrograms was averaged.

Spectral characteristics of mPFC spindles
Spindles were observed in the LFP in mPFC as 0.5-3 s oscillatory events at 8 -15 Hz with a waxing-waning temporal pattern, typically accompanied by the prominent delta-band (1-4 Hz) activity of slow-wave sleep (Fig. 1A). We detected spindle epochs automatically and registered their times of onset and offset. To account for the variable duration of spindles, we hereafter use a measure of normalized time, RST, defined as 0 at spindle onset and 1 at spindle offset (Fig. 1G). Different electrodes in mPFC registered highly similar spindle activity: spindle events showed similar time courses and aligned phase between pairs of simultaneously recorded mPFC channels ( Fig. 1A-C).
We first investigated how the spectral properties of the mPFC LFP changed over the course of individual spindle events by randomly selecting 500 spindle LFP segments from each rat (n ϭ 3000 total) and generating an average power spectrogram (Fig.  1F ). As expected, spindle onset was marked by a sharp increase in 8 -15 Hz power, which lasted for the duration of the spindle. However, the predominant frequency of the oscillation consistently increased throughout the spindle, from a mean of 10.0 Hz at RST 0.1 to 11.5 Hz at RST 0.9 ( p ϭ 6.87 ϫ 10 Ϫ141 , Wilcoxon signed-rank test). As an alternative frequency measure, we calculated instantaneous spindle frequency from the time between individual peaks and troughs in the band-pass-filtered signal (Fig. 1G). This analysis confirmed a highly significant increase in frequency, from 10.7 Hz at RST 0.1 to 12.7 Hz at RST 0.9 ( p ϭ 1.57 ϫ 10 Ϫ128 , Wilcoxon signed-rank test).
To investigate whether frequency dynamics differed between spindles of short and long duration, we generated separate spectrograms in absolute time for the 50 shortest and 50 longest spindles from each rat (n ϭ 300 total; Fig. 1D). Both spectrograms showed a clear increase in spindle-band spectral frequency from Figure 1. Characteristic spindles in a mPFC local field potential recording. A, An example 30 s trace of mPFC LFP during NREM sleep showing juxtaposed wideband signals (black) from one TRN tetrode and two mPFC tetrodes. One 6 -18 Hz band-pass-filtered mPFC tetrode signal is also shown (red). Below is a power spectrogram generated from the same period showing transient increases in power within the spindle band (8 -15 Hz, indicated by horizontal dashed white lines on the spectrogram) during identified spindle episodes, which are marked below the filtered LFP trace with black horizontal bars. B, Scatter plot of the durations of concurrent spindle events recorded on the two mPFC tetrodes shown in A. C, Cross-correlograms of spindle trough times; each line represents the mean across all mPFC-mPFC tetrode pairs in one rat. Individual spindle trough time cross-correlograms were calculated with a bin size of 2 ms and were mean normalized before averaging across tetrode pairs. D, Mean power spectrograms generated from the 50 longest and 50 shortest spindles from each rat, plotted in absolute time and time locked to the midspindle time point. The same color scale is used in both spectrograms. E, Example of how RST is calculated over a mPFC spindle. The band-pass-filtered mPFC LFP is shown, with detected spindle onset and offset times marked with vertical dashed lines. Absolute recording time is marked on the upper axis and the relative spindle time is shown on the lower axis. Spike times from a mPFC putative pyramidal cell are plotted with corresponding RST values. F, RST-normalized mean power spectrogram (n ϭ 500 spindles from 6 rats, n ϭ 3000 total). Dashed white lines indicate spindle onset and offset times. G, Mean instantaneous frequencies (dashed lines) shown over RST for 500 spindles from each rat (n ϭ 3000 total) and overall mean (solid bold line).
spindle onset to offset. The mean instantaneous frequency did not differ between long and short spindles 0.1 s after spindle onset ( p ϭ 0.477, Wilcoxon rank sum test) or 0.1 s before spindle offset ( p ϭ 0.778, Wilcoxon rank sum test). Short spindles increased on average from 11.1 to 12.4 Hz over 0.52 s, whereas long spindles increased from 10.9 to 12.3 Hz over 2.20 s, indicating that the rate of frequency increase was larger in short spindles than in long spindles and that spindle termination frequency therefore remains remarkably consistent.
Local field potentials result from the summation of synaptic and spiking activity of nearby cell populations, so we next examined the single unit activity in mPFC and TRN to determine whether the spike trains of these cells showed properties that paralleled the frequency dynamics of the mPFC LFP.
Prior research has shown strong entrainment of spike times to cycles of spindle oscillations in TRN neurons (Contreras and Steriade, 1996) and in cortical interneurons (Puig et al., 2008;Hartwich et al., 2009;Peyrache et al., 2011), so we began by examining the spike trains of these two cell types. For this particular analysis, we selected cells that increased their intraspindlefiring rates by at least 25% (5 of 17 NS cells and 17 of 20 TRN cells) to ensure sufficient samples for reliable spectral analyses of spiking during spindles. We generated power spectrograms from the spike times of these cells (Fig. 3A, four leftmost panels), which showed a similar temporal profile to the spectrograms generated The mPFC lesion is located in the prelimbic (PrL) region; the TRN lesion is located in the rostral sector of the nucleus. B, Illustration of spike waveform parameters used for cell type classification (W, spike width; P, peak amplitude; V, valley amplitude). C, Example raster plot of TRN spike bursts during a spindle epoch. ISIs are plotted above the bursts to demonstrate their accelerando-decelerando temporal structure. D, Stage 1 of mPFC unit classification: scatter plot of spike width versus peak-to-valley amplitude ratios of all mPFC units. Automatically clustered NS cells are indicated; the remaining cell population is split into Py (spike width Ն 255 s) and "unclassified" groups. E, Comparison of peak firing rates (left) and mean firing rates (right) of NS cells compared with Py cells. F, Representative examples of ISI histograms for each cell type using all spikes from the recording session. In each plot, the solid gray line indicates the mean ISI distribution for each cell type and the vertical line indicates the burst-firing threshold. G, Mean spike waveform for each cell type normalized by peak amplitude. from mPFC LFP spindles, indicating that the rhythmic-firing frequency of these cells increased from early to late RST. As an alternative measure of rhythmic firing, we used spike time autocorrelograms from the same cells (Fig. 3B), taking spikes from 200 ms windows centered on time points RST 0.1 and 0.9. This revealed a shorter lag time of the first satellite peak at RST 0.9 compared with RST 0.1, indicating an increase in rhythmic-firing frequency over the course of the spindle. These temporal dynamics of rhythmic firing in single cells were recapitulated in coherograms and cross-correlograms of spike times from simultaneously recorded pairs of TRN and NS cells (Fig. 3A,B, two rightmost panels).
To quantify the prevalence of the above effects in the TRN and NS cell populations, for each cell, we searched the spectrograms for peaks in the spindle frequency band (8 -15 Hz) at each time point from RST 0.1 to RST 0.9 in steps of 0.1 (Fig. 3C). One NS cell and two TRN cells were discarded at this point due to not having a spindle-band peak at each time point. Among TRN cells, we found an overall significant increase in rhythmic-firing frequency from RST 0.1 to RST 0.9 ( p ϭ 0.00246, n ϭ 17, paired t test). NS cells also showed a significant frequency increase ( p ϭ 0.00171. n ϭ 5, paired t test). Despite small sample sizes, these results suggest that, in both of these cell populations, the cells that fired rhythmically during spindles showed frequency dynamics that were concordant with the mPFC LFP.
The low firing rates of Py cells made them unsuitable for spectral analyses of their spike trains; however, the large number of simultaneously recorded Py-Py cell pairs permitted measurement of the level of temporally correlated network activity by calculating the correlation coefficient of their spike rates, R SC (Cohen and Kohn, 2011). We calculated the value of R SC between all pairs of simultaneously recorded mPFC Py cells (Fig. 3D, n ϭ 3270 pairs) using a 200 ms time window centered on a RST value that was shifted incrementally from Ϫ0.5 to 1.5 to capture changes in correlation over the course of the spindle. R SC appeared highest before and after the spindle, showing a progressive decrease over the spindle's time course and reaching a minimum at RST 0.65. Because uncorrelated spike rates should have a median R SC value of zero, at each time point, we tested the significance of R SC using the sign test, with Bonferroni correction to adjust the critical p-value of 0.05 for multiple observations (p crit ϭ 0.00122). R SC was significantly different from zero during the early spindle period, but not during the remainder of the spindle. R SC was also significant before the onset, and briefly after offset, of the spindle. There was a highly significant effect of RST on R SC ( p ϭ 2.12 ϫ 10 Ϫ11 , repeated-measures one-way ANOVA). This result suggests that the progression of a spindle is accompanied by a gradual decorrelation of the firing rates of mPFC pyramidal cells. Due to the small numbers of simultaneously recorded pairs of NS (n ϭ 5 pairs) and TRN cells (n ϭ 10 pairs), we were not able to perform R SC analysis on these cell types.
The dynamic characteristics of rhythmic firing in TRN and NS cells and of the correlated activity between Py cells were suggestive of widespread changes in spike timing in the thalamocortical system during spindles. This finding led us to focus subsequently on two aspects of spike timing in relation to mPFC spindles: (1) spike phase in relation to individual cycles of the oscillation and (2) firing probability over the course of the whole spindle.

Spindle phase-firing properties of prefrontal and thalamic cells
We began by characterizing the basic phase-locking properties of the cell types we recorded: for all recorded cells, we selected spikes that occurred during spindles and calculated their phase relative to the mPFC LFP oscillation (Fig. 4A,D). We thus computed the preferred phase, phase concentration (), and phase-locking significance value (Rayleigh's test) for each cell and calculated the mean preferred phase and value for each cell type (Fig. 4E).
A high proportion of cells of each subtype fired with biased phase during spindles (Fig. 4B,C). Consistent with the thalamic contribution to spindles, the highest proportion of significantly phase-locked cells was found in the TRN (all 20 cells p Ͻ 0.00001); TRN cells also had the highest mean value ( ϭ 0.404, significantly higher than PyB, n ϭ 96, and PyI cells, n ϭ 339, p Ͻ 0.05 Tukey's test; Fig. 4E), indicating that the strongest phase bias was present in these cells. mPFC NS cells (n ϭ 17) were also robustly phase-locked (71% cells, p Ͻ 0.00001; 12% cells, p Ͻ 0.001; 18% cells, p Ͻ 0.05), with a mean value of 0.330, which was higher than the other mPFC subtypes. Histograms of the preferred phases of significantly phase-locked cells revealed a concentration of preferred phases around spindle troughs in all cell types (Fig. 4F ). Among Py cells, the strongest concentration of preferred phases around spindle troughs was found in bursting cells. To ensure that our calculations of spindle phase were not compromised by using spindle activity measured from a single representative mPFC tetrode, we recalculated the preferred spin- dle phases of all mPFC cells using the LFP signals from the same tetrodes that they were recorded on (Fig. 4G). The recalculated preferred phase distributions were similar to the original distributions, suggesting that there was little divergence in spindle activity between mPFC recording sites, and therefore that a single tetrode was sufficient to represent spindle activity for all mPFC cells.
Consistent with previously published work (Contreras and Steriade, 1996;Hartwich et al., 2009;Peyrache et al., 2011), we found that spikes from the majority of TRN and cortical cells were biased toward similar phases of the spindle oscillation. However, the differences between the preferred phase distributions of the Py cell subclasses suggest that burst-firing characteristics are predictive of spindle-phase-locking behavior: cells with strong burst firing were the most likely to show a phase bias toward spindle troughs, whereas a larger proportion of regular spiking and intermediate cells preferred to fire near spindle peaks.

Perispindle firing rate characteristics of prefrontal and thalamic cells
Next, we went on to probe how the firing probabilities of cells varied over the spindle's time course. We computed firing rate histograms for each cell over RST, examples of which are shown in Figure 5A. Firing rate histograms of pyramidal cells indicated a high heterogeneity of spindle-related firing among the popula-tion; many cells showed marked increases or decreases in firing rate at particular perispindle times. For each cell, we calculated the spike rate within each spindle and within a time window 1.2-0.2 s before the start of the spindle; in this way, we calculated how many cells significantly increased or decreased their firing rates during spindles using Wilcoxon's rank sum test (Fig. 5C). There was a trend for Py cells to increase, rather than decrease, their firing rates during spindles (34% of PyB cells increased firing rate whereas 10% decreased; 30% of PyI increased whereas 7% decreased; 16% of PyRS increased whereas 11% decreased). The activity of NS cells fell principally into two stereotyped patterns, the first (35% of cells) being a gradual but marked increase in firing rate toward the midpoint of the spindle before gradually decreasing again toward the spindle offset time (Fig. 5A, first two example cells). The second firing pattern (42% of cells) was a sharp decrease in firing rate at spindle onset, followed by a gradual increase reaching a maximum at the time of spindle termination (Fig. 5A, third example cell). In contrast, TRN cells almost universally showed increased firing rates during spindles (95% of cells), with peak firing times clustered near the middle of the spindle.
For each cell type, we created perispindle-firing rate histograms of all individual cells, which were then z-score-normalized and averaged to generate mean relative firing rate curves (Fig.  5B). Comparing firing rates during spindle early (0 -0.2 s after onset) and late (0.2-0 s before offset) periods revealed increased firing rates later during spindles in all cell types (data not shown); this was significant in all three mPFC Py cell subtypes (PyB, p ϭ 0.0014, n ϭ 96; PyRS, p ϭ 0.035, n ϭ 19; PyI, p ϭ 0.0065, n ϭ 339), but not in NS (p ϭ 0.554, n ϭ 17) or TRN (p ϭ 0.314, n ϭ 20) cells (Wilcoxon signed-rank test). Among Py cells, PyRS cells tended to increase their firing rates at spindle termination more than PyB cells (p ϭ 0.0536, n ϭ 19, Kruskall-Wallis test).

Early-and late-firing cells and relationship with spindle phase-locking
The specific and varied perispindle temporal firing patterns of mPFC cells paralleled the diversity in their phase-locking to spindle oscillations, so we next investigated the possibility that perispindle-firing rate patterns might predict phase-locking behavior. Because mPFC Py cell peak firing times during spindles were nonuniformly distributed ( p ϭ 5.10 ϫ 10 Ϫ11 , n ϭ 454, Kolmogorov-Smirnov test; Fig. 5D), showing a bias toward the spindle start and end times, we suspected that separate early-and late-firing cell subnetworks might exist. To quantify the tendency to fire early or late during spindles, we calculated each unit's firing rate within an early time period, RST 0 -0.4 (R E ), and a late time period, RST 0.8 -1.2 (R L ; Fig. 6A). We ranked all Py cells (n ϭ 454) according to their R E :R L ratio, classing the upper 10% of cells as early (E) cells, and the lowest 10% as late (L) cells (n ϭ 45 cells in each group). This resulted in two groups of cells with markedly distinct perispindle firing rate profiles (Fig.  6B). PyB and PyI subtypes were similarly distributed between the two groups (PyB, 16 early, 13 late; PyI, 28 early, 27 late; Fig.   6E), however, PyRS cells showed a tendency to fall into the L group (1 early, 5 late).
To investigate whether the E and L Py cell groups were drawn from distinct mPFC subpopulations with other dissociable properties, we began by examining the overall firing rates of these two groups. To selectively measure firing rates during spindle-rich time periods, we used all time periods that were within 15 s of at least one spindle epoch. Within these time periods, we found the firing rates of E cells to be substantially lower than L cells (E mean 0.59 Hz, L mean 1.82 Hz; p ϭ 8.50 ϫ 10 Ϫ5 , Wilcoxon rank sum test; data not shown).
Next, we calculated the spindle phase-locking statistics of the cells in the early and late classes using spikes extracted from within spindle epochs. A marked difference was evident between the distributions of preferred firing phases of the two classes (p ϭ 3.99 ϫ 10 Ϫ5 , circular equal medians test; Fig. 6C): among E cells, the mean preferred phase was 5.11 radians, almost an exact phase inversion from the mean phase 2.00 radians of L cells. We also found that L cells were less strongly phase-locked to spindles, indicated by a lower value than in E cells (mean 0.277 in E cells, 0.167 in L cells; p ϭ 0.0032, Wilcoxon rank sum test; data not shown).
Because E and L cells showed different spindle phase relationships with the mPFC LFP, it followed that the temporal relationship of their spiking activity with TRN cells may also differ. We therefore selected TRN cells that increased their firing rates by at least 50% during spindles and found all possible cell pairs that included one of these TRN cells and a simultaneously recorded Py cell (n ϭ 186 pairs). We ranked the cell pairs according to the R E :R L ratio of the Py cells and classed the upper 10% of pairs as E pairs and the lower 10% as L pairs (n ϭ 19 per group). For each cell pair, we calculated spike time cross-correlograms of Py spike times relative to TRN spike times between Ϫ100 and 100 ms using spikes occurring during spindles. The distributions of cross-correlogram peak times were markedly different between early and late cells ( p ϭ 4.66 ϫ 10 Ϫ4 , Kolmogorov-Smirnov test; Fig. 6D): the mean E cell peak lag time was Ϫ44.3 ms, whereas this value was 20.5 ms in L cells.
We repeated the early versus late spindle-firing analysis on cortical NS cells, but the relatively small sample size (n ϭ 17) precluded statistical power. Nevertheless, as in the Py cell population, perispindle-firing rate histograms of individual NS cells clearly revealed cells with early and late firing patterns (Fig. 5A). When E and L subclasses were selected by taking the upper and lower 25% of R E :R L -ranked cells (n ϭ 4 per group), there was no significant difference between the mean preferred phases of E versus L ( p ϭ 0.157, circular equal medians test; data not shown). However, across the NS population, there were strong correlations between log(R E :R L ) and preferred phase (r ϭ 0.8167, p ϭ 0.0035, circular-linear Spearman rank; data not shown), and between log(R E :R L ) and (r ϭ 0.759, p ϭ 0.0075, Spearman rank; data not shown). In conclusion, we found that perispindle peak firing time is predictive of spindle phase-locking characteristics in both putative pyramidal cells and interneurons.

Intraspindle phase dynamics
Finally, we investigated a related question of whether individual TRN and NS cells, many of which showed sustained phaselocked firing during spindles, changed their phase-locking properties over the spindle's time course. To visually inspect for changes in the distribution of spike phases during spindles, we produced two-dimensional histograms of phase versus RST (Fig. 7A1). Uniquely in plots generated from TRN cells, it appeared that spike phases tended to decrease from the beginning to the end of the spindle. We calculated the mean spike phases from all TRN cells at early (RST 0 -0.2) and late (RST 0.8 -1.0) time points (Fig. 7B), which revealed that 18 of 20 cells decreased their firing phase, with an average phase shift of Ϫ0.65 radians. Testing revealed this phase shift to be significantly different from zero in TRN cells ( p ϭ 4.02 ϫ 10 Ϫ4 , n ϭ 20, circular median test), but not in any other cell class. We also found a significant decrease in TRN cells' values from early to late RST ( p ϭ 0.0047, paired t test; Fig. 7C), which was not seen in any other cell type. These findings indicate that, over the time course of a spindle, TRN spike phases precessed relative to the mPFC local field potential, whereas their phase selectivity also decreased.

Discussion
We show that neuronal activity in TRN and mPFC of naturally sleeping rats follows highly structured paths from the initiation to the termination of spindles. First, we observed a consistent increase in the intrinsic frequency of mPFC spindles, which was accompanied by increasing rhythmic firing frequencies in TRN and NS cells and a precession of TRN spike phases relative to the mPFC LFP oscillation. Second, we identified subpopulations of putative mPFC pyramidal cells with distinct perispindle peak firing times and spike-phase relationships with spindle oscillations. These findings reveal the temporally specific and dynamic nature of cortical activity during spindles, suggesting a mechanism by which cortical activity may serve apparently opposite roles during spindle initiation and termination.

Spectral dynamics of mPFC spindles and effects on thalamocortical network coordination
The robust increase in frequency of mPFC spindles from initiation to termination independent of their duration may signify that frequency is an integral determinant of a spindle's time course and/or functional impact. This is counter to one report in which spindles in human frontal and parietal cortex decreased in frequency from onset to offset (Andrillon et al., 2011). Andrillon et al. (2011) proposed that this might indicate increasing thalamic hyperpolarization altering the latency of Ca 2ϩ spike bursts in TC cells over the course of spindles (McCormick and Bal, 1997). Species or recording (particularly the use of iEEG as opposed to LFP) differences may contribute to the discrepancy between our findings, although our results showing intraspindle acceleration are consistent with the hypothesis that the depolarizing I h current in TC cells is upregulated during the late phase of spindles (Bal and McCormick, 1996). We also found that the mPFC LFP oscillation was paralleled by increases in the rhythmic firing frequencies of TRN cells and some cortical NS cells, an important corollary to the LFP effects that supports the notion of TRN being central to spindle pacemaking.
One consequence of accelerating spindle frequency is a resultant shift of relative phase between thalamic and cortical activity. Assuming that (1) the pacemaker during spindle initiation is predominantly thalamic and (2) thalamocortical volleys excite postsynaptic cortical targets with a consistent latency, an increase in thalamic driving frequency (and consequent decrease in spindle cycle duration) would necessarily shift the phase of the cortical response relative to spindle cycles. The arrival of corticothalamic feedback in the thalamus is subject to further time lags from intracortical and corticothalamic transmission and would therefore experience a further phase shift relative to the TRN-TC bursting cycle. It is thought that the cortex may aid in synchronizing spindles during the early waxing phase (Contreras and Steriade, 1996;Contreras et al., 1996;Destexhe et al., 1998) while also desynchronizing thalamic activity during the waning phase, thus contributing to spindle termination (Timofeev et al., 2001;Bonjean et al., 2011). Phase shifting of corticothalamic feedback relative to the thalamic spindle cycle may allow the cortex to serve these distinct roles during early and late stages of individual spindles. Analyses of intraspindle frequency shifts may therefore inform the nature of thalamocortical dysfunction in schizophrenia and related disorders (Ferrarelli et al., 2012).

Spindle phase locking in the thalamocortical network
A large proportion of cells showed firing during spindles that was phase selective with respect to the mPFC LFP oscillation, typically firing during the spindle's descending phase near the trough. The highest proportion of spindle-phase-locked cells was found in the TRN, consistent with these cells having a crucial role in generating spindles (Steriade et al., 1987;von Krosigk et al., 1993).
In mPFC, Py cells tended to show phase biases toward both spindle troughs and peaks, in close agreement with a study using similar mPFC recordings in naturally sleeping rats . Interestingly, in that study, the authors found that the modal spindle-firing phase was opposite in superficially and deeply located cells; the diversity of preferred phases in our data may therefore be a consequence of recordings spanning cortical layers.

Structured activity of cortical subnetworks during spindles
Although spindles differentially modulate the firing of cortical neuronal subtypes (Puig et al., 2008;Hartwich et al., 2009;Peyrache et al., 2011), the extent and nature of distinct cortical networks activated at particular times during the oscillation remains unclear. Among mPFC pyramidal neurons, the clustering of peak firing rates near the times of spindle onset and offset led us to identify distinct populations of Py cells that preferentially fired during the early or late portion of the spindle.
We postulate that the E and L Py cells form subnetworks serving distinct roles in spindle initiation and termination. E cells tended to fire spikes that were strongly entrained to a similar spindle phase to spikes from TRN cells. The spike times of L Py cells were less strongly biased by spindle phase, but their overall mean phase was opposite to that of the E cells and TRN cells, indicating an inversion of the phase of mPFC Py spiking activity over the course of the spindle. Considering the complete thalamocortical system, these results suggest a shift from TRN-mPFC synchrony during the waxing period toward asynchrony during the waning period. L cells, by firing in antiphase with the thalamic spindle cycle, might antagonize the spindle-generating mechanisms and thus serve an active spindle-terminating function.
What mechanisms might underlie the markedly different firing patterns of the E and L cell populations? E cells, which fire with a similar phase to TRN cells, could be activated via shortlatency pathways after thalamocortical spike volleys; however, the much greater phase difference between thalamic activity and L cell activity implies that these cells might be excited by polysynaptic intracortical processes perhaps involving rebound spiking after inhibition by fast-spiking interneurons. One possibility is that the two populations are spatially segregated in the cortex, with E cells located superficially and L cells located in deeper layers. This would be consistent with the more abundant thalamocortical projections to superficial layers of the mPFC (Heidbreder and Groenewegen, 2003) and with data showing that deep-layer corticothalamic cells fire out of phase with TRN cells and are most active near the end of cortical up-states under chloral hydrate anesthesia in rats .
A previous study revealed distinct fast-spiking interneuronal classes with early-and late-firing characteristics during slow oscillation up-states in rats (Puig et al., 2008). Because spindles occur within the up-states of slow oscillations, it seems plausible that the early-and late-firing NS cells we found might correspond to the two cell types described previously by Puig et al., although their recordings were primarily made in superficial frontal/motor regions and under chloral hydrate anesthesia. Although some anesthetics clearly do affect spindle properties (e.g., ketaminexylazine disrupts the waxing-and-waning pattern of spindles; Contreras and Steriade, 1996), canonical spindle features appear robust to anesthesia and we found a similar relationship of phaselocking behavior with early-and late-firing tendency in NS cells. Fast-spiking interneurons are believed to be capable of synchronously entraining the firing of large pyramidal cell populations (Hasenstaub et al., 2005), so late-firing NS cells may aid in recruiting late-firing Py cells. However, classification of mPFC cells based on extracellular recording alone is necessarily putative and we cannot rule out the possibility of some Py cells being interneurons or vice versa (Vigneswaran et al., 2011).
One possible influence on the temporal activation patterns of early-and late-firing Py and NS cells is the locus ceruleus (LC), which has diffuse norepinephrinergic projections into the cortex that are believed to strongly modulate cell excitability (McCormick et al., 1991). The firing of LC neurons is phase-locked to cycles of slow oscillation (Eschenko et al., 2012) and these neurons also show elevated firing during spindles (Aston-Jones and Bloom, 1981). Increased norepinephrinergic tone during spindles may modulate the activity or plasticity of intracortical circuits that are necessary to generate the longer-latency responses of late-firing cells.
Spindles are believed to facilitate long-term plasticity processes (Contreras et al., 1997;Sejnowski and Destexhe, 2000;Rosanova and Ulrich, 2005). Their incidence increases after hippocampus-dependent learning in both rats and humans (Gais et al., 2002;Eschenko et al., 2006) and cooccur with hippocampal sharp wave-ripple (SWR) oscillations (Siapas and Wilson, 1998;Clemens et al., 2011). It is thought that the temporal coordination of spindles and SWRs may facilitate hippocampal-PFC information integration, which is believed to underlie offline consolidation of hippocampus-dependent memory (Diekelmann and Born, 2010), although there is some doubt about how strongly SWRs are able to influence mPFC activity during spindles . To clarify the functional roles of E and L cells during spindles, it would be of interest to examine how hippocampal replay activity during SWRs affects the activity of these cells.

Conclusions
We have described two novel phenomena occurring during natural sleep spindles, both of which may make essential contributions to shaping the stereotyped waxing-waning time course of the oscillation. In a field that has been predominated by studies of spindles induced by states of anesthesia, our evidence provides insight into the dynamic network activity that occurs during natural sleep spindles. Our results extend evidence for thalamic mechanisms of spindle generation toward a model in which corticothalamic feedback contributes to spindle termination by desynchronizing thalamic activity (Timofeev et al., 2001;Bonjean et al., 2011). Furthermore, whereas recent rat studies have shown differences in spindle-firing phase in mPFC according to cell subtype or depth (Hartwich et al., 2009;Peyrache et al., 2011), our findings suggest that these distinct phase-firing patterns may be accompanied by different activation time courses during the spindle. Future work linking the spindle-associated firing of cortical populations to behavior-dependent assembly activation across thalamocortical and limbic-cortical networks may reveal the functional consequences of these mechanisms.