| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
an1,2,3
1,2,3
1,21Frankfurt Institute for Advanced Studies; 2Max Planck Institute for Brain Research, Frankfurt am Main, Germany; and 3Center for Cognitive and Neuronal Studies, Cluj-Napoca, Romania
Submitted 9 July 2007; accepted in final form 22 December 2007
 |
ABSTRACT |
|---|
|
|
|
INTRODUCTION |
|---|
|
; Steriade 2006; Steriade et al. 1993) and awake states (Canolty et al. 2006; Fries et al. 2001; Pesaran et al. 2002), and they can be organized into different frequency bands (Buzsáki and Draghun 2004). Oscillations in each band are likely to have a different underlying mechanism (Buzsáki 2006) and probably subserve a different function. Slow (0.5–1 Hz) and delta (1–4 Hz) oscillations typically engage large parts of the thalamocortical system (Steriade et al. 1993), especially during non-rapid eye movement (NREM) sleep (Steriade 2006). These oscillatory rhythms are known to correlate with hyperpolarized cellular states (Timofeev et al. 2000) during which neurons transiently engage in network bursting events (Murean and Savin 2007). Theta-band activity (4–8 Hz) is expressed in the hippocampus during locomotion or REM sleep (Maurer and McNaughton 2007). Alpha-band activity (8–12 Hz) is associated most commonly with the resting state but is also frequently expressed during behavior (Gunji et al. 2007) and various cognitive processes (Klimesch 1999; Klimesch et al. 2007), including meditation (Murata et al. 2004). Fast oscillations, with frequencies in the beta (12–30 Hz) and gamma (30–80 Hz) ranges occur during awake states and REM sleep (Steriade 2006) and, in contrast to the slow oscillations that have a more global character (Steriade et al. 1993), these fast oscillatory states are more often localized and shorter-lived (Buzsáki 2006). Importantly, the fast oscillations in the beta and gamma ranges have been shown to correlate with perception (Meador et al. 2002; Melloni et al. 2007; Rodriguez et al. 1999; Singer 1999; Tallon-Baudry et al. 1997; Vidal et al. 2006), behavioral performance (Tallon-Baudry et al. 2004), working memory (Pesaran et al. 2002), visual attention (Fries et al. 2001), motor preparation and execution (Schoffelen et al. 2005), and other cognitive processes (Lutz et al. 2004). It has been suggested that the fast, stimulus-selective gamma oscillations contribute to perceptual organization (Singer 1999), coordination of neuronal activity (Fries et al. 2007; Martinovic et al. 2007; Uhlhaas and Singer 2006), and also to the induction and maintenance of stimulus-dependent synchronization (Engel et al. 1990; Samonds and Bonds 2005). Such overwhelming evidence—that neuronal oscillations evolve on various timescales, are expressed in virtually all cortical systems, and are closely related to cognitive and executive functions—calls for robust methods of investigating their expression in neuronal activity.
Investigation of neuronal oscillations at the single-neuron level is crucial for understanding how various rhythms emerge from the interactions within large cell populations (Baker et al. 2003; Buzsáki et al. 2004). Nonetheless, estimating the degree to which a neuron follows an oscillatory rhythm is not an easy task, for several reasons. First, even when the neuronal population oscillates, as revealed by local-field potential (LFP) activity, single-cell activity frequently fails to reflect the underlying rhythm because of undersampling (Baker et al. 2003; Singer 1999). This problem can be confronted by computing spike-field coherence measures (Fries et al. 2001; Pesaran et al. 2002; Zeitler et al. 2006). Second, the firing rates of cortical neurons are relatively low (Graham and Field 2006) and thus when these neurons follow fast oscillatory rhythms they often skip oscillatory cycles. This usually leads to a reduced, often too low, number of spikes for the estimation of the cell's oscillatory behavior. Therefore, most methods estimating oscillation strength work well only for neuronal activity with high firing rates. Third, fast oscillations are often transient, whether task dependent or spontaneously occurring (Buzsáki 2006; Fries et al. 2001; Melloni et al. 2007; Rodriguez et al. 1999; Steriade 2006). Their transient nature can be attributed to functional processes or to slower cortical state changes. Because of this nonstationary expression of oscillations, methods that quantify oscillatory behavior should be applicable to short stretches of activity. If too long traces of signals (too many or too long trials) are analyzed the oscillatory episodes can be averaged out. As a consequence, and because of the often low firing rates of cortical neurons, analysis methods need to be applied that yield reliable results despite the small number of spikes. Finally, neurons can engage in multiple rhythms such that they follow fast oscillations riding on top of slow ones (Isomura et al. 2006; Steriade 2006). This can impose limitations for methods that determine the oscillation frequency and its amplitude by fitting an underlying model function, such as a Gabor function (König 1994).
Several methods have been developed for the estimation of oscillatory behavior from neuronal spike trains. Two main categories can be defined: methods based on spectral analysis of spike trains and methods based on computing auto- and cross-correlation histograms. Among the methods in the first category, probably the simplest one is the periodogram, which directly applies a fast Fourier transform (FFT) on spike trains that have been converted to sequences of zeros (no spike) and ones (spike) (Bair et al. 1994). This method suffers from bias and variance problems and is thus highly inaccurate (Jarvis and Mitra 2001). A more efficient solution is achieved by applying the FFT to spike trains that have been previously smoothed (e.g., by using a Gaussian smoothing kernel) and then windowed using specialized windows, like the Blackman window, for example (Harris 1978). Another technique to compute the spectrum of spike trains, described in Jarvis and Mitra (2001), relies on the fact that the spectrum of a so-called spike-count function (which integrates the spike counts up to a moment in time) leads to the spike spectrum (Pesaran et al. 2002). In addition, the method uses a family of smoothing functions, or "tapers," for windowing (Percival and Walden 1993). Thus it belongs to the category of so-called multitaper techniques, which apply this special type of windowing before computing the spectrum. These methods are well understood analytically (Jarvis and Mitra 2001) and have found a number of successful applications (Mitra and Pesaran 1999; Pesaran et al. 2002). Spectral techniques are relatively robust but they also have many implicit assumptions, such as stationarity, for example. Unlike methods based on correlation analysis, they are restricted to frequency and phase coherence estimation and cannot provide information on such properties as refractoriness and burstiness, for example. Importantly, methods that estimate the spectrum from spike trains do not directly quantify the strength of oscillations. For that purpose, one could, for example, use the ratio between the spectral magnitude of a particular frequency of oscillation and the average magnitude of the whole spectrum. As we subsequently see (see Comparison to other methods), such a strategy is unfortunately biased by the firing rates of neurons. In general, it is difficult to compute such a measure that is independent of firing rate, directly from the spectrum of the spike train.
The second, and perhaps most common, category of methods for characterizing the oscillatory behavior of neurons is based on auto- and cross-correlation histograms. One commonly used method fits a generalized, eight-parameter Gabor function, to the correlation histogram by an iterative Marquardt–Levenberg algorithm. The fitted parameters of the Gabor function are used to determine not only the strength of synchronization between neurons but also the frequency and strength of oscillation (König 1994). The strength of oscillation is estimated as the ratio between the size of the first satellite peak and the offset (Fig. 1 A; König 1994). The method performs well on correlograms that exhibit gamma oscillations (Brecht et al. 1999) and whose shapes match well the Gabor function, but fails to converge when frequently the Gabor assumption is violated. Moreover, it has been suggested that the presence and strength of oscillations are often overestimated by methods relying on the first satellite peak (Samonds and Bonds 2005). Overestimation usually occurs when correlograms exhibit a dip, near the central peak (due to the burst refractory period) and a satellite peak (Fig. 1A), mimicking oscillatory activity. Such cases can frequently appear, even in the total absence of oscillations, when a central peak, with refractory dips, is superimposed on a slow, correlated rate increase (Fig. 1, A and B; Brody 1999). To overcome this problem, Samonds and Bonds (2005) used normalized autocorrelation histograms (ACHs) and suggested that one should rely on the difference between the second satellite peak and second valley to quantify the strength of oscillations (Fig. 1C). They also estimated the (gamma) oscillation frequency by identifying a peak in the FFT-computed spectrum of the ACH, between 20 and 60 Hz. This method requires, in general, a relatively large amount of data (Samonds and Bonds used responses to 516 stimulus presentations, each with a duration of 2 s) and is limited when smaller amounts of data are analyzed. If only 20 stimulus presentations are given, as is common, the ACH is often noisy and the reliable estimation of the second satellite peak and second valley is no longer feasible. Moreover, oscillation frequencies are often not stationary in time such that recording a large number of trials might actually blur the oscillatory modulation of the ACH. Also, as we shall subsequently see, the difference between the secondary peak and secondary valley critically depends on the firing rate of the neuron, such that the measure does not properly quantify the oscillation strength. Another important problem is the presence of the large center peak in the ACH. Due to its narrow width and frequently associated refractory dips, it introduces a strong bias toward detecting fast oscillation frequencies in the FFT-computed spectrum of the ACH (see Fig. 3, A and B in Computation of the oscillation score). For a thorough discussion on the strengths and weaknesses of these methods, see Comparison to other methods in RESULTS.
|
|
|
|
METHODS |
|---|
|
The experiments were performed in the visual cortex of lightly anesthetized, paralyzed cats (n = 4). Anesthesia was induced with ketamine [Ketanest, Parke-Davis, 10 mg·kg–1, administered intramuscularly (im)] and xylazine (Rompun, Bayer, 2 mg·kg–1, im) and maintained with a mixture of 70% N2O and 30% O2 supplemented with halothane (0.5–1.0%). After tracheotomy, the animals were placed in a stereotactic frame. A craniotomy was performed and the skull was cemented to a metal rod. After completion of all surgical procedures, the ear and eye bars were removed, and the halothane level was reduced to 0.4–0.6%. After ensuring that the level of anesthesia was stable and sufficiently deep to prevent any vegetative reactions to somatic stimulation, the animals were paralyzed with pancuronium bromide (Pancuronium, Organon, 0.15 mg·kg–1·h–1). Glucose and electrolytes were supplemented intravenously and through a gastric catheter. The end-tidal CO2 and rectal temperature were kept in the range of 3–4% and 37–38°C, respectively.
Visual stimuli were presented binocularly on a 21-in. computer screen (Hitachi CM813ET) with 100-Hz refresh rate. To obtain binocular fusion, the optical axes of the two eyes were first determined by mapping the borders of the respective receptive fields and then aligned on the computer screen with adjustable prisms placed in front of one eye. The software for visual stimulation was a combination of custom-made programs and a stimulation tool, ActiveSTIM (http://www.ActiveSTIM.com).
The investigated neuronal activity was acquired in response to a variety of visual stimuli: 1) Sinusoidal gratings moving in 12 directions in steps of 30°. 2) Center–surround stimuli, with a small sinusoidal grating placed in the center and a larger one in the surround. The orientations of both gratings and the sizes of the surround gratings were varied, resulting in a total of 14 stimulation conditions. 3) Moving bars and gratings with rectangular changes in luminance: including single bars, two superimposed bars moving in different directions, single gratings, and two overlapping gratings (plaids). This yielded a total of eight stimulation conditions. 4) Natural movies recorded for 28 s, containing indoor and outdoor moving scenes. Three different movies with different image statistics (slow, fast moving, dark, light, etc.) were presented. In each experiment, responses were recorded for 20 stimulus presentations. Different stimulus conditions were always presented in a randomized order.
Data were recorded from area 17 of four adult cats by inserting multiple silicon-based multielectrode probes (16 channels per electrode) supplied by the Center for Neural Communication Technology at the University of Michigan (Michigan probes). Each probe consisted of four 3-mm-long shanks that were separated by 200 µm and contained four electrode contacts each (area 1,250 µm2, impedance 0.3–0.5 M
at 1,000 Hz, intercontact distance 200 µm). Signals were amplified x10,000 and filtered between 500 Hz and 3.5 kHz and between 1 and 100 Hz for extracting multiunit (MU) activity and local-field potentials (LFPs), respectively. The waveforms of detected spikes were recorded for a duration of 1.2 ms, which allowed the later application of off line spike-sorting techniques to extract single units (SUs). In the analysis we considered only data recorded between the onset and the offset of stimulation. Unless specified otherwise, the reported analyses were made on SUs.
Computation of the oscillation score
The proposed method uses the following five steps to estimate the oscillation strength from an ACH: 1) the ACH is computed, 2) then it is smoothed, 3) the central peak is removed, 4) an FFT is applied to compute the spectrum of the peakless ACH, and 5) the oscillation score is calculated. The analysis is always made for a chosen frequency band, whereby the two input variables, fmin and fmax, represent, respectively, the low and high border of the frequency band of interest.
STEP 1: COMPUTATION OF THE AUTOCORRELATION HISTOGRAM (ACH).
The first step of the analysis involves computation of the ACH. The parameter that has to be determined is the time window of the ACH (the maximum lag or maximum time shift; e.g., 80 ms in Biederlack et al. 2006). Let b be the bin size used to compute the ACH (10–3 s in our study) and fc the frequency of the correlogram, given by fc = 1/b (1,000 Hz in this study; note that fc can be different from the sampling frequency of the spikes, depending on how one chooses to bin the ACH). Let w be the size, in bins, of one flank of the ACH and W the length, in bins, of the total ACH (given by W = 2w; the full symmetric ACH including the center bin at lag 0 has 2w + 1 bins). For computational reasons (see following text), W has to be a number power of 2, so the last bin of the ACH (the +w bin) is always left out of the analysis (Fig. 2 A). The choice of w is made according to the following three criteria: 1) w x b should be sufficiently large to fit at least a few cycles of the slowest frequency of interest, fmin. We required at least three cycles on each flank of the ACH (three-peaks rule), which, e.g., for fmin = 20 Hz (50-ms period) and b = 1 ms required w = 150 bins. 2) W should be sufficiently long to allow for a good resolution of the frequency spectrum obtained by the FFT. We always required 
2 Hz per bin of the magnitude spectrum, which equaled to at least W = 500 bins for ACHs with b = 1 ms. 3) To enable the computationally efficient application of the FFT, W should be a power of 2. The three criteria for the selection of w can be expressed as
 | (1) |
represents the rounding function to the nearest smaller integer (or floor function).
|
STEP 2: SMOOTHING.
When computing ACHs on spiking data with single units and on a low number of trials, the autocorrelation histogram is often noisy (because of the low number of spikes used to estimate it). Since we seek to minimize the effect of noise and approximate the real frequency spectrum as well as possible, an option is to apply a Gaussian smoothing with a fast kernel, to increase the signal-to-noise ratio. This would interfere only with very high frequencies (usually noise) and leave the lower, biologically relevant frequencies intact. The classical Gaussian kernel G(t) is given by
 | (2) |
is the SD of the Gaussian kernel.
An example ACH smoothed by such a Gaussian kernel is shown in Fig. 2B. The relevant parameter that has to be determined is fast (for the fast smoothing kernel). Since the smoothing is essentially a low-pass filter, it is important to choose fast such that the highest frequencies of interest are not attenuated below an acceptable threshold (e.g., –3 dB). For a correlogram frequency fc = 1 kHz, choosing fast = 2 ms leads to attenuations stronger than –3 dB for frequencies >67 Hz, whereas for fast = 1 ms the same effect is achieved for frequencies >134 Hz (Fig. 2C). In our method, we impose on fast an upper bound of 2 ms (to allow for the later estimation of noise in the ACH; see Estimation of the oscillation score's confidence). On the other hand, fast is allowed to take smaller values, depending on fmax. For better precision of the smoothing, the cutoff frequency (attenuated by the smoothing with –3 dB) is chosen at 1.5 x fmax. Thus fast is given by (in bins)
 | (3) |
fast = 0.893 bin. Note that the factor of 1.5 for computing the cutoff frequency has been empirically chosen.
STEP 3: REMOVING THE CENTRAL PEAK OF THE ACH.
The most important step of the method is the removal of the central peak from the ACH. The point of the maximum height of the central peak (at lag 0) indicates solely the firing rate of the neuron, which is uninteresting information for the present analyses. When computed for a spike train of a single neuron, the width of the central peak can reflect either the degree to which the neuron exhibits bursting behavior, in which case the peak is narrow corresponding to the average burst length (DeBusk et al. 1997), or it can reflect slow changes in firing rate (Brody 1999) and/or slow oscillatory activity, in which case the peak will be much wider. For cells that often produce fast bursts of action potentials, the central peak is frequently flanked by deep troughs, the duration of which indicates the refractory period for the appearance of repeated bursting events (Gray and McCormick 1996). These troughs can also be produced by inhibitory postsynaptic potentials during an ongoing oscillation, and the distinction from refractoriness is often impossible. Importantly, even in the absence of oscillations, the troughs can be followed by "rebound" satellite peaks (Fig. 1A) that can produce a false impression of oscillatory activity (Samonds and Bonds 2005). This can be a problem for any analysis method that is not designed to distinguish between different types of satellite peaks. See Fig. 1, A and B for an illustration of how such a misleading ACH can be produced by a combination of slow changes in rate responses and the troughs of burst refractoriness. Here, a method relying on the estimation of satellite peaks, such as the one described in König (1994), might falsely report the presence of strong oscillations (Fig. 1A).
Since our method relies on the identification of the frequency with the highest magnitude in the band of interest (see step 5), the central peak of the ACH can heavily interfere with the correct estimation. The central peak introduces a very high power that contaminates the whole spectrum (Fig. 3, A and B). Moreover, in conjunction with its side troughs (that frequently occur due to burst refractoriness) it can introduce peaks in the frequency spectrum of the ACH (Fig. 3B), falsely indicating the presence of oscillations. This could impair methods that quantify the frequency of oscillations relying on the frequency spectrum of the ACH, such as the one developed by Samonds and Bonds (2005), because it induces a false estimate for the oscillation frequency.
Removal of the central peak of the ACH is not a straightforward procedure. Ideally, when oscillations are present in the ACH, the central peak should be cut at the bottom of the side troughs, which preserves the satellite peaks that reflect oscillatory activity (Fig. 3C, top). In contrast, when oscillations are not present, but instead the cell displays deep refractory troughs, false peaks in the frequency spectrum should be avoided and thus the central peak should be removed together with the troughs; i.e., the troughs should be "filled up" (Fig. 3C, middle). Unfortunately, there is no exact method for deciding whether the satellite peaks represent oscillations. It is important to mention that, even for strong oscillations, if they are transient and nonstationary, the ACH will not reflect a clear oscillatory behavior, such that only one satellite peak might be expressed. The best approach, then, would be to consider a method that can cope with a low number of experimental trials, avoiding the problem of nonstationary oscillation frequencies and that estimates the oscillation strength and frequency by considering several cycles that are expressed in the ACH. The oscillation score aims to do exactly that.
Removal of the central peak, however, can cause other problems. A low-frequency component could be artificially introduced (Fig. 3C, bottom) when the side troughs are not filled up, either because satellite peaks needed to be preserved for the oscillatory case or because the cut was too small for the nonoscillatory case. This can then lead to a false indication of low-frequency oscillations in the frequency spectrum of the ACH; thus it is important to treat this case properly by selectively looking only into frequency bands higher than the artificially introduced frequency.
A reasonable solution for removal of the central peak is achieved if we consider that the method investigates only a given frequency band at a time (see step 5). For example, if the artificially introduced low frequency (due to removal of the central peak) is lower than the lowest frequency of interest fmin, then the algorithm would not detect false slow oscillations (because in step 5, it looks only in the band fmin – fmax). The problem of removing the central peak can be formulated in terms of finding the bins of the ACH that represent the limits between which the cut should be made. To detect the cutting limits for the central peak, we first apply a slow Gaussian smoothing on the original ACH, similar to the one in step 2, but with a much larger SD slow, taking into consideration fmin
 | (4) |
For example, for a correlogram frequency of 1 kHz, and a frequency band in the beta / gamma range (fmin = 20 Hz, fmax = 40 Hz), the ACH will be smoothed with a Gaussian having slow = 8.93 ms. This would create a smoothed central peak (Fig. 3D) with an effective width of 6 x slow (
50 ms). Note that the constants in Eq. 4 were determined empirically, to achieve a reasonable compromise of cutting efficiency for ACHs with and without oscillations, in the gamma band. For lower-frequency bands (e.g., alpha), the decision on the point for cutting the central peak is less critical because there is less confound between the burst refractory period (<30 ms) and the slower oscillation periods (>30–40 ms).
After smoothing the ACH with the slow Gaussian kernel, we obtain a "slow-smoothed" ACH (Fig. 3D, thick dotted black line). Then, we need to detect the limits of the smoothed central peak. This can be achieved accurately by applying an iterative search procedure that identifies the curvature of the smoothed signal (Fig. 3D, inset). It starts at lag 0 (top of the central peak) and searches to the left (negative time lags because the ACH is symmetric) a point on the curve where the slope of the local tangent is lower than a given threshold. The local slope is given by
 | (5) |
is an infinitesimally small quantity, and scaling is used to adjust the scale of the ACH on the y-axis in units of x-axis: scaling = W/Sslow(0) (the adjustment eliminates the dependence of the slope on the firing rate of the neuron).
The formula in Eq. 5 simply represents the first-order derivative of the curve at point t. For detecting the limits of the central peak, we imposed a limit of 10° on the local slope of the slow-smoothed ACH
 | (6) |
10° (Fig. 3D, inset). The slow-smoothed ACH is used only to detect the limits for the cutting of the central peak. After detecting the limits (ACH bins) tleft and tright of the central slow-smoothed peak (Fig. 3D, inset), we go back to the original, fast-smoothed ACH (that was computed in step 2). In the latter, the ACH value at tleft is copied to all other time lags spanning the identified range (tleft...tright), thus overwriting the real central peak (Fig. 3E). Depending on the relative size of the side troughs and first satellite peak, and on fmin, this procedure can also fill up the side troughs (see examples in Fig. 3, G and H).
The advantage of using the slow-smoothed ACH (Sslow) for detecting the central peak is that, after the cutting procedure, it avoids introducing low-frequency components in the spectrum within the band of interest. Due to the time constant of the slow smoothing, the introduced frequency is smaller than fmin (Fig. 3F). Moreover, when the ACH displays strong oscillatory components, the slope of the slow Sslow will reach the threshold of 10° closer to lag 0 (due to the deep first troughs and high satellite peaks), thus preserving the legitimate satellite peaks (Fig. 3, D and E). When the ACH displays only troughs due to refractoriness, with relatively weak satellite peaks, Sslow will reach the slope threshold later (further away from the central peak), leading to an efficient removal of the central peak (Fig. 3I), and frequently of the refractory troughs as well (Fig. 3, G and H).
It is important to mention that, in addition to slow components, removal of the central peak (a highly nonlinear operation) might also introduce some nonlinear distortions in the spectrum. However, our experience showed that these distortions are negligible and do not affect the reliable estimations of the oscillation frequency and strength.
STEP 4: APPLYING THE FFT.
The frequency spectrum is computed on the fast-smoothed ACH (obtained in step 2) with the central peak (and sometimes also side troughs) removed in step 3. We consider the smooth ACH in the interval [–w...+w – 1] (a buffer having a size that is a power of 2; see Fig. 2A) and we apply an FFT. The resulting buffer is of size w and contains the magnitudes of the frequencies from 0 to fc /2 Hz. When applying the FFT, it is very important to minimize frequency leakage (Oppenheim and Schafer 1999) to increase the precision of the frequency estimate. We recommend applying a Blackman window (Harris 1978) on the smoothed ACH before computing the FFT, to minimize leakage and border effects. Because the multitaper method is just another windowing technique (Percival and Walden 1993), it can also be applied instead of the Blackman window.
It is worth noting that because of the smoothing (that attenuates very high frequencies) and because of the distribution of biological frequencies in the ACH, very high frequencies, >200 Hz, have extremely low magnitudes (Fig. 3B).
STEP 5: ESTIMATION OF THE OSCILLATION SCORE. We defined the oscillation score as the ratio between the peak magnitude in the frequency band of interest and the mean magnitude of the spectrum. Intuitively, the score can be understood as the strength of the oscillation frequency relative to all the frequencies in the spectrum. The computation of the oscillation score has three components. First, the highest magnitude in the band of interest, Mpeak, is identified from the power spectrum (Fig. 3F). The estimated oscillation frequency fosc is the frequency with magnitude Mpeak. It is important to note that the method searches for the peak magnitude only within the band of interest fmin – fmax. Recall that this is an essential aspect because removal of the central peak in step 3 can artificially introduce frequencies lower than fmin (Fig. 3, E and F). By looking only to frequencies larger than fmin, the method is safe from detecting spurious, artificial magnitude peaks.
Next, the average magnitude Mavg of the spectrum is computed by integrating the whole frequency spectrum and taking its average
 | (7) |
The average magnitude of the spectrum is reduced by removal of the central peak from the ACH. This is highly beneficial since the legitimate oscillation peaks in the spectrum will have a higher ratio to the baseline.
The oscillation score OS is then given by
 | (8) |
Finally, the estimated frequency of oscillation, fosc (Fig. 3F) is given by the frequency having the highest magnitude Mpeak in the band of interest. The strength of the oscillation is given by the oscillation score OS. We have to mention that the oscillation score measure is dimensionless because it represents the ratio of two quantities having the same units of measurement.
To illustrate that removal of the central peak of the ACH is useful, we computed the oscillation score for a data set recorded from cat visual cortex, containing 86 simultaneously recorded multi- and single units (after spike sorting). The data set exhibited strong oscillations in the beta-high band (27–30 Hz) for some experimental conditions/cells and nonoscillatory behavior for the others. We expected a bimodal distribution of oscillation scores. Indeed the oscillation score distribution was bimodal when the central peak of the ACH was removed (Fig. 3J), reflecting the weak/nonoscillatory versus strong oscillatory conditions. However, when the central peak of the ACH was kept intact, the ability to discriminate between the two conditions was lost (Fig. 3K). This effect is due to the previously described estimation problems that are introduced by the central peak of the ACH and justifies the operation of removing the peak.
Estimation of the oscillation score's confidence
A relatively low number of spikes is frequently encountered in practice (Graham and Field 2006) and thus the ACH becomes noisy and can produce the false impression of oscillations (see example in Fig. 7A). A few random peaks in the ACH might appear as an oscillatory process and produce a magnitude peak in the frequency spectrum. Moreover, if the number of spikes is very low, the baseline magnitude, Mavg, of the ACH spectrum is also low. This would boost the oscillation score to erroneously high values (see Eq. 8). Thus a method for estimating the confidence (reliability) of the oscillation score is needed.
|
 | (9) |
OSt is the estimate of SD for the oscillation score computed on individual trials, OSi is the oscillation score computed on trial i, and 
St is the average over OSi.
When the trials contain sufficient spikes and the oscillation is stimulation-stationary (equally expressed in all trials), then the estimated SD will be small. A possibility of directly quantifying how confidently the oscillation score represents the real oscillation score is by computing the coefficient of variation Cv, defined as
 | (10) |
 | (11) |
For small values of the SD relative to the mean, the CS gives a value close to 1, whereas for very poor data, with highly fluctuating oscillation scores on individual trials, the CS approaches 0. Only oscillation scores having a high confidence score should be considered as meaningful and given physiological interpretation. The confidence score is not to be confused with the statistical concept of confidence intervals. The confidence score cannot be used for significance testing because it provides only a bounded measure of how much the scores vary across individual trials. A high confidence score means that variability over trials is low, whereas a low confidence score means that variability is high. Thus the confidence score is a measure of stability of the oscillation score over individual trials.
Here we emphasize another important aspect. When the SD of the oscillation score is computed (Eq. 9), the individual oscillation scores are independently computed for each trial. This means that for each trial, a different frequency, fosc, with the peak magnitude within the band of interest (Fig. 3F), can be selected by the method. If the oscillation strength is constant but the oscillation frequency is drifting over trials, then the final oscillation score will be relatively low (because averaging the ACHs with different oscillation frequencies would smear out the oscillations). However, the confidence score will nevertheless be high (because the oscillation score was stable across individual trials). This is a consistent interpretation but it reflects only the fact that the mean oscillation frequency in the average ACH is imprecisely expressed over individual trials. To cope with such situations, one can alternatively compute a frequency confidence score by replacing in Eqs. 9 and 10 the oscillation score with the oscillation frequency. This new measure would reflect the stability of the oscillation frequency over trials. In the experimental data that we have analyzed here, the oscillation frequency was very stable across individual trials. In general, however, especially in data recorded from behavioral experiments, the oscillation frequency might change as a function of the behavioral task. In such cases, one must choose periods (or define trials appropriately) where the oscillation frequency remains relatively stable. Importantly, nonstationary oscillation frequencies are a general problem and affect most of the present methods for estimating oscillation strength.
Simulated spike trains
When estimating the range of oscillation scores for different frequency bands, we needed to have an objective estimate of the real oscillation strength in the spiking data. To this end, in addition to using recorded data from cat visual cortex, we devised an algorithm that was capable of producing realistic spike trains and realistic-like ACHs. To produce spikes, the algorithm uses an oscillatory discharge probability function that is obtained from mixing two components (one oscillatory and one background process) whose frequencies are continuously modulated by slower processes. The slow processes randomly diverge around the desired oscillation frequency, within a given limit (between a low and high boundary for the frequency), and have a given amount of memory (dependence on their past). By controlling the amount of memory and the frequency boundaries, one can obtain spike trains that oscillate more or less precisely at the frequency of interest. The oscillatory component's frequency is modulated by a relatively stable slow process with strong dependence on its past and precisely bounded in a narrow frequency band (e.g., 23–27 Hz). The background component's frequency is modulated by a highly unstable slow process with low dependence on its past and bounded only to a wide frequency range (e.g., 0–100 Hz). The mixing of the two components allows one to control the amount of oscillations independently of the firing rate of the generated spike train. The firing rate is controlled, independently of the amount of oscillations, by modulating the integral under the final oscillatory probability discharge function. In addition, we introduced realistic values (determined on recorded spiking data) for refractory periods, probability of bursting, interspike intervals, and intraburst interspike intervals. We obtained quite realistic ACHs, well resembling the ACHs computed on recorded data (Fig. A1, F and G). For details on how the artificial spike trains are produced, see the APPENDIX.
|
To characterize the degree of synchrony between LFP and spike trains we used the spike-field coherence (SFC) measure (Fries et al. 2001). To compute the SFC, first, the spike-triggered average (STA) is obtained by averaging LFP segments around each spike and, second, the power spectrum of the STA is normalized to the average power spectrum of the LFP segments used to compute the STA. The resulting measure is independent of the number of spikes and of the LFP power. We also used a second method to quantify SFC, based on multitaper spectral estimates (Pesaran et al. 2002). The coherence is computed by normalizing the cross-spectrum to the product of individual autospectra (computed on the spike train and LFP, respectively). Note that power is used instead of magnitude spectra.
|
|
RESULTS |
|---|
|
Applying the oscillation score on spiking data
We have applied the oscillation score measure on spiking data recorded from cat visual cortex. Most of the analyses have been performed on spike-sorted single-unit activity that is characterized by much lower firing rates compared with those of multiunit activity. The analyzed database contained a large number of individual units (n = 416) that were activated with various visual stimuli, ranging from drifting sinusoidal gratings and moving bars, to plaids and movies with natural scenes. We could reliably identify cells oscillating in the theta band (4–8 Hz), alpha band (8–12 Hz), and beta/gamma bands (25–33 Hz). Examples of ACHs, together with their smoothed, peakless versions, and their frequency spectra are shown in Fig. 4.
|
THE THETA SCORE.
The theta score is defined for the theta oscillation band (Maurer and McNaughton 2007), with fmin = 4 Hz and fmax = 8 Hz. In Fig. 4A, a recording site (multiunit), stimulated with a moving bar, showed robust oscillations in the theta range, with fosc = 5 Hz (Fig. 4A, bottom). Note the high magnitude of the theta oscillation in the frequency spectrum, relative to the baseline magnitude (Fig. 4A, bottom, inset), inducing a very high theta score of 60.7 (Fig. 4A, bottom). Also note that removal of the central peak of the ACH (Fig. 4A, top) introduces an artificial frequency of 2 Hz (Fig. 4A, bottom, inset), which is lower than the inferior boundary of the theta range, thus not affecting the oscillation score.
THE ALPHA SCORE.
According to Klimesch et al. (2007), we defined the alpha score as spanning frequencies between fmin = 8 Hz and fmax = 12 Hz. In Fig. 4B, we identified a single unit that had a low spike rate but exhibited alpha oscillations with a frequency of 9 Hz, yielding an alpha score of 28.8 (Fig. 4B, bottom). Note that the score is smaller than the previously shown theta score.
THE BETA SCORE.
We defined the beta score corresponding to the beta frequency band (12–30 Hz). However, because of the fuzzy, often overlapping border between the beta and gamma bands, we recommend that the beta range be split into two different subbands: beta-low (12–20 Hz) and beta-high (20–30 Hz). The literature reports are often ambiguous on the definition of borders between beta and gamma oscillations. In some reports, frequencies >20 Hz are addressed as gamma (Whittington et al. 1997), whereas in others, the range between 20 and 30 Hz is classified as beta2, and the range >30 Hz as gamma (Roopun et al. 2006). Because of our procedure of cutting the central ACH peak, we strongly recommend that beta-low and beta-high (beta2) bands be treated separately. Therefore we defined a beta-low score (with fmin = 12 Hz and fmax = 20 Hz) and a beta-high score (with fmin = 20 Hz and fmax = 30 Hz), as two separate measures for the beta range.
In Fig. 4C we identified a single unit that exhibits very strong oscillations in the beta-high frequency range (27 Hz). Note that although the oscillation is strong (Fig. 4C, top), the beta-high score reaches a value of only 19.1 (Fig. 4C, bottom), which is much lower than the observed theta and alpha scores.
THE GAMMA SCORE. For the gamma-band oscillations, we defined a gamma score in the frequency range of 30–80 Hz. As mentioned earlier, the gamma band activity has a fuzzy border with the beta-high band. Our experience has shown that the beta-high range (20–30 Hz) can be safely merged with the gamma-low band (30–50 Hz) when setting the limits fmin and fmax. Removal of the central peak will not affect frequencies between 20 and 50 Hz, if fmin is set to 20 Hz. Moreover, high oscillation frequencies (50–80 Hz) should probably be isolated from lower gamma frequencies. We recommend that one investigates these two bands (30–50 and 50–80 Hz) separately.
Interpretation of the oscillation score
The oscillation score represents the ratio between the magnitude of the oscillation frequency and the average magnitude of the spectrum. As we have previously seen, the actual oscillation score depends very much on the frequency band of interest, with low-frequency bands assuming larger scores than high-frequency bands. Indeed, it was previously reported that slow-frequency oscillations normally have a higher power than that of high-frequency ones (Freeman et al. 2000). This is partly due to the fact that slow oscillations tend to engage larger populations of neurons and persist for longer periods of time, compared with the fast, transient oscillations. It is very likely that this phenomenon is related to conduction delays (Buzsáki 2006). These tend to impair the spread of fast oscillations, but not of the slow ones. Moreover, the brain tissue is also known to have certain electrical filtering properties that allow slow oscillations to propagate more easily than fast ones, and engulf larger populations of neurons (Buzsáki 2006; Nunez and Srinivasan 2005).
In light of the evidence about the distribution of biological frequency magnitudes, we expected that oscillation scores assume higher values for low than for high frequencies (because they have higher magnitudes relative to the baseline magnitude of the whole spectrum). This was indeed the case, as indicated by the results shown in Fig. 4. The main question, then, was how should one interpret the oscillation score for a given band? Does the obtained value reflect very strong or rather weak oscillations?
Two remarks are in place here. First, one should not directly compare oscillation scores across different frequency bands (especially not across biological bands). Second, the border between the nonoscillatory and oscillatory regimes should be separately assessed for each individual band. To determine these borders, we used experimental data recorded from cat visual cortex and also simulated spike trains.
The experimental data were used to estimate the distribution of oscillation scores for the nonoscillating regime. We started with a large number of spike-sorted units (n = 725), recorded from four adult cats, in seven distinct experimental sessions, and with various types of stimuli, giving a total of 7,501 ACHs across different stimuli (see METHODS). Of these units, only 416 were different because sometimes the same unit was recorded in more than one session. We included only experimental sessions in which no oscillations were present in any of the frequency bands of interest. To identify these sessions, we visually screened the ACHs for signs of an oscillatory modulation. We could not totally eliminate very weak, stimulus-dependent oscillations, as the visual estimation is subjective. Using the large data set, we estimated, for each frequency band (theta, alpha, beta-low, beta-high, and gamma), the distribution of oscillation scores for the nonoscillating regime (Fig. 5 A, Exp. Data). As expected, the distribution of oscillation scores was wider for the theta and alpha bands. We identified the 95th percentile of each distribution and used it as an estimate of the threshold between nonoscillating and weakly oscillating regimes, for each frequency band (Fig. 5B, thick dotted line).
|
In the next step, we generated simulated data, at the other extreme, exhibiting very strong oscillations (Fig. A1G). Again, we computed the distribution of oscillation scores (Fig. 5A, Simulation Osc.) and, this time, used the 5th percentile, on the lower side of the distribution, to estimate the threshold between moderate and very strong oscillations (Fig. 5B, thick black line between the transition and oscillating areas). As expected, low-frequency oscillations yielded much larger scores than did high-frequency oscillations. We found a good match between the range of scores obtained with simulated and experimental data exhibiting oscillations (Fig. 4). On average, experimental data with strong oscillations produced scores in the range of 40–90 for the theta band and scores in the range of 10–20 for the beta-high/gamma bands (data not shown).
The area with scores between the nonoscillating and strongly oscillating regimes can be considered as a fuzzy, transition regime (Fig. 5B, Transition). Here, there was a continuum of monotonically increasing oscillation scores, from weak to strong oscillations.
We emphasize that the distributions in Fig. 5A are empirically determined and the thresholds of 95 and 5% are guidelines that delineate the borders between different oscillating regimes. They are not to be considered as direct thresholds for statistical significance testing. In practice, the oscillation strength varies continuously from weak oscillations to very strong ones and one should devise appropriate statistical tests for the exact question being asked. For example, one could compare the oscillation scores in two different experimental conditions and test the significance of differences, considering the appropriate type of test, depending on the shape of score distributions, correcting for multiple comparisons, and so forth.
Performance of the oscillation score
Our previous analyses of experimental and simulated spike trains indicate that the oscillation score indeed reflected, to a high degree, the oscillation strength. However, we wanted to have another, independent quantification of the meaningfulness of the oscillation score. For that purpose, we relied on the well-known relationship between oscillations in the LFP and oscillations in spiking activity. It has been reported that neuronal activity can be either weakly related to nearby LFP activity (Kreiman et al. 2006) or strongly related, in which case the spike-field coherence is high (Fries et al. 2001; Pesaran et al. 2002). For strongly oscillating LFP activity, in a given frequency band, it is expected that neurons will also strongly oscillate. When spike-field coherence is high in the given frequency band (Fries et al. 2001), the respective spike trains exhibit precise phase locking with the respective LFP oscillations and thus also exhibit a strong oscillatory modulation. We tested whether, for strongly oscillating LFP activity, the oscillation score of individual neurons correlated with spike-field coherence.
We have analyzed simultaneously recorded spiking and LFP activity from cat area 17. We selected a session with strong oscillations of LFPs in the beta-high band (27 Hz). We first quantified the strength of LFP oscillations in a manner similar to the oscillation score: we computed the FFT spectrum on the LFP signal and divided the peak magnitude, in the beta-high band, to the baseline magnitude of the whole spectrum. We called this measure "LFP oscillation score" and it revealed strong and reliable LFP oscillations in all stimulation conditions (14 stimuli; see METHODS) and for all 32 electrodes (Fig. 6 A). In the second step, we spike-sorted the units recorded on the 32 electrodes and obtained 86 single units. For each unit, we considered in the analysis only the ACHs for which the oscillation score could be reliably determined (484 ACHs; see Confidence estimation). The oscillation scores for these units, shown in Fig. 6B, revealed a bimodal distribution, with a group of nonoscillating or weakly oscillating responses (OS < 10), and a separate group of strongly oscillating responses (OS > 10). Some units did not exhibit oscillations for any of the 14 stimulation conditions.
|
Next, using the method described in Fries et al. (2001), we computed the spike-field coherence between single units and the corresponding LFP signal, in the beta-high band (20–30 Hz). We found that spike-field coherence was highly correlated to the spike oscillation score (r = 0.92; r2 = 0.85; Fig. 6D). When coherence was computed with a second method, based on multitaper spectral estimates (Pesaran et al. 2002), similar results were obtained (r = 0.86; r2 = 0.74; data not shown). Because the LFPs exhibited strong oscillations, a correlation between the real strength of oscillations in the spike train and coherence was expected (when the coherence between spikes and a strongly oscillating LFP is high, it follows that the spikes also exhibit strong oscillatory behavior). The observed correlation between the spike oscillation score and coherence indicates that the former is a meaningful and reliable measure of the strength of neuronal oscillations. Also, it is important to note that the spike oscillation score was computed after applying the nonlinear transformation of removing the central peak from the ACH and computing a ratio between spectral magnitudes. The fact that the spike oscillation score correlated strongly with the spike-field coherence is a good indication that these procedures do not impair the estimation of the real oscillation strength.
Confidence estimation
So far, we have not estimated the degree to which an oscillation score, for a given ACH, can be trusted. For single units that exhibit very few spikes, it might happen that the ACH (computed on all trials of a given condition) has false peaks, due to chance. There is probably no perfect way to estimate whether such peaks occur by chance or represent legitimate oscillation peaks. However, when multiple trials with the same stimulation condition are available, one can require that most trials exhibit roughly the same ACH structure (stimulation-condition stationary structure).
We developed a measure for quantifying the confidence with which one can trust the oscillation score, called the confidence score (see METHODS). We have to mention here that the confidence score is not related to the concept of confidence intervals, from statistics. Rather, it reflects the inverse of the variation, thus stability. The measure is bound between 0 (no confidence) and 1 (full confidence) and is based on the estimation of the coefficient of variation for the oscillation score, on a set of recorded trials. We analyzed a mixed data set with 86 single- and multiunits, some exhibiting strong oscillations in the beta/gamma band (score >10; n = 53), some others showing very weak or no oscillations at all (score 10; n = 33).
Figure 7 A shows as an example a single unit with a very low firing rate. The analyses revealed an oscillation frequency fosc = 19.5 Hz (beta band) and an oscillation score OS = 10.4, indicating relatively strong oscillations (see also Fig. 5B). However, the shape of the ACH does not appear reliable. Indeed, the confidence score was only CS = 0.384, indicating that one should not trust the oscillation score OS. Figure 7B, by contrast, depicts a single unit with an oscillation frequency fosc = 29.3 Hz, having an oscillation score OS = 15.7 (very strong oscillations; see also Fig. 5B) and a high confidence score CS = 0.85.
Finally, we investigated whether there was a dependence between the number of available spikes for computing the ACH and the confidence score. Most methods relying on correlation analysis cannot reliably estimate the oscillation frequency and oscillation strength when the total number of spikes is very low, and thus they need to rely on multiunits (König 1994) or consider many experimental trials in the analysis (Samond and Bonds 2005). We investigated how the confidence score CS depends on the number of spikes that were used to compute each of the 1,204 ACHs (86 units x 14 stimulation conditions). Figure 7C revealed a surprising result. For spike counts up to about 200, the confidence score increased linearly with the logarithm of the spike count (note that the horizontal axis, representing the number of spikes used to compute the ACHs, is shown on a logarithmic scale). Then, the confidence saturated around CS = 0.8, occasionally reaching higher values (0.9), as the spike count increased further. There are two important conclusions that stem from these results. First, oscillation scores computed on ACHs that have a confidence score CS <0.65 should not be trusted. The reason is that, by adding more spikes to the ACH, the confidence would rapidly increase, showing that the oscillation score tends to become more stable across trials when more spikes are available. One cannot know whether the estimation of the oscillation score is reliable, when CS <0.65, because it is relatively unstable over individual trials. Second, the confidence with which one can trust the oscillation score converges exponentially fast with the number of spikes (in this case for spike counts up to 200). This means that even a small number of spikes might suffice to yield a confident estimation of the oscillation score. In our data set, 200 spikes represented firing rates of about 2.5 Hz (20 trials, each of 4-s duration). Thus the oscillation score measure lends itself to the analysis of single-unit neuronal activity and is also applicable to cases where the number of available experimental trials is small.
Comparison to other methods
We have seen that the oscillation score performs remarkably well on experimentally recorded data. However, it is very important to analyze the relation between the method proposed here and other existing methods both qualitatively and quantitatively. For this purpose, we implemented three additional methods: the generalized Gabor fitting (GGF), described in König (1994); the secondary peak–valley difference (SPVD), described in Samonds and Bonds (2005); and spectral techniques applied directly on the spike train. We applied each method on three artificially generated spiking data sets (see APPENDIX) having trials with a length of 30 s. Each of the three data sets had a different baseline firing rate: 10, 27, and 50 Hz, respectively. The oscillation frequency was fixed at 25 Hz (beta-high band). In addition, at a fixed firing rate, we independently and systematically varied the oscillation strength (Fig. 8 A) by using a parameter (o, described in the APPENDIX) bounded between 0 (no oscillation) and 1 (very strong oscillation). For each estimation of oscillation strength, only a single trial was used, to push the methods to their limit. We computed 100 estimations for each pair of firing rate–oscillation strength and computed the mean, SD, and coefficient of variation of these estimations.
|
and used in many studies, approximates the shape of the ACH with a generalized Gabor function having eight parameters. The function is iteratively fit to the correlogram by a Marquardt–Levenberg algorithm that tries to minimize the fitting error. Since the method relies on an iterative search for a minimum, it suffers from the problem of local minima. The larger the number of parameters, the worse the problem becomes. To compensate for this, we used 1,000 initial guesses for the parameters, also using an informed technique that takes into account the mean of the ACH, the size of the central peak, and so forth. After a successful fitting, the strength of oscillation is computed as the ratio between the size of the first satellite peak (relative to the baseline of the ACH) and the baseline of the ACH (see Fig. 1A; Brecht et al. 1999; Konig 1994). We identified and summarized a set of qualitative problems of the GGF method. The first and most important one is related to the poor convergence since the fitting-error landscape has many local minima. As a consequence, one must use many initial conditions. Choosing these initial conditions, however, is not an easy task since initial conditions can depend on such parameters as the frequency of interest and the baseline of the correlogram, for instance. In addition, the large number of initial conditions (1,000 in our case) can render the method very slow. Second, the GGF algorithm is not straightforward to implement since it requires advanced minimization techniques such as Marquardt–Levenberg or Newton. Third, as mentioned earlier, the method does not fit correlograms that depart from the Gabor model. This was quite frequently the case for the data that we analyzed, yielding particularly large fitting errors when oscillations were weak, or when rate covariation or slower oscillations were present. Also, the method cannot cope with multiple oscillation frequencies since there is only one main frequency of the Gabor. Enriching the model to incorporate multiple oscillation frequencies would require additional parameters and would pose even more difficulties to the fitting algorithm. Fourth, the choice of the correlation window is very important when rate covariation is present. In contrast to the oscillation score, which automatically adjusts the window as a function of the frequency of interest, the GGF requires the user to choose the window and this is not exactly straightforward. In Fig. 8B we present a case where the choice for an overly large window of an ACH exhibiting slow rate covariation/slow oscillations leads to the underestimation of the size of the first satellite peak and thus the oscillation strength. For a large data set, with many ACHs to fit, an automated estimation method based on the GGF cannot easily decide on the optimal size for the ACH window. Fifth, as discussed before, the GGF method might have problems when refractoriness and rate covariation are simultaneously present, producing false oscillation peaks. Finally, we also noticed that for cases without or with weak oscillations the GGF method produces highly distorted distributions for the oscillation strengths, because the ACHs that are not fitted well, frequently have an estimated oscillation strength of 0 (Fig. 8C; notice the very high peak at 0). Such distorted distributions can pose serious problems to parametric statistical methods that are not robust against the violation of certain assumptions, such as unimodality. In contrast, the oscillation score provides comparatively smooth distributions that can easily be used in statistical tests (Fig. 8D).
We next estimated the "quantitative" behavior of the GGF by applying it to the three artificially generated data sets having various firing rates and systematically varied oscillation strengths. For high firing rates (27 and 50 Hz), the GGF method performed relatively well. However, for the case of very weak/no oscillations, even at high firing rates, the GGF method became very imprecise (Fig. 8A, Coeff. Variation). This phenomenon was due to the fact the GGF model can only poorly fit the ACHs with weak or no oscillations. Very frequently, the strength of oscillation, for these cases, was 0. In some cases, the oscillation strength was also different from 0, having small values. This relation critically depended on the optimality of automatic fitting, leading to a distorted distribution of oscillation strengths, similar to the one shown before in Fig. 8C. In addition, for any oscillation strength, at a given firing rate, the coefficient of variation for the GGF method was larger than its counterpart computed for the oscillation score (Fig. 8A, GGF and OS).
For a low firing rate of 10 Hz, which produced very noisy ACHs (due to the length of the spike train of only 30 s), the GGF model became very unreliable. The coefficient of variation was very large at all oscillation strengths (Fig. 8A, GGF). Although the expected value for the oscillation strength was relatively good, there was a high fluctuation on individual spike trains, even for the case of strong oscillations, indicating poor fittings.
SECONDARY PEAK–VALLEY DIFFERENCE (SPVD).
A different method for computing the strength of oscillations, proposed in Samonds and Bonds (2005), relies on the difference between the second satellite peak and the second valley of the normalized ACH. The method is relatively simple; the most complex step involves only the computation of the normalized ACH, which also removes shift-predictor components (for details see Samonds and Bonds 2005). The method suffers most when the secondary peak and valley cannot be easily identified, such as in very noisy ACHs (when few spikes are available). Quantitatively, we measured the performance of the method on the three artificial data sets. As expected, the SPVD method performed poorly for low firing rates (10 Hz), where the ACH was very noisy. In that case, it failed to distinguish between strong and weak oscillations (Fig. 8A, SPVD, 10 Hz; note the SDs). For higher firing rates (27 and 50 Hz) the method performed increasingly well. The coefficients of variation were relatively low for all firing rates and oscillation strengths. However, there was an already obvious and very important problem with the SPVD method: the oscillation strength estimations highly depended on the firing rate (Fig. 8A, SPVD). The problem stems from the fact that the difference between the secondary peak and secondary valley does not reflect in any way the baseline of the ACH. The same difference could be "riding" on top of either a small or a large baseline. Since the strength of oscillations should reflect how strongly the oscillatory component is expressed relative to the other components, the SPVD method provides highly inaccurate results.
DIRECT SPECTRAL TECHNIQUES.
Spectral techniques are frequently used to compute the spectrum of spike trains. However, they do not provide a direct quantification of the strength of oscillations. A measure similar to the oscillation score has to be derived for that purpose. We implemented two techniques to compute the spike magnitude spectrum: the multitaper technique (Jarvis and Mitra 2001) and the direct FFT with a Blackman window. After computing the magnitude spectrum, we quantified the strength of oscillations by applying exactly the same procedure as that for the oscillation score (see METHODS, Computation of the oscillation score, step 5). For that reason this category of methods will be denoted by "oscillation score on raw spikes" (OSRS). All the other measures described in conjunction with the oscillation score also apply to the OSRS (confidence score, etc.).
We first implemented the multitaper technique (Percival and Walden 1993) adapted to spike trains, following the method described by Jarvis and Mitra (2001). The method uses a set of orthogonal tapers (here we used discrete prolate spheroidal sequences) to produce a smooth spectrum. The larger the number of tapers, the smoother the estimate and the lower the variance, however, the peak reflecting the oscillation frequency becomes broader and more imprecise. The number of tapers used to compute spike spectra is usually in the range of nine (Pesaran et al. 2002). In our case, such a large number of tapers yielded smooth but very broad peaks in the spectrum at the oscillation frequency. This created imprecision in identifying the peak and also created a spread of power that reduced the size of the peak. For this reason, we used only three tapers and a sliding window of 1024 ms for analysis. Qualitatively, the multitaper method was the most difficult to implement. Also, it was slower compared with the SPVD and the oscillation score. Although the multitaper technique was suggested to yield high performance in estimating the spectrum (Pesaran et al. 2002), our own experience on the well-controlled simulated spike trains showed a rather poor performance of the multitaper-based OSRS. First, we have identified that the modulation of the oscillation frequency in the model spike trains (see the APPENDIX) produced strong violations of stationarity (and this is also the case, in general, for recorded data). For that reason, a proper size for the sliding window needs to be chosen. Windows that are too small yield insufficient spikes for the estimation, whereas windows that are too large could suffer from the nonstationarity of the data (stationarity is a very important assumption when applying the FFT). Quantitatively, the performance of the OSRS for low firing rates was rather poor (Fig. 8A, OSRS, 10 Hz). In general, increasing the length of the spike train does not improve performance because the spectral estimation is local to the sliding window and the noise, spread in all frequency bands, does not completely average out. Increasing the size of the sliding window is not an option either, for the reasons described earlier.
For higher firing rates, the OSRS had better performance, but showed the same firing rate dependence as that of the SPVD method (Fig. 8A, OSRS). To understand the phenomenon, we computed separately the magnitude of the DC component, the magnitude of the peak at the oscillation frequency, and the integral of the rest of the spectrum. The important finding was that, although the DC and the peak scaled with the firing rate, the rest of the magnitude spectrum scaled much slower. This effect, due to the noise of the spectrum estimation (noise that does not get averaged out across different sliding windows), leads to the faster scaling of the peak with respect to the integral of the spectrum and thus the rate dependence of the OSRS. We next investigated whether this scaling was related to the fact that we used the magnitude spectrum instead of the power spectrum for computing the OSRS. According to the Wiener–Khinchin theorem, the FFT of the autocorrelation function gives the power spectrum. Thus we hypothesized that the difference between the oscillation score and OSRS might stem from the fact that we used the magnitude instead of the power spectrum to compute the OSRS. We next computed the power spectrum on the raw spikes by using the same multitaper technique and again computed the OSRS. The results were qualitatively the same, exhibiting strong rate dependence.
In addition to the multitaper technique, we also applied a more classical method, that is, the FFT directly on the spike train, with a sliding window of 1,024 ms and Blackman windowing. We simply computed the FFT, either directly on the spike train or on a smoothed signal, obtained from the spike train using a Gaussian kernel with SD of 2 ms. Both methods yielded similar results. Surprisingly, the spectrum estimation with this simplified method was much more precise than that for the case of the multitaper. The rate dependence was still present and strong, although discriminability increased dramatically (data not shown) even for low firing rates. The SDs were very small and so were the coefficients of variation. Qualitatively, the performance curves looked very similar to those in Fig. 8A (OSRS), with the difference that SDs were much smaller and thus the estimations more precise. However, both the OSRSs, based on direct FFT and on multitaper spectral techniques, were heavily biased by the firing rate.
We conclude that the rate dependence of the OSRS stems from the highly nonlinear and nonuniform scaling of different components of the spectrum with the firing rate. The effect, for the particular artificial data sets that we analyzed, is due to the accumulating noisy spectral components that are spread in all frequency bands and that do not scale with rate as fast as the DC and the oscillatory component. It is expected that this effect is also very pronounced for experimentally recorded spike trains that exhibit highly nonstationary, transient frequency components, especially in awake preparations.
THE OSCILLATION SCORE. In addition to the three different methods, we also computed the oscillation score on the three artificial data sets. The oscillation score was the fastest method, and by far the simplest to implement (except for applying the FFT directly on spikes). It also yielded the best results. For high firing rates (27 and 50 Hz) the estimations of the strength of oscillations were very precise (see Fig. 8A, OS, coefficients of variation) and completely independent of the firing rate.
For the low firing rate of 10 Hz, when the ACH was particularly noisy, the oscillation score (OS) still performed remarkably well. For the case without oscillations, the OS had a tendency to overestimate the strength of oscillations (Fig. 8A, OS, 10 Hz). The scores obtained were still in the safe range of the very weak/nonoscillating regime (see Fig. 5B at 25 Hz). This problem was described before as being introduced by the noise in the ACH, which produces a slight overestimation of the strength of oscillations (see Fig. 7A for an extreme case). However, with the confidence score at hand, one can easily detect these cases and carefully judge how to interpret the score. Again, we mention that the overestimation is still in the "safe" range/boundary of the very weak oscillations. For medium and strong oscillations and a firing rate of only 10 Hz, the OS had a tendency to slightly underestimate the oscillation strength. However, the same was the case with the GGF and this was probably due to the signal-to-noise ratio that was low, due to the small number of spikes. It is worth noting that the discriminability between weak/no-oscillations and strong oscillations was excellent, even for the case with low firing rates and highly noisy ACHs.
|
|
DISCUSSION |
|---|
|
) and theoretically (Brunel and Hakim 1999). The second relies on preferential locking of neurons, induced by their individual frequency preferences (Hutcheon and Yarom 2000; Murean and Savin 2007). The latter process is also called membrane resonance and is known to be voltage modulated (Puil et al. 1994). It is very likely that neuronal oscillations result from a combination of these mechanisms, dependent on the frequency of oscillations. It is also likely that different types of oscillation have a variety of distinct functional roles in the brain. Direct evidence for the functional role of oscillations was put forward in very recent studies. In Montgomery and Buzsáki (2007) gamma oscillations were shown to dynamically couple hippocampal CA3 and CA1 regions. Lakatos et al. (2007) have shown that oscillatory activity can dramatically modulate multisensory interaction in the primary auditory cortex, providing a "context" within which sensory information is integrated. Their results also suggested a complex patterning of different oscillatory rhythms and this was also found to be the case in Isomura et al. (2006) who showed that slow oscillations can coordinate the integration/segregation of activity in the hippocampus and entorhinal cortex. Schaefer et al. (2006) further suggested that neuronal oscillations can enhance stimulus discrimination by ensuring a precise control over the timing of action potentials. The timing of spikes, relative to the cycle of the ongoing gamma oscillation, was also proposed as a flexible neuronal code (Fries et al. 2007) and it was suggested that the windows of opportunity created by the oscillation cycles may promote coherent coordination of neuronal populations in the temporal domain (Fries 2005). There is now also evidence from studies on humans that neuronal oscillations constrain the timing of firing of cortical and hippocampal neurons, across different frequency bands (Jacobs et al. 2007). Taken together, these recent studies and many others as well (e.g., Engel et al. 1990; Fries et al. 2001; Klimesch et al. 2007; Lutz et al. 2004; Martinovic et al. 2007; Meador et al. 2002; Melloni et al. 2007; Murata et al. 2004; Pesaran et al. 2002; Rodriguez et al. 1999; Samonds and Bonds 2005; Schoffelen et al. 2005; Singer 1999; Tallon-Baudry et al. 1997, 2004; Uhlhaas and Singer 2006; Vidal et al. 2006) suggest that neuronal oscillations do matter in the brain. The method that we have proposed here could provide a useful tool for investigating the functional role of neuronal oscillations because it allows the assessment of oscillatory patterns in the spiking activity of single neurons. The oscillation score gives a reliable estimate of the strength of oscillations in a given frequency band. The described method is simple and fast. First, the original ACH is smoothed with a fast kernel, to remove the noise. Then a slower Gaussian kernel is applied to detect the extension of the central peak, which is subsequently removed from the fast-smoothed ACH. An FFT estimates the spectrum of the peakless smoothed ACH, and next, the oscillation score is computed as the ratio between the maximum frequency magnitude in the band of interest and the baseline magnitude of the entire frequency spectrum. As we have seen, the obtained oscillation score depends on the frequency band, being normally higher for lower frequencies. Moreover, when the LFPs show strong oscillations, the score is highly correlated with the spike-field coherence. Finally, the method is applicable to units that have very low firing rates (i.e., spike-sorted single units) and converges to a reliable, confident oscillation score, exponentially fast with the number of spikes used to compute the ACH.
In the following, a few important aspects need to be discussed. First, we shall dwell upon the relation of our method to other similar techniques for estimating the oscillation strength in neurons. Second, we will scrutinize the correct applicability and the weaknesses of the oscillation score method, proposing a few solutions. Finally, we will discuss some general observations that stem from our own experience on analyzing large experimental data sets.
Relation to other methods
The oscillation score falls into the category of methods that rely on the correlation analysis of neuronal spike trains. Thus, it is not easily comparable to methods that are based on spectral estimations directly from spike trains (Bair et al. 1994) or from taper-smoothed spike trains (Jarvis and Mitra 2001). However, it can be readily compared with methods that estimate oscillation strength from auto- or cross-correlograms (König 1994; Samonds and Bonds 2005). The oscillation score is different in many respects from these methods. First, the presented method does not assume any underlying model that needs to be fitted to the correlogram (like the Gabor model; König 1994). This fact allows it to cope properly with a variety of correlogram shapes, and it is, thus, not sensitive to any kind of initial conditions (since it does not fit a function, initial conditions are not defined). We have to mention, however, that the oscillation score relies exclusively on autocorrelograms and is not applicable to cross-correlograms. As a result, our method is not useful for the quantification of neuronal synchrony between pairs of neurons, but only for the quantification of oscillation strength (note that synchrony can also be expressed independently of oscillations, depending on the size of the involved population; Bhattacharya 2001; Brecht et al. 1998). For synchronization analysis other methods are available and perform reasonably well (König 1994; Womelsdorf et al. 2007).
Second, our method does not rely on the estimation of peak sizes in the ACH. As a consequence, even when the number of available spikes is small, it can reliably estimate the oscillation strength from the frequency spectrum of the ACH. Other methods need a considerable amount of data to properly identify and estimate the sizes of peaks and troughs in the correlogram (Samonds and Bonds 2005). Moreover, we removed the central peak, and sometimes also its side troughs, so as not to induce false oscillatory peaks in the frequency spectrum.
Third, the oscillation score can cope with multiple simultaneous oscillation frequencies and it can correctly estimate their relative expression in the spiking activity. The reason is that it searches selectively in a desired frequency band to identify the oscillation frequency. When such a frequency is identified, the method normalizes its magnitude to the baseline magnitude of the frequency spectrum. Thus, an estimation of the expression of the identified frequency, relative to other frequencies, is obtained. Such estimations are difficult with methods that rely on the size of satellite peaks in the correlograms.
Finally, the qualitative and quantitative comparison of different methods (Fig. 8) provided very insightful results. The generalized Gabor fitting (König 1994) performs well when the fitting is good. However, the quality of the fit depends on many factors rendering an automated procedure very difficult. One must be very careful in making sure that the fitting was made properly because it can dramatically influence the estimations of the strength of oscillations. The secondary peak–valley difference (Samonds and Bonds 2005) method performs very poorly for low firing rates and is also highly biased by the firing rate, leading to false conclusions about the strength of oscillations. Spectral techniques (Bair et al. 1994; Jarvis and Mitra 2001; Pesaran et al. 2002) are sensitive to the spectral smoothing that is applied and also to the size of the sliding window. Violations of stationarity can dramatically affect the spectral estimations. They are also biased by the firing rate since the baseline of the spectrum does not scale with the firing rate as fast as the peak at the oscillation frequency (at least in the model spike trains that we used). We conclude that the best is to compute the spectrum on the ACH, not directly on spikes. In addition, with the central peak of the ACH removed, one can obtain very accurate results, such as those provided by the oscillation score. It is also worth mentioning that the techniques based on ACHs can take advantage of longer traces of data. Even for much lower firing rates than those we used for testing (10–50 Hz), if one has a sufficient number of trials or sufficiently long trials, the ACHs can become smooth enough to provide very accurate estimations (see Fig. 7 and comments). This is, in general, not the case with spectral techniques because the constant noise in the sliding window, spread across all frequency bands, did not get averaged out with longer traces of data, at least for the case of spiking data investigated here. The comparison of different methods, presented in Fig. 8A, strongly suggests that methods based on correlation analysis are still very effective tools that can outperform more modern, direct spectral techniques, at least on some tasks, such as the estimation of the oscillation strength in spike trains. The oscillation score is, in particular, a hybrid method relying on both correlation analysis and spectral techniques, and thus, it can take advantage of both.
Applicability and weaknesses
The only input parameters that are supplied to the method are fmin and fmax—the lower and higher bounds that define the frequency band of interest. These parameters are also the most important since they determine how much information is removed from the vicinity of the central peak of the ACH and also the band within which the oscillation frequency is identified. An incorrect setting for the two parameters could lead to underestimations of the oscillation strength and also to the detection of false oscillation frequencies. However, one should keep in mind that the oscillation score has been designed to look at a single frequency band at a time! If biologically relevant frequency bands are separately used as frequency limits, then the method gives accurate and clean results. These bands are: theta (4–8 Hz), alpha (8–12 Hz), beta-low (12–20 Hz), beta-high or beta2 (20–30 Hz), gamma-low (30–50 Hz), and gamma-high (50–80 Hz). Also note that the frequency band of interest is implicitly required for other methods as well (e.g., for the methods relying on correlation analysis, one needs to choose the size of the correlation window as a function of the frequency of interest).
Another potential problem of the method can occur if one looks in a band higher than the actual oscillation frequency. Very rarely, if the oscillation frequency has a high power relative to the baseline, and the ACH exhibits high levels of noise, the frequency spectrum can display strong harmonics having higher frequencies, at multiples of the true oscillation frequency. By looking into bands higher than the oscillation frequency, one might actually quantify, by mistake, the harmonics of the real frequency. This problem can be easily resolved by calculating the confidence score, which, in our analyses was usually very low (due to noise) for ACHs that produced strong harmonics. Moreover, one should pay attention to oscillation frequencies with a high score, detected at exactly the lower limit of the interval of interest (fmin). This might indicate the presence of a lower frequency of oscillation spreading to higher frequencies in the spectrum due to leakage that occurs at the border. In such cases one should inspect the ACH and the frequency spectrum and, if necessary, select a lower frequency band for inspection.
As a general rule, it is safer to inspect the frequency spectrum and ACH before drawing definitive conclusions on the interpretation of the oscillation score. Moreover, one has to take into account the confidence score before giving any meaningful interpretation to the value of the oscillation score. When the confidence score exceeds the value of 0.65, one is likely to be on the safe side. The software package we provide as support (see APPENDIX) allows one to have access to all relevant variables: the ACH, the smoothed ACH, the peakless smoothed ACH, the frequency spectrum, and so forth.
General observations and concluding remarks
Before concluding, we share a few important remarks. First, we have observed that oscillating single-unit activity yields smoother ACHs and higher oscillation scores than those of the corresponding multiunit activity (originating from the same electrode). Although this result is somewhat counterintuitive (because multiunits have more spikes), this difference can be explained by keeping in mind that cells usually engage in oscillatory behavior when stimulation is optimal (Engel et al. 1990), and thus in such conditions, single units have quite elevated rates, too. Moreover, individual cells that contribute to the multiunit activity, like pyramidal cells and interneurons, could have slightly different oscillation properties. This leads to a less well defined oscillatory behavior in the multiunit average. Our method is also very robust with single-unit activity and therefore one should rely on sorted, single-unit data, whenever possible.
Second, since the oscillation score is a global measure for the ACH that directly characterizes its frequency distributions, the reliability of the oscillation score across individual trials also reflects the reliability, and stationarity, of the ACH. In principle, one could use the confidence score, as defined in this study, to quantify how reliable the ACHs are, from the perspective of stationary responses to the stimulus.
Third, other nonbiological bands can be considered, but one should choose very narrow bands, 5–10 Hz in width. We successfully quantified the oscillation frequencies and oscillation strengths for very slow oscillations (2–3 Hz) and very fast ones (90–100 Hz). However, we used narrow bands of inspection and we also relied on the frequency spectrum to guide this choice.
Fourth, as we have seen, the oscillation score can be well estimated on a relatively small number of spikes. This allows one to take into consideration only a few trials and single-unit activity when estimating the oscillatory behavior of the cell. Ideally, one should use the minimum amount of trials, for which the confidence score reaches a maximum, to avoid the problems of nonstationary oscillation frequencies caused by changes in long-term excitability, cortical state, and so forth. Increasing the number of trials might actually impair the detection of oscillations, if their frequency/expression is drifting in time. The oscillation score together with the confidence score allow for a flexible decision on how to select the experimental data for analysis.
In conclusion, our method behaves very well on single-unit neuronal data. It converges to a confident estimation exponentially-fast with the number of available spikes. We are not aware of any other method that is so robust, relative to the number of spikes required to achieve accurate estimations. The oscillation score, with its various flavors (theta score, alpha score, beta score, and gamma score), should prove a very useful tool for characterizing the oscillatory responses of individual neurons and for helping to advance our knowledge on the very complex, interesting, and fascinating realm of neuronal oscillations.
|
|
APPENDIX |
|---|
|
Open source implementations for the computation of the oscillation score and estimation of the confidence are freely available for download at: http://www.raulmuresan.ro/sources/oscore. The algorithms have been implemented in C++, Delphi, and MATLAB and are also available as Windows Dynamic Link Libraries (DLL).
Algorithms for computing the oscillation score and the confidence score
Generating artificial spike trains
The spike train model tries to mimic basic properties of real neurons such as firing rate, oscillation frequency, and spike-timing properties such as intraburst interspike intervals (ISIs) or refractoriness. These properties were quantitatively determined from the recordings used throughout this study. We used two distinctive control processes to generate the artificial spikes: a global discharge probability function ps(t) (Fig. A1A) and a finer temporal process that controls exact spike timing (refractoriness, bursting, etc.; Fig. A1C).
The spike discharge probability function ps(t) should have the following properties: first, it should allow the spike train to exhibit a preferred oscillation frequency for transient periods of time; second, it should be able to control the stability and strength of the desired oscillation; and third, it should allow some control over the firing rates. The amount of oscillation is controlled by the mixing factor o (Eq. A1) of two discharge probabilities po(t) and pb(t), which correspond to one oscillatory process and a background process, respectively. To obtain a transient oscillation, we start from a sine probability function and we modulate its frequency (Fig. A1A) with a random process fo(t) that takes into account the past values of the frequency and thus has memory and offers stability (Eq. A4). The importance of the history fades in time exponentially, with a decay time constant 
. The function rand[–1,1](t) adds noise, with values in the interval [–1, 1], to the frequency, with a strength factor m. With and m we control how fast the oscillation frequency changes in time. Boundaries can be imposed on the frequency range, such that the frequency of po(t), and thus oscillation preference, becomes stable in a certain range (Fig. A1B, fo). The background activity pb(t) is generated in a similar fashion but, in this case, the stability of the process is lowered and oscillation frequency fb(t) varies in a much broader frequency interval (Fig. A1B, fb). The spikes are generated only in regions where ps(t) is >0 (Fig. A1A), and are concentrated toward the peaks of ps(t) (where the probability of producing spikes is higher; Fig. A1C). The number of generated spikes depends on the positive area of ps(t). With the offset (Fig. A1A) we can control how well the spikes are aligned to the desired oscillation [the peaks of ps(t) approach a probability of 1], whereas with the area of the positive ps(t) we can influence how many spikes are generated overall and thus the firing rate. An interesting case is when ps(t) is <0 most of the time, having high positive values for only short periods of time [when the amplitude of ps(t) is very high while the offset is also very large; see Fig. A1A]. This situation results in very precise firing patterns. In the end, with ps(t) we can control the periodicity of the spike trains in a manner that realistically mimics the preferred oscillation frequency of recorded neurons
 | (A1) |
 | (A2) |
 | (A3) |
 | (A4) |
At smaller timescales the model controls the timings between the spikes and the burstiness of the generated data. Experimentally recorded spike trains frequently exhibit bursts with duration between 3 and 17 ms (Fig. A1D). The simplest way to model the burst probability, pburst(t), is to set it as a fraction of ps(t) (Eq. A5). Once a discharge is initiated this probability function is used to determine whether a tonic spike or a burst is generated (Fig. A1C).
To obtain realistic spike trains, we first computed the distribution of burst lengths (Fig. A1D) and intraburst ISIs (Fig. A1E) by using experimentally recorded spike trains. In the simulated data, if a burst occurs then its length is set according to the measured distribution of burst lengths (Fig. A1D). The intervals between the spikes in the burst are set based on the measured intraburst ISI distribution (Fig. A1E). Immediately after a tonic spike, or a burst, a refractory period (rtonic and rburst, respectively, in Fig. A1C) prevents any discharge from being initiated. These refractory periods are uniformly distributed between 5–10 and 10–20 ms, respectively. Refractory periods, intraburst ISIs, burst lengths, and a burst probability constant, k = 0.8, induce a biologically plausible center-peak shape in the ACH (Fig. A1G), which is lost if k is set to 0
 | (A5) |
In Fig. A1F we show the ACH of an experimentally recorded neuron from cat area 17, whereas in Fig. A1G we present the ACH of a corresponding artificial spike train. The spike train was generated to match the basic properties of the neuron in Fig. A1F: the oscillation frequency (29 Hz; see also Fig. A1B), discharge rate (7.5 Hz, real; 8.1 Hz, artificial), and the bursting parameters (Fig. A1, D and E). The two ACHs, for the real and simulated spike trains, look remarkably similar.
|
|
GRANTS |
|---|
|
iei Cercet
rii i Tineretului (MECT) / Unitatea Executiv pentru Finaarea Învmântului Superior i a Cercetrii 
tiinifice Universitaire; and Partnerships Program Neurobot Contract 11039/18.09.2007, financed by MECT/Autoritatea Naional pentru Cercetare tiinific and a Deutsche Forschungsgemeinschaft grant NI 708/2 - 1.
|
|
ACKNOWLEDGMENTS |
|---|
|
ca for useful comments on earlier versions of the manuscript and interesting discussions. We also thank A. Furcoane and P. Yikström for support during revision of the manuscript. |
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: R. C. Murean, Center for Cognitive and Neuronal Studies, Str. Saturn 24-26, 400504 Cluj-Napoca, Romania (E-mail: contact{at}raulmuresan.ro)
|
|
REFERENCES |
|---|
|
Baker SN, Pinches EM, Lemon RN. Synchronization in monkey cortex during a precision grip task. II. Effect of oscillatory activity on corticospinal output. J Neurophysiol 89: 1941–1953, 2003.
Bhattacharya J. Reduced degree of long-range phase synchrony in pathological human brain. Acta Neurobiol Exp (Warsz) 61: 309–318, 2001.[Medline]
Biederlack J, Castelo-Branco M, Neuenschwander S, Wheeler DW, Singer W, Nikoli D. Brightness induction: rate enhancement and neuronal synchronization as complementary codes. Neuron 52: 1073–1083, 2006.[CrossRef][Web of Science][Medline]
Brecht M, Singer W, Engel AK. Correlation analysis of corticotectal interactions in the cat visual system. J Neurophysiol 79: 2394–2407, 1998.
Brecht M, Singer W, Engel AK. Patterns of synchronization in the superior colliculus of anesthetized cats. J Neurosci 19: 3567–3579, 1999.
Brody CD. Correlations without synchrony. Neural Comput 11: 1537–1551, 1999.[CrossRef][Web of Science][Medline]
Brunel N, Hakim V. Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput 11: 1621–1671, 1999.[CrossRef][Web of Science][Medline]
Buzsáki G. Rhythms of the Brain. Oxford, UK: Oxford Univ. Press, 2006.
Buzsáki G, Draguhn A. Neuronal oscillations in cortical networks. Science 304: 1926–1929, 2004.
Buzsáki G, Geisler C, Henze DA, Wang XJ. Interneuron diversity series: circuit complexity and axon wiring economy of cortical interneurons. Trends Neurosci 27: 186–193, 2004.[CrossRef][Web of Science][Medline]
Canolty RT, Edwards E, Dalal SS, Soltani M, Nagarajan SS, Kirsch HE, Berger MS, Barbaro NM, Knight RT. High gamma power is phase-locked to theta oscillations in human neocortex. Science 313: 1626–1628, 2006.
DeBusk BC, DeBruyn EJ, Snider RK, Kabara JF, Bonds AB. Stimulus-dependent modulation of spike burst length in cat striate cortical cells. J Neurophysiol 78: 199–213, 1997.
Engel AK, König P, Gray CM, Singer W. Stimulus-dependent neuronal oscillations in cat visual cortex: inter-columnar interaction as determined by cross-correlation analysis. Eur J Neurosci 2: 588–606, 1990.[CrossRef][Web of Science][Medline]
Freeman WJ, Rogers LJ, Holmes MD, Silbergeld DL. Spatial spectral analysis of human electrocorticograms including the alpha and gamma bands. J Neurosci Methods 95: 111–121, 2000.[CrossRef][Web of Science][Medline]
Fries P. A mechanism for cognitive dynamics: neuronal communication through neuronal coherence, Trends Cogn Sci 9: 474–480, 2005.[CrossRef][Web of Science][Medline]
Fries P, Nikoli D, Singer W. The gamma cycle. Trends Neurosci 30: 309–316, 2007.[CrossRef][Web of Science][Medline]
Fries P, Reynolds JH, Rorie AE, Desimone R. Modulation of oscillatory neuronal synchronization by selective visual attention. Science 291: 1560–1563, 2001.
Graham DJ, Field DJ. Sparse coding in the neocortex. In: Evolution of Nervous Systems, edited by Kaas JH, Krubitzer LA. Oxford, UK: Elsevier/Academic Press, 2006, vol. 1, p. 181–187.
Gray CM, McCormick DA. Chattering cells: superficial pyramidal neurons contributing to the generation of synchronous oscillations in the visual cortex. Science 274: 109–113, 1996.
Gunji A, Ishii R, Chau W, Kakigi R, Pantev C. Rhythmic brain activities related to singing in humans. Neuroimage 34: 426–434, 2007.[CrossRef][Web of Science][Medline]
Harris FJ. On the use of Windows for harmonic analysis with the discrete Fourier transform. Proc IEEE 66: 51–83, 1978.[CrossRef]
Hutcheon B, Yarom Y. Resonance, oscillation and the intrinsic frequency preference of neurons. Trends Neurosci 23: 216–222, 2000.[CrossRef][Web of Science][Medline]
Isomura Y, Sirota A, Özen S, Montgomery S, Mizuseki K, Henze DA, Buzsáki G. Integration and segregation of activity in entorhinal-hippocampal subregions by neocortical slow oscillations. Neuron 52: 871–882, 2006.[CrossRef][Web of Science][Medline]
Jacobs J, Kahana MJ, Ekstrom AD, Fried I. Brain oscillations control timing of single-neuron activity in humans. J Neurosci 27: 3839–3844, 2007.
Jarvis MR, Mitra PP. Sampling properties of the spectrum and coherency of sequences of action potentials. Neural Comput 13: 717–749, 2001.[CrossRef][Web of Science][Medline]
Klimesch W. EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res Rev 29: 169–195, 1999.[CrossRef][Medline]
Klimesch W, Sauseng P, Hanslmayr S. EEG alpha oscillations: the inhibition-timing hypothesis. Brain Res Rev 53: 63–88, 2007.[CrossRef][Medline]
König P. A method for the quantification of synchrony and oscillatory properties of neuronal activity. J Neurosci Methods 54: 31–37, 1994.[CrossRef][Web of Science][Medline]
Kreiman G, Hung CP, Kraskov A, Quiroga RQ, Poggio T, DiCarlo JJ. Object selectivity of local field potentials and spikes in the macaque inferior temporal cortex. Neuron 49: 433–445, 2006.[CrossRef][Web of Science][Medline]
Lakatos P, Chen C, O'Connell M, Mills A, Schroeder C. Neuronal oscillations and multisensory interaction in primary auditory cortex. Neuron 53: 279–292, 2007.[CrossRef][Web of Science][Medline]
Lutz A, Greischar LL, Rawlings NB, Ricard M, Davidson RJ. Long-term meditators self-induce high-amplitude gamma synchrony during mental practice. Proc Natl Acad Sci USA 101: 16369–16373, 2004.
Martinovic J, Gruber T, Müller MM. Induced gamma band responses predict recognition delays during object identification. J Cogn Neurosci 19: 921–934, 2007.[CrossRef][Web of Science][Medline]
Maurer AP, McNaughton BL. Network and intrinsic cellular mechanisms underlying theta phase precession of hippocampal neurons. Trends Neurosci 30: 325–333, 2007.[CrossRef][Web of Science][Medline]
Meador KJ, Ray PG, Echauz JR, Loring DW, Vachtsevanos GJ. Gamma coherence and conscious perception. Neurology 59: 847–854, 2002.
Melloni L, Molina C, Pena M, Torres D, Singer W, Rodriguez E. Synchronization of neural activity across cortical areas correlates with conscious perception. J Neurosci 27: 2858–2865, 2007.
Mitra PP, Pesaran B. Analysis of dynamic brain imaging data. Biophys J 76: 691–708, 1999.[Web of Science][Medline]
Montgomery SM, Buzsáki G. Gamma oscillations dynamically couple hippocampal CA3 and CA1 regions during memory task performance. Proc Natl Acad Sci USA 104: 14495–14500, 2007.
Murata T, Takahashi T, Hamada T, Omori M, Kosaka H, Yoshida H, Wada Y. Individual trait anxiety levels characterizing the properties of zen meditation. Neuropsychobiology 50: 189–194, 2004.[CrossRef][Web of Science][Medline]
Murean RC, Savin C. Resonance or integration? Self-sustained dynamics and excitability of neural microcircuits. J Neurophysiol 97: 1911–1930, 2007.
Nunez PL, Srinivasan R. Electric Fields of the Brain: The Neurophysics of EEG (2nd ed.). New York: Oxford Univ. Press, 2005.
Oppenheim AV, Schafer RW. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice Hall, 1999.
Percival DB, Walden AT. Spectral Analysis for Physical Applications. Cambridge, UK: Cambridge Univ. Press, 1993.
Pesaran B, Pezaris JS, Sahani M, Mitra PP, Andersen RA. Temporal structure in neuronal activity during working memory in macaque parietal cortex. Nat Neurosci 5: 805–811, 2002.[CrossRef][Web of Science][Medline]
Puil E, Meiri H, Yarom Y. Resonant behavior and frequency preferences of thalamic neurons. J Neurophysiol 71: 575–582, 1994.
Rodriguez E, Lachaux JP, Martinerie J, Renault B, Varela FJ. Perception's shadow: long-distance synchronization of human brain activity. Nature 397: 430–433, 1999.[CrossRef][Medline]
Roopun AK, Middleton SJ, Cunningham MO, LeBeau FEN, Bibbig A, Whittington MA, Traub RD. A beta2-frequency (20–30 Hz) oscillation in nonsynaptic networks of somatosensory cortex. Proc Natl Acad Sci USA 103: 15646–15650, 2006.
Samonds JM, Bonds AB. Gamma oscillation maintains stimulus structure-dependent synchronization in cat visual cortex. J Neurophysiol 93: 223–236, 2005.
Schaefer AT, Angelo K, Spors H, Margrie TW. Neuronal oscillations enhance stimulus discrimination by ensuring action potential precision. PLoS Biol 4: e163, 2006.[CrossRef][Medline]
Schoffelen JM, Osstenveld R, Fries P. Neuronal coherence as a mechanism of effective corticospinal interaction. Science 308: 111–113, 2005.
Singer W. Neuronal synchrony: a versatile code for the definition of relations? Neuron 24: 49–65, 1999.[CrossRef][Web of Science][Medline]
Steriade M. Grouping of brain rhythms in corticothalamic systems. Neuroscience 137: 1087–1106, 2006.[CrossRef][Web of Science][Medline]
Steriade M, McCormick DA, Sejnowski TJ. Thalamocortical oscillations in the sleeping and aroused brain. Science 262: 679–685, 1993.
Tallon-Baudry C, Bertrand O, Delpuech C, Pernier J. Oscillatory gamma-band (30–70 Hz) activity induced by a visual search task in humans. J Neurosci 17: 722–734, 1997.
Tallon-Baudry C, Mandon S, Freiwald WA, Kreiter AK. Oscillatory synchrony in the monkey temporal lobe correlates with performance in a visual short-term memory task. Cereb Cortex 14: 713–720, 2004.
Timofeev I, Grenier F, Bazhenov M, Sejnowski TJ, Steriade M. Origin of slow cortical oscillations in deafferented cortical slabs. Cereb Cortex 10: 1185–1199, 2000.
Uhlhaas PJ, Singer W. Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron 52: 155–168, 2006.[CrossRef][Web of Science][Medline]
Vidal JR, Chaumon M, O'Regan JK, Tallon-Baudry C. Visual grouping and the focusing of attention induce gamma-band oscillations at different frequencies in human magnetoencephalogram signals. J Cogn Neurosci 18: 1850–1862, 2006.[CrossRef][Web of Science][Medline]
Whittington MA, Stanford IM, Colling SB, Jefferys JG, Traub RD. Spatiotemporal patterns of gamma frequency oscillations tetanically induced in the rat hippocampal slice. J Physiol 502: 591–607, 1997.
Whittington MA, Traub RD. Interneuron diversity series: inhibitory interneurons and network oscillations in vitro. Trends Neurosci 26: 676–682, 2003.[CrossRef][Web of Science][Medline]
Womelsdorf T, Schoffelen JM, Oostenveld R, Singer W, Desimone R, Engel AK, Fries P. Modulation of neuronal interactions through neuronal synchronization. Science 316: 1609–1612, 2007.
Zeitler M, Fries P, Gielen S. Assessing neuronal coherence with single-unit, multiunit, and local field potentials. Neural Comput 18: 2256–2281, 2006.[CrossRef][Web of Science][Medline]
This article has been cited by other articles:
|
|
 |
|
  K. Mirpour and H. Esteky State-Dependent Effects of Stimulus Presentation Duration on the Temporal Dynamics of Neural Responses in the Inferotemporal Cortex of Macaque Monkeys J Neurophysiol, September 1, 2009; 102(3): 1790 - 1800. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
TOP
ABSTRACT
INTRODUCTION
