Quantifying Spike Train Oscillations: Biases, Distortions and Solutions

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 nonhuman 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. estimationWe introduce a novel objective measure, the "modulation index", which overcomes these biases, and enables reliable detection of oscillations from spike trains and a direct estima-tion 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.

In this manuscript, we expose major biases and distortions which arise from the quantification of neuronal spike train oscillations. These, previously neglected, biases hinder the comparison of oscillations across brain regions, neuronal types and behavioral states, leading inevitably to severe misinterpretation of experimental results. We demonstrate the biases computationally, formulate them analytically and validate their appearance and magnitude in an experimental dataset recorded from Parkinsonian nonhuman primates. detectionNext, following a formulation of the distortions, we introduce a novel objective measure, the "modulation index", which overcomes these biases, and enables a reliable de-tection of oscillations from spike trains and a direct estimation of the oscillation magnitude. The modulation index is validated on the same experimental data demonstrating the unbiased detection of beta oscillation in the globus pallidus during Parkinsonism. The manuscript provides a solid infrastructure for oscillation analysis Which benefits multiple neuroscience fields ranging from basic science to clinical studies, moreover its results may be expanded to encompass additional fields in biology Which require the spectral analysis of point process data. This is a PLOS Computational Biology Methods paper

Neuronal oscillations are typically classified into a range of frequencies including the delta (1—4HZ), theta (4—8Hz), alpha (8—12HZ), beta (12—30 Hz), and gamma (30—80Hz) bands [3]. Enhanced expression of specific oscillatory frequencies or oscillations in a broader band is perceived as indicative of different normal functions, such as the enhanced gamma preceding movement [1] and different pathological conditions such as the enhanced beta associated with Parkinsonism [4]. Analysis of the power spectrum is a common method for identifying enhanced (or reduced) oscillations in neuronal data, and is widely used on a variety of brain signals spanning multiple orders of magnitude, such as electroencephalograms (EEG), local field potentials (LFP), multiunit activity (MUA), single unit spike trains and cellular membrane potentials [5—8].

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

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.

The spiking activity of oscillatory neurons can be modeled as an inhomogeneous Poisson process Whose rate Mt) is modulated by a cosine function (Fig 1A):

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 (fo) 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, p(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 ath 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.

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 nonhuman primates (NHPs) following injections of 1-meth-yl-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—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 SNRZ 5, increased with the firing rate of the neurons (Fig 2C). This implies that in real experimental data, nai've usage of the power spectrum results in a biased detection of oscillatory activity that can easily lead to misinterpretation of the experimental dataset.

The values of each condition were obtained from 1000 simulated Poisson spike trains generated from a single oscillatory (fo = 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 (fo) (Eq 1) the power spectrum is: (see methods: power spectrum of an inhomogeneous Poisson process with an oscillatory rate function) The term sian‘T) 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(2f0 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-fo = 0), and then sinc(0) equals 1 for every T. Consequently, the formulation may be simplified for the case of the base frequency (f: f0); i.e., the peak power, to: The magnitude of the peak power and its relation to the baseline power are dependent on multiple factors; namely r0, T, m. Thus, a measure that is independent of subjective properties is required Which we term the modulation index (161). This measure can be extracted from the simplified equation of the peak power (Eq 4)

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 SPT(f;£ f0) approaches r0, 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 mean firing rate (f0) and the total recording time (T), the modulation index is extracted, and the estimated rate function 01(0) can be fully described (Fig 4A, bottom).

For each of the 1000 generated spike trains we calculated the mean peak power, corrected it for comparison with Welch’ 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 (1% = 0.23 :I: 0.13 mean 1 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 (R2 = 0.02, p >0.1), unlike the results obtained using the power spectrum (Fig 2B).

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 (mp) 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 (mp = 0.35 :I: 0.12

This factor must be set prior to the experiment itself to avoid a case of failed detection due to the experiment’ s time constraint (Fig 3F) rather than an actual lack of oscillations. The analogous analytic SNR is defined as: And for large enough T: The expected analytic SNR and the numerical estimated 87V} are similar, so we can replace the analytic SNR With the desired S/N\R, 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 NWZ-ns is the number of Windows used. From this equation, the required time for significance of the 87V} is: We calculated the required recording time for a window length of 1 second for various firing rates, for a 87V} of 1 to 7 while the modulation was fixed at 0.25 (Fig 6A), and for a fixed 87V} 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.

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.

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.

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

Estimating the required recording duration. Analytic results forthe required recording duration, (A) over different firing rates and for increasing SNR, with modulation of 0.25 and (B) over different firing rates and for various modulations, calculated for a SNR of 5.

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.

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.

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

Several methods for the estimation of oscillatory activity in spike trains have been suggested based on the autocorrelation function [10,21] and analysis of the power spectrum [7]. As the spectrum may be defined as the Fourier transform of the auto-cor-relation function (Wiener-Khinchin theorem), these two groups of methods are biased. The bias could be visualized by the autocorrelation function of the spike trains as well as in its power spectrum (SlA—SlF Fig). As in the power spectrum, the SNR of the autocorrelation 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 autocorrelation function, while during low f1r-ing rate, the oscillation cannot readily be seen (81G Fig). In order to overcome the bias in the autocorrelation, one of the methods defines an oscillatory score that is less sensitive to the f1r-ing 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.

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.

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

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.

In all spectral calculations, the 1000 bin segments (see methods: Efiect 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 (At) 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.

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 )L (t). A portion of the process realization in the time interval [0 T] is given by the finite summa-tion: Where t,- 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 1‘,- is an independent random variable With the probability density function: The Fourier transform of pT(t) is: The power spectral density can be calculated: 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 Mt) in the interval [0 T].

For a rate function defined as: The expected number of spikes Within a period T is: The Fourier transform of the instantaneous rate in the time interval [0 T] for a large T is: Resulting in the power spectrum of the inhomogeneous Poisson process:

To do so, Welch's power needs to be multiplied first by the sampling frequency, and then divided by 2. Then, a correction due to the lower spectral resolution of Welch's power needs to be made: first multiply the power by the number of windows, and then subtract from it the baseline power multiplied by the number of windows minus 1. Finally, a correction due to the windowing needs to be applied: the correction term is the mean of the window multiplied by the length of the window, and then divided by the power of the window. This is expressed in the final formulation: Where F5 is the sampling frequency, Nwms is the number of Windows used, w is the Window, and wl is the length of the Window.

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 (fp) is calculated using the recorded spike train. (2) An oscillatory rate function is generated of length T and baseline rate Which is set to fp. The modulation of the rate function (mp) is set to the analytically calculated modulation index (161) of the spike train. (3) From this rate function, multiple new Poisson-like spike trains are generated, including the given deviations from the Poisson model. The modulation indeX measure of each spike train is estimated (161,), and the mean modulation indeX is calculated (161—5). (4) The mean simulated modulation indeX is compared to the original modulation indeX of the spike train. If the simulated modulation indeX is different from the original modulation index (:67, 7E :61), increase (in the case of an underestimation) or decrease (in the case of an overestimation) the modulation of the rate function (mp), and repeat step 3. If the simulated modulation indeX is similar to the original modulation indeX (161—5 2 :61) Within the required error boundaries, the Poisson equivalent modulation indeX is the modulation of the rate function (mp).

The correction process was explicitly implemented to compensate for the underestimation of the modulation indeX arising from the refractory period: (1) The firing rate of the spike train assuming no refractory period is calculated using the properties of the recorded spike train: the number of spikes (Nspikes), the recording duration (T) and the refractory period (Tref), Which could be extracted from the interspike interval (181) histogram: (2) The oscillatory rate function is generated of length T using a baseline rate (fp) and the modulation index (1%) calculated from the spike train. (3) From this rate function, 100 neW Poisson spike trains are generated using a refractory period of Tref, and their mean modulation indeX is calculated (161—5). (4) If the simulated modulation indeX is smaller than the original modulation index (2675 < 161), increase the modulation of the rate function (mp), and repeat step 3. If the simulated modulation indeX is similar to the original modulation indeX (261—5 2 m), the corrected modulation indeX is the modulation of the rate function (11% p ). Effect of the window size on the SNR

When the frequency of the oscillation is stable and tightly locked to a specific frequency, a longer window will yield a more precise power spectrum (Fig 7A). However, when the frequency jitters, using a window with a resolution smaller than the jitter size will not lead to a more precise spectrum, but rather to a spread of the power over the entire frequency jitter range (Fig 7B). Real neurons are not perfect oscillators in a single precise frequency, but instead are oscillatory within a jittered frequency range, which in the GPi of the Parkinsonian NHP is about 1 Hz (Fig 7C), justifying a window size of around the jitter size (a window length of 1000 leads to a 1Hz resolution). The length of the window affects the SNR, such that the SNR grows

sp/s. (D-F) The autocorrelation function of the same spike trains as in A-C. (G) The SNR of the peak in the autocorrelation as a function of the mean firing rates. The SNR is defined as the number of STDs from the mean of the autocorrelation function to the first peak near the zero point. The circles indicate the SNR of the peak shown in examples D-F, and the solid line indicates the fitted linear function. (H-]) The first order 181 histograms of the same spike trains shown in AC. The dashed vertical lines indicate the 1/f0 ISI.

We thank Yosef Pinhasi for help With the initial formulation, Moshe Abeles, Dorin Yael and Edward Stein for helpful comments, Anan Moran, Yaara Erez and Hadass Tischler for the experimental data.

Performed the experiments: AM IBG. Analyzed the data: AM IBG. Wrote the paper: AM IBG.

Appears in 54 sentences as: (2) Spike train (1) spike train (31) Spike trains (1) spike trains (21)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- However, the stochastic nature of neuronal activity leads to severe biases in the estimation of these oscillations in single unit spike trains .Page 1, “Abstract”
- Different biological and experimental factors cause the spike train to differentially reflect its underlying oscillatory rate function.Page 1, “Abstract”
- estimationWe introduce a novel objective measure, the "modulation index", which overcomes these biases, and enables reliable detection of oscillations from spike trains and a direct estima-tion of the oscillation magnitude.Page 1, “Abstract”
- 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.Page 1, “Abstract”
- In this manuscript, we expose major biases and distortions which arise from the quantification of neuronal spike train oscillations.Page 1, “Author Summary”
- detectionNext, following a formulation of the distortions, we introduce a novel objective measure, the "modulation index", which overcomes these biases, and enables a reliable de-tection of oscillations from spike trains and a direct estimation of the oscillation magnitude.Page 1, “Author Summary”
- Analysis of the power spectrum is a common method for identifying enhanced (or reduced) oscillations in neuronal data, and is widely used on a variety of brain signals spanning multiple orders of magnitude, such as electroencephalograms (EEG), local field potentials (LFP), multiunit activity (MUA), single unit spike trains and cellular membrane potentials [5—8].Page 2, “Introduction”
- 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.Page 2, “Introduction”
- These neuronal spike trains may be viewed as a stochastic point process where a discrete event represents each action potential [9].Page 2, “Introduction”
- 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.Page 2, “Introduction”
- Spectral analysis of an experimentally recorded spike train under this assumption is thus assumed to reflect the oscillations of the underlying instantaneous rate directly.Page 2, “Introduction”

See all papers in *April 2015* that mention spike train.

See all papers in *PLOS Comp. Biol.* that mention spike train.

Back to top.

Appears in 49 sentences as: firing rate (42) firing rates (21)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- 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 nonhuman primates.Page 1, “Abstract”
- The generation of each spike within the train is assumed to be dependent on an underlying instantaneous firing rate .Page 2, “Introduction”
- 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].Page 2, “Introduction”
- 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.Page 2, “Introduction”
- This model is termed the inhomogeneous Poisson process when the instantaneous firing rate changes over time.Page 2, “Introduction”
- 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].Page 2, “Introduction”
- where r0 is the baseline firing rate , 0 g m g 1 is the modulation index, and f0 is the oscillation frequency.Page 3, “Results”
- The SNR of these simulated neurons varies linearly as a function of the base firing rate of the neuron (Fig 1H).Page 3, “Results”
- 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.Page 3, “Results”
- Simulated neurons can illustrate the potential firing rate induced bias in the assessment of spectral activity.Page 3, “Results”
- 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 nonhuman primates (NHPs) following injections of 1-meth-yl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP) which led to the appearance of Parkinsonian symptoms.Page 3, “Results”

See all papers in *April 2015* that mention firing rate.

See all papers in *PLOS Comp. Biol.* that mention firing rate.

Back to top.

Appears in 14 sentences as: Poisson process (13) Poisson process: (1)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- This model is termed the inhomogeneous Poisson process when the instantaneous firing rate changes over time.Page 2, “Introduction”
- The spiking activity of oscillatory neurons can be modeled as an inhomogeneous Poisson process Whose rate Mt) is modulated by a cosine function (Fig 1A):Page 3, “Results”
- The power spectrum of an inhomogeneous spike train over a period (T) is (see methods: power spectrum of a finite inhomogeneous Poisson process ):Page 5, “Results”
- In the specific case of cosine rate modulation over a base frequency (fo) (Eq 1) the power spectrum is: (see methods: power spectrum of an inhomogeneous Poisson process with an oscillatory rate function)Page 5, “Results”
- Refractoriness prevents the neuron from firing two successive spikes within a short interval, and thus the spike train is never a true Poisson process .Page 9, “Results”
- Furthermore, we presented adaptations to the modulation index to account for deviations from the Poisson process assumptions.Page 11, “Discussion”
- Explicitly, we show its usage in an inhomogeneous Poisson process with an absolute refractory period, which dramatically alters the spectrum of real neurons.Page 11, “Discussion”
- 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 .Page 14, “Discussion”
- The spike trains were modeled as an inhomogeneous Poisson process , Whose rate is modulated by a 12 Hz cosine function.Page 15, “Spike train simulations”
- Power spectrum of a finite inhomogeneous Poisson processPage 15, “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 )L (t).Page 15, “Power spectrum of a finite inhomogeneous Poisson process”

See all papers in *April 2015* that mention Poisson process.

See all papers in *PLOS Comp. Biol.* that mention Poisson process.

Back to top.

Appears in 7 sentences as: experimental data (7)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- 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 nonhuman primates.Page 1, “Abstract”
- The modulation index is validated on the same experimental data demonstrating the unbiased detection of beta oscillation in the globus pallidus during Parkinsonism.Page 2, “Author Summary”
- Finally, we derive a solution for the evaluation of the actual recording duration required for the detection of spike train oscillations in experimental data .Page 3, “Introduction”
- This implies that in real experimental data , nai've usage of the power spectrum results in a biased detection of oscillatory activity that can easily lead to misinterpretation of the experimental dataset.Page 3, “Results”
- 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.Page 11, “Discussion”
- 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.Page 11, “Discussion”
- We thank Yosef Pinhasi for help With the initial formulation, Moshe Abeles, Dorin Yael and Edward Stein for helpful comments, Anan Moran, Yaara Erez and Hadass Tischler for the experimental data .Page 20, “Acknowledgments”

See all papers in *April 2015* that mention experimental data.

See all papers in *PLOS Comp. Biol.* that mention experimental data.

Back to top.

Appears in 6 sentences as: Poisson model (6)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- However, real neurons deviate from the Poisson model , primarily as a result of the refractory period.Page 9, “Results”
- In real neurons, there are deviations from the Poisson model .Page 14, “Discussion”
- 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.Page 14, “Discussion”
- Correction of the modulation index for deviations from the Poisson modelPage 18, “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.Page 18, “Correction of the modulation index for deviations from the Poisson model”
- (3) From this rate function, multiple new Poisson-like spike trains are generated, including the given deviations from the Poisson model .Page 18, “Correction of the modulation index for deviations from the Poisson model”

See all papers in *April 2015* that mention Poisson model.

See all papers in *PLOS Comp. Biol.* that mention Poisson model.

Back to top.

Appears in 4 sentences as: brain region (1) brain regions (3)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- 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.Page 1, “Abstract”
- These, previously neglected, biases hinder the comparison of oscillations across brain regions , neuronal types and behavioral states, leading inevitably to severe misinterpretation of experimental results.Page 1, “Author Summary”
- 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.Page 11, “Discussion”
- 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.Page 14, “Discussion”

See all papers in *April 2015* that mention brain regions.

See all papers in *PLOS Comp. Biol.* that mention brain regions.

Back to top.

Appears in 3 sentences as: neuronal activity (3)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- Estimation of the power spectrum is a common method for identifying oscillatory changes in neuronal activity .Page 1, “Abstract”
- However, the stochastic nature of neuronal activity leads to severe biases in the estimation of these oscillations in single unit spike trains.Page 1, “Abstract”
- Classically, the oscillatory nature of the neuronal activity is assessed using spectral estimation methods.Page 3, “Results”

See all papers in *April 2015* that mention neuronal activity.

See all papers in *PLOS Comp. Biol.* that mention neuronal activity.

Back to top.

Appears in 3 sentences as: r0 (3)

In *Quantifying Spike Train Oscillations: Biases, Distortions and Solutions*

- where r0 is the baseline firing rate, 0 g m g 1 is the modulation index, and f0 is the oscillation frequency.Page 3, “Results”
- The magnitude of the peak power and its relation to the baseline power are dependent on multiple factors; namely r0 , T, m. Thus, a measure that is independent of subjective properties is required Which we term the modulation index (161).Page 7, “Results”
- However, when there is no underlying oscillation in that frequency, the result will tend to be zero, as the value of SPT(f;£ f0) approaches r0 , as shown in Eq 5.Page 7, “Results”

See all papers in *April 2015* that mention r0.

See all papers in *PLOS Comp. Biol.* that mention r0.

Back to top.