Abstract
Neuronal oscillations in the gamma frequency range have been reported in many cortical areas, but the role they play in cortical processing remains unclear. We tested a recently proposed hypothesis that the intensity of sensory input is coded in the timing of action potentials relative to the phase of gamma oscillations, thus converting amplitude information to a temporal code. We recorded spikes and local field potential (LFP) from secondary somatosensory (SII) cortex in awake monkeys while presenting a vibratory stimulus at different amplitudes. We developed a novel technique based on matching pursuit to study the interaction between the highly transient gamma oscillations and spikes with high time–frequency resolution. We found that spikes were weakly coupled to LFP oscillations in the gamma frequency range (40–80 Hz), and strongly coupled to oscillations in higher gamma frequencies. However, the phase relationship of neither low-gamma nor high-gamma oscillations changed with stimulus intensity, even with a 10-fold increase. We conclude that, in SII, gamma oscillations are synchronized with spikes, but their phase does not vary with stimulus intensity. Furthermore, high-gamma oscillations (>60 Hz) appear to be closely linked to the occurrence of action potentials, suggesting that LFP high-gamma power could be a sensitive index of the population firing rate near the microelectrode.
Introduction
Gamma oscillations (40–80 Hz) have been commonly observed in several cortical areas during a variety of cognitive tasks (Roelfsema et al., 1997; Miltner et al., 1999; Fries et al., 2001; Pesaran et al., 2002; Schoffelen et al., 2005; Womelsdorf et al., 2007). Recent invasive EEG studies on epileptic patients have shown further that functional activation of cortex increases power even at frequencies much higher than the traditional gamma range (called “high-gamma”) (Crone et al., 2006, and references therein). However, the role of gamma or these high-gamma oscillations in cortical processing remains unknown.
Here, we test a specific functional role of gamma oscillations in the context of stimulus intensity coding, which was proposed recently (Fries et al., 2007). Some investigators have proposed that gamma oscillations could provide a temporal framework for the firing of neurons, such that information is coded in the timing of spikes relative to the ongoing gamma cycle (Buzsáki and Chrobak, 1995; Buzsáki and Draguhn, 2004; Mann and Paulsen, 2007). Fries et al. (2007) suggested that the gamma cycle could serve as a fundamental computation mechanism, whereby excitatory input to a pyramidal cell is converted into a temporal code relative to the gamma cycle, with stronger inputs leading to earlier responses (this is referred to in this report as the “gamma phase-coding” hypothesis). An analogous concept of phase coding has been documented in the hippocampus with low-frequency theta (4–10 Hz) oscillations (O'Keefe and Recce, 1993; Huxter et al., 2003; Buzsáki and Draguhn, 2004), where the location of the animal in its environment is represented in the phase of the action potentials with respect to theta oscillations. However, this prediction has not been directly tested for cortical gamma oscillations.
We recorded single-unit activity and local field potential (LFP) simultaneously from secondary somatosensory cortex (SII) in awake macaque monkeys while presenting vibratory tactile stimuli at different amplitudes, and we tested whether an increase in the excitatory drive caused by an increase in stimulus amplitude changes the phase relationship between spikes and gamma or high-gamma oscillations. Because these oscillations are highly transient, traditional methods such as short time Fourier transform (STFT) or wavelet transform are not well suited toward studying this phase relationship. We developed a new technique based on the matching pursuit (MP) algorithm (Mallat and Zhang, 1993), which iteratively decomposes a signal into a linear expansion of waveforms selected from a large over-complete dictionary, with waveforms chosen to best match the local signal structures. This allowed a rigorous comparison of LFP signal components associated with spikes under varying stimulus intensity, although these signal components overlapped with each other in both spectral and temporal domains.
We found that neither the gamma nor high-gamma oscillations that were coupled with spikes changed their phase relationship even with a 10-fold increase in stimulus intensity. Thus, stimulus intensity is not coded in the phase of gamma (or high-gamma) oscillations in SII cortex.
Materials and Methods
Animals and animal behavior.
Single-unit recordings were obtained from SII cortex in two macaque monkeys (Macaca mulatta, one female, 4–6 kg). The animals were trained to sit quietly in a chair during the experiment. Water rewards were given at pseudorandom intervals during the experiment to keep the animals awake and alert. All procedures and experimental protocols complied with the guidelines of the Johns Hopkins University Animal Care and Use Committee and the National Institutes of Health Guide for the Care and Use of Laboratory Animals.
Experimental design.
A sinusoidal stimulus was delivered perpendicular to the skin through a probe on the distal pad of the second or third digits (D2 or D3) by a vibratory stimulator (Mini-shaker; Brüel and Kjær model 4810). A minority of the neurons (34 of 203, or ∼17%) had receptive fields on the palm and for those neurons the stimulator was positioned on the palm instead of the finger pads. Because no quantitative difference could be found between the properties of the neurons with receptive fields on the digits and palm, the responses were pooled. The stimulus was presented for 1 s, with an interstimulus interval of 1.2 s. Three different stimulus frequencies (50, 100, and 200 Hz) and four different amplitudes (in the ratio 1:2:5:10 and denoted by G1, G2, G5, and G10 throughout this study) were used, with 50 trials per frequency and amplitude combination. Stimuli were presented in pseudorandom order. Stimulus amplitude at the lowest amplitude (G1) was approximately five times higher than the perceptual threshold (on average ∼5, 3, and 0.8 μm at 50, 100, and 200 Hz, respectively, so that robust firing rates and gamma activity were recorded at all four stimulus amplitudes). Mini-shaker movements were monitored using an accelerometer (Type 8636C5; Kistler Instruments) with a range of ±5 g. The Fourier analysis of the acceleration revealed a single peak at the stimulus frequency, and the corresponding displacement amplitudes were computed by double integration of the acceleration.
Recordings.
Stainless-steel recording cylinders were placed surgically in each hemisphere over SII using aseptic procedures under sodium pentobarbital anesthesia (25 mg/kg, i.v., initial dose plus 5–15 mg/kg/h, i.v.). Neural recordings were done at least 1 week after surgery to allow sufficient time for recovery and eliminate any effect of barbiturates on neural oscillations. Neural activity was recorded using seven platinum-iridium extracellular microelectrodes spaced 400 μm apart and driven by a Reitboeck microdrive (Mountcastle et al., 1991). The signals from the seven channels were amplified by a headstage amplifier (10×) and divided into two streams for the collection of LFP and spikes, respectively. One stream of the incoming signal was amplified (100×) and filtered (0.3–300 Hz, 6 dB/octave) using Grass amplifiers (model 15 LT with 12A54 Quad amplifiers; Grass-Telefactor/Astro-Med). The sampling rate was 5 kHz. The second input stream was bandpass filtered (500–1000 Hz) and amplified (10 to 100×), and spikes were isolated using a window amplitude discriminator. Only neurons for which the action potentials were well isolated from noise were selected for analysis. As a control, to ensure that the observed activity was not an electrical artifact, we repeated the experiment for several neurons after taking the stimulator off the finger. For such trials, all stimulus-induced activity disappeared.
Two hundred and three neurons were recorded from three hemispheres (55 and 71 neurons from hemispheres 1 and 2, respectively, of monkey 1, and 77 neurons from one hemisphere in monkey 2).
Stimulus frequency.
Most neurons responded at all three stimulus frequencies (50, 100, and 200 Hz), although the response was typically much stronger for the stimulus frequency of 50 Hz. This is expected because slowly adapting type 1 (SA1) and rapidly adapting (RA) receptors found in the periphery respond more strongly to stimulus frequencies at or below ∼50 Hz. In this study, only the results for a stimulus frequency of 50 Hz are shown (the results for stimulus frequencies of 100 and 200 Hz are qualitatively similar; data for 100 Hz are also discussed for control studies in the supplemental material, available at www.jneurosci.org).
Time–frequency analysis.
Time–frequency analysis was performed using the MP algorithm (Mallat and Zhang, 1993). MP is an iterative procedure to decompose a signal as a linear combination of members of a specified family of functions gγn, which are usually chosen to be sine modulated Gaussians, i.e., Gabor functions or “Gabor atoms,” because Gabor atoms give the best compromise between frequency and time resolution. In this algorithm, a large overcomplete dictionary of Gabor atoms is first created. In the first iteration the atom gγ0 which best describes the signal f(t) is chosen from the dictionary and its projection onto the signal is subtracted from it. The procedure is repeated iteratively with the residual replacing the signal. Thus, during each of the subsequent iterations, the waveform gγn is matched to the signal residue Rnf, which is the residue left after subtracting the results of previous iterations. Here, 〈x, y〉 = ʃx(t)y(t)dt denotes the inner product of x and y, and Ω denotes the dictionary.
Each Gabor atom can be expressed as follows: Here K(γ) is a normalization parameter and γ = {s,u,ω,φ} stands for the set of parameters of the functions in the dictionary.
Different sampling strategies to sample the continuous parameters {s,u,ω} result in different types of dictionaries. For a signal of length N, a dyadic dictionary with the following values of {s,u,ω}can be used: Following Mallat and Zhang (1993), to further improve the resolution, we used a dyadic dictionary for approximate atom localization and then refined the decomposition with a Newton search on a finer grid depending on the scale [log2(s)] of the atom. In this modification, the following values of {s,u,ω} were used:
This resulted in improved time resolution for small-scale atoms and improved frequency resolution for large-scale atoms. For example, atoms with a small scale (j ≤ 3) responsible for the negative peak of the “sharp transient” (see Fig. 3C) had the resolution equal to sampling interval (0.2 ms in this study).
We also included the Dirac δ function [δ(t), which, for a discrete signal is 1 at t = 0 and zero elsewhere], as well as Fourier atoms (pure sinusoids) along with the Gabor atoms. These atoms were useful for representing highly localized time and frequency components of the signal.
Time–frequency plots were obtained by calculating the Wigner distribution of individual atoms and taking the weighted sum (Mallat and Zhang, 1993). Here, M denotes the number of iterations, and denotes the Wigner distribution of the atom g (here, * denotes the complex conjugate but it is redundant in this case because g is real). This approach is especially useful because it eliminates all cross terms of the actual Wigner distribution of the signal.
The spike-triggered average (STA) was computed by first taking LFP segments between −51.2 and 51 ms centered on spikes (n = 512, at 0.2 ms resolution) and taking their mean (spikes were taken from different time intervals as described below in Analysis periods). For the STA under each stimulus condition, we fitted 100 atoms.
For the computation of time–frequency plots (see Fig. 1), the LFP was down-sampled by a factor of 5. The MP was performed on signals of length n = 2048 bins (−0.523–1.524 s at 1 ms resolution, where zero denotes the time of stimulus onset), yielding a 2048 × 2048 array of time–frequency values (with a time resolution of 1 ms and frequency resolution of 500/2048 Hz, where 500 Hz is the Nyquist frequency after down-sampling). For each signal, we fitted 500 atoms so that high-frequency atoms (generally of a much lower energy) could also be selected: typically this fit accounted for >99.9% of the signal energy. These time–frequency plots were used for the construction of the spike-triggered time frequency average (STTFA) as shown in Figure 4.
Noise and stimulus artifact removal.
With matching pursuits, line noise (60 Hz and harmonics) is represented by long atoms (spread along the time axis) concentrated at ∼60 Hz and its harmonics. Similarly, a sinusoidal stimulus artifact is represented by long atoms concentrated at the stimulus frequency and its harmonics. Such atoms are excluded from the analysis to eliminate these artifacts. This procedure is more powerful than traditional filtering because at least some energy at 60 Hz attributable to a physiological origin is preserved. For example, a short burst of activity, which MP represents as an atom localized in time but spread in frequency, will not be excluded from analysis.
For each neuron, noisy trials were eliminated before analysis by visual inspection and window discrimination of the filtered LFP. For each stimulus frequency and amplitude combination, between 25 and 50 trials were analyzed (on average, approximately three trials of 50 were deleted).
Selection of neurons.
Because the detailed behavioral state of the monkeys was not controlled, we performed a series of statistical tests to ensure that the results were not biased by a change in the level of alertness/attention of the monkey over the duration of recording. First, we tested whether the firing rates as well as the power in the gamma range (details of computation are described below) were stable over time. We performed a regression analysis between the firing rates versus trial number in the baseline as well as during the stimulus period. Similarly, a regression analysis was performed between gamma power and trial number. All neurons for which the first order regression coefficient (either in the rate versus trial or the power versus trial case) was significantly different from zero (p < 0.05, with Bonferroni correction for multiple comparisons) were discarded (n = 50). For these 50 discarded neurons, an equal number of positive and negative regression values were observed (25 instances of each). This selection criterion was very conservative because at low firing rates random fluctuations often resulted in significant regression coefficients.
Next, for every neuron, we divided each set of trials into three groups (early trials, middle trials, late trials) and performed a three-way ANOVA on the rate values as well as the power values, both in the baseline and the stimulus period. All neurons that showed significant differences in firing rate or gamma power in the three intervals (p < 0.05, with Bonferroni correction) were discarded. An additional 12 neurons were discarded from this analysis.
Our rejection criteria were extremely conservative, but they ensured that firing rates and gamma power were stable for the entire recording duration for the remaining 141 neurons. The results shown in this study were obtained from this subset of 141 neurons, although similar results were obtained when the analyses were repeated on the full set of 203 neurons.
Neural classification.
From visual inspection of the raster plots as well as a comparison of the firing rates in the stimulus period versus baseline, the neurons were divided into three categories: excited (78 neurons), inhibited (21), and not-driven (42) (see Fig. 1A). Excited neurons had statistically significant increases in firing rate for the first 200 ms after stimulus onset at stimulus amplitude of G10, compared with baseline (t test, p < 0.05). Neurons were categorized as inhibited if they showed an overall decrease in firing rate relative to baseline while the stimulus was indented in the skin. The remaining neurons were classified as not- driven. The not-driven neurons were included in analysis because, although the firing rates did not change, we found significant changes in gamma, high-gamma, and β (16–24 Hz) power in the LFP after stimulus onset (see Fig. 1B).
Analysis periods.
We studied the spike–LFP correlations at three different time intervals relative to stimulus onset. A closer look at the firing rates and gamma power for excited neurons revealed that the stimulus induced a sharp transient increase in firing rates and gamma activity 50–200 ms after stimulus onset (supplemental Fig. 1, available at www.jneurosci.org as supplemental material). This period is called the “early stimulus period” (esp). To study the spike–LFP correlation under conditions far away from any transients, the analyses were also performed at longer latencies, called the “late stimulus period” (lsp; 250–500 ms after stimulus onset). The baseline period was chosen between −200 to −50 ms so that it was symmetrical around stimulus onset with the early stimulus period, although several other periods (chosen between the interval −350–0 ms and of duration 150–350 ms) were also used, and very similar results were obtained.
Frequency bands.
We used the gamma frequency range (sometimes referred to in this report as “low-gamma” for emphasis) between 40 and 60 Hz and high-gamma frequency range between 60 and 150 Hz. However, the choice of these ranges was not critical. For example, all analyses for the gamma band were repeated for 40–80 Hz, and for the high-gamma band were repeated for 80–150 Hz as well as for 60–200 Hz, and very similar results were obtained. In general, the lower frequency limit of a frequency band is more critical because the LFP power decreases with increasing frequency (a typical “1/f” falloff) (Buzsáki and Draguhn, 2004; Henrie and Shapley, 2005); thus, the total power in any frequency band is dominated by the power at lower frequencies within the band. For example, total power between 60 and 100 Hz is very similar to the power between 60 and 150 Hz, or even between 60 and 200 Hz. In this regard, the high-gamma frequency band used in this study is similar to the “gamma” frequency range in the notation used by Bauer et al. (2006), who used a range between 60 and 95 Hz. Although the gamma range extended beyond ∼200 Hz, the power at these higher frequencies was often too low to be distinguished from noise. The results are thus shown only up to 175 Hz.
Construction of STTFA and D.
The STTFA was computed directly from the time frequency plot obtained from the MP algorithm. It is analogous to the spike-triggered average, but instead of taking a one-dimensional (1D) segment from the signal centered on each spike and averaging over all segments, we took a two-dimensional segment from the time–frequency plot of the signal centered on the spike (50 ms either side of the spike) and took the average of those 2D segments. This is illustrated in supplemental Figure 2 (available at www.jneurosci.org as supplemental material). The STTFA derived from the baseline, early stimulus period and late stimulus period were denoted by STTFAbl, STTFAesp, and STTFAlsp, respectively (supplemental Fig. 2B, only the baseline and early stimulus periods are shown for clarity, available at www.jneurosci.org as supplemental material).
For the STA, positive and negative non-phase-locked components in the signal cancel each other in the average, giving only the part of the signal that is phase locked to the spikes. The STTFA, however, is derived from the time–frequency plot, which is always positive. To compensate for the time–frequency components that are not locked to the spikes, we randomized the spike times (supplemental Fig. 2A, magenta spike raster, available at www.jneurosci.org as supplemental material) and recomputed the STTFA of the randomized spikes (the randomized spike times were uniformly distributed in the analysis interval). This average yields the randomized STTFA (rSTTFA), which represents the average power present in the signal that is not phase locked to spikes.
To isolate the spike components, we subtracted the randomized STTFA from the STTFA, and defined the “normalized STTFA,” denoted by D: Here, x denotes the time period (baseline, early stimulus period, or late stimulus period) from which the spikes are taken for the computation of the STTFA and randomized STTFA. The statistics D(bl) and D(esp) are shown in supplemental Figure 2C (available at www.jneurosci.org as supplemental material). The statistics are noisy in this figure because they are computed from a single neuron.
This method was tested for consistency as well as for possible sources of bias, such as the randomization strategy to get the rSTTFA, stimulus phase locking and choice of interval size (supplemental Discussion 1, available at www.jneurosci.org as supplemental material). The results shown in the study could not be explained by any of these factors.
Results
The responses from single neurons and LFP signals were collected simultaneously from the secondary somatosensory cortex of three hemispheres of two awake behaving monkeys while vibratory stimuli at three different frequencies (50, 100, and 200 Hz) and four different intensities (relative amplitude 1:2:5:10) were presented. From a visual inspection of the raster plots, as well as comparison of the firing rates in the stimulus period versus baseline, the neurons were divided into three categories: excited (78 neurons), inhibited (21), and not-driven (42) (Fig. 1A) (for details, see Materials and Methods). All results shown are for a stimulus frequency of 50 Hz; the results for other stimulus frequencies were similar.
Time–frequency analysis
Time–frequency analysis of the LFP showed a consistent increase in power in the gamma and high-gamma frequency ranges as well as a decrease in power in the β range (16–24 Hz) after stimulus onset (Fig. 1B). The results shown in Figure 1B were derived by calculating the relative change in power (in decibels) in each time–frequency bin with respect to baseline power, which is defined as the mean power 200–50 ms before stimulus onset. The different rows in this figure show the time–frequency power changes at the four stimulus amplitudes, denoted by G1, G2, G5, and G10. The plots show mean changes in power for the three populations of excited, inhibited, and not-driven neurons; the plots for individual neurons were similar. As shown in Figure 1B, the increase in gamma power and suppression in β power was more prominent at high stimulus amplitudes. This modulation of power was observed even in inhibited or not-driven neurons, for which the firing rates decreased or did not change (middle and right columns).
Spike-triggered average
To study the correlation between spikes and LFP, we computed the STA of the LFP for the entire population of recorded neurons (Fig. 2A). The STA was computed from the spikes in the baseline period (gray line) or from spikes evoked in the early stimulus period (50–200 ms after stimulus onset; black line). As expected, the STA had a large negative component for all neuron types (Fig. 2A, denoted by 2), which is most likely caused by a decrease of Na+ concentration at the time of the spike in the extracellular space near the neuron. This negativity was followed by a positive deflection in the LFP (Fig 2A, denoted by 3), presumably caused by hyperpolarization of the neuron after the spike. The STA constructed from the spikes in the baseline period had a greater negative deflection, and appeared to be phase locked to the negative edge of a β oscillation (Fig. 2B). We also observed a small positive peak (Fig 2A, denoted by 1) preceding the negativity, possibly because of capacitive currents (Gold et al., 2006). When computed separately, the shapes of the STA were similar for the three neuron categories (see Fig. 4A, top).
Decomposition of STA using MP
Here, we illustrate the MP algorithm by decomposing the mean STA of the entire population during the baseline period (Fig. 2A, gray line). The MP algorithm iteratively decomposes a signal into a linear combination of functions taken from a large overcomplete dictionary, which best explain the shape of the signal (see Materials and Methods). We used a dictionary composed of Gabor, Fourier, and δ functions, referred to as “atoms” in this study. Figure 2B shows the first three atoms with center positions close to t = 0 (which is the time of occurrence of the spike as recorded by the data collection system), together with the measured STA and the estimated STA as reconstructed by adding these three atoms (black dotted trace). The three atoms are listed, in the legend, in the order of appearance in the iterations of the MP algorithm, i.e., in the order of decreasing energy. The first (highest-energy) atom is a signal in the β range (blue trace), suggesting that the spikes were indeed phase locked to β oscillations in the baseline period. The second atom (magenta trace) is a Gabor function with zero modulation frequency (ω = 0 in Eq. 2), i.e., a Gaussian function, which captured the sharp transient associated with the action potential. The third atom is a sinusoidal signal in the high-gamma frequency range (red trace).
This signal decomposition is critical for testing gamma phase coding, because the phase-coding mechanism necessarily requires an oscillatory signal in the gamma frequency range. Traditional methods such as STFT decompose any signal into a linear sum of windowed sinusoids, so that a decomposition based on STFT will necessarily generate a series of oscillatory components even if they are not actually present in the signal. These oscillatory components in this case are meaningless: the same signal could equally well be represented by a different set of basis functions that are not sinusoidal. It is mathematically correct to decompose a signal into a sum of sinusoids, but it is inappropriate to assign physiological significance to these individual sinusoidal components. For example, the sharp transient component, shown by the magenta trace, will instead be represented by a sum of oscillatory signals at several different frequencies if STFT is used for signal decomposition. MP decomposition allows the comparison of only the oscillatory signals (such as the red trace) in different frequency bands as a function of stimulus amplitude, as described in the next section.
Local components of the STA
We decomposed the STA of each neuron separately, obtaining 100 atoms per STA. In general, the center positions and frequencies of these atoms varied considerably for each neuron. Based on the center frequencies of the atoms, the STA was decomposed into several “local components.” Figure 3A shows the raw mean STAs (for the entire population of neurons) during baseline (dotted black line) and early stimulus period (four solid lines in different shades of gray for the four stimulus intensities). We reconstructed the signal from atoms with center frequencies between 10 and 40 Hz to isolate the “low-frequency component” (mainly β rhythm at ∼20 Hz). Figure 3B shows the mean low-frequency component taken by averaging over low-frequency components of individual neurons. As shown in Figure 3B, the spikes were indeed phase locked to the negative trough of the β rhythm, which decreased in magnitude with increasing stimulus intensity.
For most STAs, there was an atom similar to the magenta trace shown in Figure 2B, and some atoms with center frequencies >200 Hz. Supplemental Figure 3 (available at www.jneurosci.org as supplemental material) shows the histogram of the center positions (u as defined in Eq. 2) and scale [log2(s), s as defined in Eq. 2] of the highest energy atom at frequencies f = 0 and f > 200 Hz for each STA [where f = ω/(2π), ω as defined in Eq. 2]. A large proportion of atoms with frequency f = 0 were centered near the time of occurrence of the action potential and had very small scales (≤3). These atoms were similar to the magenta trace as shown in Figure 2B, which captured the negative deflection of the action potential. Similarly, the atoms with f > 200 Hz centered near the time of occurrence of the action potential captured the other two positive transients (peaks 1 and 3 in Fig. 2A). Together, these atoms represented the sharp transients typically associated with the action potential itself. We reconstructed the “sharp transient” component by adding all atoms with a center frequency of either zero or >200 Hz (Fig. 3C). Similarly, we added the atoms with center frequencies between 40 and 60 Hz (Fig. 3D) and between 60 and 150 Hz (Fig. 3E) to generate low-gamma and high-gamma components, respectively. The high-gamma component was much stronger than the low-gamma component (Fig. 3D,E, same scale used for comparison).
This decomposition is somewhat arbitrary; these local components need not originate from separate biophysical sources because MP decomposition does not isolate independent components. Nonetheless, this approach dissociates the “spike-like” sharp transient, which has energy spread over a broad frequency range (including the gamma and high-gamma ranges) from other oscillatory components in the low-gamma or high-gamma range while preserving the local structure of the signal. Importantly, we found that neither the sharp transient nor the gamma or high-gamma components varied with stimulus intensity. This is in contradiction to the phase-coding hypothesis, which predicts a rightward shift in the low-gamma or high-gamma component with increasing intensity. Similar results were obtained when a different range was used for low-gamma (40–80 Hz) or high-gamma frequencies (80–150 or 60–200 Hz).
We were concerned that our results could be biased by the dyadic structure of the dictionary used for MP decomposition (see Materials and Methods). To test whether the dyadic grid affected our results, we repeated the analysis after moving the STA signal relative to the grid. The details are described in supplemental Figure 4 (available at www.jneurosci.org as supplemental material) and the corresponding text. We show that the absence of a rightward shift of the STA components with increasing stimulus intensity cannot be attributed to a bias in the MP grid structure.
These results demonstrate that simple stimulus features such as intensity do not modify the different local structures of the LFP (Fig. 3B–E) associated with an action potential. However, from these plots the complete temporal and spectral characteristics of these local structures are not clear. Because these components have a narrow time support, their energy is spread over a broad frequency range. To observe all the components associated with a spike in both spectral and temporal domains, we computed the spike-triggered average of the time–frequency spectrum of the LFP (denoted by STTFA and described in detail in Materials and Methods). Here, instead of averaging over segments of the LFP centered on each spike, we averaged the two-dimensional time–frequency energy plot of the LFP around each spike.
Time–frequency analysis of the STA
The STTFA represents the average time–frequency spectrum of the LFP that is associated with a spike (see Materials and Methods). To isolate the spectral-temporal components specifically associated with a spike, we normalized the STTFA (denoted by D, defined in Eq. 4) using a randomization technique (see Materials and Methods) (supplemental Fig. 2, available at www.jneurosci.org as supplemental material).
Figure 4A shows the average STA of the three neuron types during the baseline period (top), as well as the corresponding normalized STTFA in the baseline period, i.e., D(bl) (bottom). We observed moderate low-gamma and strong high-gamma power near the origin, suggesting that spikes were tightly coupled to time-localized activity in gamma and particularly in high-gamma ranges.
Figure 4B shows the normalized STTFA during the early stimulus period, i.e., D(esp), for the four different stimulus amplitudes. We observed that the low- and high-gamma power associated with the spikes remained similar for the different stimulus amplitudes for all three neuron types. Consistent with earlier results, neither the power of these low- and high-gamma oscillations nor their position with respect to the action potential appeared to vary with changes in stimulus amplitude.
The method of construction of STTFA eliminates possible sources of bias because of the use of a dyadic dictionary. As described, STTFA is constructed by first computing the time–frequency plot of signal energy and then extracting 2D segments (centered on the spikes) from this plot. This is different from first taking the 1D segments (i.e., the STA) from the signal and then computing their time–frequency plots. The latter method may be biased because of biases in the dyadic dictionary, but in the former method MP is performed only once, and the 2D segments that are taken out are centered on spikes (which are not regularly spaced and do not coincide with the dyadic grid). Therefore, it is unlikely that these 2D segments have any systematic bias. Furthermore, in a previous study (Ray et al., 2003), we compared the time–frequency spectrum generated by the dyadic dictionary with a time–frequency spectrum generated by a different (stochastic) dictionary, and showed that although the atomic decompositions were different with the two dictionaries, the time–frequency spectra were very similar.
We performed detailed statistical tests by computing the time of peak negativity in the low-gamma (Fig. 3D) and high-gamma component (Fig. 3E) of individual STAs. The peak negativity time was defined as the time at which the signal reached its minimum value in the range between ± 10 ms (other ranges such as ±15 or ±20 ms were also used and similar results were obtained). Figure 5A shows the mean and SDs of the low-gamma peak negativity times for the three neuron types during baseline period and the early stimulus period. The peak negativity of the low-gamma component of the STAs during the stimulus period was not different from baseline for any of the four intensities (p > 0.05, t tests; n = 78, 21, and 42 for the three neuron types). The histogram in Figure 5B shows when these peak negativity times occurred in the neural population (the histograms constructed separately for the three neuron types were similar; hence, the data were pooled together). The histogram shows that the distribution of peak negativity times was fairly uniform (with peaks at the extremes of the analysis window), which is not surprising given the low magnitude of the low-gamma signal. Scatter plots of the peak negativity times at different stimulus intensities also did not show any systematic relationship between peak negativities and stimulus amplitude for any neuron type (t tests, p > 0.05) (data not shown).
The results were similar when the analysis was performed on the time–frequency STTFA plots shown in Figure 4. The time (in the interval ±10 ms relative to the spike) at which the power in the STTFA between 40 and 60 Hz was maximal is shown in Figure 5C. The histogram in Figure 5D shows when these peak power times occurred in the neural population. The peak time did not change with stimulus amplitude (p > 0.05, t tests; n = 78, 21, and 42 for the three neuron types). These results are at odds with the gamma phase-coding hypothesis that predicts well-defined peaks in the histograms of negativity times and power times (Fig. 5B,D), and a rightward shift in these peaks with an increase in stimulus amplitude.
Similar statistical tests performed on the high-gamma range (60–150 Hz) showed no systematic change in the peak negativity times (Fig. 6A) or peak power times (Fig. 6C) with stimulus intensity (p > 0.05, t tests; n = 78, 21, and 42 for the three neuron types). The histograms of peak negativity times (Fig. 6B) revealed a peak at 0 ms and a smaller peak at ∼8 ms, corresponding to the two negative peaks in Figure 3E (during the baseline period or at low stimulus intensities, the variance was higher because the STAs were constructed from fewer spikes and, hence, were more noisy). Similar results were obtained from the histogram of peak power times (Fig. 6D). Thus, spikes and high-gamma oscillatory signals showed a specific phase relationship, but there was no change in this phase relationship with stimulus intensity.
Similar results were obtained when the analysis was performed at the late stimulus period (250–500 ms after stimulus onset), in which the variation in firing rates was much lower (data not shown). In supplemental Discussion 2 (available at www.jneurosci.org as supplemental material), we discuss additional tests that were performed to test the phase-coding hypothesis. The results are consistent with the main findings of this study.
Discussion
We studied the neural correlates of gamma oscillations and their possible functional role in cortical processing. In particular, we tested a recently proposed hypothesis that excitatory input to cortical pyramidal cells is converted to a phase code with respect to LFP gamma oscillations, such that more excited cells fire earlier in the gamma cycle (Fries et al., 2007). We presented vibratory stimuli at different amplitudes to the finger pads of awake behaving monkeys, recorded the responses (spikes and LFPs) of neurons in SII cortex, and used a novel technique based on matching pursuits to study the spike–LFP correlation as a function of stimulus amplitude. We observed weak low-gamma (40–60 Hz) and strong high-gamma (60–150 Hz) oscillations coupled to spikes. However, neither the phase nor the power of these oscillations depended on stimulus amplitude. The results do not support the gamma phase-coding hypothesis, at least under our experimental conditions.
Our experimental conditions were well suited to test the gamma phase-coding hypothesis. The phase coding scheme cannot be operational at the earliest levels of sensory processing (Fries et al., 2007), but SII cortex is at an advanced stage in somatosensory processing (Hsiao et al., 2002). Furthermore, neurons in SII show increases in neuronal synchronization with selective attention (Steinmetz et al., 2000) and there is evidence of gamma oscillations in this area during attentional processing (Bauer et al., 2006). The stimulus amplitudes that we used varied by a factor of 10 and clearly increased the excitatory drive reaching SII (Fig. 1), without any saturation effects. Furthermore, the MP algorithm is capable of detecting both rhythmic components with high-frequency resolution as well as transients with high temporal resolution (down to ∼0.2 ms in this study).
Matching pursuit versus other time–frequency methods
Several methods have been used to study the correlation between spikes and LFP in spectral and temporal domains. In most of these methods, such as multitapering (Mitra and Pesaran, 1999; Jarvis and Mitra, 2001) and empirical mode decomposition (EMD) (Liang et al., 2005), the signal is decomposed into functions that are localized in frequency but extended in time (at least on the order of 100 ms). However, gamma oscillations, especially at higher frequencies, often appear for a short duration but may have frequency content in a wide range. Thus, methods based on short-time Fourier transforms or Hilbert transforms (used in EMD) fail to capture these transients with a single function, instead decomposing the signal into a set of functions that may not be physiological. For example, a sharp transient (similar to the magenta trace in Fig. 2B) will be represented by a series of “oscillations” at different frequencies. Although wavelet transforms eliminate some of these limitations by decomposing the signal into functions that have high frequency resolution at low frequencies and high temporal resolution at higher frequencies, they create a biased division of the time–frequency space. In matching pursuit, by choosing a large dictionary of Gabor atoms, we get fewer a priori limitations on decomposition and more free parameters than other methods and are able to detect local patterns in the signal with the best possible compromise between time and frequency resolution. Furthermore, the technique is better adapted to capture the complex phenomena studied here than methods with rigid decompositions of the time–frequency space because the division of the time–frequency space in MP depends on the structure of the LFP signal itself.
“Bleed-through” of a spike in the LFP
One important concern in spike–LFP interaction analysis is that low-frequency components of the action potential may bleed through in the LFP. One method of avoiding bleed-through is to record spikes and LFP from two different microelectrodes. However, this method is only approximate, because low-frequency components of an action potential are physiological entities, not recording artifact, and their contribution in the gamma range must be dealt with head on. The use of matching pursuit, especially testing the phase-coding hypothesis on oscillatory components, is crucial in addressing the issue. In supplemental Discussion 3 (available at www.jneurosci.org as supplemental material), we discuss these issues in detail and also perform cross-STTFA analysis from two microelectrodes. The results are consistent with the findings of this study.
Band-limited versus broad-band gamma
We divided the gamma frequency bands into low and high ranges only for completeness, because the upper limit of the traditional gamma band is ∼80 Hz. However, Figure 1B shows that the increase in power after stimulus onset is fairly uniform in a broad frequency range, starting at around 40 Hz and up to 200 Hz and beyond. This broad-band increase in power, with a typical peak in the high-gamma range, has been observed in a number of invasive EEG studies on epileptic patients in a variety of functional domains, including sensorimotor (Crone et al., 1998; Pfurtscheller et al., 2003, Miller et al., 2007), auditory (Crone et al., 2001a; Edwards et al., 2005), visual (Lachaux et al., 2005; Tanji et al., 2005), language (Crone et al., 2001b; Sinai et al., 2005), attention (Tallon-Baudry et al., 2005; Ray et al., 2008), and working memory (Mainy et al., 2007). A similar broad-band increase in power with an increase in stimulus contrast was also reported in a single-unit neurophysiology study in primary visual cortex by Henrie and Shapley (2005). These authors reported an increase in power in a broad frequency band between 25 Hz to well beyond 200 Hz, although the maximum increase in power was at 25–90 Hz.
There are several reasons why high-gamma oscillations above ∼80 Hz have not been extensively reported in LFP studies. First, high frequencies carry considerably less power than lower frequencies (Buzsáki and Draguhn, 2004; Henrie and Shapley, 2005), with a typical power falloff of 1/f. Consequently, high-quality recordings as well as sophisticated signal processing tools are required to distinguish oscillations in high-gamma frequencies from noise. Furthermore, low-gamma oscillations appear to last for a longer duration than high-gamma oscillations (supplemental Fig. 1, available at www.jneurosci.org as supplemental material). Traditional time–frequency analysis methods such as STFT have long temporal integration windows (typically 250 ms or longer), and thus yield the average power over a long temporal window. The power of the highly transient high-gamma oscillations is consequently much lower than the more sustained low-gamma oscillations when averaged over a long period.
Our results suggest that the power in the high-gamma range is tightly coupled to the occurrence of action potentials, and therefore could be a neural correlate of population firing rate. These results are in agreement with a recent study by Liu and Newsome (2006) that showed that LFP power in the gamma range and above were tightly coupled to multiunit activity near the microelectrode tip. In their study, the maximum modulation was in the high-gamma range and very similar to the range used by us (60–150 Hz). Indeed, if each action potential is associated with a burst of energy in the high-gamma frequency range (Fig. 4), the overall LFP high-gamma power would be tightly correlated with the total number of action potentials generated by the neurons near the microelectrode tip, i.e., the multiunit firing rate.
Functional role of gamma oscillations
Our results do not support the hypothesis that stimulus strength is encoded in cortex by the position of action potentials relative to the phase of population oscillations in the gamma (or high-gamma) range. Our results do not, however, contradict the view that gamma oscillations may play an important role in neural processing. Strong experimental evidence has been obtained for such a role in various paradigms, including selective attention (Fries et al., 2001), corticospinal interaction (Schoffelen et al., 2005), sensorimotor integration (Roelfsema et al., 1997; Womelsdorf et al., 2006), movement planning/selection (Grammont and Riehle, 2003; Scherberger et al., 2005), working memory (Pesaran et al., 2002), and associative learning (Miltner et al., 1999). It is possible that gamma oscillations have functional importance in the implementation of such higher order processes, but are not used to code lower level stimulus features such as stimulus amplitude. Because our experimental task did not manipulate higher processes such as attention, it may not be surprising that no modulation in spike–LFP correlation in the gamma range was observed. Furthermore, gamma oscillations may play a role in neural processing through other mechanisms, such as changing the response gain of neurons (Azouz and Gray, 2000, 2003; Shu et al., 2003; Hasenstaub et al., 2005) or controlling the timing and probability of action potential generation in pyramidal cells (Hasenstaub et al., 2005). Thus, gamma oscillations may well have a functional role in neural processing although the phase-coding hypothesis could not be corroborated in our experiment.
The techniques developed in this study provide a framework for studying the fine temporal interaction of gamma oscillations with neuronal spiking. We show that simple stimulus features, such as amplitude, do not modulate these spike–LFP interactions. A study of these interactions in the spectral and temporal domain under experimental conditions that modulate higher functions such as attention or memory may lead to a deeper understanding of the functional role of gamma oscillations.
Footnotes
-
This work was supported by National Institutes of Health Grants R01-NS18787 (S.S.H.), R01-NS40596 (N.E.C., E.N., P.J.F.), 5R01EY016281-02 (E.N.), R01-NS34086 (S.S.H.), and NS43188-01A1 (E.N., S.S.H.). We thank J. H. R. Maunsell and V. Stuphorn for their helpful comments, P. J. Fitzgerald for his extensive help in data collection, and J. Killebrew and J. F. Dammann for technical support.
- Correspondence should be addressed to Supratim Ray, Department of Neurobiology, Harvard Medical School and the Howard Hughes Medical Institute, 201 Goldenson, 220 Longwood Avenue, Boston, MA 02115. Supratim_Ray{at}hms.harvard.edu