| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
INNOVATIVE METHODOLOGY
1Center for Neural Computation, 2Department of Physiology—Hadassah Medical School; 3Department of Statistics, 4 Roland Center for Neurodegenerative Diseases, The Hebrew University, Jerusalem; and 5 Gonda Brain Research Center, Bar-Ilan University, Ramat-Gan, Israel
Submitted 18 January 2005; accepted in final form 2 January 2006
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; 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.
|
; 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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 nr 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 pn = k(nr+1–n) · p; n 
nr, where n is the number of bins after the preceding spike and nr is the length of the refractory period. An absolute refractory period is defined by 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 fosc and amplitude posc. In each time bin, the term posc · sin(2
foscn · 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 pcorr. 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 nth bin was pcorr higher if a spike had occurred in the spike train representing the common input in the nth 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 (106 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): log10(
1) + [norm_inv() · log10e]/
L, where 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
The spectrum of the spiking activity of GP neurons in both normal (Fig. 1A) 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).
|
). 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. 3A). 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. 4A) 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).
|
|
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. 3B) but rather at the time where the closest spike to time T has occurred (Fig. 3C). 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, Tr, which changes randomly in a small range. In this study, we used Tr 
[150,200] ms. Prior to shuffling of any segment, we draw Tr out of a uniform distribution U(150,200). Then the closest spike to Tr (soft division) is identified to set the limits of the segment, and the actual shuffling is done within that segment.
An example of the original PSD and the PSD of the locally shuffled spike train is shown in Fig. 5A. 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. 5A, bottom).
|
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) |
) > Sshuf (), 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 5B 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 (posc, see METHODS). The shuffling method is more efficient in identifying oscillations of smaller amplitude in our simulated data.
|
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. 7A. 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. 7B). 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. 8A). 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) |
|
|
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. 7E 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. 7F). 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. 7D).
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 9A (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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
ms. That is
 | (A1) |
 | (A2) |
 | (A3) |
, therefore
 | (A4) |
 | (A5) |
(i) = 1 for i = 0, and 0 otherwise.
The Fourier transform of Rx(i – k) is Sx()e–jk. Therefore the Fourier transform of both sides (if we take Rx(i) for negative indexes), yields
 | (A6) |
 | (A7) |
.
|
|
appendix B: Variance of the compensation |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
to be
 | (B1) |
Proof: the variance of a quotient of two random variables can be approximated as follows (Mood et al. 1974)
 | (B2) |
) = Sorg() and Y() = Sshuf() 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) |
 | (B4) |
 | (B5) |
. Therefore the variance of the compensation for all frequencies is constant.
|
appendix C: Cross-correlation function between two processes with absolute refractory periods |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
ms and a zero mean. The empirical cross-correlation function between them is defined as
 | (C1) |
 | (C2) |
Because the summation is over all values of the autocorrelation functions, we obtain for T >> t, m
 | (C3) |
). 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: M. Rivlin-Etzion (E-mail: michriv2{at}alice.nc.huji.ac.il)
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: Auto-spectrum of...appendix B: Variance of...appendix C: Cross-correlation...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Bair W, Koch C, Newsome W, and Britten K. Power spectrum analysis of bursting cells in area MT in the behaving monkey. J Neurosci 14: 2870–2892, 1994.[Abstract]
Bar-Gad I, Ritov Y, and Bergman H. The neuronal refractory period causes a short-term peak in the autocorrelation function. J Neurosci Methods 104: 155–163, 2001a.[CrossRef][ISI][Medline]
Bar-Gad I, Ritov Y, Vaadia E, and Bergman H. Failure in identification of overlapping spikes from multiple neuron activity causes artificial correlations. J Neurosci Methods 107: 1–13, 2001b.[CrossRef][ISI][Medline]
Bergman H, Wichmann T, Karmon B, and DeLong MR. The primate subthalamic nucleus. II. Neuronal activity in the MPTP model of Parkinsonism. J Neurophysiol 72: 507–520, 1994.
Bloomfield P. Fourier Analysis of Time Series: An Introduction. New York: Wiley, 1976
Brillinger DR. The identification of point processes systems. Ann Probab 3: 909–929, 1975.
Brillinger DR. Time Series: Data Analysis and Theory. San Francisco: Holden-Day, 1981.
Brown P. Oscillatory nature of human basal ganglia activity: relationship to the pathophysiology of Parkinson's disease. Mov Disord 18: 357–363, 2003.[CrossRef][ISI][Medline]
Cox DR and Isham V. Point Processes. London: Chapman and Hall, 1980.
Edwards B, Wakefield G, and Powers NL. The spectral shaping of neural discharges by refractory effects. J Acoust Soc Am 93: 3353–3364, 1993.[CrossRef][ISI][Medline]
Filion M and Tremblay L. Abnormal spontaneous activity of globus pallidus neurons in monkeys with MPTP-induced Parkinsonism. Brain Res 547: 142–151, 1991.[CrossRef][ISI][Medline]
Franklin J and Bair W. The effect of a refractory period on the power spectrum of neuronal discharge. SIAM J Appl Math 55: 1074–1093, 1995.[CrossRef]
Gauss R and Seifert R. Pacemaker oscillations in heart and brain: a key role for hyperpolarization-activated cation channels. Chronobiol Int 17: 453–469, 2000.[CrossRef][ISI][Medline]
Gray CM. Synchronous oscillations in neuronal systems: mechanisms and functions. J of Computational Neuroscience 1: 11–38, 1994.
Halliday DM, Rosenberg JR, Amjad AM, Breeze P, Conway BA, and Farmer SF. A framework for the analysis of mixed time series/point process data–theory and application to the study of physiological tremor, single motor unit discharges and electromyograms. Prog Biophys Mol Biol 64: 237–278, 1995.[CrossRef][ISI][Medline]
Heimer G, Bar-Gad I, Goldberg JA, and Bergman H. Dopamine replacement therapy reverses abnormal synchronization of pallidal neurons in the 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine primate model of Parkinsonism. J Neurosci 22: 7850–7855, 2002.
Hurtado JM, Gray CM, Tamas LB, and Sigvardt KA. Dynamics of tremor-related oscillations in the human globus pallidus: a single case study. Proc Natl Acad Sci USA 96: 1674–1679, 1999.
Hurtado JM, Rubchinsky LL, Sigvardt KA, Wheelock VL, and Pappas CT. Temporal evolution of oscillations and synchrony in GPi/muscle pairs in Parkinson's disease. J Neurophysiol 93: 1569–1584, 2005.
Hutchison WD, Lozano AM, Tasker RR, Lang AE, and Dostrovsky JO. Identification and characterization of neurons with tremor- frequency activity in human globus pallidus. Exp Brain Res 113: 557–563, 1997.[CrossRef][ISI][Medline]
Jarvis MR and Mitra PP. Sampling properties of the spectrum and coherency of sequences of action potentials. Neural Comp 13: 717–749, 2001.
Jenkins GM and Watts DG. Spectral Analysis and Its Applications. San Francisco,CA: Holden-Day, 1968.
Kuffler SW, Fitzhugh R, and Barlow HB. Maintained activity in the cat's retina in light and darkness. J Gen Physiol 40: 683–702, 1957.
Lenz FA, Tasker RR, Kwan HC, Schnider S, Kwong R, Murayama Y, Dostrovsky JO, and Murphy JT. Single unit analysis of the human ventral thalamic nuclear group: correlation of thalamic "tremor cells" with the 3–6 Hz component of Parkinsonian tremor. J Neurosci 8: 754–764, 1988.[Abstract]
