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) (, see Methods for details). All results shown are for a stimulus frequency of 50 Hz; the results for other stimulus frequencies were similar.
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 beta range (16−24 Hz) following stimulus onset (). The results shown in were derived by calculating the relative change in power (in dB) in each time-frequency bin with respect to baseline power, which is defined as the mean power 200 to 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 , the increase in gamma power and suppression in beta 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).
To study the correlation between spikes and LFP, we computed the spike-triggered average (STA) of the LFP for the entire population of recorded neurons (). 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 (denoted by ‘2’ in ), which is most likely caused by a decrease of Na+
concentration at the time of the spike in the extra-cellular space near the neuron. This negativity was followed by a positive deflection in the LFP (denoted by ‘3’), presumably caused by hyper-polarization of the neuron following 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 beta oscillation (). We also observed a small positive peak (denoted by ‘1’) preceding the negativity, possibly due to capacitive currents (Gold et al., 2006
). When computed separately, the shapes of the STA were similar for the three neuron categories (, upper panel).
Figure 2 Spike-triggered average (STA) analysis. A) The mean STA of all recorded neurons. The STA is constructed from spikes taken from the baseline period (shown in gray) or early stimulus period (50 − 200 ms, black). B) Atomic decomposition of the STA (more ...)
Decomposition of STA using MP
Here we illustrate the MP algorithm by decomposing the mean STA of the entire population during the baseline period (gray line in ). 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 Methods). We used a dictionary composed of Gabor, Fourier and delta functions, referred to as ‘atoms’ in this paper. shows the first 3 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 beta range (blue trace), suggesting that the spikes were indeed phase-locked to beta oscillations in the baseline period. The second atom (magenta trace) is a Gabor function with zero modulation frequency (ω = 0 in equation 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’. 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 beta rhythm at ~20 Hz). shows the mean ‘low frequency component’ taken by averaging over low frequency components of individual neurons. As shown in , the spikes were indeed phase-locked to the negative trough of the beta rhythm, which decreased in magnitude with increasing stimulus intensity.
For most STAs, there was an atom similar to the magenta trace shown in , and some atoms with center frequencies greater than 200 Hz. Supplementary Figure 3
shows the histogram of the center positions (‘u’ as defined in equation 2
) and scale (log2
(s), s as defined in equation 2
) of the highest energy atom at frequencies f=0 and f>200 Hz for each STA (where f=ω/(2π), ω as defined in equation 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 , 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 ). Taken 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 greater than 200 Hz (). Similarly, we added the atoms with center frequencies between 40−60 Hz () and between 60−150 Hz () to generate ‘low-gamma’ and ‘high-gamma’ components, respectively. The high-gamma component was much stronger than the low-gamma component (note that are on the same scale for comparison).
This decomposition is somewhat arbitrary – these local components need not originate from separate biophysical sources, since 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 range) 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 Hz 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 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 Supplementary Figure 4
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 (as shown in ) associated with an action potential. However, from these plots the complete temporal and spectral characteristics of these local structures are not clear. Since 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 the Methods section). Here instead of averaging over segments of the LFP centered on each spike, we averaged the 2-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 Methods). To isolate the spectral-temporal components specifically associated with a spike, we normalized the STTFA (denoted by D, defined in equation 4
) using a randomization technique (see Methods and Supplementary Figure 2
shows the average STA of the three neuron types during the baseline period (upper panel), as well as the corresponding normalized STTFA in the baseline period i.e, D(bl) (lower panel). 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.
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 due to the use of a dyadic dictionary. As described, STTFA is constructed by first computing the time-frequency plot of signal energy and then extracting 2-D segments (centered on the spikes) from this plot. This is different than first taking the 1-D segments (i.e., the STA) from the signal and then computing their time-frequency plots. The latter method may be biased due to biases in the dyadic dictionary, but in the former method MP is performed only once, and the 2-D segments that are taken out are centered on spikes (which are not regularly spaced and don't coincide with the dyadic grid). Therefore, it is unlikely that these 2-D segments have any systematic bias. Further, 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 even though 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 () and high-gamma component () 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). shows the mean and standard deviations 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 shows when these peak negativity times occurred in the neural population (the histograms constructed separately for the three neuron types were similar and 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).
Figure 5 Changes in low-gamma phase with stimulus amplitude. A) Time of the peak negativity in the low-gamma component (as described in ) for the three different neuron types at baseline and four different stimulus intensities. B) Histogram of times at (more ...)
The results were similar when the analysis was performed on the time-frequency STTFA plots shown in . The time (in the interval ±10 ms relative to the spike) at which the power in the STTFA between 40−60 Hz was maximal is shown in . The histogram in 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 (), 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 () or peak power times () with stimulus intensity (p>0.05, t-tests, n=78, 21 and 42 for the three neuron types). The histograms of peak negativity times () revealed a peak at 0 ms and a smaller peak at about 8 ms, corresponding to the two negative peaks in (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 (). Thus, spikes and high-gamma oscillatory signals showed a specific phase-relationship, but there was no change in this phase relationship with stimulus intensity.
Changes in high-gamma phase with stimulus amplitude. Panel configurations are similar to .
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 Supplementary Discussion 2
, we discuss additional tests that were performed to test the phase-coding hypothesis. The results are consistent with the main findings of this report.