## Abstract

Spectral analysis of neuronal spike trains is an important tool in understanding the characteristics of neuronal activity by providing insights into normal and pathological periodic oscillatory phenomena. However, the refractory period creates high-frequency modulations in spike-train firing rate because any rise in the discharge rate causes a descent in subsequent time bins, leading to multifaceted modifications in the structure of the spectrum. Thus the power spectrum of the spiking activity (autospectrum) displays elevated energy in high frequencies relative to the lower frequencies. The spectral distortion is more dominant in neurons with high firing rates and long refractory periods and can lead to reduced identification of low-frequency oscillations (such as the 5- to 10-Hz burst oscillations typical of Parkinsonian basal ganglia and thalamus). We propose a compensation process that uses shuffling of interspike intervals (ISIs) for reliable identification of oscillations in the entire frequency range. This compensation is further improved by local shuffling, which preserves the slow changes in the discharge rate that may be lost in global shuffling. Cross-spectra of pairs of neurons are similarly distorted regardless of their correlation level. Consequently, identification of low-frequency synchronous oscillations, even for two neurons recorded by a single electrode, is improved by ISI shuffling. The ISI local shuffling is computed with confidence limits that are based on the first-order statistics of the spike trains, thus providing a reliable estimation of auto- and cross-spectra of spike trains and making it an optimal tool for physiological studies of oscillatory neuronal phenomena.

## INTRODUCTION

Spectral analysis is widely used in many science and engineering fields (Miller and Sigvardt 1998; Percival and Walden 1993). Spectral methods are best suited for systems that display rhythmic behavior such as the nervous system, which tends to fire in periodic oscillations in many normal and pathological states (Brown 2003; Gauss and Seifert 2000; Gray 1994; McCormick 2002). One key example is the low frequency modulation of firing activity that arises in Parkinson's disease (PD) (Hutchison et al. 1997; Nini et al. 1995). These neuronal oscillations, which can be found in the globus pallidus (GP), thalamus, and other basal-ganglia-related structures, may be important clues to the understanding of the clinical symptoms of Parkinsonism, especially the rest tremor.

The spike train of a neuron is a point process frequently modeled as a Poisson process. The power spectrum of a Poisson process exhibits equal power in all frequencies. Nevertheless, the distribution of the power spectrum of a typical spike train is not homogeneous. The power at high frequencies (>100 Hz) is close to that of a Poisson process with the same average firing probability (Bair et al. 1994; Franklin and Bair 1995). However, the power in the low frequencies is significantly smaller than expected, and a trough can be seen in the spectrum of the neuron (Fig. 1). This distortion in the spectrum of the neuronal firing makes it difficult to identify any oscillations in the low frequencies. As a result, Parkinsonian pallidal cells that oscillate in the range of 5–10 Hz may not be classified correctly as oscillatory unless the oscillation power is high enough to overcome the distortion. This issue has been neglected in many previous studies (e.g., Raz et al. 2000, 2001), which have only depicted the lower range (<30 Hz) of the spectrum.

The spectral distortion phenomenon has been described previously in cortical area MT (Bair et al. 1994; Franklin and Bair 1995) and in the auditory system (Edwards et al. 1993). It occurs because the Poissonian model does not take the refractory period, which is an important property of the neuron, into consideration. The refractoriness prevents the neuron from firing two successive spikes within a very short interval. Hence a spike train is never a true Poisson process and does not contain two spikes that occur in a range smaller than the refractory period. In this paper, we overcome these problems by applying the ISI shuffling method. We demonstrate the basic process for compensation of the refractory period and for the detection of oscillatory neurons. Then we present several refinements that enable better assessment of the frequency domain properties of spike trains.

## METHODS

### Electrophysiological recordings

Electrophysiological examples were obtained from extracellular recordings made in a vervet monkey (*monkey C, Cercopithecus aethiops*) and a cynomolgus monkey (*monkey R, Macaca fascicularis*). Recordings were made in the GP before and after treatment with 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP), which induced Parkinsonism. Standard recording and analysis methods were used and are detailed elsewhere (Heimer et al. 2002). Briefly, units were detected and isolated in real time using a template matching algorithm (MSD, AlphaOmega Engineering). The isolation quality was evaluated on-line according to the signal-to-noise ratio and the stability of the spike waveform. Recording quality and stability were further evaluated off-line, and neurons were included in this study only if they were well isolated according to their ISI distribution and their firing rate was stable.

### Spike train simulations

The spike train of a simulated neuron was modeled as an orderly point process (i.e., only 1 event can occur in each bin) (Brillinger 1975; Cox and Isham 1980). It was assumed to be quasi-stationary and to satisfy a cross-mixing condition (i.e., differential changes in discharge rate that are widely spaced in time are statistically independent) (see Brillinger 1975). The simulated spike trains were modeled as a renewal process with a refractory period and a constant firing probability thereafter (Bar-Gad et al. 2001a, b). In this model, neurons had a constant firing probability (*p*) for each discrete time bin (Δ*t*) except that after a spike occurred, the neuron entered a refractory period (with a length of *n _{r}* time bins) in which its probability of firing was smaller than the steady-state probability. We used a refractory period represented by the exponential function

*p*=

_{n}*k*

^{(nr}

^{+1−n)}·

*p*;

*n*≤

*n*, where

_{r}*n*is the number of bins after the preceding spike and

*n*is the length of the refractory period. An absolute refractory period is defined by

_{r}*k = 0*leading to zero firing probability for the entire period; a relative refractory period is defined by 0 <

*k*< 1. This model of the refractory period simplifies the simulations and the mathematical analysis and was thus used instead of the more realistic gamma distribution (Kuffler et al. 1957; Stein 1965). For oscillation modeling, we added a sinusoidal modulation of the firing rate with frequency

*f*

_{osc}and amplitude

*p*

_{osc}. In each time bin, the term

*p*

_{osc}· sin(2π

*f*

_{osc}

*n*· Δ

*t*) was added to the firing probability. The sine term produced oscillatory spike trains that contained bursts that started and ended slowly, like those seen in human patients with Parkinsonian tremor (Zirh et al. 1998), as opposed to bursts with a quick start that slowly decay.

To examine the spectral effects of the refractory period on the cross-spectrum, we studied a pair of neurons with a common input of strength *p*_{corr}. To simulate the common input, a spike train representing the hidden common source was generated as described in the preceding text. Then two other spike trains were generated in a similar manner, but in both of them, the probability of having a spike in the *n*th bin was *p*_{corr} higher if a spike had occurred in the spike train representing the common input in the *n*th bin. To simulate two neurons that were recorded from the same electrode leading to incomplete spike identification due to the shadowing effect (Bar-Gad et al. 2001b), we generated both spike trains as described in the preceding text and then deleted simultaneous spikes or spikes that occurred 1 ms apart.

Throughout both the electrophysiological recordings and the simulations, the default time bin (Δ*t*) was 1 ms. The length of the simulation was in the range of the average duration of the electrophysiological recordings (10^{6} bins ≈ 17 min).

### Spectral calculations

The power spectral density (PSD), cross-spectral density (CSD), and the coherence functions were estimated based on Welch's method (Halliday et al. 1995; Welch 1967). In all spectral calculations, the windows used were nonoverlapping Hanning windows with a length of 4,096; in each window, a discrete Fourier transform (DFT) was performed. The sampling frequency for spike detections was 1,000 Hz. Therefore the resolution of all spectral functions was ∼0.25 Hz and the maximal frequency was 500 Hz. All spectral calculations were performed using standard Matlab 6.5 (The MathWorks, Natick, MA) code. Other methods for calculating the power spectral density are available (Percival and Walden 1993) and yield qualitatively identical results.

In all our neurons, simulated and real, the PSDs were flat in the high frequencies (>100 Hz) range and converged to the expected PSD of a Poisson process with the same average firing rate. Therefore we based our initial confidence level—which determines the frequencies in which the neuron tends to fire more than expected by the Poisson process—on the high frequencies range. Because the spectra of all neurons in this study were flat over the last 100–500 Hz, all the figures present the spectrum ≤300Hz, and the initial confidence levels were constructed based on the last 10% of the bins. The means ± SD of the PSD values (in logarithmic scale) in the range of 270–300 Hz were calculated, and assuming a normal distribution (according to the central limit theorem), a confidence level of 99% (normalized to the total number of bins) was constructed. Another possible asymptotic 100α% confidence level for a Poisson process was obtained by (Halliday et al. 1995): log_{10}(*P̂*_{1}) + [norm_inv(α) · log_{10}*e*]/√*L*, where *P̂*_{1} is the estimated firing rate, *L* is the number of subrecords that the data were divided into in the spectral calculation and norm_inv(α) is the inverse of the normal cumulative distribution function at the corresponding probabilities in α. We tested this method on our data and the results were almost identical to those based on the last 10% bins.

We only calculated confidence levels for the absolute value of the complex function of the CSD. Because the CSD was also flat in the high frequencies (100–300 Hz), the confidence level was constructed in a manner similar to the construction of the confidence level of the PSD. In results, we show why this choice of confidence levels is not helpful for both PSD and CSD, and we propose alternative confidence levels.

The coherence function is also commonly used to explore the spectral relations between two processes. It ranges from 0 to 1, and a significant value indicates linear relations between the two processes at a particular frequency. A 100α% confidence level for the coherence function is given by 1 − (1 − α)^{1/(}^{L}^{−1)} where *L* is the number of subrecords the data were divided into when calculating the coherence (Bloomfield 1976; Brillinger 1981; Rigas 1996a). The confidence limits and sensitivity of this and other methods are compared with the ISI shuffling methods in results.

All the functions used for both the simulation and analysis are Matlab 6.5 (Mathworks, Natick, MA) compatible and may be found at http://basalganglia.huji.ac.il/links.htm.

## RESULTS

### Distorted spectrum of spike trains

The spectrum of the spiking activity of GP neurons in both normal (Fig. 1*A*) and Parkinsonian (Fig. 1, *B* and *C*) primates reflects an uneven distribution of power in different frequencies. The power is stereotypically lower in the range of 3–100 Hz, whereas it is stable in the higher (>100 Hz) frequencies. The low-frequency trough can be found in all pallidal neurons, including neurons that tend to fire spikes periodically. As a result, it is often difficult to identify the expected peak in the 4- to 14-Hz range in the Parkinsonian state (Fig. 1, *B* and *C*). The significance of the 10-Hz oscillations of the neurons shown in Fig. 1, *B* and *C,* is not obvious according to our initial *P* = 0.01 confidence levels (that are based on the mean and variability at 270–300 Hz, see methods). Moreover, the figures demonstrate that calculating the confidence levels based on a different frequency range (e.g., 3–30 Hz) (Raz et al. 2000, 2001) may lead to other errors in estimations of the oscillatory activity. The recorded neurons demonstrating the distortion in the spectrum are well isolated, and the refractory period is clearly evident in their ISI histogram (Fig. 1, *D–F*). Indeed, the refractory period is the main source for the spectral distortion. To illustrate the effect of the refractory period, Fig. 2 shows the PSD and ISI histograms of four simulated neurons with a pure Poisson process (Fig. 2, *A* and *E*) versus three examples of simulated neurons with different refractory periods and oscillatory activity (Fig. 2, *B–D, F–H*).

An intuitive explanation for the spectral distortion of spike trains can be given by examining the autocorrelation function that also suffers from the effects of refractoriness. The autocorrelation function of the high-frequency discharge GP cells reveals a short term peak followed by a small trough due to the refractory period (Bar-Gad et al. 2001a). For long offsets, the autocorrelation seems to be flat, but in fact, it contains many small peaks that are similar to the short-term peak but are less dominant. They appear because any rise in the autocorrelation function above its expected value will cause a drop from the expected value in the following time bins due to the refractory period. Therefore the autocorrelation function of a neuronal spike train with a refractory period keeps rising and falling from its expected value very rapidly. The size of these fluctuations in the autocorrelation function apparently decreases with the offset from zero lag because of the accumulation of variance in the timing of the following spikes. Nevertheless, as a result of these small fluctuations, the power of the high frequencies of the neural autocorrelation function is high compared with the low frequencies. The uneven power distribution creates a trough in the low frequencies of the auto-spectrum function that is the Fourier transform of the autocorrelation function. The mathematical analysis of this phenomenon can be found in (Bair et al. 1994; Edwards et al. 1993; Franklin and Bair 1995). A simpler explanation is given in appendix a for an absolute refractory period.

### ISI shuffling

An oscillatory spike train is a nonrenewal process because the oscillations are usually generated by periodic bursts of spikes. Therefore to detect significant oscillations, we consider as a null hypothesis a renewal process with the same distribution of ISIs. Unlike the original spike train that may be oscillatory, the ISIs of the new process do not depend on the previous ISIs or on time. The renewal process is achieved by using the ISI shuffling method (Perkel et al. 1967; Tam et al. 1988). ISI shuffling uses the first-order ISIs (i.e., the time differences between adjacent spikes) as the building blocks of the new spike train. The ISIs are therefore randomly permuted, and a new spike train is generated based on the shuffled ISIs (Fig. 3*A*). The PSD of the new spike train is determined solely by the first-order ISIs of the original spike train. All higher-order effects (i.e., the time difference between 2 spikes that are separated by 1 spike or more) are abolished by the shuffling process. Hence comparison of the original PSD to the shuffled one enables the detection of patterns such as periodic burst oscillations that are generated by higher-order ISIs. Repeating the shuffling process *N* times and averaging the results provides a less noisy estimator of the contribution of the refractory period to the PSD. The PSD of the shuffled train is shown for a simulated (Fig. 4*A*) and for real neurons (*B*). As expected, the PSD of the shuffled train does not contain the peak around 10 Hz. Yet unlike the simulated neurons, the original and the shuffled PSDs of the real neuron do not fully converge in the range of the trough. This difference was found to be correlated with higher power at low frequencies that probably arises from the persistent slow changes in the neural firing rate (Wichmann et al. 2002). We can overcome the effect of slow changes in the firing rate by performing local shuffling (see following text).

### Local shuffling

Because slow changes are apparent in the neural firing rate of real neurons, our initial null hypothesis of a renewal process with the same distribution of ISIs is not appropriate. We therefore modified our null hypothesis into a renewal process with the same local distribution of ISIs. The ISI distribution of this renewal process does not depend on the previous ISIs, but it changes slowly over time. This hypothesis can be examined by using the local shuffling method. In the global shuffling method described in the preceding text, we permute all ISIs recorded throughout the recording period, thus abolishing temporal changes in firing rate. In the local shuffling method, the spike train is divided into segments of length *T*, and the ISIs in each segment are shuffled within that segment. Consequently, any changes in the firing rate that occur over a period longer than *T* will be preserved in the locally shuffled train.

Two additional procedures further improve the outcome of the local shuffling: first, “soft” rather than “hard” division of the spike train: the actual partition is not done every time *T* (Fig. 3*B*) but rather at the time where the closest spike to time *T* has occurred (Fig. 3*C*). This way, there are no additional spikes in the new spike train, and the ISI histogram is identical to the original histogram.

Second, random rather than fixed segment duration: once *T* is set (for soft or hard division) the shuffling within the segment is done randomly, but the spikes that define the borders of the segments of the local shuffling are constant for a given spike train, and will appear in all *N* random repetitions. These constant spikes can lead to weak oscillations of 1/*T* Hz in the mean shuffled PSD. To avoid this, *T* should be defined as a random variable, *T _{r}*, which changes randomly in a small range. In this study, we used

*T*∈ [150,200] ms. Prior to shuffling of any segment, we draw

_{r}*T*out of a uniform distribution

_{r}*U*(150,200). Then the closest spike to

*T*(soft division) is identified to set the limits of the segment, and the actual shuffling is done within that segment.

_{r}An example of the original PSD and the PSD of the locally shuffled spike train is shown in Fig. 5*A.* The convergence was improved as compared with the PSD of the globally shuffled spike train. The improvement emerged not only in the range of 3–50 Hz but also in the range of 0–3 Hz. This is because the local shuffling maintains the high energy of slow changes in the firing rate that appear in these frequencies (Fig. 5*A*, *bottom*).

### Spectral compensation and confidence levels

Here we present a statistical test for rejection of the null hypothesis: the confidence level for determining whether the spike train is oscillatory can be constructed for a compensated spectrum. The goal of the compensation is to overcome the distortion in the spectrum and equalize the level of detection of low-frequency oscillations to that of the high frequencies. Spectral compensation may be achieved by subtracting the shuffled PSD from the original one, or by dividing them. An advantage of compensation by division is that the variance of the compensated term is the same for all frequencies, and therefore the confidence limits are more reliable (see appendix b). Hence the recommended spectral compensation for the refractoriness of the neuron in each frequency ω is (1) where *S*_{org} is the original PSD and *S*_{shuf} is the PSD of the shuffled train (or the mean of the shuffled trains' PSDs, for *N* > 1). Note that if *S*_{org} = *S*_{shuf} the original PSD is merely the outcome of the first-order ISIs resulting in a compensated PSD of 1 for all frequencies. If *S*_{org}(ω) > *S*_{shuf} (ω), then the original spike train has a power in the frequency ω that is beyond the expected power from its first-order ISIs, and the compensation will be >1 for that ω. Therefore a possible confidence level for the compensation term is constructed by its distance from the expected value 1. The distance is determined by the SD (i.e., the *P* value) of the compensated function in the high frequencies (here we again used the 270- to 300-Hz range), assuming a normal distribution. Figure 5*B* demonstrates the compensated function and its confidence level that is based on the high frequencies. Although the compensation for real neurons is less accurate than the compensation for the simulated neurons (not shown), the significance of the low-frequency oscillations—that might be neglected without the shuffling compensation—is often maintained with the shuffling procedure.

Other potential confidence levels that use bootstrap procedure can be constructed based on a χ^{2} assumption (Jarvis and Mitra 2001; Percival and Walden 1993). However, these confidence levels require many repetitions of the shuffling process and therefore are more costly in terms of calculation. To establish the possible differences between the suggested methods, we compare the performance of the ISI shuffling method to Halliday's method (Halliday et al. 1995). Figure 6 demonstrates the percentage of the detection of the oscillations in each method as a function of oscillation strength (*p*_{osc}, see methods). The shuffling method is more efficient in identifying oscillations of smaller amplitude in our simulated data.

### Distortion and compensation of cross-spectrum and the coherence function

The cross-spectrum of two neurons suffers from the same distortion described in the preceding text for the autospectrum even if the neurons are independent. The size of this distortion depends on the refractory period of both cells. The reason is similar to the one already given regarding the PSD but relates to the cross-correlation function: the expected cross-correlation of two independent cells is flat, but again, once the cross-correlation at time *t* is higher than expected, the probability distribution of the value of the following bin changes accordingly and its expected value is lower than the original expected value (see appendix c for two independent neurons with an absolute refractory period). Therefore the cross-correlation function has small and rapid fluctuations in it. This leads to the relatively low power in the low frequencies and therefore to a trough in the low frequencies of the cross-spectrum functions.

The CSD of two simulated oscillatory cells is shown in Fig. 7*A.* These neurons have no common input, but they oscillate in the same frequency; therefore the CSD contains a peak in the common frequency. The CSD of two oscillatory neurons that have a nonoscillatory common input in addition to an oscillatory correlation has the same shape, but the power over the spectrum increases due to the nonoscillatory correlation (Fig. 7*B*). The reduction in the energy of low frequencies in the CSD is very common in pairs of recorded neurons (Fig. 8, *A* and *B*) and can be observed even if they were nonsimultaneously recorded (Fig. 8*A*). A reliable detection of real oscillatory correlations can be based, again, on the ISI shuffling method that was described above. The shuffling is done in both spike trains, and the compensation in each frequency ω is (2) where *C*_{org} is the original CSD and *C*_{shuf} is the CSD of the shuffled trains (or the mean CSDs, for *N* > 1). The significant correlation frequencies are those found above the confidence level constructed by a normal distribution of mean 1 and a SD that is based on the 270- to 300-Hz band (last 10% of the bins). The global compensation in the CSD, unlike PSD compensation, tends to be accurate and less affected by slow changes (Fig. 8, *C* and *D*). However, local shuffling may be performed to achieve better results.

Note that this method can detect nonoscillatory correlations in addition to the oscillatory correlations (Fig. 7, *C* and *D*), but the former are better identified by the conventional cross-correlation (time domain) function, which is more sensitive to correlations of this kind and is more informative about their nature.

Typically, the coherence function is the cross-spectrum normalized by the autospectra. Therefore the coherence function is not influenced by the refractory period, as can be seen in Figs. 7*E* and 8, *E* and *F.* Nonetheless, the simulated pair with the common input has nonoscillatory correlations and its coherence function tends to cross the confidence level in the high frequencies rather than in the low frequencies (except for the 10-Hz oscillatory correlations; Fig. 7*F*). The coherence of two neurons with a very high firing rate will cross the confidence level more vigorously. This results from the refractoriness that reduces the correlations over the low frequencies and is also demonstrated in the compensation term (Fig. 7*D*).

### Distortion and compensation of cross-spectrum of two neurons recorded from the same electrode

When two or more neurons are recorded from a single electrode, spikes that occurred within a short time interval overlap each other and therefore are hard to identify (Lewicki 1998). This identification failure is termed the “shadowing effect” and can produce various artifacts in the auto- and cross-correlation of the recorded neurons (Bar-Gad et al. 2001b). As a result, the cross-spectrum of two spike trains recorded from the same electrode is distorted not only by the refractory period, but from the shadowing effect as well. Figure 9*A* (black curve) illustrates the cross-spectrum between two simulated spike trains that were recorded from the same electrode (see methods): the trough in the low frequencies still exists, but it becomes narrower. In addition, there is a slow decay in the power of the high frequencies. When applying the shuffling method to both spike trains, the compensation term is far from being optimal (Fig. 9, *A* and *B,* gray curve). To obtain better compensation, the characteristics of the shadowing effect must be preserved during the shuffling. First, any simultaneous spikes must be prevented (as well as adjacent spikes; depend on the shadowing duration). Second, longer ISIs should tend to occur together because each time we lose a spike in both spike trains due to the shadowing effect, an ISI that is longer than twice the refractory period (at least one refractory period before the lost spike and at least one refractory period after the lost spike) is produced in each spike train. These long ISIs tend to occur at parallel times in both spike trains and cause the enhanced correlation. ISI shuffling with the shadowing limitations is possible if the shuffling is done step by step and in parallel for both spike trains. In each step, one ISI of one of the neurons is chosen randomly from the ISIs that have not yet been chosen. This ISI is concatenated to the end of the built train (the part of the shuffled train that was already built for that neuron). If the last spike that was added occurs in the same bin or one bin apart (the shadowing duration in our case) from a spike in the other built train, the shadowing effect must take place. We therefore remove the last ISIs out of both built trains, and we choose (of the ISIs left for each neuron) for each built train an ISI that is longer than twice the refractory period. The results of the step-by-step shuffling are shown in Fig. 9, *A* and *B* (red curve), and indicate that step-by-step shuffling should be performed for the study of synchronous oscillations of neurons recorded by the same microelectrode.

## DISCUSSION

This manuscript discusses the distortion in the spectrum of spike trains that occurs due to the refractory period of neurons. The effect of the refractory period has been previously examined with respect to spectral shape (Bair et al. 1994; Edwards et al. 1993; Franklin and Bair 1995). However, these studies do not address the difficulties in identifying oscillatory phenomena created by the spectral shape of real spike trains. We present a method to overcome these difficulties using ISI shuffling of the spike train. The shuffling of the train may be performed in a local manner so that it overcomes the additional distortions induced by slow changes in the neural firing rate. The shuffled train has exactly the same first-order properties of refractoriness as the original spike train, and therefore the spectra of both trains have the same structure. The periodic properties are revealed by division of the original by the shuffled spectra. The division process ensures equal distribution of power across all frequencies, providing reliable ways to establish confidence limits.

One example of the importance of detecting oscillatory neural activity is the study of neural activity in PD, where basal ganglia and motor thalamic cells tend to fire periodically in the frequency of the tremor (Bergman et al. 1994; Filion and Tremblay 1991; Hurtado et al. 1999, 2005; Hutchison et al. 1997; Lenz et al. 1988; Levy et al. 2002; Nini et al. 1995). The method presented here enables reliable detection of auto- and cross-oscillations in the Parkinsonian state. The effect of the refractory period is crucial in neurons with high discharge rate because the distortion in their spectrum is larger in magnitude. As a result, the method is extremely useful when analyzing neurons with tonic high firing rates or neurons that tend to increase their firing rates under certain circumstances such as a stimulus or behavioral task.

The method provides confidence limits based on the first-order statistics of the spike train. Therefore our method is practical for detecting burst oscillations in all frequencies for all kinds of spike trains, regardless of the discharge rate or the length of the refractory period. This naturally includes spike trains of well-isolated neurons where the refractory period is dominant and the distortion is large as well as multi-unit recordings where there is typically little or no effect of the refractory period and thus the spectrum is not distorted.

The shuffled spike train does not contain any structure of the original spike train that is inherent to second and higher order ISIs. All structures that are generated by the first-order ISI are preserved in the shuffled train and therefore are abolished in the compensated function. Usually, periodic oscillations in the spike train of neurons with burst oscillations are caused by higher-order ISIs; therefore the compensation process suggested here is useful for detecting them. However, this is not always correct. For example, in the case of a regular neuron that fires single spikes periodically (Ahissar and Vaadia 1990), the oscillations are caused by the first-order ISIs, and the PSD of the shuffled train will be identical to the original PSD. The shuffling method suggested here is therefore not recommended for these cases.

This paper sheds light on another common phenomenon in spectral studies of neuronal activity—the excess energy in low (1–3 Hz) regimes. These frequencies often “leak” to the frequencies of interest and may mask significant oscillations in the 5- to 10-Hz region. This occurs due to the ongoing slow changes in the firing rate, which do not exist in the globally shuffled spike train. The larger distortion might be surprising because the ISIs are the same in both spike trains; therefore the expectation is to have a similar distortion in both spectra. However, because the firing rate is not constant, the original spike train tends to have segments of short ISIs that contribute to the large distortion and other segments of longer ISIs that are less effective. Because the PSD function is the average of a nonlinear function that is calculated for each segment, a stronger distortion in the nonequally distributed segments is obtained.

Unlike global shuffling, the local shuffling depends on the parameter *T* that is set by the researcher. Naturally, *T* should be smaller than the size of spectral window (4,096 ms in our case). Moreover, *T* should be as small as possible so the changes in the firing rate are better preserved. On the other hand, as *T* gets smaller, the local shuffling loses its power because the number of ISIs in each segment is close to one, and the shuffling becomes useless. The *T* chosen here was ∼200 ms and seemed to be satisfactory for most of our pallidal cells with mean discharge rate ∼50 Hz. Nevertheless, because the changes in the firing rate are probably variable themselves, there is no way to choose *T* so that the compensation in the spectrum will be perfect.

To summarize, the ISI shuffling method is an important tool for compensating for the distortion in the spectrum of neuronal spike trains and therefore provides reliable confidence limits for spectral functions. It enables the detection of periodic oscillations in spike trains that could have been ignored otherwise, and thus it is recommended for the study of neural oscillations in normal and pathological states.

## APPENDIX A: Auto-spectrum of a point-process with an absolute refractory period

Let *x*_{1}, *x*_{2},… be a process with an absolute refractory period of τ ms. That is (A1) Let *R _{x}*(

*i*) be the autocorrelation function of the process.

*A*(

_{x}*i*) is the event in which no spike occurred at the

*i*th bin. Then (A2) where

*A*is the complementary event in which a spike did occur at the

_{x}^{c}(i)*i*th bin.

We obtain (A3) Recall that the events *A _{x}*(

*i*),

*A*(

_{x}*j*) are disjoint for 0< |

*i*−

*j*| ≤ τ, therefore (A4) and the last equation can be extended to (A5) for δ(

*i*) = 1 for i = 0, and 0 otherwise.

The Fourier transform of *R _{x}*(

*i*−

*k*) is

*S*(ω)

_{x}*e*

^{−j}

^{ω}

*k. Therefore the Fourier transform of both sides (if we take R*(

_{x}*i*) for negative indexes), yields (A6) Because we enforced the autocorrelation function to be 0 for a negative index and because the autocorrelation is a symmetric function, the real Fourier transform is twice the real part of what we got (A7) Figure A1 illustrates the results for different firing probabilities

*p*and refractory periods τ.

## appendix B: Variance of the compensation

If we take the compensation in each frequency ω to be (B1) where *S*_{org} is the original PSD and *S*_{shuf} is the PSD of the shuffled train (or the mean PSDs, for *N* > 1), the variance of the compensation term is the same for all frequencies.

Proof: the variance of a quotient of two random variables can be approximated as follows (Mood et al. 1974) (B2) Let *X*(ω) = *S*_{org}(ω) and *Y*(ω) = *S*_{shuf}(ω) for any ω.

Due to the random shuffling, *X*(ω) and *Y*(ω) are independent and therefore cov(*X, Y*) = 0.

The estimate of the PSD has a distribution which is analogous to a χ^{2} distribution in each ω with the same degrees of freedom (see Brillinger 1981, p. 164). Because the ratio of the SD to the mean for such a random variable depends solely on the degrees of freedom, it does not depend on ω (B3) Therefore (B4) and the same is true for *Y* (with const_{Y}).

If *E*(*X*) = *E*(*Y*), we have (B5) and the term is a constant that does not depend on ω.

Therefore the variance of the compensation for all frequencies is constant.

## appendix C: Cross-correlation function between two processes with absolute refractory periods

Let *X* and *Y* be two independent processes with an absolute refractory period of τ ms and a zero mean. The empirical cross-correlation function between them is defined as (C1) The covariance between the empirical cross-correlation function at time *t* and time *t* + *m* is (C2) where *R _{X}*(

*u*),

*R*(

_{Y}*u*) are the autocorrelation functions of

*X*and

*Y,*respectively.

Because the summation is over all values of the autocorrelation functions, we obtain for *T* ≫ *t, m* (C3) Recall that the autocorrelation function of each of the processes is positive at *u*= 0, negative along the refractory period, and then becomes positive again and finally converges to 0 and is symmetric around that value (Bar-Gad et al. 2001a). Therefore the covariance between the cross-correlation function at different time bins is negative if the distance between them is on the same order of magnitude as the refractory period, which yields a periodic structure of the empirical cross-correlation function.

For further details and clarifications, see (Jenkins and Watts 1968; Rigas 1996b).

## GRANTS

This study was partly supported by a Center of Excellence grant administered by the Israel Science Foundation, by the German-Israel Binational Foundation and the BMBF Israel-Germany collaboration in medical research and by a “Fighting against Parkinson” grant administrated by the Netherlands Friends of the Hebrew University. Y. Ritov was supported in part by an ISF grant.

## Acknowledgments

We thank N. Parush, I. Nelken, and R. Kass for helpful discussions, J. A. Goldberg for useful suggestions and for sharing data with us, G. Morris and Y. Prut for comments on earlier versions of this 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 © 2006 by the American Physiological Society