Abstract
Neuronal oscillations allow for temporal segmentation of neuronal spikes. Interdependent oscillators can integrate multiple layers of information. We examined phase–phase coupling of theta and gamma oscillators in the CA1 region of rat hippocampus during maze exploration and rapid eye movement sleep. Hippocampal theta waves were asymmetric, and estimation of the spatial position of the animal was improved by identifying the waveform-based phase of spiking, compared to traditional methods used for phase estimation. Using the waveform-based theta phase, three distinct gamma bands were identified: slow gammaS (gammaS; 30–50 Hz), midfrequency gammaM (gammaM; 50–90 Hz), and fast gammaF (gammaF; 90–150 Hz or epsilon band). The amplitude of each sub-band was modulated by the theta phase. In addition, we found reliable phase–phase coupling between theta and both gammaS and gammaM but not gammaF oscillators. We suggest that cross-frequency phase coupling can support multiple time-scale control of neuronal spikes within and across structures.
Introduction
Neuronal oscillations are a natural mechanism for forming cell assemblies (Hebb, 1949) since most brain oscillations are inhibition based, and inhibition can segregate spike train messages (Buzsáki, 2010). The cerebral cortex generates multitudes of oscillations at different frequencies, but how the various rhythms interact with each other is not well understood. A well-studied mechanism is cross-frequency coupling. As described first in the hippocampus, the phase of theta oscillations biases the amplitude of the gamma waves [phase–amplitude (P–A) coupling or “nested” oscillations] (Buzsáki et al., 1983, 2003; Soltesz and Deschênes, 1993; Bragin et al., 1995; Lisman and Idiart, 1995; Mormann et al., 2005; Canolty et al., 2006; Demiralp et al., 2007; Sirota et al., 2008; Tort et al., 2008, 2009, 2010; Sauseng et al., 2009; Wulff et al., 2009). The principle of cross-frequency P–A coupling generalizes to many frequency bands (Kandel and Buzsáki, 1997; Penttonen et al., 1999; Schroeder and Lakatos, 2009; Canolty and Knight, 2010); for example, phase modulation of gamma waves by α (Palva et al., 2005), delta (Lakatos et al., 2005), slow (Steriade et al., 1993; Sanchez-Vives and McCormick, 2000; Isomura et al., 2006; Mukovski et al., 2007), and ultraslow (Leopold et al., 2003) oscillations has been observed in multiple neocortical structures (for review, see Buzsáki, 2006; Jensen and Colgin, 2007; Schroeder and Lakatos, 2009; Fell and Axmacher, 2011).
While P–A coupling can reflect spatiotemporal organization of cell assemblies (Lisman and Idiart, 1995; Varela et al., 2001; Buzsáki, 2010), its mechanisms are not well understood. It was hypothesized that P–A coupling between theta and gamma rhythms occurs because perisomatic basket cells contribute to both rhythms by firing theta rhythmic trains of action potentials at gamma frequency (Buzsáki et al., 1983; Bragin et al., 1995; Traub et al., 1996; Kamondi et al., 1998; Penttonen et al., 1998; Gloveli et al., 2005; Rotstein et al., 2005). This hypothesis predicts that temporal coordination of both rhythms by the basket cells also introduces a phase–phase (P–P) coupling relationship (or phase synchronization) (Tass et al., 1998) between theta and gamma oscillations. The coding advantages of cross-frequency P–P coupling have been illustrated by computational models (Lisman and Idiart, 1995; Varela et al., 2001; Fell and Axmacher, 2011) and experiments on humans (Palva et al., 2005; Sauseng et al., 2009). The precise temporal coordination of neuronal spikes at multiple time scales by cross-frequency phase–phase coupling may be beneficial for both transferring information and spike timing-dependent plasticity (Gerstner et al., 1996; Markram et al., 1997; Fell and Axmacher, 2011). In this report, we describe P–P coupling between hippocampal theta and different sub-bands of gamma rhythms.
Materials and Methods
Detailed information of the various methods used in this paper has been published previously (Diba and Buzsáki, 2008; Mizuseki et al., 2009). Briefly, six male Long–Evans rats (250–400 g) were implanted with four- or eight-shank silicon probes in the dorsal hippocampus. Individual shanks (200 μm apart) had eight staggered recording sites (20 μm vertical distance). The shanks were aligned to the septotemporal axis of the hippocampus, located 3.5 mm posterior to bregma and 2.5 mm from midline. An additional silicon probe was also implanted in the entorhinal cortex in three rats. Two stainless-steel screws were located above the cerebellum and used as ground and indifferent electrodes during recordings. After completion of the experiments, the animals were deeply anesthetized and perfused through the heart with 0.9% saline followed by 10% formaline solution. All protocols were approved by the institutional Animal Care and Use Committee of Rutgers University.
Data collection and analysis.
A week after recovery from surgery, electrophysiological recordings were performed under two different behavioral states. For the active running state, animals were required to run back and forth for a 30 μl water reward on both ends of an elevated linear track (250 × 7 cm). For the sleep state, neuronal activity was recorded during sleep in the home cage before and after the task. During the recording session, the signals were amplified (1000×), bandpass filtered (1 to 5000 Hz), and acquired continuously at 20 kHz on a 128-channel DataMax system (16-bit resolution; RC Electronics). The wideband signal was downsampled to 1250 Hz and used as the local field potential (LFP) for further analysis. Throughout this paper, LFP positive polarity is upward. For off-line spike sorting, the wideband signal was high-pass filtered (800 Hz), and isolated waveforms were resampled.
Neurophysiological and behavioral data were explored using NeuroScope (http://neuroscope.sourceforge.net) (Hazan et al., 2006). Spikes were sorted in two steps, first automatically, using a custom made software, KlustaKwik (http://klustakwik.sourceforge.net) (Harris et al., 2000), followed by manual adjustment of the clusters (using Klusters software package; http://klusters.sourceforge.net) (Hazan et al., 2006). Only units with clear refractory periods and well-defined cluster boundaries were included in the analyses (Harris et al., 2000). Additional details of recording, unit separation, and cell-type classifications have been described previously by Csicsvari et al. (1999b) and Mizuseki et al. (2009).
Wavelet analysis.
For wavelet analysis, the discrete wavelet transform (65 levels, 1–150 or 1–300 Hz) was computed by using a MATLAB wavelet software package provided by C. Torrence and G. Compo (http://paos.colorado.edu/research/wavelets/software.html). Each band (i.e., frequency) of the wavelet transform was individually normalized by the mean and SD of wavelet power during each session.
Multitaper FFT.
For spectral analysis of oscillatory patterns, we used a modified version of the multitaper FFT MATLAB package (Mitra and Pesaran, 1999). When whitening was applied, a low-order autoregressive model was used (Mitra and Pesaran, 1999). We use an FFT window length of 1024 points and a 75% overlap over frequencies ranging from 1 to 600 Hz.
Slow-wave and rapid eye movement (REM) sleep were separated as described previously (Sirota et al., 2008; Mizuseki et al., 2009). Briefly, theta segments were first automatically detected using the ratio between the power of the theta band (4 to 12 Hz) and the power of nearby bands (1–4 and 12–20 Hz), followed by manual adjustment with the aid of visual inspection of the whitened power spectra using a low-order autoregressive model (Mitra and Pesaran, 1999) and raw traces.
Waveform-based estimation of theta phase and phase precession.
To examine the waveform of CA1 theta oscillations, we first applied a broad bandpass filter (1–80 Hz) and determined local minima and maxima in the theta range. Local minima and maxima were then used to determine the duration of the ascending and descending part of the theta cycle to calculate the asymmetry index. For the determination of theta waveshape phases, we calculated the median LFP over all CA1 recording channels that were sufficiently correlated (threshold, r > 0.7) to avoid inclusion of channels with strong phase shifts. From this median LFP signal, local minima and maxima were determined to extract the waveform phase by linear interpolation between them, with minima equaling 0° and maxima equaling 180°. With this approach, independent of the length of the ascending and descending parts of the theta cycle, each part always covers 180°. Other phase estimates and phase precession analyses were done as described previously (Schmidt et al., 2009).
To refine the waveform-based estimation of theta phase, the filtered signal (1 to 60 Hz) was used to find peaks and troughs as well as ascending and descending zero crossings of theta waves. Peaks (180°) of theta waves were identified as local maxima and troughs (0°) as local minima. Ascending points (90°) were identified as the zero crossings of the signal between the trough and peak, and descending points (270°) were identified as the zero crossings of the signal between peak and trough. The waveform-based theta phase was obtained by interpolating phase values between these phase quadrants.
Phase–phase plots.
Phase-locking statistics were used to estimate the significance of the phase covariance between theta and gamma signals without the interference of amplitude correlations (Lachaux et al., 1999). Phase–phase locking between theta and gamma oscillations was determined from the instantaneous phase series of the filtered theta and gamma bands. Statistical significance was estimated by temporal shuffling of the signals relative to each other. To construct shuffled phases, the detected gamma phases were shifted between 1 and 200 ms. Phase–phase plots were then constructed using the shifted gamma phases and the nonshifted waveform-based theta phases. This procedure was repeated 1000 times to compute mean and standard deviation for each bin of the phase–phase plots.
We also analyzed the data by filtering the signal between 1 and 28 Hz (instead of 1 and 60 Hz), and the results were largely the same (data not shown). One caveat of the filtering method is that large gamma waves near the trough or peak of theta may bias the detection of theta troughs and peaks, despite the filtering. To overcome this potential caveat, five controls were performed. (1) An altered signal was used. An altered version of the original signal was used to compute the waveform-based theta phases for the phase–phase plots. In this procedure, the original LFP signal was filtered and separated into three bands (1–25, 25–120, and 120–500 Hz). The 25–120 Hz band filtered signal was shifted randomly (100 ± 10 ms) and summed with the 1–25 and 120–500 Hz filtered derivatives of the original signal. The peak, trough, and zero crossings of the theta waves were then determined from this altered trace with a random relationship between gamma and theta waves. The goal of the procedure was to exclude the possibility that large amplitude gamma waves coinciding with the trough or peak of theta would bias theta phase detection. In the analysis, waveform-based theta phase of this altered signal was calculated and used to calculate the phase–phase plots with the original gamma phases. (2) Signal A was compared with Signal B. The phase–phase plots in this case were calculated using the gamma phases from one shank of the silicon probe and the waveform-based theta phases obtained from a shank at least 600 μm away. (3) The waveform-based theta phase was computed from the mean of the signals from all shanks, except from the one used for the calculation of gamma phases. (4) Phase–phase plots were calculated with the waveform-based theta phases calculated from the LFP recorded in Layer 3 of the entorhinal cortex against gamma phases detected in CA1. The rationale of this approach was that while theta waves between the two structures are highly coherent, coherence between gamma waves are low (Chrobak and Buzsáki, 1998). (5) Theta phases were obtained from the wavelet band between 6 and 10 Hz with maximum power.
n:m phase–phase locking.
Phase locking between two oscillators (Oscillator n and Oscillator m) occur in an n:m ratio when there are m cycles of the “driven” oscillator for every n “stimuli” (Guevara and Glass, 1982). In n:m phase locking, one observes m events associated with the driven cycle occurring at n different times or phases in the “stimulus” cycle (Tass et al., 1998). Although both theta and gamma oscillations co-occur spontaneously, and thus there is no external driving stimulus, their relationship can be analyzed similarly to forced coupling since the n:m pattern is a repeating sequence periodic in time. The “mean resultant length” of the distribution of the difference between phases [Δphasen:m(t) = n * theta phases(t) − m * gamma phases(t)] was calculated for different n:m ratios, e.g., 1:1, 1:2, 1:3, 1:4, 1:5, etc. When a distribution of n:m combination deviates from a uniform distribution (high radial distance valued, with r = 0 for uniform and r = 1 for a perfect unimodal distribution; Rayleigh test for uniformity), it means that the phase difference (Δphasesn:m) oscillates around some constant value. These distributions are regarded to show synchronization between the signals. This measure is independent of the amplitude of the signals (Tass et al., 1998) and related to “Arnold tongues” or synchronization tongues (Guevara and Glass, 1982). These tongues are regarded as regions of synchronization (García-Alvarez et al., 2008).
Cell classification.
Principal cells and interneurons were separated based on their autocorrelograms, waveforms, and mean firing rates (Csicsvari et al., 1999b, Mizuseki et al. 2009). Additionally, monosynaptic interactions between pairs of units were detected using a nonparametric significance test based on jittering of spike trains. A significant peak with a latency between 1 and 5 ms was considered to be due to an excitatory connection while a trough between 1 and 5 ms for at least 1 ms was considered indicative of an inhibitory connection (Fujisawa et al., 2008).
Spike removal.
Spike removal was performed on the wideband signal (1–6 kHz; 20 kHz sampling rate) separately for pyramidal cells and interneurons. Mean waveforms of single spikes, first spikes of bursts (bursts were defined as spikes series <6 ms interspike intervals), and subsequent spikes within bursts were calculated. The mean waveform was then subtracted from the raw signal at the corresponding spike times (window, 1.9 ms before and 2.9 ms after spike trough). As the mean waveform differs from the exact shape of each spike, at every spike time a smooth version of each remaining waveform was subtracted (window, 0.2 ms before and 0.25 ms after spike time). Note that while the “despiking” method can illustrate the contribution of spikes to the LFP spectrum, it does not provide a reliable quantitative measure because it depends strongly on both the number and amplitude of the detected units. The numerous more distant neurons whose spikes remain below the clustering threshold also contribute to the power spectrum, but their contribution cannot be assessed by this approach. Furthermore, the contribution of the spike to the power spectrum also depends on the definition of the spike duration. However, since spike afterpotentials (i.e., nonsynaptic events) coincide in time with synaptic, network-mediated events, quantification of their contribution to the spectral power requires more complex methods.
Results
LFPs and unit firing were recorded in the hippocampal CA1 pyramidal layer. Recordings were performed while the animal ran on a linear maze (250 cm long; hereafter, theta periods during behavioral tasks are referred to as RUN). We also recorded during sleep in the animal's home cage, and each session contained at least one REM episode (Diba and Buzsáki, 2008; Mizuseki et al., 2009).
Asymmetry of theta waves
Theta oscillations in the LFP are usually studied after bandpass filtering between 4 and 12 Hz, generating a largely sinusoidal signal. However, the hippocampal theta rhythm deviates significantly from a harmonic (sinusoidal) oscillation. Theta waves are asymmetric with an inverse saw-tooth appearance in the CA1 region (Buzsáki et al., 1985). To quantify the waveform asymmetry of CA1 theta oscillations, a broad bandpass filter (1–80 Hz) was applied, followed by the detection of local minima and maxima in the theta period range (Fig. 1A). Theta waveform asymmetry was quantified by the “asymmetry index,” defined as the ratio of the duration of the ascending part and the duration of the descending part, expressed on a logarithmic scale (asymmetry index of zero corresponds to a symmetric theta shape). The mean asymmetry index was negative (Fig. 1B) and consistent across rats, indicating shorter ascending and longer descending parts of theta, and the index was significantly more negative during RUN than REM (RUN, −0.27 ± 0.0004; REM, −0.13 ± 0.0004; p < 0.001; t test). Figure 1C illustrates why the waveform-based phase yields different phase estimates than the Hilbert (narrow bandpass, 4–12 Hz) phase. The estimate based on the Hilbert transformation reflects equalized time within theta cycles. In contrast, the waveform-based phase squeezes and stretches time commensurate with waveform asymmetries. During the ascending part of the theta cycle, the phase changes faster (in time), while during the descending part it changes slower (Fig. 1C). A practical consequence of the asymmetry is that the trough of the theta cycle occurs later than in the bandpass-filtered (4–12 Hz) theta wave (Siapas et al., 2005). Furthermore, estimation of phase from Hilbert also assumes a “time-warped” signal, i.e., that the theta waveshape (symmetry) remains the same regardless of its frequency. This is not the case since wave asymmetry increases with theta frequency (Buzsáki et al., 1985).
Temporal coding properties of spikes relative to theta phase (“phase precession”) (O'Keefe and Recce, 1993) are typically examined using the phase of the narrow bandpass-filtered LFP (4–12 Hz; Fig. 1D). To illustrate the relation of phase precession and the theta waveform, we assigned spikes to the ascending and descending phases of theta waves, respectively (Fig. 1D). Spikes that occur on the ascending part of the theta cycle represent largely the early part of the place field. The spike density showed a skewed phase distribution relative to the phase of the CA1 LFP (Fig. 1E). Importantly, theta asymmetry matched the asymmetry observed in the firing-rate histogram (Geisler et al., 2010).
To examine the importance of the theta waveform-based phase estimation method, phase precession of spikes was calculated in relation to both the bandpass phase and waveform-based phase. For each cell, spike phases were correlated with the relative position of the rat in the place field (Fig. 1F). In the resulting distributions of correlation coefficients, pooled across animals, the phase-position correlation was significantly more negative (reflecting a more linear phase precession) for waveform-based phases than for those derived from the bandpass-filtered LFP (Fig. 1F; one-sided Kolmogorov–Smirnov test; p = 3.4 × 10−15). These differences were equally robust if phase precession was calculated on single trials (p = 2.7 × 10−38) (Schmidt et al., 2009). The effect of the theta waveform on phase precession indicates the computational significance of the theta asymmetry for the place representation of CA1 pyramidal cells.
These observations underscore the importance of incorporating waveform information into phase estimates of LFP waves (Buzsáki et al., 1985). Since theta–gamma coupling has been suggested to be involved in neuronal computation (Bragin et al., 1995; Lisman and Idiart, 1995; Jensen and Lisman, 1996; Chrobak and Buzsáki, 1998; Canolty et al., 2006; Tort et al., 2008, 2010), we reassessed the relationship between theta and gamma oscillations, using waveform-based estimation of phase. To improve theta phase estimation further, the peak, trough, and zero (baseline) crossings to both ascending and descending directions were determined for each theta wave (1–60 Hz, Butterworth filter), and each theta quadrant represented one quarter (90°) of the theta period (Siapas et al., 2005). Within the quadrants, the phase was linearly interpolated.
Theta phase modulation of multiple gamma frequency oscillations
While “gamma” rhythm originally referred to a band of 35–85 Hz (Bressler and Freeman, 1980), over the past several years it has been expanded to cover higher frequencies (90–140 Hz) known as “fast” or “high” gamma (Fig. 2A) (Csicsvari et al., 1999a; Canolty et al., 2006; Sullivan et al., 2011). The term “epsilon” was also suggested for this higher-frequency band (Freeman, 2007). These respective frequency bands are locked to different phases of the theta cycle (Fig. 2) (Csicsvari et al., 1999a; Canolty et al., 2006; Colgin et al., 2009).
Analyzing the relationship between gamma power (amplitude of the complex Hilbert transform) and theta phase confirmed that the largest gamma power was associated with the descending phase of theta oscillations in the CA1 pyramidal layer during both RUN and REM (Fig. 2C) (Buzsáki et al., 2003; Csicsvari et al., 2003; Colgin et al., 2009). In contrast, the maximum power of fast gamma (epsilon, ϵ) occurred at the trough of the theta cycle during RUN and shifted to the peak during REM, largely reflecting the maximum firing probability of pyramidal neurons (Fig. 2B) (Montgomery et al., 2008). These findings indicate that the mechanisms by which the different gamma bands are modulated by theta oscillations are state dependent.
To explore the relationship between theta and gamma waves further, the wavelet distribution in the 30–150 Hz band was calculated and related to the waveform-based theta phase (Fig. 3A,B). Visual inspection of these plots identified three distinct gamma bands associated with different phases of theta waves during RUN. We refer to these frequency bands as slow gamma (gammaS; peak power, 40.5 Hz), midfrequency gamma (gammaM; peak power, 60.6 Hz), and fast gamma (gammaF, or ε band; peak power, 118.9 Hz) (Csicsvari et al., 1999a; Canolty et al., 2006; Sullivan et al., 2011) oscillations. GammaM power was largest near the peak (200.1°), and gammaS on the descending phase (251.1°), while gammaF was more prominent near the trough (340°) of the theta cycle. The span of the gamma sub-bands were state dependent because during REM, the peak frequency of gammaM was significantly higher (90.71 Hz) compared to that during RUN (Fig. 3A; p < 0.001; n = 11 sessions in 3 rats; t test), but its theta phase preference remained the same (195.4°) as during RUN. The preferred theta phase of gammaF (ε band) during REM shifted to near the peak (199.4°), giving a continuous band at this theta phase from 50 to 150 Hz (Figs. 2C, 3A). Analysis of the distribution of gamma oscillation episodes within individual theta cycles showed that the majority of theta waves (41.1 ± 1.5%) contained both gammaS and gammaM, and 37.3 ± 0.7% contained gammaM only, whereas only a small fraction of theta cycles (12 ± 0.2%) was dominated by only gammaS (Fig. 3B). The remaining 9.6 ± 0.9% of theta waves did not contain detectable gamma oscillations.
In the time domain, the troughs of the filtered gamma cycles (30–50 Hz for gammaS and 50–90 Hz for gammaM) were used as a reference to generate trough-triggered average LFP (1–1250 Hz) and average theta waves (4–12 Hz; Fig. 3C). In another approach, the theta-filtered (4–12 Hz) LFP was cross-correlated with the power of either for gammaS or gammaM (Fig. 3D). These analyses confirmed that gammaM peaked significantly earlier than gammaS during both RUN and REM (p < 0.01; t test; n = 11 sessions in 3 rats).
Sublayer specificity of gamma oscillations
We found that gammaS and gammaM show distinct preferred depths of the CA1 pyramidal layer (Fig. 4). The vertical span (140 μm) and the precise distribution of the eight recording sites on the probe shanks (20 μm vertical steps) allowed for the determination of the relative depths within the CA1 pyramidal cell layer (Mizuseki et al., 2011). The site with the largest power in the >200 Hz band was regarded as the middle of the pyramidal layer and served as the reference depth for comparison across electrodes shanks and rats. The largest power of gammaM was present in the deep part of the pyramidal layer (i.e., adjacent to the stratum oriens), whereas the power of gammaS was largest in the superficial part, adjacent to the stratum radiatum. GammaF (ε) occupied a narrow zone in the middle of the CA1 pyramidal layer (Fig. 4A,B). We expanded these observations from data recorded by 16-site linear probes (100 μm steps). To identify the dipoles more precisely, the current source densities (CSDs) of the LFPs were calculated. The depth profile of the CSDs of the filtered sub-bands allowed us to predict the potential anatomical sources of these sinks and sources. As expected, gammaS had the largest sink in the mid-stratum radiatum, indicating that this low-frequency gamma oscillation is largely transferred by the CA3 afferents (Bragin et al., 1995; Csicsvari et al., 2003). In contrast, the largest CSD of gammaM was observed in the stratum lacunosum moleculare, indicating the entorhinal origin of this band (Fig. 4C,D) (Charpak et al., 1995; Colgin et al., 2009). The fast (gammaF, ε) gamma waves were confined to the pyramidal cell layer. The theta phase-based distribution of wavelet power at different depths confirmed that these three gamma bands are generated by different inputs and mechanisms (Fig. 4D).
Spike contribution to gammaF (ε)
Because the preferred phase of gammaF power largely corresponded to the maximum spiking of neurons, we examined the possible spectral contamination of the spectral power by spikes (Montgomery et al., 2008; Ray et al., 2008; Whittingstall and Logothetis, 2009; Quilichini et al., 2010; Ray and Maunsell, 2011). The mean waveform of each clustered putative pyramidal cell and interneuron was subtracted from the individual actions potentials of the corresponding neuron in the wideband signal. Comparison of power spectra of the “whitened” (1/f slope removed) LFP showed that spikes of both pyramidal cells and interneurons began to contribute significantly to spectral power above 100 Hz, i.e., in the gammaF (ε) band (Fig. 5). We should note that in these analyses, only spike waveforms of the clustered minority neurons were deleted. The much larger but not identified spiking neuron population surrounding the electrodes likely further contributes to the power. Also, removing wider waveforms, containing also the spike afterpotentials, would likely affect the lower frequency bands as well.
Theta–gamma cross-frequency phase–phase coupling
P–P synchronization between two oscillators can be determined by correlating the instantaneous phase values of the oscillations (Tass et al., 1998). Plotting the waveform-based phase of theta oscillations against the phase of gamma oscillation (filtered 30–90 Hz) during RUN revealed a correlation in every animal tested (Fig. 6A). To examine the statistical reliability of P–P coupling, the filtered gamma and theta traces were randomly shuffled, and the resultant pseudocorrelation (Fig. 6A, second panel) was subtracted from the original theta–gamma correlation matrix (Fig. 6B). The diagonal stripes in the phase–phase plots indicate a statistically reliable relationship between the respective phases of theta and gamma oscillations during both RUN and REM (Fig. 6B).
Assuming two phase-locked sinusoid oscillators with an integer relationship between them (e.g., 1:5), multiple identical P–P coupling magnitudes are expected (e.g., 72, 144, 216, 288, and 360° of the slow rhythm). Instead, the magnitudes of phase coupling were not perfectly distributed within the theta cycle (as reflected by the different heat colors of the stripes), but phase coupling was stronger on the descending phase of theta. This observation is in support of previous findings that theta phase also modulates the amplitude of gamma power (Bragin et al., 1995).
To exclude the possibility that phase coupling occurred spuriously due to the detection of an occasional prominent gamma wave at the trough or peak of waveform-based theta cycle, additional control analyses were performed (Fig. 6C). In the first control analysis, the filtered gamma band (30–90 Hz) was subtracted from the LFP, and the phase values of the gamma oscillations were determined. The gamma signal was shifted randomly and added back to the gamma-subtracted LFP. The trough, peak, and zero crossings of the theta waves were then determined from this altered signal (Fig. 6Ca, altered signal). The goal of this analysis was to reduce the possibility of spurious phase–phase coupling due to biasing the detection of theta troughs/peaks by contaminated gamma cycles. In the second control analysis, the theta phase and gamma phase were measured from the furthest shanks of the silicon probe, respectively (600–1200 μm distance; Fig. 6Cb). Third, gamma was detected from one site, whereas the theta phase was calculated from the sum of LFPs recorded at the sites of the remaining shanks (Fig. 6Cc). The rationale for the second and third control analyses was that gamma amplitude varies significantly with distance (Sirota et al., 2008). Fourth, gamma phase was determined from the CA1 pyramidal layer, whereas theta phase was obtained from another probe in Layer 3 of the entorhinal cortex (Fig. 6Cd). Finally, the phase of both theta and gamma oscillations was determined by the wavelet method (Fig. 6Ce). The shuffled derivatives were subtracted from each of constructed data. Despite these control manipulations, phase coupling between theta and gamma oscillations remained significant in each case (n = 11 sessions in 3 rats). Examination of the relationship between theta and gammaF (ε band) showed no reliable P–P relationship between the two frequency bands (Fig. 6D).
While the phase–phase plots in Figure 6 show prominent diagonal stripes, several stripes are discontinuous and have different phase angles, an indication of the presence of distinct gamma bands with different phase–phase relationship to the waveform-based theta cycle. To explore this qualitative observation quantitatively, we computed radial distance values (R) of the circular distribution from the difference between theta and gamma phases for different n:m relationships (Fig. 7A) (Tass et al., 1998). These plots yielded two large discrete peaks at 5 and 9 n:m ratios, corresponding to possible 5 gammaS and 9 gammaM cycles within a theta period. The relative size but not the position of the peaks depended on the filtered gamma band used for the determination of the instantaneous phase of gamma oscillations (Fig. 7A), reminiscent of “Arnold tongues” (Guevara and Glass, 1982). The existence of two discrete peaks provided additional support for the presence of relatively separate gammaM and gammaS bands within the theta cycle with different theta phase preferences (Fig. 7).
The distribution of the ratios of gamma and theta frequencies, calculated from the peak frequency of the power in the 30–90 Hz band and the reciprocal value of the theta period is shown in Figure 7B. Fitting Gaussians to the distribution revealed peaks at 4.9 ± 0.3 and 8.7 ± 0.6 during RUN and 5.4 ± 0.1 and 9.5 ± 0.3 during REM, supporting the phase ratio-based calculations. The distributions during RUN and REM were not significantly different (p > 0.05, Kolmogorov–Smirnov test) with a median of 7.3 (RUN) and 7.8 (REM). Overall, these findings show that the previously inferred seven gamma cycles per theta cycle (Bragin et al., 1995; Lisman and Idiart, 1995; Chrobak and Buzsáki, 1998) can be more quantitatively described as a theta phase-ordered sequence of gammaS and gammaM cycles nesting in the theta wave, corresponding to 1:5 and 1:9 P–P coupling, respectively.
A prediction of the cross-frequency P–P coupling is that the frequencies of theta and gamma oscillations are correlated (Bragin et al., 1995). The mean frequencies of gammaM and gammaS depended on the period (1/f) of the theta waves, as shown in an example animal during REM (Fig. 8A). The theta frequency versus gamma frequency relationship was clearer with gammaM and more pronounced during REM than RUN (data not shown). Taking advantage of the theta phase segregation of gammaS and gammaM, the gamma power at particular phases (between 125 and 216° and between 216 and 288°, respectively) of each theta cycle was determined and plotted against the frequency (1/duration) of the theta cycle (Fig. 8B). These resulted in significant correlations but with different slopes for RUN and REM (Fig. 8B). In a related analysis, the duration of each theta cycle was plotted against the mean duration of nested gammaS (30–50 Hz) and gammaM (50–90 Hz) waves (Fig. 8C). This analysis again showed a positive relationship between the frequency of theta and gamma rhythms (Fig. 8C; RUN, R2 = 0.46, p < 0.01; REM, R2 = 0.42, p < 0.01). The wave duration analysis showed significantly different theta–gamma coupling ratios in RUN and REM (Fig. 8C; p = 0.013; ANCOVA).
Firing patterns of pyramidal cells and interneurons during gamma oscillations
Since the observed LFP theta–gamma relationship arises from the coordinated timing of membrane potential oscillations at theta and gamma frequencies with consequent phase coupling of their spikes (Buzsáki et al., 1983; Soltesz and Deschênes, 1993; Bragin et al., 1995; Penttonen et al., 1998; Whittington and Traub, 2003; Hájos and Paulsen, 2009), we examined the spike versus LFP phase relationship of pyramidal cells and interneurons. For these analyses, the LFP segments were sampled with the spike centered on the sampling window (Fig. 9). As expected, wavelet maps of spike-triggered LFP samples also showed the temporal separation of gammaS, gammaM, and gammaF (ε) in the theta cycle (Fig. 9B,E). For these analyses, LFP patterns were sampled from the neighboring shank of the probe, relative to shank from which the unit was recorded from to reduce the contribution of the spike to the LFP.
A fraction of pyramidal cells (21.9%; n = 109) were significantly phase locked to gamma oscillations (30–90 Hz). However, when the distinct sub-bands were used for the examination of phase locking, 17.7% of all pyramidal cells (n = 88) were phase locked to gammaS, and 36.5% (n = 182) to gammaM. The increased number of phase-locked cells to gammaM is likely explained by the excess phase jitter present in the wider band. Most pyramidal cells (75.1%; n = 374) were phase locked to the filtered gammaF (90–150 Hz, ε) waves. It was reported previously that spikes of pyramidal cells show phase preference to either the trough (0 to 30° and 240 to 360°; Fig. 9A, gammaT neurons) or the rising phase (30 to 240°; gammaR neurons) of the gamma waves (Senior et al., 2008; Mizuseki et al., 2011). Spike-triggered wavelet maps for the gammaT and gammaR groups show that both groups contributed similarly to gammaS and gammaM (Fig. 9C).
Wavelet maps of the LFP triggered by spikes from putative interneurons also confirmed the theta phase separation of gammaS and gammaM oscillations (Fig. 9E). While we could not reliably classify interneuron classes (Freund and Buzsáki, 1996; Klausberger and Somogyi, 2008), gamma phase preference of interneurons also showed bimodality (Fig. 9D) (Senior et al., 2008). Spike-sampled LFPs of these two groups showed very similar wavelet maps (data not shown). The majority (87.7%; n = 53) of putative interneurons was significantly phase locked to gamma oscillations (30–90 Hz). When examined separately for the different sub-bands, 88.6% were significantly phase locked to gammaS and 90% to gammaM. Most interneurons (93.6%; n = 57) were phase locked to the filtered gammaF (90–140 Hz, ε) waves. The large percentage of interneurons phase locked to multiple bands reflects that most interneurons were phase locked to two or all three patterns. The small numbers of interneurons precluded their further quantitative groupings. By examining each interneuron separately, we observed that nine interneurons were phase locked exclusively to gammaS and two exclusively to gammaM oscillations. Using these subsets, the generated wavelet maps of spike-sampled LFP showed distinct profiles (Fig. 9F).
Discussion
Our findings show an orderly relationship between gamma oscillations and the asymmetric theta waves. GammaS and gammaM oscillations have an integer relationship to the theta rhythm (5:1 and 9:1, respectively, in the waking rat) and, in addition, the power of both gamma rhythms is phase modulated by theta waves. GammaS and gammaM rhythms are distinct patterns from the faster (>90 Hz; ε band) oscillations, which is strongly affected by the spike waveforms of nearby active neurons.
Multiple gamma bands, multiple mechanisms?
What defines the identity of a rhythm? LFP frequency bands have been traditionally and artificially defined by frequency border criteria rather than mechanisms. The term “gamma oscillation” has been suggested to reflect a single mechanism in the olfactory bulb for a frequency range between 35 and 85 Hz (Bressler and Freeman, 1980). Subsequently, the terms “fast gamma” (Csicsvari et al., 1999a) and “high gamma” (Canolty et al., 2006) have been coined to describe a frequency band between 90 and 140 Hz in the cortex (or even up to 600 Hz) (Gaona et al., 2011). Our findings suggest that these two general bands are generated by different mechanisms, and the exact borders of these sub-bands vary with brain state. Whereas the “traditional” gamma oscillations showed P–P coupling with the theta rhythm, no such relationship was found with fast gamma (ε) oscillations. A further complication is that although bona fide bursts of oscillations in the 90–140 Hz band have been described in the CA1 pyramidal layer (Csicsvari et al., 1999a; Colgin et al., 2009; Sullivan et al., 2011), both previous works and the present findings indicate that a large part of the power in this high-frequency band derives from spikes and spike afterpotentials (Montgomery et al., 2008; Ray et al., 2008; Quilichini et al., 2010; Ray and Maunsell, 2011). For these reasons, it may be justified to distinguish the 30–90 Hz and 90–140 Hz bands by different names and refer to the latter as the ε band (Freeman, 2007).
Within the gamma band proper, two sub-bands (gammaS, 30–50 Hz; gammaM, 50–90 Hz) could be distinguished by their P–A coupling to theta waves and layer distribution in the CA1 region. The mechanisms of gammaM and gammaS are not known, and it needs to be established whether and how they are related to the different models of gamma generation (Leung, 1982; Whittington et al., 1995; Traub et al., 1996; Wang and Buzsáki, 1996; Traub et al., 2004; Tiesinga and Sejnowski, 2009). The presence of gammaS in the intact hippocampus is in good agreement with previous studies (Bragin et al., 1995; Hormuzdi et al., 2001; Buhl et al., 2003;Buzsáki et al., 2003; Csicsvari et al., 2003; Tort et al., 2008, 2009, 2010; Colgin et al., 2009; Wulff et al., 2009; Korotkova et al., 2010; Chen et al., 2011). However, in most of these studies, the “high” gamma has been defined from 50 to 150 Hz or higher, thus mixing gammaM and the ε band (but see Tort et al., 2008). On the basis of the present findings, it remains ambiguous whether the behavioral and computational mechanisms attributed to the 50–150 Hz band in those papers can be explained by gammaM oscillations or the ε band. Since the ε band is typically contaminated by spiking activity (Ray and Maunsell, 2011), in future studies, it will be important to examine the selective behavioral correlates and computational roles of the distinct gamma and ε bands.
Phase–amplitude and phase–phase coupling between theta and gamma oscillations
Numerous experiments have reported nested gamma waves in theta cycles (P–A coupling) (Buzsáki et al., 1983, 2003; Soltesz and Deschênes, 1993; Bragin et al., 1995; Mormann et al., 2005; Canolty et al., 2006; Demiralp et al., 2007; Tort et al., 2008, 2010; Sauseng et al., 2009). The computational advantages of cross-frequency P–P coupling (or n:m phase synchronization) (Tass et al., 1998) have been emphasized repeatedly by theoretical models (Lisman and Idiart, 1995; Varela et al., 2001; Fell and Axmacher, 2011), and P–P coupling between theta and gamma rhythms has been observed previously in scalp-recorded EEG (Sauseng et al., 2009). In principle, P–A and P–P coupling can occur independently between two (or more) rhythms, provided that the oscillators are generated by independent mechanisms and linked by some specific conduit (Fell and Axmacher, 2011). On the other hand, if the two oscillators emerge from the same neuronal substrate, they likely share mechanisms, e.g., contribution of perisomatic inhibition to both gamma and theta oscillations (Buzsáki et al., 1983; Bragin et al., 1995; Gloveli et al., 2005; Rotstein et al., 2005). In this case, P–P and P–A coupling may reflect two aspects of a single mechanism. We hypothesize that, in brain networks, P–P and P–A coupling most often co-occur and suggest that the reasons why P–P coupling is rarely reported, relative to the ubiquitous P–A coupling, include the lack of appropriate methods to estimate the waveform-based phase of nonharmonic oscillators and properly segregate frequency bands (Lachaux et al., 1999; Mureşan et al., 2008; Tort et al., 2009; Burns et al., 2011). The importance of the proper estimation of phase is amply illustrated by our observations that relating the waveform-based theta phase of spikes to the current position of the rat provided a stronger relationship than phase estimation based on the assumption that theta is a sinusoidal oscillation.
In computational models, identical gamma cycles with an equal number of spikes in each cycle are distributed across the entire theta cycle to support a multiitem working memory buffer (Lisman and Idiart, 1995; Jensen and Lisman, 1996). Each gamma cycle contains a discrete item (or position in space), and approximately seven gamma cycles can store 7 ± 2 sequential items. As the rat moves forward (or as memory content moves forward), a new item occupies the first gamma cycle, the contents of the cycles shift, and the content of the last gamma cycle in the previous theta period is deleted, analogous to a serial first-in, first-out shift register in electronics. Although the needed multiple nested gamma cycles are present in the theta cycle in vivo (Bragin et al., 1995), the gamma periods vary both in magnitude and duration within the theta cycle. The largest amplitude gammaM cycles are dominant after the peak of the theta cycle, whereas the strongest gammaS occurs on the descending phase. In contrast, the maximum probability of pyramidal spikes in the waking rat is concentrated at the trough of the theta, largely coinciding with the strongest power of gammaF (ε band). The implication of these findings is that if gamma period-associated cell assemblies represent separate items (or positions), the items are not equally weighted. This may be advantageous because the asymmetric nature of the gamma power distribution and unit firing in the waveform-based phase space of the theta cycle may facilitate strengthening of time-forward connections compared with backward connections, an important organizational property of episodic memory (Kahana, 1996).
Footnotes
This work was supported by National Institute of Health Grants NS34994 and MH54671, NSF, the Human Frontier Science Program, the James S. McDonnell Foundation, the Pew Charitable Trust Latin American Fellows Program in the Biomedical Sciences (M.A.B.), and Bundesministerium für Bildung und Forschung (BMBF) Grants 01GQ0410 and 01GQ0972.
- Correspondence should be addressed to György Buzsáki, Center for Molecular and Behavioral, Neuroscience, Rutgers University, 197 University Avenue, Newark, NJ 07102. buzsaki{at}axon.rutgers.edu