Abstract

Spontaneous reactivation of previously stored patterns of neural activity occurs in hippocampus and neocortex during non-rapid eye movement (NREM) sleep. Notable features of the neocortical local field potential during NREM sleep are high-amplitude, low-frequency thalamocortical oscillations including K-complexes, low-voltage spindles, and high-voltage spindles. Using combined neuronal ensemble and local field potential recordings, we show that prefrontal stored-trace reactivation is correlated with the density of down-to-up state transitions of the population of simultaneously recorded cells, as well as K-complexes and low-voltage spindles in the local field potential. This result strengthens the connection between reactivation and learning, as these same NREM sleep features have been correlated with memory. Although memory trace reactivation is correlated with low-voltage spindles, it is not correlated with high-voltage spindles, indicating that despite their similar frequency characteristics, these two oscillations serve different functions.

Euston et al. (2007) used two statistical techniques to demonstrate reactivation in rat medial prefrontal cortex (mPFC). “Explained variance” (EV) is based on the similarity between the firing rate correlation matrices in the ensemble of simultaneously recorded neurons during behavior and subsequent sleep (Kudrimoti et al., 1999). The other technique, template matching, identifies multineuronal firing sequences during waking behavior that recur during postbehavioral sleep (Louie and Wilson, 2001; Tatsuno et al., 2006). EV is a relatively robust measure, largely independent of both the temporal order in which neuronal patterns appear and the speed of playback. Template matching captures the temporal order of cell-firing sequence replay, but is highly sensitive to both permutations of event order and fluctuations in playback timescale.

Reactivation occurs predominantly during non-rapid eye movement (NREM) sleep (Kudrimoti et al., 1999), which is characterized by high-amplitude, low-frequency oscillations of local field potentials (LFPs) (Steriade, 2003). Although NREM sleep in rodents is not generally divided into different stages [but see Jarosiewicz et al. (2002) and Jarosiewicz and Skaggs (2004)], distinctive electrographic features reflecting different underlying network dynamics can be identified. K-complexes, down-to-up state transitions (Steriade et al., 1993a,b), and the slow oscillation are different manifestations of the same neural population event (Amzica and Steriade, 1997, 1998) and are identifiable both in the LFP and in the sum of single spikes (Fig. 1a,b). Additionally, despite occupying similar and sometimes overlapping frequency bands, two different kinds of sleep spindles can be differentiated (Kandel and Buzsáki, 1997). The amplitude of one is three to five times larger than the other (Buzsáki, 1991), and thus we refer to them as high-voltage spindles (HVSs: 7–8 Hz) and low-voltage spindles (LVSs: 6–20, primarily 10–20 Hz).

We investigated how reactivation strength in the mPFC is related to neural population events observed in the LFP and the sum of individual neuron activity. The motivation for this effort was twofold: first, to determine when, during sleep, memory replay occurs and, second, to strengthen the link between memory-trace reactivation and memory consolidation. While reactivation has not yet been correlated with increased learning, several electrographic features of NREM sleep have been (Gais et al., 2002; Clemens et al., 2005; Eschenko et al., 2006; Marshall et al., 2006). A correlation between these same features and stored-trace reactivation strength would support the hypothesized connection between learning and sleep reactivation.

Materials and Methods

Subjects.

Three male Brown Norway/Fisher 344 hybrid rats (7–9 months old at the time of surgery, 350–400 g) were housed in Plexiglas home cages, maintained on a reversed, 12:12 h light–dark cycle. Training and recording sessions occurred during the dark portion of this cycle.

Data acquisition.

Neural recordings were obtained via a chronically implanted “hyperdrive” consisting of 12 independently manipulatable tetrodes. Each tetrode consisted of four polyimide-coated, nichrome wires (diameter 14 μm, Kanthal Palm Coast) twisted together (McNaughton et al., 1983; Recce and O'Keefe, 1989; Wilson and McNaughton, 1994). Hyperdrive construction was as described by Gothard et al. (1996) except that tetrodes were inserted into silica tubing (65 μm inner diameter, 125 μm outer diameter, Polymicro Technologies) for added rigidity, and secured with cyanoacrylate glue after insertion into the drive. Approximately 4–5 mm of wire protruded, ensuring that only the tetrodes penetrated the brain and not the silica sheath. During recording sessions, the hyperdrive was connected to a unity-gain headstage (HS-54, Neuralynx) that enabled low-noise transmission of neural data to the recording system. The headstage also contained an array of light-emitting diodes (LEDs) that could be detected by an overhead camera, enabling tracking of the position of the rat on the maze and during sleep at 60 frames/s (30 frames/s interleaved) with a pixel resolution of 640 × 480. All data were recorded using a Neuralynx Cheetah recording system running in combination with the computer. Single-unit data from each tetrode were amplified, filtered between 0.6 and 6 kHz, and digitized at 32 kHz. Video spatial resolution was ∼3 pixels/cm. Twelve LFP traces were collected, one from a single channel on each tetrode. The LFP traces were amplified and filtered between 1 and 475 Hz with a 1 or 2 kHz sampling rate.

Surgery and electrode placement.

National Institutes of Health guidelines and the University of Arizona Institutional Animal Care and Use Committee-approved protocols were followed for all surgical and behavioral procedures. Each rat was anesthetized with isoflurane (1–2% by volume in oxygen at a flow rate of 1.5–2.5 L/min), placed in a stereotaxic holder, and injected with penicillin G (30,000 U per hindleg, i.m.). The skull was cleared of skin and fascia, and craniotomies were opened for two stimulating electrodes targeting the medial forebrain bundle (MFB) [rat 1: −3.25 mm anteroposterior (AP), 1.65 mm mediolateral (ML) bilaterally, 8.5 mm ventral from dura; rat 2: −2.5 mm AP, 1.8 mm ML right, 8.5 mm ventral from dura and −4.0 mm AP, 1.5 mm ML right, 8.2 mm ventral from dura; rat 3: −2.5 mm AP, 1.8 mm ML right, 8.5 mm ventral from dura and −4.0 mm AP, 1.5 mm ML right, 8.2 mm ventral from dura] and a hyperdrive. Each stimulating electrode consisted of two Teflon-coated, stainless-steel wires (coated diameter 0.0045 inches, A-M Systems) twisted together with ∼0.5 mm of insulation removed from one tip. The tips were separated by ∼0.5 mm in the long axis. In rats 1 and 2, the hyperdrive was centered over the left mPFC at 3.0 (rat 1) or 2.9 (rat 2) mm AP, 1.3 mm ML, and angled at 9.5° toward the midline. In rat 3, the hyperdrive actually consisted of two groups of electrodes, only half of which went into the frontal cortex. The other six electrodes were positioned over the hippocampus, but the data from these electrodes are not reported here. The six prefrontal tetrodes were centered over the left mPFC at 3.0 mm AP, 1.3 mm ML, and angled at 9.5° toward the midline. It should be noted that the bundle of cannulae holding the electrodes, and therefore the electrodes themselves, extended ∼1.4–2 mm in the anterior–posterior dimension, so that cells were sampled from a relatively wide region of the medial frontal cortex. Medial–lateral extent varied between 0.5 and 1.4 mm, depending on the geometry of the cannula bundle. Rats were returned to ad libitum feeding and allowed to recover for 3–4 d after surgery. Single units were recorded with respect to a reference electrode positioned deep in the mPFC. The tetrodes were moved gradually from 2000 to ∼4000 μm from the brain surface, keeping all electrodes at the same depth as nearly as possible. For rats 1 and 2, a full set of behavior (see below) was obtained at a given depth, and then all of the electrodes were moved down 80–120 μm to ensure a fresh group of cells. After that, tetrodes were moved only as necessary to obtain good recordings until the next group turning event. This procedure was repeated until electrodes reached a depth of 4000 μm (rat 1) or 3200 μm (rat 2). In rat 3 the prefrontal electrodes were lowered to depths between 1800 and 2400 μm and held constant throughout the experiment except for occasional movements to obtain more cells. All tetrode positioning was done after a given recording session, to allow the tetrodes at least 18 h to stabilize before the next recording session.

After the conclusion of the experiments, a histological analysis was performed to confirm that all electrodes were located in the medial precentral, anterior cingulate, and prelimbic cortices [nomenclature from Krettek and Price (1977)]. A description of the analysis and pictures of the histological sections can be found in the supplemental material (supplemental Fig. S1, available at www.jneurosci.org).

Reward.

Following stimulation electrode placement, MFB stimulation was used as reinforcing reward [for a review and training techniques, see Olds and Milner (1954), Olds and Fobes (1981), and Liebman and Cooper (1989)]. All stimulation used two wires, though not necessarily the two twisted wires in one stimulating electrode. The choice of electrodes was determined empirically based on the rat's response. Despite the bilaterally implanted stimulation electrodes in rat 1, all stimulation was delivered exclusively on right side of the brain, contralateral to the hyperdrive. A range of stimulation parameters was explored using an operant conditioning chamber equipped to deliver MFB stimulation when the rat performed a nose poke. Stimulus parameters were selected based on the minimum net current needed to sustain repeated self-stimulation. The final selected MFB stimulation consisted of a train of 400-μs-wide, 70–100 μA, diphasic current pulses, delivered at 150 Hz for 320–370 ms.

Behavioral procedures.

The behavioral task has been described in detail previously (Euston and McNaughton, 2006; Euston et al., 2007) but is briefly described here. All rats were trained to find reward at one of the eight equally spaced reward zones on the edge of a 1.2 m circular arena. The correct zone was indicated by a flashing LED. Upon reaching the goal zone, reward was delivered, and after a brief delay (0.5–1.0 s), the next zone was cued by a flashing LED. Onset of the LED was also accompanied by a brief 4 kHz tone. This training continued until each rat made direct trajectories to reward locations. A training session lasted 50–60 min.

Rats were initially trained to run to randomly chosen reward zones, but were subsequently switched to a specific sequence of zones. After a rat completed a sequence three times with guidance from LED cues (a “cued” block of sequences), a 5 s delay was inserted between the nonspatial, audio cue and the illumination of the cue light, providing time for the rat to move to the next reward location without the aid of the visual cue. Given the typical running speed of a rat, the vast majority of cue-delay trials in well learned sequences were completed without the LED and are hence referred to as “noncued” trials. After the rat completed the noncued sequence three times, audio and visual cues were presented simultaneously, again, starting another cued block. Blocks of three complete traversals of the sequence alternated between cued and noncued throughout the duration of the recording session.

Rats were run in one of two different kinds of sequences: one eight elements long and the other six elements long. The eight element sequence contained two repeated segments in the shape of a “V,” followed by an alternating choice between two zones. This sequence was the same as that used by Euston and McNaughton (2006), who gave a thorough description. The other sequence had six elements with no repeated segments. For each sequence, rats were trained until they reached asymptotic performance. Asymptotic performance was defined in the same way as in Bower et al. (2005); the animal completed 70% of the cue-delay trials with sufficient continuity of movement and accuracy. Asymptotic performance was usually reached within 3 d, using two training sessions per day. For rats 1 and 2, electrodes were then pushed down to acquire new cells, and a new sequence was initiated. The new sequence was created by flipping (i.e., a mirror image) and rotating the original sequence so as to create a sequence novel to the rat. The flip ensured that the order of turns was reversed, while the rotation ensured that a different configuration of places was rewarded. Rats 1 and 2 were presented with a simple six-element nonrepeating sequence for 2 d before starting each new repeat-element sequence. Rat 3 was trained on the same eight-element “V” sequence for all days reported here. All data presented in the present paper were gathered on days when the rat had reached asymptotic performance. The rats were run daily for 1–4 months.

At the start of the experiment, the rat was brought into the recording room and placed in a towel-lined bowl in the center of the recording arena. The headstage was attached to the hyperdrive. The signal gain and threshold for each electrode were manually set on the basis of visual inspection, data recording was initiated, and the rat was allowed to rest for 20–40 min in the towel-lined bowl (pretask sleep, S1). The exact duration of the rest period was not systematically controlled. Following this rest period, the rat was moved from the bowl to the arena. The task phase (M1) began when the rat was placed on the behavioral platform. The rat was then returned to the towel-lined bowl and allowed to rest for another 20–40 min (posttask sleep, S2). Rats then ran another task phase (M2) and were given a final rest period (posttask sleep, S3).

Data preprocessing.

Artifact arising from MFB stimulation was removed from the single-unit data by deleting all spikes recorded during the stimulation and those occurring within 20 ms of stimulation termination (i.e., the black-out window was stimulus duration + 20 ms). The LFP was not analyzed during behavioral epochs. Individual units were isolated offline (for details, see supplemental material, available at www.jneurosci.org); the end result was of a collection of timestamps associated with each detected action potential from a given cell. The LFP traces were low-pass filtered at 100 Hz and downsampled to 200 Hz.

Determination of motionless periods.

Position during each video frame was extracted by fitting a circle to the ring of LEDs on the headstage in each video frame. The positions were then smoothed by convolution of both x and y position data with a normalized Gaussian function with standard deviation of 120 video frames (∼1.8 s). After smoothing, the instantaneous velocity was found by taking the difference in position between subsequent video frames. An epoch during which the headstage velocity dropped below 0.78 pixels/s for >2 min was judged to be a period of motionlessness. This threshold was set based on visual examination of the velocity profiles to exclude frame-by-frame jitter and subtle readjustments of body position. In practice, each 30–40 min sleep epoch contained several motionless periods with lengths varying between 2 and ∼30 min (mean = 18 min, standard deviation = 8.8 min).

Sleep analysis window.

For each “sleep” epoch the periods during which the animal was motionless for at least 2 min (sleep blocks) were extracted. For sleep blocks that were longer than 2 min, a sliding window was established. The 2 min window started at the beginning of the sleep block and moved forward in 2 min increments so that none of the analysis windows overlapped. At each of these steps, all of the analysis parameters were calculated for that time window. Thus, every increment of the window has a corresponding data point. Experimental sessions lacking any sleep blocks were excluded from the analysis. For each data point, the EV and the number of template matches were also calculated over the same time period.

Digital filtering.

Fourth-order Chebyshev type 1 filters implemented using the Matlab Signal Processing Toolbox (The MathWorks) were used for the digital filtering of both the local field potentials and the sum of the binned spike trains.

Detection of down states.

Down states are commonly recognized by the absence or near absence of spiking in the population of cells (Sirota and Buzsáki, 2005; Hoffman et al., 2007; Luczak et al., 2007). Consequently, down states were defined as times when all of the individual cells were silent for at least 100 ms. To find these times, all of the spike trains from populations of 50–150 simultaneously recorded cells within one session were binned into 50 ms (20 Hz) nonoverlapping bins, and all of the binned spike trains were summed. The incidences when this summed signal was zero for two or more consecutive points were counted as down states.

Detection of HVSs.

Although easily distinguished by eye, HVS epochs are difficult to detect automatically due to the presence of LVSs in the same frequency band (LVS power is mostly in the 10–20 Hz range but can be as low as 6 Hz). During a bout of HVSs, a large number of cells are progressively recruited and become phase locked to the oscillation. Because of this powerful oscillatory entrainment, HVS epochs are associated with high-amplitude deflections in the LFP. The depth of modulation in HVSs is greater than the depth of modulation in LVSs (Kandel and Buzsáki, 1997), and this was used to distinguish the two events. HVS epochs were identified by locating those periods of time when the dominant power in the LFP was in the spindle range and when the depth of signal modulation in the spindle range (as measured by the variance) was high.

Detection was accomplished by implementing the following procedure: (1) Two filtered versions of the LFP were made, one representing the K-complex frequency (2–6 Hz) and one representing the entire spindle range (6–20 Hz). (2) A 250 ms sliding window was used to evaluate the standard deviation of the spindle-filtered LFP over time. Starting at the beginning of the sleep block, the window was advanced 1/10 of its full extent (25 ms) until the end of the window reached the end of the LFP trace. In each increment of the 250 ms window, the standard deviation of the LFP filtered in the spindle frequency range was calculated over that window. The output of this step is referred to as the “standard deviation signal.” (3) To isolate the incidences of high variance, the points in the standard deviation signal that exceeded two times the standard deviation of that signal were selected. (4) To eliminate false identification of K-complexes as HVSs, the group of selected points was restricted to those falling during a time when the power in the spindle band was dominant. This was insured by averaging the windowed and filtered versions of the LFP over 500 ms and finding the times during which the magnitude of the power in the K-complex band was <80% of the magnitude of the power in the spindle band. (5) This process was performed independently for each of the LFP traces. For each of these, the number of peaks in a sleep block passing the criteria described above was summed. To quantify the duration of HVSs in each 2 min sleep block, an “HVS score” was computed by averaging the number of peaks identified in all of the LFP traces from all of the recording electrodes during that sleep block. The HVS score does not indicate the number of continuous bouts of HVSs, but rather the amount of time in each 2 min window that was occupied by HVS oscillations.

The bin width and power threshold values were subjective. However, these values were chosen to provide the optimal detection of sleep events and were therefore not truly arbitrary. Importantly, assignment of these values was made independently from the reactivation results and consistently applied across sessions and between subjects. A change in filtering parameters or thresholds would thus be expected to influence the number of HVS (or LVS or K-complex) epochs detected, but should not directly influence the relationship between sleep parameters and strength of reactivation, an assessment of which was the primary goal of this paper.

Detection of K-complexes.

K-complexes and LVSs both occur during times when the sum of single-cell activity fluctuates between up states and down states. As can be seen in Figure 1, K-complexes and down-to-up state transitions cooccur and are in fact different reflections of the same phenomenon (Amzica and Steriade, 1997, 2002). The detection of K-complexes in the LFP was performed using the same steps that were used for HVS detection with a few exceptions. First, in step 2 the absolute change in amplitude of the signal, rather than the standard deviation, was used to identify an event. This was because the K-complex is a high-amplitude, biphasic event, the full extent of which generally fits within one window. A window with a large peak-to-peak difference, therefore, is very likely to be centered on a K-complex. In each window, the absolute difference between the maximum amplitude and minimum amplitude was calculated, and thus we refer to the output signal as the “difference signal.” Second, the difference signal was smoothed with a 10-point moving average filter, and a peak-finding algorithm was used to identify the time points corresponding to local maxima. Those peaks that were greater than two times the standard deviation of the difference signal and that also occurred during times of high K-complex power were selected (the ratio of spindle power to K-complex power was <80%). The number of peaks falling during each sleep block was summed and, again, the number of peaks averaged over all of the LFPs from all of the recording electrodes was used as a measure of the number of K-complexes. Unlike the HVS score, the K-complex score represents the number of events, and is not a measure of time.

K-complexes and down states. A, B, The LFP (A) and the sum spikes from 135 single neurons (B) from a single experimental session are averaged on the down state to up state transitions. Before averaging, the LFP is normalized to the maximum value recorded on that electrode during that recording session. The red dot marks the transition point. Standard deviations are in green. In A, one-hundred fifty-five 2 s sections of one normalized LFP trace are centered on the down–up transitions identified in the sum of spikes. The characteristic shape of the K-complex is clearly revealed in the mean activity. The mean of the total activity (50 ms bins) in B is from the same 135 cells from which the down states were identified. C, D, During an epoch marked by up-and-down states (taken from a single session), the down states tend to recur at ∼1 Hz, as is indicated by the distribution of intervals (C), which peaks at ∼1 s; however, the distribution of intervals between down states is approximately exponential as shown on a logarithmic scale (D), indicating that down state occurrence is not a true oscillation, but essentially a random process.

Detection of LVSs.

The LVS detection method was the same as the HVS detection method, with two notable differences. First, in step 2, instead of the standard deviation, the power in the spindle band was used. Points exceeding two times the standard deviation of the binned spindle power were selected. Second, in step 4, the peaks were required to fall during times when the K-complex frequency was the dominant power in the signal. This effectively excluded all of the points that had already been designated as HVSs and required that identified incidences of LVSs cooccur with K-complexes. As before, the number of peaks occurring during each sleep block was averaged over all of the LFP traces from each of the tetrodes, and this was used to represent the amount of time occupied by LVSs.

Spindle separation.

To confirm that LVSs and HVSs are different populations of events, oscillatory activity in the spindle band of the LFP was further examined with respect to the dominant frequency and the depth of modulation, the two most distinguishing characteristics of LVSs and HVSs. For each rat, a session was selected in which a substantial amount of both LVSs and HVSs was identified. An LFP trace from one of the recording electrodes (randomly chosen) was binned (250 ms). The power and variance in the entire spindle band (6–20 Hz) as well as the power in the HVS (6–10 Hz) and LVS (10–20 Hz, capturing most of the LVS power) bands were calculated for each bin. The peaks of the spindle band, defined as the bins in which the power exceeded two times the standard deviation of the power in the 6–20 Hz band, were selected. The HVS power, the LVS power, and the variance corresponding to these peaks were used to characterize the activity in the spindle band.

Regression analysis.

The data points for each of the measures and for the EV for the two sleep sessions following behavioral epochs (S2 and S3) were grouped and expressed as a z-score. The data points for each measure and for the EV were then collapsed across recording sessions, and each measure was linearly fit to the EV using a least-squares method. Fits were considered significant for a p value <0.05.

Statistical comparison of sleep features before and after behavior.

To determine whether there was a significant difference in the incidence of down states, K-complexes, LVSs, and HVSs between the first (S1) and the last (S3) sleep sessions, a two-sample t test was performed to test the hypothesis that the means were the same. For each session, the score for each event (i.e., the number of down states) in S1 was divided by the number of 2 min windows in S1. Likewise, the score for each event in S3 was divided by the number of 2 min windows in S3. These normalized scores were pooled over all sessions. The equal mean hypothesis was tested at the 5% level.

Explained variance.

For the first sleep epoch, the behavioral epochs, and each of the sleep blocks pulled from the second and third sleep epochs, spike trains were binned into T nonoverlapping intervals of 50 ms (for comments on bin size, see supplemental material, available at www.jneurosci.org), producing sequences of spike counts xi(t). The normalized correlation C between each pair (ij) of spike trains was computed using the following equation:
derived from Pearson's correlation coefficient, where μ is the mean and σ is the standard deviation of x(t). Thus the same population of cells produced matrices of firing rate correlations, corresponding to the pretask sleep, the task, and each of the posttask sleep blocks. The correlation matrix is symmetric, so the lower off-diagonal elements are extracted, producing a vector of all pairwise cell correlations. The separate correlation matrices are entered into the subsequent EV calculations, described next.

Any measure of reactivation must show that activity patterns induced during the task reappear during the subsequent sleep. In addition, activity patterns that existed during the first sleep epoch should be discounted because these were evidently not induced by the task experience. The EV measure is the square of the partial correlation of task and posttask sleep cell-pair correlations, partialing out any preexisting correlations during pretask sleep:
where T, Sa, and Sb represent task, pretask sleep, and posttask sleep periods (Kudrimoti et al., 1999). EV reflects the proportion of variability in the Sb sleep correlations that can be accounted for by task correlations after accounting for Sa sleep correlations. Any pretask sleep replay of task activity (e.g., from a previous session) or, indeed, any restriction on the allowable correlations due to connectivity constraints would tend to reduce the EV measure.

Template matching.

Each template was an N × M matrix representing the firing rate of N cells (rows) over M time bins (columns). Six or eight templates were generated, one for each segment of the sequence. A segment was defined as the time between arrival at one reward zone and arrival at the next reward zone, excluding times during which stimulation was delivered. Each row of the template comprised the spike counts from one cell within a series of 100 ms bins covering one segment, averaged over all repetitions of the sequence. Repetitions were first screened to eliminate segments that took inordinately long times (i.e., segments during which the rat was off-task). For each segment, any repetitions during which the traversal time exceeded four times the distance of the quartile from the median were excluded. For the remaining repetitions, each segment was scaled so that traversal time equaled the median time. The spikes from these scaled repetitions were then averaged in 100 ms bins. Note that median times were different for the different segments, meaning that the number of template bins varied between 13 and 42.

The bin width used in template matching was chosen for both behavioral and analytical reasons. Based on observation, a bin width of 100 ms effectively captures the rate of change of cellular firing during behavior. With smaller bins, a given cell's firing rate may not change significantly from one bin to the next. With longer bins, temporal sequences across cells are no longer discernable. Bin width is also constrained by the need to account for temporal compression of replay. For a 100 ms template bin width at six times compression (as used in this paper), the target bin was 16.6 ms. In practice, 16 ms is near the lower limit of bin sizes, which yields strong template correlations for this data. Finally, the 100 ms bin width has proven effective for detecting reactivation in a previously published paper (Euston et al., 2007).

Template matching was based on a subset of cells that were active during all phases of the experiment and that exhibited task-related firing rate changes. Any cells with <50 spikes in any given period (pretask sleep, task, or posttask sleep) were excluded. Two additional criteria were applied to determine whether a cell should be included in the template-matching analysis. The first was a firing rate criterion. Cells in which the average spike count exceeded five spikes during a given template (i.e., the sum of all spike counts within that cell's row of the template exceeded five) were judged to have a sufficient spike count for that template. The second criterion was that the cell show strong task-related activity. Toward this end, each row of each template was first smoothed using a five-point (i.e., five 100 ms bins) moving average. This removed high-frequency random fluctuation of firing activity but left task-related slow modulation. The coefficient of variation (standard deviation divided by the mean, computed on the 13–42 values in a given row) was then computed for each cell that exceeded the firing rate criterion. A cell was judged to have significant task-related firing in a given template if its coefficient of variation exceeded 0.40. To be included in the template-matching analysis, neurons needed to pass both the firing rate and coefficient of variation criteria on at least half of the templates (three in the case of six templates and four in the case of eight templates). This requirement resulted in a significant reduction in the number of cells; typically between one-third and two-thirds of the original number of cells were included in the template-matching analysis. In addition, it was also confirmed that the template-matching analysis was robust to variations of the threshold.

The spike activity of all included cells from a recording session is stored in an N × T spike matrix Q (McNaughton, 1998), where rows correspond to N cells and columns to T discrete bins. The bin contents represent the number of spikes each cell fired during each time bin. In the limit of infinitely small bins, each row becomes a spike raster. As described above, the template matrix is an averaged firing pattern during each segment of the sequential task. A single template, X, is an N × M matrix, where M corresponds to the number of bins in the template. A target matrix Y, with the same dimensions as the template, is also selected from the Q matrix. In matrix form, both template matrix X and target matrix Y are represented as follows:
The template-matching method seeks to calculate the similarity of these two matrices (Louie and Wilson, 2001). Based on our previous study on different measures of template match strength (Tatsuno et al., 2006), we used the standardized Pearson correlation coefficient measure. In this measure, each row of the X and Y matrices is standardized to zero-mean and unit variance by subtracting its row mean, xn¯ and yn¯, and dividing by the row standard deviation, σx,n and σy,n, respectively. xn¯, yn¯, σx,n, and σy,n are defined as follows:
This normalization transforms the elements xnm and ynm to z-score variables wnm and znm through the following:
By construction, mean firing rate differences among different rows are fully suppressed with this normalization. Therefore, this measure is sensitive to the firing order relationships among different neurons. The standardized Pearson measure is then defined using the two-dimensional Pearson correlation coefficient, COR, for the normalized W and Z matrices:
where the mean w̄ and z̄ are calculated as follows:
In a previous study (Euston et al., 2007), we determined that sleep replay occurred approximately six times faster than the speed of firing during behavior. Therefore, in this study, we modified the bin width of the target matrix to reflect a compression factor of 6. Examples of the template-matching results can be found in the supplemental material (supplemental Fig. S4, available at www.jneurosci.org).

To assess a template “match,” the z-score for each match was computed using a sampling distribution estimated from shuffled templates. For this analysis, the column vectors of the template matrix were randomly shuffled, and 100 shuffled templates were generated for each template matrix. The template match was then calculated for each shuffled template, and a sampling distribution was generated for each time point. The raw template match value at each time point was then converted to a z-score, and the z-score >4 was marked as a template “match.” The number of template matches found for each of the templates during each 2 min window was summed. This number, the total number of template matches, was used as an alternative to explained variance to quantify the magnitude of reactivation.

Results

For these experiments, three rats were chronically implanted with 6 or 12 moveable tetrodes in the prefrontal cortex. These tetrodes were used to record simultaneously the single-unit firing and the local field potentials generated by the brain during active performance of a sequence running task. Based on the histologically reconstructed electrode tracks and the depths to which electrodes were lowered, all electrodes were confirmed to be located in the most ventral portion of medial precentral, anterior cingulate, or prelimbic cortex.

Each recording session consisted of an initial 20–40 min sleep epoch (S1), a 50 min behavioral epoch in which the animal ran a trained sequence task (M1), an intermediate 20–40 min sleep epoch (S2), another 50 min behavioral epoch (M2), and a final 20–40 min sleep epoch (S3). The experimental protocol was slightly modified for the third rat; the final sleep epoch was extended to ∼120 min. Here, however, we are only using the first 40 min of this extended sleep period to make the results between animals more comparable.

The analysis was restricted to times when the rat was either sleeping or resting motionlessly. Video-tracker data were used to identify the times in the sleep epochs during which the animal was not moving. Previous analyses showed that the frequency of occurrence of REM sleep during the first 40 min was extremely rare, and so the analyses focused on aspects of NREM sleep.

Characterization of sleep features in the LFP and single-unit activity

In a plot of the LFP during a sleep session, it is clear that motion (wakefulness) is marked by low-amplitude, high-frequency activity. During motionlessness, on the other hand, the LFP is marked by high-amplitude, low-frequency activity. Likewise, motionlessness in the sum of individual spikes is indicated by times when all cells stop firing and the sum of spikes goes to zero (down states) (Fig. 2a), a phenomenon that is not observed during waking. Closer visual examination of the LFP during motionlessness reveals that this state is marked by high-frequency, low-amplitude oscillations punctuated by high-amplitude, low-frequency oscillations (Fig. 2b). Again, these states are not only visible in the LFP, they are also apparent in the sum of single-unit activity.

Examples of LFP and the sum of single-cell activity during periods of motion and motionlessness. A, During the full sleep epoch, the animal alternated between times of waking motion and quiet motionlessness. Motionless periods are indicated by a gray background. In the top panel is one LFP trace, in the middle panel is the absolute value of the total velocity, and the bottom panel is the sum of single-cell activity (124 cells, 30 ms bins). Wakefulness is marked by high-frequency, low-amplitude oscillations in the LFP and consistent spiking in the single units (no bins contain 0 spikes). Motionlessness, on the other hand, is marked by periods of high-amplitude, low-frequency oscillations in the LFP and multiple down states (number of spikes goes to 0 in several bins) in the single-unit activity. B, When a motionless period is examined on a finer timescale, it is clear that the system is sometimes in a state similar to that seen during motion. However, that state is punctuated by two additional low-frequency, high-amplitude “states.” These are HVSs (left arrow) and K-complexes/LVSs (right arrow). Arrows mark the approximate middle of these events. Bouts of HVSs and K-complexes/LVSs are visible not only in the LFP (top), but also in the sum of single-cell activity (bottom).

The bouts of high-amplitude, low-frequency oscillations can be further differentiated into periods of high-voltage spindling and periods marked by K-complexes and low-voltage spindles (Figs. 2b, 3). This differentiation is based on frequency characteristics in the LFP, entrainment of spikes in the single-unit activity, and the presence or absence of down/up states. In Figure 3 the activity of a group of cells recorded simultaneously is compared with the LFP from the same electrodes during HVSs (Fig. 3a) and K-complexes/LVSs (Fig. 3b). It is easily demonstrated that HVSs correspond not only to a distinctive 7–8 Hz oscillation in the LFP, but also to a 7–8 Hz oscillation in the sum of spikes. Although the activity frequently reaches zero during the oscillation, down states, as defined above, do not occur during HVSs.

Example of HVS epoch and down state/K-complex/LVS epoch. A, HVS epochs occur during periods of low-frequency, high-amplitude activity in the LFP (top). These oscillations are also identifiable in the total spike activity, which oscillates at the same frequency as the LFP during HVSs (124 cells, 30 ms bins, middle). The red bars indicate the onset of the HVS epoch. In the spectrogram of the LFP (bottom, hot colors = high power), HVSs are marked by a strong peak at 7–8 Hz and smaller peaks at several harmonics. B, K-complex/LVS epochs also occur during periods of low-frequency, high-amplitude activity in the LFP (top). These periods correspond to periods of up state/down state fluctuations in the total spike activity (124 cells, 30 ms bins, middle). A clear example LVS preceded by a K-complex is marked by the arrow in the top panel. Each K-complex in the LFP is matched by a network down state in the total spike activity (red bars). In the spectrogram of the LFP (bottom), the K-complexes correspond to high power in the 2–6 Hz range, while LVSs correspond to a peak in the 10–20 Hz range. Although the LVS example shown here shows power in a frequency range distinct from HVSs, in fact many LVS episodes show power in frequencies ranging from 6 to 20 Hz, thus largely overlapping with HVSs. Spectrograms were computed over 1 s Hamming windows with 95% overlap between the windows and evaluated at 210 frequency points. Shown here are normalized versions; hot colors = high power.

The LFP frequencies associated with LVSs tend to be higher and more distributed (10–20 Hz) (Fig. 3b), although we have observed them to be as slow as 6 Hz. Unlike HVSs, incidences of LVSs are frequently coupled with K-complexes in the LFP and down/up states in the sum of spikes. Another important difference between HVSs and LVSs is that LVS oscillations do not recruit the entire population of cells such that they can be identified in the single-unit activity (Fig. 3b). The fact that HVS oscillations entrain a much larger fraction of the cells relative to LVSs presumably accounts for the much larger LFP peak-to-peak amplitude seen in HVSs. This demonstrates the extent to which the population of cells becomes entrained to the HVS oscillation.

The different qualities of HVSs and LVSs are further illustrated in Figure 4. Peaks in the spindle range (6–20 Hz) have a bimodal distribution in both frequency and depth of modulation. This clustering shows that HVSs and LVSs represent two distinct populations of spindles, the functional implications of which should be evaluated separately.

Spindle distributions. For one session from each rat (rows), the LFP trace was binned (250 ms) and the power and variance of each bin were calculated. Bins where the power exceeded two times the standard deviation of the power in the 6–20 Hz band were selected and further analyzed. In the first column, the 10–20 Hz power of the selected bins is plotted against the 7–8 Hz power of the same bins. The two clusters correspond to HVSs (green) and LVSs (red). In the second column, the depth of modulation (variance) in the LFP for the same bins is plotted against the depth of modulation in the sum of individual spikes over the same time period. Clusters form in both frequency and in depth of modulation, indicating that the two different types of spindles are drawn from different populations. These results are typical for sleep sessions with both HVSs and LVSs.

K-complexes are strongly associated with down state to up state transitions (Figs. 1, 3b) and thus, like HVSs, are identifiable in the single-unit activity. During those periods when K-complexes interspersed with low-voltage spindles are observed in the LFP, the sum of spikes fluctuates slowly between complete silence and strong activity (down/up states) (Fig. 5). In fact, the distribution of the total spike activity transitions from unimodal to bimodal at the onset of K-complex activity in the LFP.

Distribution of population activity states during waking and motionlessness. A, B, Dividing the LFP (A) into periods without K-complexes (left of the red bar) and periods with K-complexes (right of the red bar) is equivalent to dividing the total spike activity (135 cells, 30 ms bins, B) into periods without down states (left of the red bar) and periods of fluctuation between down states and up states (right of the red bar). C, This difference is reflected in the distribution of the total spike activity; the distribution (taken from the data shown in B) is unimodal when there are no down states (left) and bimodal during periods of down state/up state fluctuations (right).

Relationship between sleep characteristics and memory reactivation

For analysis, the motionless periods were subdivided into 2 min nonoverlapping windows, each of which defines a data point. “Sleep” epochs without any motionless periods were not analyzed. For each window, four sleep-characterization parameters were calculated: the number of down states, the number of K-complexes, the amount of LVSs, and the amount of HVSs. The EV (taking into account S1) over each window was also determined so that each four-dimensional data point was matched to an EV value. A variable number of data points were computed for each session depending on the amount of time the animal remained motionless during that session. For the first rat, 29 sessions and 356 data points were analyzed, for the second rat, 3 sessions and 72 data points were analyzed, and for the third rat, 10 sessions and 196 data points were analyzed. On average, 98 cells (range: 66–142), 74 cells (range: 69–78), and 63 (range: 43–76) cells were recorded per session in rats 1, 2, and 3, respectively. To compare results across sessions, the EV measure and each of the four characterization parameters for each session were expressed as z-scores (means and standard deviations of the un-normalized parameters for each of the three rats can be found in Table 1). These parameters were subsequently evaluated for their ability to predict the magnitude of the EV using a least-squares regression analysis. For the first two rats, we also used a template-matching method as an alternative measure of memory reactivation. The number of template matches occurring during each time window was tallied, and a regression analysis for the characterization parameters was performed in the same way as for the EV. Because the template-matching method is extremely computationally intensive and time consuming, the analysis was only run on a subset of the sessions for two of the rats.

In Figure 6 the linear regression plots for the individual parameters are shown for each of the three rats. For all three of the rats, the number of down states, the number of K-complexes, and the amount of LVSs in a time window were significantly positively correlated with the EV value (p < 0.05). In contrast, in all three rats, the amount of HVSs was significantly, but weakly, inversely correlated with EV.

Linear regression plots for EV. In each plot, the data points for all sleep epochs and all recording sessions are compiled. The z-scores of the explained variance, EV, are on the y-axis; the z-score of the number of down states, the number of K-complex events, the number of LVS events, and the number of HVS events are on the x-axis. All correlations are significant (p < 0.05). The results for rat 1 include data from 33 d and a total of 359 data points. The EV ranged between 0.0 and 0.27. The results for rat 2 include data from 3 d and total of 72 data points. The EV value ranged between 0.0 and 0.19. The results for rat 3 include data from 12 d and 583 data points. The EV value ranged between 0.0 and 0.35. For all three rats, the EV was significantly positively correlated with the down states, K-complexes, and LVSs. HVSs are anticorrelated with EV, signifying that the two different kinds of spindles have different functional implications.

Figure 7 shows the linear regression plots for each of the characterization parameters and the number of template matches in S3 for seven sessions recorded in the first rat (first column), and for three sessions recorded in the second rat (second column). These results are consistent with the EV results for all sessions; the number of down states and K-complexes and the amount of LVSs are positively correlated with the number of template matches. The amount of HVSs is weakly inversely correlated with the number of template matches in rat 2 and uncorrelated with the number of template matches in rat 1. Thus, both measures of reactivation produce similar results. That this should be the case is not surprising, as the two measures are strongly correlated with each other (Fig. 7).

Linear regression plots for template matching. Data for rat 1 are shown in column 1. In each plot, all data points from the third sleep epoch of seven recording sessions (52 data points) are compiled. The z-scores of the template matches are plotted on the y-axis; the z-scores of the down states, K-complexes, LVSs, HVSs, and EV are on the x-axis. The number of template matches ranged from 8 to 1325. A similar set of plots for rat 2 are shown in column 2. In each plot, all data points from the third sleep epoch of three recording sessions (41 data points) are compiled. The number of template matches ranged from 41 to 399. With the exception of the template match/HVS relationship for the first rat (marked with an asterisk), all correlations were significant (p < 0.05). The relationships between the template matching and the various events are similar to the relationships between the EV and those events, and the EV and template-matching measures are also highly correlated. This allows for a consistent interpretation of the correlations (or lack thereof) between the sleep events and memory replay.

Sleep features in prebehavioral and postbehavioral sleep

A t test was used to determine whether the incidence of the different sleep characteristics was different between S1 and S3. As Table 2 shows, the number of down states, the number of K-complexes (except in rat 2), and the amount of LVSs were significantly greater in S3 than in S1. In contrast, the amount of HVSs was not significantly different between S1 and S3. This suggests that while down states, K-complexes, and LVSs are affected by the behavior, HVSs are relatively unaffected. Thus, K-complexes and LVSs, which are associated with reactivation, are much more numerous after a task known to induce memory replay, but HVSs, which are actually negatively associated with reactivation, do not increase in frequency.

Relationship of sleep features to each other

The four characterization parameters are not linearly independent of each other, as shown in Table 3. As it has been previously suggested that down states and K-complexes are different manifestations of the same event (Amzica and Steriade, 1997, 2002), it is not surprising that these two measures are strongly correlated with each other. Because down states and K-complexes were extracted using completely independent methods (one based on LFP, the other on the sum of spikes from all cells), this correlation, as well as the analysis shown in Figure 1, stands as confirmation of the observations of Amzica and Steriade (1997). Our methods may have strengthened correlations between K-complex activity and LVSs and weakened those between K-complexes and HVSs because the presence of K-complex activity was used to discriminate between HVSs and LVSs. Hence, correlations between characterization parameters are inconclusive with respect to the relationship between K-complexes, HVSs, and LVSs. This fact does not detract, however, from our primary findings of a relationship between EV and K-complexes and LVSs. The similar results for down states, LVSs, and K-complexes confirm that the methods used to extract these parameters, though computationally different, are, at least, internally consistent and, therefore, more likely to be reliable.

Discussion

These results indicate that the reactivation of memory traces occurs preferentially during periods with a high density of down states and associated K-complexes and low-voltage spindles. This conclusion is supported by the positive correlation of the EV, a measure of reactivation strength, with the number of down states, K-complexes, and time spent in LVSs. These three variables were also correlated with the number of template matches, an alternative measure of memory reactivation that is computed using entirely different principles from the correlation-based EV measure. In contrast, the amount of time spent in HVSs was weakly inversely correlated with the magnitude of EV and the number of template matches.

In recent years a great deal of interest has been focused on NREM sleep as more evidence emerges to support its role in learning and memory. Long considered a quiescent state, it has now been demonstrated that during NREM sleep cortical neurons are actually quite active, although sensory inputs are dramatically diminished. Some of this activity can be attributed to reactivation of previously stored memory traces (Wilson and McNaughton, 1994; Qin et al., 1997; Hoffman and McNaughton, 2002; Euston et al., 2007; Ji and Wilson, 2007), a process that is thought to serve memory consolidation (Marr, 1971; McClelland et al., 1995; Buzsáki, 1998). Large-amplitude, low-frequency corticothalamic oscillations of the LFP are the hallmark of this sleep state (Steriade, 2003). The oscillations themselves, especially spindles, have been hypothesized to play an important role in memory (Steriade and Timofeev, 2003), a theory that is supported by a growing body of evidence (Gais et al., 2002; Eschenko et al., 2006). Our results reveal a correlation between memories representing recent experience in the prefrontal cortex and network down states, K-complexes, and low-voltage spindling. Notably, replay is negatively correlated with high-voltage spindling. These results shed light on which brain states are associated with memory replay and, when combined with other findings, suggest a correlation between memory replay and memory consolidation.

One of the important contributions of the present findings is to sharpen the distinction between what we call high-voltage and low-voltage spindles. Two distinct types of spindles were described by Contreras and Steriade (1996) and Contreras et al. (1997), a waxing and waning type and an exclusively waning type; in our terminology these correspond to HVSs and LVSs, respectively. The authors described that LVSs were grouped by the slow oscillation and were preceded by a depth-positive EEG wave (K-complex). HVSs gradually waxed by entraining more neurons and then desynchronized the same way. The maximum amplitude of the oscillation was dependent on the level of background activity. The authors reported that, although both types of spindles occur in natural sleep, under ketamine–xylazine anesthesia only LVSs are observed. Under barbiturate anesthesia, on the other hand, spontaneous LVSs are not observed, but spontaneous HVSs are. This difference was ascribed to the different levels of corticothalamic network activity observed under the two drugs (barbiturate = low activity, ketamine–xylazine = high activity). The conclusion of the authors was that the cortex is very important in shaping and synchronizing thalamically generated spindles.

Both LVSs and HVSs are thalamically generated oscillations; damage to the reticular nucleus of the thalamus abolishes both phenomena (Steriade et al., 1985; Buzsáki et al., 1988). The frequency ranges of HVSs and LVSs are also close together and sometimes overlap; HVSs predominantly fall within a narrow 7–8 Hz frequency band with several prominent harmonics, whereas LVSs cover a wider and higher frequency range, ∼10–20 Hz, but may be as low as 6 Hz. The current source density profiles for HVSs and LVSs are also similar (Kandel and Buzsáki, 1997). Despite these similarities, however, the two oscillations can be qualitatively separated based on waveform shape, amplitude, duration, and the extent to which they entrain the population of cells. Figure 4 demonstrates that LVSs and HVSs can also be quantitatively separated. This study adds a possible functional distinction between the two oscillations as well: LVSs are associated with memory trace replay and possibly memory consolidation, while HVSs play some other yet-to-be-determined role.

There are several reasons to believe the cortical spindling plays an important role in memory. Unfortunately, however, the distinction between HVSs and LVSs is not always clear in the literature. Delineating the two spindle types is quite important, as our results show that the two oscillations have different implications concerning reactivation and memory. Therefore, the classification of LVSs and HVSs must be inferred where it is not explicitly stated.

In humans, increased spindle density in the 12–15 Hz band is correlated with learning a declarative memory task (Gais et al., 2002), and also with retention of verbal memories (Clemens et al., 2005). In rats, a corresponding increase in spindle power was observed after learning, and after relevant recall of a remote memory (Eschenko et al., 2006). Also in rats, stimulation of synapses that had been previously potentiated by tetanizing stimulation more reliably evoked 7–14 Hz spindles than stimulation of unpotentiated synapses (Werk et al., 2005).

Spindles in the neocortex (both 7–14 Hz and 12–18 Hz), which we have shown to be correlated with neocortical reactivation, have also been temporally correlated with hippocampal sharp waves (Siapas and Wilson, 1998; Sirota et al., 2003). Hippocampal sharp waves have also been correlated with cortical transitions from down to up states and with the slow oscillation (Battaglia et al., 2004). These correlations are interesting not only because hippocampal reactivation occurs preferentially during sharp waves (Wilson and McNaughton, 1994; Kudrimoti et al., 1999), but also because suppression of sharp waves has been shown to impair spatial memory (Girardeau et al., 2009). These results, together, reinforce the connection between neocortical reactivation, hippocampal reactivation, and learning, consistent with consolidation theories (McClelland et al., 1995). In humans, induction of the slow oscillation by transcranial stimulation increased the number of sleep spindles (8–12 Hz) after learning a paired-associate task (Marshall et al., 2006). Remarkably, stimulation of the slow oscillation also significantly enhanced memory retention. Thus, it appears that LVSs and the slow oscillation (and therefore the K-complex) are involved in the neocortical–hippocampal dialogue in a manner that facilitates learning.

Our data are consistent with a growing body of literature suggesting that HVSs may be associated with memory deficits. HVS activity in rats increases with age (Radek et al., 1994), while memory is known to decline (Burke and Barnes, 2006). In aged rats, density of HVSs correlates with passive avoidance retention deficits (Aldinio et al., 1985) and spatial memory impairment (Radek et al., 1994). Alternatively, these results may simply reflect a decrease in the time spent in LVSs, which, based on the above discussion, might also lead to memory deficits. Rates of LVSs were not reported in these studies.

The results presented here are correlative and cannot be considered a proof; however, they are consistent with the postulate that stored-trace reactivation contributes to the establishment of long-term memory. Our data show that LVSs, down-to-up state transitions, and K-complexes are positively correlated with reactivation. A separate body of literature has shown that these same sleep features are predictive of learning and memory (Gais et al., 2002; Clemens et al., 2005; Eschenko et al., 2006; Marshall et al., 2006). Our findings thus provide possible link between these two domains.

Approximately 50% of the inter-down state intervals occur within the 1–2 s range, which is also approximately the range of LVS duration. Therefore, although further analysis is required, it appears likely that reactivation actually occurs during the low-voltage spindle oscillations themselves.

Our results, in addition to the studies cited above, indicate that the cortical activity state is a determining factor in when memory replay can occur. This conclusion is in agreement with the conclusion of Marshall et al. (2006) that the slow oscillation/K-complexes/network down state is causally involved in memory replay and consolidation. The reason for this causal link is unknown, but it is reasonable to speculate that the complete silencing of neocortical activity during down states facilitates the retrieval of weakly biased attractors corresponding to the most recent memories. According to the idea that the sleeping cortex passes through complex sequences of states reflecting learned or innate asymmetries in its weight matrix, this ongoing activity may prevent the network from converging on recent memory states. Briefly silencing the activity may therefore enable the network to reset to the strongest attractor states, which would typically be the most recently established ones.

Note added in proof.

While this manuscript was under review, an article with relevant content was published (Peyrache et al., 2009).

Footnotes

This work was supported by Grant NS020331 from the National Institute of Neurological Disorders and Stroke, by Grant MH046823 from the National Institute of Mental Health, and by the Alberta Heritage Foundation for Medical Research. We thank P. Lipa, R. Roop, G. Van Acker, J. Dees, L. Ash, A. Uprety, S. Van Rhoads, A. Egurrola, J. Wang, and J. Meltzer for help with data collection and analysis.

Correspondence should be addressed to Dr. Bruce L. McNaughton,
Canadian Centre for Behavioural Neuroscience, The University of Lethbridge, 4401 University Drive West, Lethbridge, AB T1K 3M4, Canada.bruce.mcnaughton{at}uleth.ca