Quantifying Spike Train Oscillations: Biases, Distortions and Solutions

Estimation of the power spectrum is a common method for identifying oscillatory changes in neuronal activity. However, the stochastic nature of neuronal activity leads to severe biases in the estimation of these oscillations in single unit spike trains. Different biological and experimental factors cause the spike train to differentially reflect its underlying oscillatory rate function. We analyzed the effect of factors, such as the mean firing rate and the recording duration, on the detectability of oscillations and their significance, and tested these theoretical results on experimental data recorded in Parkinsonian non-human primates. The effect of these factors is dramatic, such that in some conditions, the detection of existing oscillations is impossible. Moreover, these biases impede the comparison of oscillations across brain regions, neuronal types, behavioral states and separate recordings with different underlying parameters, and lead inevitably to a gross misinterpretation of experimental results. We introduce a novel objective measure, the "modulation index", which overcomes these biases, and enables reliable detection of oscillations from spike trains and a direct estimation of the oscillation magnitude. The modulation index detects a high percentage of oscillations over a wide range of parameters, compared to classical spectral analysis methods, and enables an unbiased comparison between spike trains recorded from different neurons and using different experimental protocols.

The time of occurrence of action potentials emitted by a single neuron; i.e., single unit spike trains, are a major source of neurophysiological data stemming from both intracellular and extracellular recordings. These neuronal spike trains may be viewed as a stochastic point process where a discrete event represents each action potential [9]. The generation of each spike within the train is assumed to be dependent on an underlying instantaneous firing rate. The resultant point process reflects the originating rhythm only partially since the individual events are stochastically generated from the rate function [10]. Thus, despite an underlying oscillatory firing rate, in most cases the neuron will skip a large portion of the oscillation cycle or even entire cycles [11]. The most simplistic statistical spike train model assumes that the generation of each spike is dependent solely on the underlying instantaneous firing rate, and is independent of all other previous spikes. This model is termed the inhomogeneous Poisson process when the instantaneous firing rate changes over time. Spectral analysis of an experimentally recorded spike train under this assumption is thus assumed to reflect the oscillations of the underlying instantaneous rate directly. However, different properties of the neuron or the experiment can cause the spike train to reflect its underlying oscillatory firing rate well or poorly, and hence can facilitate or impede the detection of an underlying oscillation [7,10].
In this paper, we address the problem of oscillation detection in spike trains. First, we quantify the influence of different biological and experimental factors on the detection of the spectral peak and its significance. We then analytically derive a solution for the modulation index of an oscillatory Poissonian neuron and present a novel objective measure that enables reliable detection of oscillations in spike trains. We investigate the oscillations in experimentally recorded neurons from Parkinsonian primates, and compare these experimental results to our numerical and analytic findings. Finally, we derive a solution for the evaluation of the actual recording duration required for the detection of spike train oscillations in experimental data.

Results
The spiking activity of oscillatory neurons can be modeled as an inhomogeneous Poisson process whose rate λ(t) is modulated by a cosine function (Fig 1A): where r 0 is the baseline firing rate, 0 m 1 is the modulation index, and f 0 is the oscillation frequency. Classically, the oscillatory nature of the neuronal activity is assessed using spectral estimation methods. The power spectrum of this rate function enables the identification of the oscillation base frequency (f 0 ) which appears as a clear peak in that frequency, while all the other frequencies have no power ( Fig 1B). However, estimation of the power spectrum of the Poisson generated spike train, ρ(t), from this rate function (Fig 1C) results in a reduced peak at the base frequency and increased power in all the other frequencies, such that detection of the base frequency is not straightforward (Fig 1D). The power spectrum can be normalized to reflect the signal-to-noise ratio (SNR) in standard deviations of the power in each frequency relative to the mean power in the 100-500 Hz frequency band of the spike train which serves as a noise baseline for different values of the baseline rate ( Fig 1E-1G). The power spectrum of this generated spike train presents an increased peak at f 0 relative to the baseline. The SNR of these simulated neurons varies linearly as a function of the base firing rate of the neuron ( Fig  1H). As a result, the detection of significant oscillations crossing a specific SNR threshold is not possible for a neuron with a low baseline firing rate (Fig 1E), compared to neurons with higher firing rates (Fig 1F-1G), which have a higher SNR and are therefore identified as oscillatory. Simulated neurons can illustrate the potential firing rate induced bias in the assessment of spectral activity. To examine the effect of the mean firing rate on the detection of oscillations in experimentally recorded data, we calculated the power spectrum of globus pallidus internus (GPi) neurons recorded in four non-human primates (NHPs) following injections of 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) which led to the appearance of Parkinsonian symptoms. The firing pattern of a large proportion of the GPi neurons in the Parkinsonian NHP resembles Poissonian firing, which is modulated in an oscillatory manner at a base frequency in the beta band (10-15 Hz) [12][13][14] (Fig 2A). In the experimental dataset (229 neurons), which consisted of neurons with firing rate in the wide range of 20 to 180 spikes/s and a positive corrected SNR value, the SNR of the highest peak within the beta band varied considerably but exhibited a linear relation with the firing rate ( Fig 2B). Thus, the fraction of neurons identified as oscillatory; i.e., neurons with a SNR!5, increased with the firing rate of the neurons ( Fig 2C). This implies that in real experimental data, naïve usage of the power spectrum results in a biased detection of oscillatory activity that can easily lead to misinterpretation of the experimental dataset.
Simulations were next used to quantify the effect of different biological (mean firing rate, oscillation's modulation magnitude) and experimental factors (recording duration) on the SNR. The values of each condition were obtained from 1000 simulated Poisson spike trains generated from a single oscillatory (f 0 = 12 Hz) rate function. When the duration was fixed, and the SNR was estimated over a range of firing rates and modulations, the simulations showed that the SNR had a linear relation to the baseline firing rate (Fig 3A), and a squared relation to the modulation index ( Fig 3B). When the modulation index was fixed, and the SNR was estimated over varying firing rates and durations, the simulations indicated a linear relation of the SNR to the firing rate (Fig 3C), and a square root relation to the total length of the recording ( Fig 3D). These relations suggest a dependence on these factors that impacts the detection of significant oscillatory activity (Fig 3E-3F). Increases in the firing rate, modulation index ( Fig 3E) and recording duration ( Fig 3F) result in increased detection of the spectral peak. This dependence on biological and experimental variables thus shows that the ability to objectively detect a peak in the power spectrum is limited. This dependence hinders the comparison of different behavioral states or brain areas and leaves them prone to biases. The formulation of an objective method of measuring oscillations is thus a necessity to enable unbiased comparisons of spike trains arising from different biological and experimental sources.
The power spectrum of an inhomogeneous spike train over a period (T) is (see methods: power spectrum of a finite inhomogeneous Poisson process): In the specific case of cosine rate modulation over a base frequency (f 0 ) (Eq 1) the power spectrum is: (see methods: power spectrum of an inhomogeneous Poisson process with an  The term sinc(fT) decays as 1/T. Thus, for large enough values of T, most of the sinc terms will decay (e.g., for oscillations at 12 Hz, sinc(2 f 0 T) is negligible within a few seconds of recordings), except for the case when the frequency term inside the sinc is zero (i.e. f-f 0 = 0), and then sinc(0) equals 1 for every T. Consequently, the formulation may be simplified for the case of the base frequency (f = f 0 ); i.e., the peak power, to: And for all other frequencies (f 6 ¼ f 0 ); i.e., the baseline power, to: The magnitude of the peak power and its relation to the baseline power are dependent on multiple factors; namely r 0 , T, m. Thus, a measure that is independent of subjective properties is required which we term the modulation index (m). This measure can be extracted from the simplified equation of the peak power (Eq 4) This equation enables the extraction of the modulation index for any frequency. When there is a real underlying oscillation in that frequency, the outcome of the equation will be the described modulation index. However, when there is no underlying oscillation in that frequency, the result will tend to be zero, as the value of S ρ T (f 6 ¼ f 0 ) approaches r 0 , as shown in Eq 5. Through this measure, we can reconstruct the rate function of an inhomogeneous Poissonian oscillatory spike train (Fig 4A, top). First, the power spectrum is estimated, and the peak power is extracted (Fig 4A, middle). Then, by using the estimated peak power (Ŝ r T f ¼ f 0 ð Þ), the mean firing rate (r 0 ) and the total recording time (T), the modulation index is extracted, and the estimated rate function (lðtÞ) can be fully described (Fig 4A, bottom). We calculated the modulation index for the simulation described above. For each of the 1000 generated spike trains we calculated the mean peak power, corrected it for comparison with Welch 0 s estimator, which is a common method for estimating the power spectrum, (see methods: correction for the power estimated by Welch's method), and used it for the calculation of the modulation index. Comparison of the estimated modulation and the original modulation demonstrated that the modulation index was constant across a wide range of parameters over the simulated data ( Fig 4B). The significance level for detecting oscillatory neurons using the modulation index depends on the firing rate, and is defined as the mean result for the modulation index of zero at a specific firing rate + 2 standard deviations, as revealed by the simulations. The results of the detection of oscillatory neurons using the modulation index indicate that it is dependent on the baseline rate, such that the detection is better for higher firing rates ( Fig 4C). For example, when the modulation is 0.4, and the required detection probability is 80%, the firing rate should be higher than 15 spikes/s. However, for all firing rates, the detection using the modulation index was better than the detection using the spectral peak ( Fig 4D). We applied the modulation index to the recordings from the GPi of the Parkinsonian NHP ( Fig 4E) and obtained a constant modulation index (m ¼ 0:23 AE 0:13 mean ± SD). In our experimental dataset, 98 neurons were modulated; i.e., had a modulation index greater than zero, and 78 of these were significant, according to the aforementioned significance test. The results for all neurons revealed that the oscillation modulation was independent of the firing rate (R 2 = 0.02, p >0.1), unlike the results obtained using the power spectrum ( Fig 2B).
The modulation index was derived above for an inhomogeneous Poisson spike train where the power spectrum exhibits equal power in all frequencies except for the base frequency. However, real neurons deviate from the Poisson model, primarily as a result of the refractory period. Refractoriness prevents the neuron from firing two successive spikes within a short interval, and thus the spike train is never a true Poisson process. In the case of an oscillatory spike train, refractoriness distorts the oscillation, and the modulation appears to be smaller than the real modulation of its rate function (Fig 5A). This occurs due to the larger effect of the refractory period around the peak of the oscillation as more spikes are expected at that time. This underestimation of the modulation index is more prominent for high firing rates due to the increased effect of the refractory period ( Fig 5B). The modulation index can however be corrected to accommodate for the refractory period (see methods). The correction procedure reconstructs the modulation index (m p ) of the rate function that will result in the analytically calculated modulation index, which is smaller than expected as a result of the refractory period, and thus extracting the "true" modulation of the driving rate function (Fig 5C). We applied this correction to GPi spike trains. This resulted in higher modulation index values (m p ¼ 0:35 AE 0:12 mean ± SD) that are independent from the firing rate (R 2 = 0.0005, p >0.5) (Fig 5D). The percentage of significantly modulated neurons (97/229, 42%) was significantly larger than the percentage of neurons that were identified as significantly oscillatory by the SNR measure (64/229, 28%) (x 2 test, p < 0.01).
The detection of oscillatory activity depends on a critical experimental factor-the total recording duration (T). This factor must be set prior to the experiment itself to avoid a case of failed detection due to the experiment 0 s time constraint (Fig 3F) rather than an actual lack of oscillations. The analogous analytic SNR is defined as: And for f = f 0 : And for large enough T: The expected analytic SNR and the numerical estimated d SNR are similar, so we can replace the analytic SNR with the desired d SNR, such that: This equation is true when no windows are used. When using Welch's method, the SNR must be multiplied by the root of the number of windows used (see methods: effect of the window size on the SNR), which is defined as: where wl is the window length, and N wins is the number of windows used.
From this equation, the required time for significance of the d SNR is: We calculated the required recording time for a window length of 1 second for various firing rates, for a d SNR of 1 to 7 while the modulation was fixed at 0.25 (Fig 6A), and for a fixed d SNR of 5 for various modulations (Fig 6B). High firing rates as well as high modulation index values will result in a shortening of the required recording duration. In addition, setting the significance of the SNR to lower values will shorten the required recording duration. Thus, for example, for a GPi neuron with a mean firing rate of 75 spike/s and an estimated modulation index of 0.25, less than half a minute of recording will suffice to detect a spectral peak with SNR of 5. However, in order to detect a spectral peak with SNR of 7 for the same neuron, about 1 minute of recording is necessary.

Discussion
Classical oscillation estimation methods for neuronal data represent the magnitude of oscillatory activity in terms of elevated peaks in a specific frequency, as compared to the baseline defined over other frequencies. In the first part of this study we quantified the dependence of the magnitude of the spike train oscillation, as revealed by power spectrum estimation, on different biological and experimental factors. We showed that in both simulated and experimental data, the estimated magnitude of the oscillation depends highly on the mean firing rate of the spike train. We used simulations to quantify this dependence and investigated the influence of these factors on the detectability of oscillatory spike trains.
In the second part of this study, we introduced a novel method for assessing the oscillatory nature of a spike train by calculating the modulation index of the oscillation. The modulation index can be estimated on the basis of the mean firing rate, the total recording duration, and the magnitude of the peak in the power spectrum, and provides an objective measure of the oscillation magnitude. We applied this measure to the simulated spike trains, and showed that the measure produces a correct estimation over a wide range of biologically plausible parameters. The application of this measure to experimental data recorded from GPi neurons in the Parkinsonian NHP revealed that the modulation index is independent of the firing rate. We introduced the corrections that need to be applied to the estimated power revealed by Welch's method, to compare the index to the analytically calculated power. Furthermore, we presented adaptations to the modulation index to account for deviations from the Poisson process assumptions. Explicitly, we show its usage in an inhomogeneous Poisson process with an absolute refractory period, which dramatically alters the spectrum of real neurons. Finally, we proposed a practical method for estimating the recording time required to accurately detect oscillations in neurophysiological experiments.
The effect of the firing rate on the size of the spectral peak is dramatic. In the mammalian brain, the range of firing rates between brain areas and neuronal types is considerable, and even within a specific brain region and a specific neuronal type, the heterogeneity of firing rates within the population is large. For example, the firing rates of the GPi neurons in the Parkinsonian NHP used in our dataset ranged from 25 spikes/s to over 100 spikes/s. Furthermore, the rate distribution can vary between behavioral states: in the normal state, firing rates increase when attention is directed toward a stimulus [15] and changes in the firing rates of hippocampal neurons occur as the result of experience [16]. Pathological conditions influence firing rates as well, as shown in the firing rate changes throughout the basal ganglia over the course of Parkinson's disease [17,18]. These differences consequently lead to a situation in which identification of oscillations using the spectral peak may be difficult in some conditions, and moreover, bias the comparison across different brain areas, neuron types and behavioral states. Nonetheless, most studies do not take different firing rates into account when inferring their spectral results. For example, previous studies of the Parkinsonian primate have shown that there is a larger fraction of oscillating neurons as well as an increased magnitude of oscillations in the GPi relative to the globus pallidus externus (GPe) [14,19]. These studies neglected to take the substantially different firing rates between the two brain areas during parkinsonism into account: the mean firing rate of GPe neurons in the MPTP treated NHP is 45 spikes/s, whereas the mean firing rate of the GPi neurons is 75 spikes/s [19]. According to our simulation results (Fig 3E), the detection of significant oscillations in neurons with a realistic modulation index of 0.25 and a firing rate of 45 spikes/s is less than 20%, whereas when the firing rate is 75 spikes/s the detection is about 80%. Thus, the conclusions relating to different oscillatory activities may be derived from the firing rates themselves and not from the oscillatory nature of the neurons.
The firing rate effect on the spectral peak is also a major caveat in the interpretation of the oscillatory nature of MUA. Previous studies have reported that MUA appears to be more oscillatory than single-unit activity (SUA) [14,20]. However, this phenomenon is affected to a large extent by the higher firing rate of the MUA that contains more than one spike train, such that the comparison between the spectral peak of SUA and MUA is problematic and may lead to a misinterpretation of the results.
An additional factor that biases the spectral peak estimation considerably is not dependent on the neuronal properties but rather on an experimental property: the duration of the recording. As the total recording duration increases, the SNR will increase. Yet, most studies are not aware of this bias when comparing the results of spectral analyses of recordings with different durations. Moreover, in some situations, the duration of the recording can be extremely short, such as recordings in the operation room which can only yield a few seconds of recording [20], and thus are not long enough for the detection of a spectral peak despite an underlying oscillatory rate function. Even worse, transient changes in the oscillatory activity of a neuron in response to a behavioral task might prevent the detection of a spectral peak which is unique to the transient period. Such transient changes occur for instance in the gamma frequency in relation to movement [1].
Thus, a measure quantifying oscillations is needed that goes beyond the identification of significant spectral peaks. Several methods for the estimation of oscillatory activity in spike trains have been suggested based on the auto-correlation function [10,21] and analysis of the power spectrum [7]. As the spectrum may be defined as the Fourier transform of the auto-correlation function (Wiener-Khinchin theorem), these two groups of methods are biased. The bias could be visualized by the auto-correlation function of the spike trains as well as in its power spectrum (S1A-S1F Fig). As in the power spectrum, the SNR of the auto-correlation function is dependent on the firing rate, such that when the firing rate is higher, the oscillatory nature of the spike train is more evident in the auto-correlation function, while during low firing rate, the oscillation cannot readily be seen (S1G Fig). In order to overcome the bias in the auto-correlation, one of the methods defines an oscillatory score that is less sensitive to the firing rate [10] but still depends on the frequency band, where higher bands yield higher scores. Other methods have dealt with the problem of finite recording durations by applying corrections to the confidence limits of the spectral estimations [7] but have not handled the different firing rates explicitly.
The major drawback in the standard methods of assessing spike train oscillations were addressed in consecutive steps within this manuscript. Initially, we quantified the factors that bias the magnitude of the spectral peak. Next, we introduced a measure that overcomes these factors, and leads to a reliable detection of oscillations and a direct estimation of the strength of oscillations: the modulation index. This method is simple and fast, and can be customized to accommodate the size of the window and the spectral smoothing applied prior to the spectral estimation. Finally, a practical method for estimating the required recording duration was proposed, based on the mean firing rate, the evaluated modulation index and the desirable SNR. In order to evaluate the mean firing rate and the modulation index, a few short preliminary recordings can be performed, and these two parameters can be grossly assessed. The significance levels for the SNR have to be chosen, and then, the required recording duration may be estimated. In other situations, where the recording duration is fixed, and cannot be adjusted, the experimenter should be aware of the limitations of detecting oscillations from the recordings and use our analytic results to estimate the fraction of unidentified oscillatory neurons.
Throughout most of this study we assumed that the spiking activity of single neurons follows a Poisson distribution [22][23][24], i.e. they fire stochastically and the probability of generating a spike at a certain time depends solely on the underlying rate function. In real neurons, there are deviations from the Poisson model. The most prominent deviation, the refractory period, was addressed in this study by correcting the estimation of the modulation index. The correction procedure is based on information extracted from the given spike train, and on the underlying rate function, including the deviation from the Poisson process. The process numerically finds the original modulation index of the rate function that results in the modulation index calculated analytically from the recorded data. For the case of refractory period, we modeled the rate function with an absolute refractory period following a spike. However, in some unique cases, when the refractory period causes the oscillation frequency peak to be too small, which results in a modulation index of zero, the correction procedure will not be able to reveal the original modulation. In this situation, a shuffling procedure on the first order ISI's could compensate for the distortion of the spectrum caused by the refractoriness and lead to an increased accuracy of the peak detection [25]. The general rate function model could be expanded to accommodate for other deviations from the Poisson model, such as relative refractory period and bursting activity; i.e., an elevated firing probability immediately after the refractory period.
In conclusion, the modulation index can provide an objective quantification of spike train oscillations, and thus an unbiased comparison across brain regions, behavioral states and separate recordings with different recording lengths. Using this quantification method can expand our understanding of neural oscillations, and their role in normal and pathological states.
All the code required for the estimation of the modulation index and the required recording time is available as custom MATLAB code (compatible with versions 2007B-2014A) on our website: http://neurint.ls.biu.ac.il/software/

Ethics statement
The neuronal recordings were taken from four Cynomologus monkeys (Macaca Fascicularis), that underwent 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) injections leading to a stable Parkinsonian state. The monkeys' water and food consumption and weight were checked daily and their health was monitored by a veterinarian. All procedures were in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals, Bar-Ilan University Guidelines for the Use and Care of Laboratory Animals in Research and the recommendations of the Weatherall Report. All procedures were approved and supervised by the Institutional Animal Care and Use Committee (IACUC). These procedures were approved by the National Committee for Experiments on Laboratory Animals at the Ministry of Health (permit number BIU150605). Full details of the experimental protocol appear elsewhere [14,26].

Spike train simulations
The spike trains were simulated with a resolution of Δt = 1 ms. The spike trains were modeled as an inhomogeneous Poisson process, whose rate is modulated by a 12 Hz cosine function. Spike trains with refractory period were modeled by setting the instantaneous firing rate to zero for 2 ms following a spike.

Power spectrum
The power spectrum was estimated based on Welch's method [27] for spectral estimation using non-overlapping segments. In all spectral calculations, the 1000 bin segments (see methods: Effect of the window size on the SNR) were windowed using a Hamming window. The power spectral densities of each segment were calculated, and were then averaged. The bin size (Δt) of both the experimental and simulated data was 1 ms, leading to a spectral resolution of 1 Hz and a maximal frequency of 500 Hz.

Power spectrum of a finite inhomogeneous Poisson process
The following is based on the general formulation by Pinhasi and Lurie [28]: Let's consider a Poisson process of occurring delta functions, with a non-constant and time dependent rate λ (t). A portion of the process realization in the time interval [0 T] is given by the finite summation: where t i are the successive firing instances, k(t) is a random variable that counts the number of delta functions occurring during the time interval T, with an expected value of: The event appearance t i is an independent random variable with the probability density function: The Fourier transform of ρ T (t) is: The power spectral density can be calculated: In a Poisson process: s k 2 T ð Þ ¼ k T ð Þ, and thus, the second moment is: Inserting this in Eq 17 yields: The statistical average: Leading to the power spectrum formulation as: where the first term in brackets is the statistical average of the number of events over T, and the second term is the power spectrum of the instantaneous rate function λ(t) in the interval [0 T].

Power spectrum of an inhomogeneous Poisson process with an oscillatory rate function
For a rate function defined as: The expected number of spikes within a period T is:

Correction of the modulation index for deviations from the Poisson model
The deviation of the spike train from the Poisson model results in an incorrect estimation of the modulation index. A correction process could compensate for this: (1) The firing rate of the spike train assuming a true Poisson process (r p ) is calculated using the recorded spike train.
(2) An oscillatory rate function is generated of length T and baseline rate which is set tô r p . The modulation of the rate function (m p ) is set to the analytically calculated modulation index (mÞ of the spike train.