## Abstract

We present a method that estimates the strength of neuronal oscillations at the cellular level, relying on autocorrelation histograms computed on spike trains. The method delivers a number, termed *oscillation score*, that estimates the degree to which a neuron is oscillating in a given frequency band. Moreover, it can also reliably identify the oscillation frequency and strength in the given band, independently of the oscillation in other frequency bands, and thus it can handle superimposed oscillations on multiple scales (*theta*, *alpha*, *beta*, *gamma*, etc.). The method is relatively simple and fast. It can cope with a low number of spikes, converging exponentially fast with the number of spikes, to a stable estimation of the oscillation strength. It thus lends itself to the analysis of spike-sorted single-unit activity from electrophysiological recordings. We show that the method performs well on experimental data recorded from cat visual cortex and also compares favorably to other methods. In addition, we provide a measure, termed *confidence score*, that determines the stability of the oscillation score estimate over trials.

## INTRODUCTION

Neuronal oscillations are believed to play an important role in brain function. They are expressed during both sleep (Buzsáki 2006; 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 (Mureşan 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. 1*A*), 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. 1*C*). 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.

Here, we propose a novel method for estimating the strength of oscillations that relies on ACHs. The method combines the advantages of both spectral methods and correlogram-based techniques because it relies on the estimation of the frequency spectrum from ACHs. It is simple, fast, and works well with noisy ACHs. Thus it can be used when the number of available experimental trials is small and even with single units. The method first uses a fast Gaussian kernel to smooth the ACH and to remove high-frequency noise. Then, a slow Gaussian kernel is applied, for the sole purpose of detecting the boundaries of the central peak. Using this information, the central peak is efficiently removed from the buffer containing the fast smoothing, which is then subjected to FFT. Eventually, the *oscillation score* is computed as the ratio between the highest-frequency magnitude within the band of interest and the baseline (average) magnitude of the spectrum. The method is robust against noisy ACHs, it minimizes the false detection of oscillations, and allows one to quantify the strength of oscillations in any frequency band (e.g., *theta*, *alpha*, *beta*, *gamma*). We also propose a measure to estimate the degree of confidence with which the oscillation score can be computed from a given ACH.

## METHODS

### Experimental procedures, stimulation, and recording

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% N_{2}O and 30% O_{2} 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 CO_{2} 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 μm^{2}, impedance 0.3–0.5 MΩ at 1,000 Hz, intercontact distance 200 μm). Signals were amplified ×10,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, *f*_{min} and *f*_{max}, 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 *f _{c}* the frequency of the correlogram, given by

*f*= 1/

_{c}*b*(1,000 Hz in this study; note that

*f*can be different from the sampling frequency of the spikes, depending on how one chooses to bin the ACH). Let

_{c}*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*= 2

*w*; the full symmetric ACH including the center bin at lag 0 has 2

*w*+ 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*×

*b*should be sufficiently large to fit at least a few cycles of the slowest frequency of interest,

*f*

_{min}. We required at least three cycles on each flank of the ACH (three-peaks rule), which, e.g., for

*f*

_{min}= 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) where ⌊ ⌋ represents the rounding function to the nearest smaller integer (or floor function).

For example, consider an interval of frequencies in the *beta-*/*gamma-band* (*f*_{min} = 20 Hz; *f*_{max} = 40 Hz) and a bin *b* = 1 ms (*f _{c}* = 1 kHz). The minimum

*w*according to the three-peaks rule should be 150 bins, whereas according to the criterion of having a minimum FFT precision of 2 Hz per bin, it should be 250 ms (

*W*= 500 ms). The formula in

*Eq. 1*takes the maximum closest estimate that is a power of 2. It follows that the ACH should be computed using time lags from −256 to +256 ms (

*w*= 256 ms). The algorithm will use the bins from −256 to +255, that is,

*W*= 512 bins to compute the FFT in step 4.

##### 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) where σ is the SD of the Gaussian kernel.

An example ACH smoothed by such a Gaussian kernel is shown in Fig. 2*B*. 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 *f _{c}* = 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. 2

*C*). 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

*f*

_{max}. For better precision of the smoothing, the cutoff frequency (attenuated by the smoothing with −3 dB) is chosen at 1.5 ×

*f*

_{max}. Thus σ

_{fast}is given by (in bins) (3) For example, for a correlogram frequency

*f*= 1,000 Hz and a maximum frequency of interest

_{c}*f*

_{max}= 100 Hz, the smoothing should be made with a Gaussian kernel having σ

_{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. 1*A*) 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. 1*A*).

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. 3*B*), 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. 3*C*, *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. 3*C*, *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. 3*C*, *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 *f*_{min}, then the algorithm would not detect false slow oscillations (because in step 5, it looks only in the band *f*_{min} − *f*_{max}). 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 *f*_{min} (4)

For example, for a correlogram frequency of 1 kHz, and a frequency band in the *beta* / *gamma* range (*f*_{min} = 20 Hz, *f*_{max} = 40 Hz), the ACH will be smoothed with a Gaussian having σ_{slow} = 8.93 ms. This would create a smoothed central peak (Fig. 3*D*) with an effective width of ≥6 × σ_{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. 3*D*, 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. 3*D*, *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) where *slope*(*t*) is the slope of the slow-smoothed ACH at lag *t* (lag 0 corresponds to the central peak), *S _{slow}* represents the ACH smoothed with the slow Gaussian kernel, ε 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*/

*S*

_{slow}(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) where *t _{left}* is the time lag (on the left of the central peak) where the slope of the slow-smoothed ACH is ≤10° (Fig. 3

*D*,

*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) *t _{left}* and

*t*of the central slow-smoothed peak (Fig. 3

_{right}*D*,

*inset*), we go back to the original, fast-smoothed ACH (that was computed in step 2). In the latter, the ACH value at

*t*is copied to all other time lags spanning the identified range (

_{left}*t*…

_{left}*t*), thus overwriting the real central peak (Fig. 3

_{right}*E*). Depending on the relative size of the side troughs and first satellite peak, and on

*f*

_{min}, 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 (*S _{slow}*) 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

*f*

_{min}(Fig. 3

*F*). Moreover, when the ACH displays strong oscillatory components, the slope of the slow

*S*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,

_{slow}*D*and

*E*). When the ACH displays only troughs due to refractoriness, with relatively weak satellite peaks,

*S*will reach the slope threshold later (further away from the central peak), leading to an efficient removal of the central peak (Fig. 3

_{slow}*I*), 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. 2*A*) and we apply an FFT. The resulting buffer is of size *w* and contains the magnitudes of the frequencies from 0 to *f _{c}* /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. 3*B*).

##### 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, *M _{peak}*, is identified from the power spectrum (Fig. 3

*F*). The estimated oscillation frequency

*f*is the frequency with magnitude

_{osc}*M*. It is important to note that the method searches for the peak magnitude only within the band of interest

_{peak}*f*

_{min}−

*f*

_{max}. Recall that this is an essential aspect because removal of the central peak in step 3 can artificially introduce frequencies lower than

*f*

_{min}(Fig. 3,

*E*and

*F*). By looking only to frequencies larger than

*f*

_{min}, the method is safe from detecting spurious, artificial magnitude peaks.

Next, the average magnitude *M _{avg}* of the spectrum is computed by integrating the whole frequency spectrum and taking its average (7) where

*Magnitude*(

*f*) is the estimated magnitude of frequency

*f*in the FFT-computed spectrum.

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 *O _{S}* is then given by (8)

Finally, the estimated frequency of oscillation, *f _{osc}* (Fig. 3

*F*) is given by the frequency having the highest magnitude

*M*in the band of interest. The strength of the oscillation is given by the oscillation score

_{peak}*O*. 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.

_{S}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. 3*J*), 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. 3*K*). 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. 7*A*). 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, *M _{avg}*, 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.

We used two measures to quantify the degree of confidence of the oscillation score. Let us assume that there are *N* available trials (when only one trial is available, the confidence cannot be estimated). We can then estimate the SD of the oscillation scores, computed on individual trials *i*, by (9) where *σ _{OSt}* is the estimate of SD for the oscillation score computed on individual trials,

*O*is the oscillation score computed on trial

_{Si}*i*, and

*Ō*is the average over

_{St}*O*.

_{Si}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 *C _{v}*, defined as (10) The

*C*estimates the ratio between the SD and the mean, and thus it is an unbounded measure. To get a measure between 0 (no confidence) and 1 (high confidence) we defined a measure that we called

_{v}*confidence score*(

*C*) as (11)

_{S}For small values of the SD relative to the mean, the *C _{S}* gives a value close to 1, whereas for very poor data, with highly fluctuating oscillation scores on individual trials, the

*C*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.

_{S}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, *f _{osc}*, with the peak magnitude within the band of interest (Fig. 3

*F*), 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.

### Spike-field coherence measure

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

We investigated the performance of the described method on neuronal spike trains, both recorded from cat visual cortex, and simulated. We first estimated the ability of the method to quantify the strength of oscillations in different frequency bands. Then, we worked out the proper interpretation of the oscillation score, which depended on the investigated frequency band. Furthermore, we studied how the oscillation score correlated with the oscillation of local-field potentials (LFPs) when the spike-field coherence was high, thus validating the oscillation score as a meaningful measure. Finally, we investigated how the confidence of the score related to the stationary expression of neuronal oscillations in different experimental trials and also how it correlated to the number of available spikes.

### 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.

Our method investigates a given frequency range, defined by the boundary limits *f*_{min} and *f*_{max}. We strongly recommend that the choice for the boundary frequencies closely matches the boundaries of relevant biological bands (*theta*, *alpha*, *beta*, *gamma*), for the following reason: To safely remove the central peak, without contaminating the frequency spectrum in the band of interest, the method uses a smoothing procedure that is slower than the period of the slowest frequency of interest, *f*_{min} (see *Eq. 4*). If *f*_{min} ≪ *f*_{max}, and if the ACH exhibits robust high-frequency oscillations, then the algorithm for removing the central peak (see methods, *Computation of the oscillation score*, step 3) might also remove legitimate satellite peaks (because it tries to cut a portion around the central peak that is larger than the period of *f*_{min}). For example, choosing *f*_{min} = 5 Hz and *f*_{max} = 100 Hz determines an effective cut around the central peak of −100 to +100 ms, thus removing satellite peaks (<50 ms) for potential *beta* and *gamma* oscillations. Had the frequency boundaries been chosen properly, spanning only one biological band (*theta*), *f*_{min} = 5 Hz and *f*_{max} = 8 Hz, the cutting would not have interfered with the highest possible frequency in the band: 8 Hz (which has a period of oscillation of 125 ms > 100 ms). We have determined that for each biological band, the method safely removes the central peak, without interfering with legitimate satellite peaks of the oscillations in the given band. We next defined oscillation scores for each specific band—*theta*, *alpha*, *beta*, and *gamma*—and named the scores accordingly: *1*) *theta score*, *2*) *alpha score*, *3*) *beta score*, and *4*) *gamma score*, respectively.

##### THE THETA SCORE.

The *theta score* is defined for the *theta* oscillation band (Maurer and McNaughton 2007), with *f*_{min} = 4 Hz and *f*_{max} = 8 Hz. In Fig. 4*A*, a recording site (multiunit), stimulated with a moving bar, showed robust oscillations in the *theta* range, with *f _{osc}* = 5 Hz (Fig. 4

*A*,

*bottom*). Note the high magnitude of the

*theta*oscillation in the frequency spectrum, relative to the baseline magnitude (Fig. 4

*A*,

*bottom, inset*), inducing a very high

*theta score*of 60.7 (Fig. 4

*A*,

*bottom*). Also note that removal of the central peak of the ACH (Fig. 4

*A*,

*top*) introduces an artificial frequency of 2 Hz (Fig. 4

*A*,

*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 *f*_{min} = 8 Hz and *f*_{max} = 12 Hz. In Fig. 4*B*, 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. 4*B*, *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 *f*_{min} = 12 Hz and *f*_{max} = 20 Hz) and a *beta-high score* (with *f*_{min} = 20 Hz and *f*_{max} = 30 Hz), as two separate measures for the *beta* range.

In Fig. 4*C* 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. 4*C*, *top*), the *beta-high score* reaches a value of only 19.1 (Fig. 4*C*, *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 *f*_{min} and *f*_{max}. Removal of the central peak will not affect frequencies between 20 and 50 Hz, if *f*_{min} 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. 5*B*, thick dotted line).

To estimate the range of oscillation scores for strong oscillations, we needed an objective quantification of oscillation strength. Therefore, we simulated oscillating spike trains (see methods and appendix) in which we could manipulate the strength of oscillations. We first produced spike trains without any oscillations, computed the distribution of oscillation scores (Fig. 5*A*, Simulation No Osc.), and then we identified the 95th percentile of each distribution. This provided an absolute threshold for the nonoscillating state (Fig. 5*B*, thick black line, delimiting the dark gray area). Note that the distributions for the simulated data matched reasonably well those for experimentally recorded data. In the latter case, the nonoscillating regime scores were slightly higher because, as mentioned previously, in some of the 7,501 ACHs, very weak oscillatory activity could not be subjectively distinguished from nonoscillating cases.

In the next step, we generated simulated data, at the other extreme, exhibiting very strong oscillations (Fig. A1*G*). Again, we computed the distribution of oscillation scores (Fig. 5*A*, 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. 5*B*, 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. 5*B*, Transition). Here, there was a continuum of monotonically increasing oscillation scores, from weak to strong oscillations.

We emphasize that the distributions in Fig. 5*A* 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. 6*B*, revealed a bimodal distribution, with a group of nonoscillating or weakly oscillating responses (*O _{S}* < 10), and a separate group of strongly oscillating responses (

*O*> 10). Some units did not exhibit oscillations for any of the 14 stimulation conditions.

_{S}To examine the relationship between the LFP and spike oscillations, we first computed the oscillation scores for units and for the LFPs recorded from the same electrode as the corresponding unit. We found that the oscillation scores for LFPs and single units were only weakly correlated (*r* = 0.21; *r*^{2} = 0.04). One reason might have been the bimodal distribution of spike oscillation scores. However, even if only the group exhibiting strong oscillations was considered, the correlation between spike and LFP oscillation scores remained low (Fig. 6, *B* and *C*, portion on the right of the dotted line). Note that the LFP oscillation scores were all high.

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; *r*^{2} = 0.85; Fig. 6*D*). When coherence was computed with a second method, based on multitaper spectral estimates (Pesaran et al. 2002), similar results were obtained (*r* = 0.86; *r*^{2} = 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 *f _{osc}* = 19.5 Hz (

*beta*band) and an oscillation score

*O*= 10.4, indicating relatively strong oscillations (see also Fig. 5

_{S}*B*). However, the shape of the ACH does not appear reliable. Indeed, the confidence score was only

*C*= 0.384, indicating that one should not trust the oscillation score

_{S}*O*. Figure 7

_{S}*B*, by contrast, depicts a single unit with an oscillation frequency

*f*= 29.3 Hz, having an oscillation score

_{osc}*O*= 15.7 (very strong oscillations; see also Fig. 5

_{S}*B*) and a high confidence score

*C*= 0.85.

_{S}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 *C _{S}* depends on the number of spikes that were used to compute each of the 1,204 ACHs (86 units × 14 stimulation conditions). Figure 7

*C*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

*C*= 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

_{S}*C*<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

_{S}*C*<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.

_{S}### 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.

##### GENERALIZED GABOR FITTING.

The GGF method, described in König (1994) 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. 1*A*; 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. 8*B* 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. 8*C*; 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. 8*D*).

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. 8*A*, 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. 8*C*. 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. 8*A*, 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. 8*A*, 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. 8*A*, 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. 8*A*, 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. 8*A*, 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. 8*A*, 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. 8*A* (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. 8*A*, 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. 8*A*, OS, 10 Hz). The scores obtained were still in the safe range of the very weak/nonoscillating regime (see Fig. 5*B* 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. 7*A* 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

Neuronal oscillations and their synchronization are receiving increasingly greater attention because there is now ample evidence that they might play important functional roles in the brain. At a mechanistic level, neuronal oscillations are likely to originate from two distinct processes. The first involves the interaction between excitatory and inhibitory cells and has been studied extensively, both experimentally (Whittington and Traub 2003) and theoretically (Brunel and Hakim 1999). The second relies on preferential locking of neurons, induced by their individual frequency preferences (Hutcheon and Yarom 2000; Mureşan 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. 8*A*, 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 *f*_{min} and *f*_{max}—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 (*f*_{min}). 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

### Software package

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

Variables:

fmin—the lower bound for the frequency band of interest;

fmax—the lower bound for the frequency band of interest;

fc—correlogram frequency in Hz;

w—the width of one flank of the ACH;

W—the width of the buffers = 2*w;

ACH—the original autocorrelation histogram;

sgm

_{fast}—the standard deviation of the fast Gaussian kernel;S

_{fast}—fast smoothed ACH;sgm

_{slow}—the standard deviation of the slow Gaussian kernel;S

_{slow}—slow smoothed ACH;slope—the slope of S

_{slow}at a given offset;t

_{left}—the index on the left, for cutting the central peak;S-P

_{fast}—fast smoothed, peakless ACHS

_{fft}—the spectrum of S-P_{fast};M

_{peak}—the peak magnitude in the band of interest;fosc—the oscillation frequency;

M

_{avg}—the baseline magnitude of the spectrum;Os—the oscillation score;

Ost—array holding the oscillation scores for each individual trial;

sgm_Ost—the standard deviation computed on the array Ost;

mean_Ost—mean of the Ost array;

Cv—coefficient of variation;

Cs—confidence score;

//1. Computation of the auto-correlation histogram (ACH)

w := pow(2,floor(max(log2 (3*fc/fmin), log2 (fc/4)))+1);

W := 2*w;

ACH := Compute_Autocorrelation(neuron,condition,trials);

//2. Smoothing

sgm

_{fast}:= min(2,134/(1.5*fmax))*fc/1000.0;S

_{fast}:= Gaussian_Smooth(ACH,sgm_{fast});

//3. Removing the central peak of the ACE

sgm

_{slow}:= 2*134/(1.5*fmin)*fc/1000.0;S

_{slow}:= Gaussian Smooth (ACH,sgm_{slow});for i:=0 down to −w+1 do //identify the limits of the peak

slope:=(S

_{slow}[i]−S_{slow}[i−1])*W/S_{slow}[0];if (slope <= tan (PI*10/180)) then

t

_{left}:= i;break;

endif;

endfor;

if (i=0) then t

_{left}:= 0; //index not found//Effective removal of the peak:

S-P

_{fast}:= S_{fast};for i:=t

_{left}+1 to abs (t_{left})−1 do //the ACH is symmetric;S-P

_{fast}[i] := S_{fast}[t_{left}]; //overwrite the peak;

endfor;

//4. Applying the FFT

S

_{FFT}:= FFT(S-P_{fast}); //the spectrum of the peakless fast-smoothed ACH;

//5. Estimation of the oscillation score

fosc := fmin;

M

_{peak}:= S_{FFT}[HzToBins(fmin)]; //assume the highest magnitude is at fminfor i:=HzToBins(fmin) to HzToBins(fmax) do //identify fosc

if (S

_{FFT}[i] > M_{peak}) thenM

_{peak}:= S_{FFT}[i];fosc := BinsToHz(i);

endif;

endfor;

M

_{avg}:= Compute_Average(S_{FFT});Os := M

_{peak}/M_{avg}; //the final oscillation score

//Estimation of the confidence score

for i:=1 to n_trials do Ost[i] := Oscillation_Score(i);

sgm_Ost := STDEV(Ost);

mean_Ost := Compute_Average(Ost);

Cv := sgm_Ost / mean_Ost;

Cs := 1/(1+Cv);

### 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 *p _{s}*(

*t*) (Fig. A1

*A*) and a finer temporal process that controls exact spike timing (refractoriness, bursting, etc.; Fig. A1

*C*).

The spike discharge probability function *p _{s}*(

*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*. A

*1*) of two discharge probabilities

*p*(

_{o}*t*) and

*p*(

_{b}*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. A1

*A*) with a random process

*f*(

_{o}*t*) that takes into account the past values of the frequency and thus has memory and offers stability (

*Eq*. A

*4*). 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

*p*(

_{o}*t*), and thus oscillation preference, becomes stable in a certain range (Fig. A1

*B*,

*f*). The background activity

_{o}*p*

_{b}*(t*) is generated in a similar fashion but, in this case, the stability of the process is lowered and oscillation frequency

*f*(

_{b}*t*) varies in a much broader frequency interval (Fig. A1

*B*,

*f*). The spikes are generated only in regions where

_{b}*p*(

_{s}*t*) is >0 (Fig. A1

*A*), and are concentrated toward the peaks of

*p*(

_{s}*t*) (where the probability of producing spikes is higher; Fig. A1

*C*). The number of generated spikes depends on the positive area of

*p*(

_{s}*t*). With the offset (Fig. A1

*A*) we can control how well the spikes are aligned to the desired oscillation [the peaks of

*p*(

_{s}*t*) approach a probability of 1], whereas with the area of the positive

*p*(

_{s}*t*) we can influence how many spikes are generated overall and thus the firing rate. An interesting case is when

*p*(

_{s}*t*) is <0 most of the time, having high positive values for only short periods of time [when the amplitude of

*p*(

_{s}*t*) is very high while the offset is also very large; see Fig. A1

*A*]. This situation results in very precise firing patterns. In the end, with

*p*(

_{s}*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. A1*D*). The simplest way to model the burst probability, *p _{burst}*(

*t*), is to set it as a fraction of

*p*(

_{s}*t*) (

*Eq*. A

*5*). Once a discharge is initiated this probability function is used to determine whether a tonic spike or a burst is generated (Fig. A1

*C*).

To obtain realistic spike trains, we first computed the distribution of burst lengths (Fig. A1*D*) and intraburst ISIs (Fig. A1*E*) 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. A1*D*). The intervals between the spikes in the burst are set based on the measured intraburst ISI distribution (Fig. A1*E*). Immediately after a tonic spike, or a burst, a refractory period (*r _{tonic}* and

*r*, respectively, in Fig. A1

_{burst}*C*) 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. A1

*G*), which is lost if

*k*is set to 0 (A5)

In Fig. A1*F* we show the ACH of an experimentally recorded neuron from cat area 17, whereas in Fig. A1*G* 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. A1*F*: the oscillation frequency (∼29 Hz; see also Fig. A1*B*), 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

This work was supported by the Hertie Foundation, three grants of the Romanian Government: Human Resources Program RP-5/2007 Contract 1/01.10.2007, Ideas Program ID_48/2007 Contract 204/01.10.2007, both financed by Ministerul Educaţiei Cercetǎrii şi Tineretului (MECT) / Unitatea Executivǎ pentru Finaţarea Învǎţǎmântului Superior şi a Cercetǎrii Ştiinţifice Universitaire; and Partnerships Program Neurobot Contract 11039/18.09.2007, financed by MECT/Autoritatea Naţionalǎ pentru Cercetare Ştiinţificǎ and a Deutsche Forschungsgemeinschaft grant NI 708/2 - 1.

## Acknowledgments

The authors thank R. V. Florian and S. P. Paş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

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society