## Abstract

The coherence between oscillatory activity in local field potentials (LFPs) and single neuron action potentials, or spikes, has been suggested as a neural substrate for the representation of information. The power spectrum of a spike-triggered average (STA) is commonly used to estimate spike field coherence (SFC). However, when a finite number of spikes is used to construct the STA, the coherence estimator is biased. We introduce here a correction for the bias imposed by the limited number of spikes available in experimental conditions. In addition, we present an alternative method for estimating SFC from an STA by using a filter bank approach. This method is shown to be more appropriate in some analyses, such as comparing coherence across frequency bands. The proposed bias correction is a linear transformation derived from an idealized model of spike–field interaction but is shown to hold in more realistic settings. Uncorrected and corrected SFC estimates from both estimation methods are compared across multiple simulated spike–field models and experimentally collected data. The bias correction was shown to reduce the bias of the estimators, but add variance. However, the corrected estimates had a reduced or unchanged mean squared error in the majority of conditions evaluated. The bias correction provides an effective way to reduce bias in an SFC estimator without increasing the mean squared error.

## INTRODUCTION

The relationship between local field potential (LFP) oscillations and single neuron action potentials, or spikes, has become a topic of increasing interest. The significance of coherence between spikes and LFPs has been shown to be important for auditory and visual sensory processing (Eggermont and Smith 1995; Fries et al. 1997; Gray and Singer 1989), motor tasks (Courtemanche et al. 2002; Donoghue et al. 1998; Murthy and Fetz 1996), memory representation (Harris et al. 2002; Lee et al. 2005; O'Keefe and Recce 1993; Quyen et al. 2008), attention (Fries et al. 2001), and neurological disorders (Foffani et al. 2007; Goldberg et al. 2004). For example, neurons activated by an attended stimulus are likely to be synchronized with gamma band (30–70 Hz) LFP oscillations (Fries et al. 2008). This is traditionally demonstrated through a measure called *spike field coherence* (SFC), which gives the strength of coupling between spike times and the phase of LFP at any given frequency. However, due to limitations in the finite amount of data available to measure this coherence, methods developed to study this phenomenon can only estimate the coherence. Therefore a better understanding of the quality of the estimation is needed if biological inferences are to be drawn from this measure.

Although it is easily describable as a qualitative physiological phenomenon, the exact quantitative definition of SFC is a matter of debate. No single definition can capture all the aspects of the spike–field relationship. Therefore the best we can do is to select an aspect of physiological importance and model it to satisfactorily define a quantity that represents this aspect. The strength of phase coupling between two continuous signals is traditionally inferred through the quantity called *coherence*, defined as the squared magnitude of the cross-spectrum divided by the product of the autospectra (Bendat and Piersol 1971). The same definition of coherence has also been used by some to describe spike field coupling, although we do not consider this definition here. Instead we examine estimators of alternatively defined quantities, based on the spike-triggered average (STA), which also represent spike–field coupling strength. This type of SFC estimator emphasizes different physiological information than estimators that begin in the frequency domain. For example, an STA-based estimator weights its value equally by spike, whereas traditional coherence weights all sections of the time series equally, regardless of spike density. It is not our intent to claim that STA-based methods are superior to others, only that they deserve close examination of the quality of estimation if they are to be used scientifically. Although the quantities estimated by the STA-based methods presented here are not mathematically equivalent to the traditional definition of coherence, they are also referred to as coherence or SFC hereafter.

The STA was originally constructed to measure postsynaptic potentials of dorsal roots (Mendell and Henneman 1968, 1971), but has since been applied to LFP. For each spike, an LFP segment of fixed duration is taken that is centered on the spike. All spike-centered segments are added together and divided by the number of spikes. This produces an average of the LFP happening within some window of each spike. One existing method for estimating coherence from the STA uses a normalized power spectrum of the STA to estimate the coherence (Fries et al. 1997). Another method, introduced here, uses a filter bank approach to estimate SFC from the amplitude at the center of the STA. However, in any estimate, there are two sources of potential error: bias and variability.

To define the bias problem, we must identify these sources of error. In the theoretical case, one would have an infinite amount of data from which to measure the asymptotic (population) coherence. In addition, these data would be second-order stationary (Brillinger 1981) in that the coherence does not change over time. However, since an infinite number of spikes is never available, a more realistic scenario is to have a finite subset of these data from which to estimate the sample coherence. In practice it is common to perform single trial estimation averaging over all available spikes. However, this creates a condition in which it is unknown whether the error in this single trial estimate comes from bias or variability. To investigate these sources of error, one can estimate the coherence of simulated data over multiple trials using a fixed number of spikes in each trial. In this way, a single trial estimate is defined as the quantification of an STA constructed from a given number of spikes *n*. A multitrial estimate of the coherence can thus be indexed by two factors: the number of trials *m* and the number of spikes in a trial *n*. In the analyses presented here, we choose *m* to be large enough (*M*) that we have an accurate measure of the expected value and variance of the coherence for a given *n*. We then use this framework to examine the bias and variance associated with the expected value of the coherence as the number of spikes is varied.

The bias is defined as the difference between the expected value of the coherence estimate for a given number of spikes and the asymptotic value of the coherence. The asymptotic coherence *c* is equal to the estimated coherence *ĉ*, as the number of spikes *n* approaches infinity
*f*. Because an infinite number of spikes never exists in practical situations, the asymptotic coherence is approximated as *ĉ _{N}*, where

*N*denotes a number of spikes large enough so that the estimated coherence is a good approximation of the asymptotic coherence. This will be quantified in methods. These two terms—

*ĉ*and

_{N}*c*—are used interchangeably for the remainder of this article and are both referred to as the asymptotic coherence.

The fact that the bias of the coherence estimator arises from the limitation of the finite number of spikes used in the STA can be best understood by the following thought experiment. Consider an experiment in which a total number of spikes *N* is recorded simultaneously with an LFP. Furthermore, the spikes and LFP are independent of one another and the asymptotic coherence is therefore zero. To evaluate an estimate of coherence, one could perform a single trial estimate using all of the spikes. Alternatively, one could select a subset of spikes *n* multiple times, to obtain an expected estimate of the coherence over the multiple trials. If there were no bias, the expected value of the estimated coherence and the asymptotic coherence would be the same and equal to zero. However, there is a difference, and thus a bias, because the power of the STA is dependent on the number of spikes. A single coherence estimate using a large number of spikes *N* will always be less than or equal to the expected value of the coherence using a subset of *n* spikes (*n* < *N*) for each trial. The case of zero coherence is a specific example of the general behavior exhibited by the STA in cases of nonzero coherence, as will be demonstrated in the following text.

We propose a correction for this bias that enables an estimation of SFC that is not dependent on the number of spikes used in the STA and, in methods, derive this correction analytically. To test this bias correction approach, we apply it to both the power spectrum based (PS-SFC) and filter bank based (FB-SFC) spike–field coherence methods identified earlier and described in methods. We use different models of spike distributions and LFP to evaluate the effect of the bias correction on the estimate of the SFC. We find that the correction reduces the effect of finite number of spikes (the bias) on the estimate, but adds variance. However the mean squared error (MSE) of the estimates is reduced after correction in the majority of cases simulated. We further show that these results are consistent with the effect of the bias correction on experimentally collected data.

## METHODS

To assess the effect of our bias correction, we evaluated the bias, variance, and MSE of the uncorrected and corrected PS-SFC and FB-SFC estimation methods using simulated and experimentally collected data. As stated in the introduction, bias is defined as the difference between the expected value of the coherence estimate for a given number of spikes and the asymptotic coherence. Variance (*Var*) is defined as the average squared deviation of each single trial estimate from the expected value of the estimate for the same conditions. MSE is defined as the sum of the bias squared and the variance. These sample statistics were calculated over a large number of trials *M*, to accurately approximate the population statistics they represent. Each is a function of frequency of oscillation *f* and indexed by the number of spikes *n* used in each trial *m*. Denoting 〈 … 〉* _{m}* as the arithmetic mean and

*E*[ … ] as the expected value operators

### Simulations

The LFP signal used in the simulations was represented as either a 50 Hz sine wave or as a real LFP signal filtered to an approximate 50 Hz wave. In the latter case, the raw (0.5–500 Hz) signal was digitized and filtered to 45–55 Hz. See *Surgery, data collection, and spike sorting* for details. Filtering was performed using a zero-phase (MATLAB filtfilt.m) finite impulse response filter with a Hamming window and order 2,000. Both LFP representations had a sampling rate of 2 kHz. The oscillatory signal was used as a reference for placing spikes according to the different spike distributions, to control the coherence strength at 50 Hz. The simulations could have been repeated at any other frequency band, but to reduce redundant results this was the only point on the frequency spectrum to be examined in simulations. A full investigation of how coherence in one band affects coherence in another is outside the scope of this presentation.

Spikes were distributed around a selected phase of the 50 Hz oscillation by one of two distributions: either a binary distribution or a Gaussian distribution. These distributions were centered on an arbitrary phase of the oscillation called the preferred phase. The binary distribution is described by mixing coherent spikes with incoherent spikes, so that each spike was either coherent or not coherent. Coherent spikes occurred at exactly the preferred phase of the oscillation, whereas the incoherent spikes occurred as a random Poisson variable, independent of the LFP. The parameter that controlled the coherence strength was the ratio of coherent spikes to total spikes. Specifically, the simulated coherence was given by
*p* is the probability of choosing *n _{X}* coherent spikes from the total number of spikes. The Gaussian distribution is a more realistic model, where individual spikes were neither fully coherent nor fully incoherent. Instead, all spikes were randomly assigned a time of occurrence, but according to a Gaussian distribution. The mean of this distribution occurred at the preferred phase of the oscillation and the SD was the parameter that controlled coherence strength. Therefore a narrower distribution would result in a higher coherence. This relationship between SD of the spike distribution σ and the simulated coherence γ is a Gaussian-weighted sum of sine waves, given by

For each trial, a 10 min section of signal was generated or selected at random, depending on the LFP representation. A phase was chosen randomly to be the preferred phase of the 50 Hz oscillation. Each spike then required another random selection to determine which cycle in the 10 min section it would be placed. Lastly, the exact spike time was obtained according to the distribution being used and referenced to the preferred phase in the chosen cycle. Once 2,000 spikes were placed to represent the total population, the estimation methods were applied. These methods used ≤1,000 of these spikes to create an STA, estimating the SFC after the addition of each new spike. The bias correction was also applied to the estimates after each new spike. The end result of a single trial was four estimates (uncorrected and corrected PS-SFC and FB-SFC) of coherence from 1 through 1,000 spikes. We did not show results beyond 100 spikes because the effect of the bias correction was less pronounced and it became difficult to visualize differences in estimation statistics beyond this point. For each of the four simulations, *M* = 1,000 trials were performed to accurately determine bias and variance.

### Analysis of experimentally collected data

To demonstrate efficacy of the bias correction on real data, several neurons were recorded simultaneously during a 1 h interval for each of three rats. From initial examination of the broad coherence spectrum, it was determined that these hippocampal neurons were most coherent with the theta band (centered on 6 Hz). Therefore we restricted our analysis to this frequency band. It would technically require an infinite number of spikes to perfectly determine the asymptotic coherence; however, a very good approximation was determined by the following method. First, the uncorrected coherence was estimated from all available spikes *N* in the time period. We wanted to ensure that this quantity would change very little with any amount of additional spikes (*N* + 1, *N* + 2, … , *N* + *k*), assuming stationary coherence. For a given SFC estimate using *N* spikes, adding more spikes will change the estimate by a maximum of 1/*N*. This is true in general because the maximum rate of decrease in a coherence estimate as the number of spikes increases is given by the case of independent spikes and LFP (zero asymptotic coherence). Here, the expected coherence estimate at a given number of spikes *n* is equal to 1/*n*; thus this is the maximum possible amount this estimate can change as *n* increases. Since an estimate of any nonzero asymptotic coherence will not change more than for the case of zero asymptotic coherence, the upper bound on the amount of uncertainty in a coherence estimate is given by 1/*n*. In examining the data, we chose only the neurons that displayed uncertainties of <1% of the asymptotic coherence. This left 9 neurons, which had asymptotic coherence values ranging from 0.005 to 0.1. Once these values were determined, estimation methods were evaluated in a manner similar to that of simulation. Coherence was estimated for each number of randomly chosen spikes *n*, ranging from 1 to 100. The spikes chosen did not necessarily occur in succession, but were selected with uniform probability from the entire recording. For each value of *n*, 1,000 trials were performed and the mean and variance of the estimated coherence were calculated before and after bias correction.

### Estimation methods

The PS-SFC was estimated in a manner similar to that reported by others (Fries et al. 2001). First, an STA was constructed from the raw signal with a window 150 ms wide around the spike. The power spectral density was estimated by use of a multitaper approach with two orthogonal tapers specified from discrete prolate spheroidal sequences (MATLAB pmtm.m). The number of tapers was kept to a minimum to isolate the power in the band surrounding the 50 Hz point on the spectrum. With these parameters the frequency resolution was 10 Hz. This resolution implied a fractional bandwidth of 0.2 at 50 Hz, the same bandwidth as that of the filter used for simulation. For each LFP segment used to construct the STA, the power spectrum was estimated using the above-cited technique. These spectra were then averaged together to obtain the average spectrum of the LFP segments. Finally, the normalized spectrum was obtained by dividing the power spectrum of the STA by the average spectrum of the LFP. This normalized power spectrum of the STA is the definition of SFC for this method

The FB-SFC method, introduced here, uses a filter bank to measure the amplitude of each frequency component in the STA spectrum. Instead of constructing the STA from the raw LFP, the signal was first filtered into a narrow band (45–55 Hz), as described earlier. In general, this technique would require multiple filters to estimate power over a spectrum of interest, but since simulated spikes were correlated only to the 50 Hz oscillation, only one filter was required. An STA was constructed from the filtered signal in the same way as that done for the raw signal. For each LFP segment, the absolute peak heights were found and averaged together to obtain the average peak height for the segment. Each segment was then normalized by dividing it by its average peak height. The STA was rectified and the peak closest to the center was found. Each segment was individually normalized so there was no need to normalize the central peak. Since the signal bandwidth is narrow, the STA is approximately sinusoidal, with *Power* ≅ *Peak*^{2}/2. Therefore the square of the normalized peak amplitude is the definition of SFC for the FB-SFC method and is analogous to the power ratio from the PS-SFC definition

### Derivation of bias correction

The PS-SFC estimation method requires that the variance (power) must be estimated for each frequency component in the STA. However, this derivation requires the assumption that the power of a signal can be measured exactly. Fortunately, the bias correction is a robust procedure in that this assumption does not have to be met for the correction to perform well. This is demonstrated in the results, where the bias is reduced when the power is only estimated. Although this derivation is based on the variance of the STA, it applies equally well to the FB-SFC definition of coherence, since ideally the filtered STA has variance equal to the peak amplitude squared.

The variance in an STA can be thought of as coming from two independent sources: the variance of the signal that is time-locked to the spikes and the variance of the uncorrelated or random noise present in the LFP. For a given frequency, the *j*th LFP segment in an STA can be represented as
*x* is the coherent signal component and *y* is the uncorrelated noise component. A stationary coherence is assumed, in that the coherent component of each LFP segment is the same. The parameter we seek to estimate is the asymptotic coherence *c*. In this formalism, we can define asymptotic coherence as the fraction of variance in the time-locked signal *x*, contained in the variance of the total signal *x* + *y*
*A*, constructed from *n* spikes, is defined as
*Eq. 3*. The bias in this estimator becomes apparent as the terms are expanded; *x* and *y* are independent variables by definition so the variance of their sum is equal to the sum of their variances
*Eq. 5* as the variance of a sampling distribution of means. Since *y* has no dependence on the spike times, *y*_{1}, *y*_{2}, … , *y _{n}* are independent. Therefore

*Eq. 2*we know that the second term in

*Eq. 7*can be expressed in terms of

*c*

*Eq. 7*, we have

*c*gives the final form of the bias correction

*n*spikes and the asymptotic coherence. Of course, in a practical application of this estimation method, there are inevitable deviations from the assumptions. This causes variability between the estimated and asymptotic coherence. However, our results demonstrate that when coherence is estimated over multiple trials, the trial-averaged estimate satisfies this relationship, indicating a reduction in bias.

### Surgery, data collection, and spike sorting

Experimentally collected data were recorded from rats chronically implanted with tetrode bundles. All procedures were approved by the Institutional Animal Care and Use Committee at Drexel University and followed National Institutes of Health guidelines. Three male Long Evans rats (300–350 g) were given atropine (0.05 mg/kg), anesthetized with isoflurane (2%), and placed in a standard stereotactic apparatus for electrode implantation. A craniotomy was performed and a bundle of eight tetrodes were slowly lowered into the CA3 region of the dorsal hippocampus and sealed in place with dental acrylic. Coordinates from bregma were 3.5 mm posterior, 2.0 lateral, and 3.5 mm ventral. Tetrodes consisted of four HVF-coated tungsten wires, 12 μm in diameter (California Fine Wire, Grover Beach, CA), twisted together. Recordings were performed while rats were awake and freely moving. Signals were filtered to a wide band (0.8 Hz to 6 kHz) and amplified at a gain of ×600 using a 31 channel wireless headstage (Triangle Biosystems, Durham, NC). These signals were then sampled at 40 kHz and stored on a computer for off-line analysis. All further filtering was performed digitally with zero-phase filters (Matlab, filtfilt.m) to preserve timing between spikes and LFP. LFP data were obtained by low-pass filtering the raw signal at 500 Hz and down-sampling to 2 kHz. To discriminate single units, the raw signal was high-pass filtered at 400 Hz and then 800 μs segments were extracted wherever negative peaks exceeded 5SD of peak heights. From these waveforms, the energy and the projections on the first three principal component axes were used as features for automatic cluster cutting (MClust-3.5A, A. D. Redish et al., University of Minnesota, and KlustaKwik-1.7, K. D. Harris, Rutgers University). Clusters were then manually adjusted on the basis of waveform shape, autocorrelation, and cross-correlation between clusters.

## RESULTS

Both the PS-SFC and the FB-SFC estimation methods were applied to the simulated data sets described in methods and our experimentally collected data. The simulated data sets were divided into categories based on the LFP model and the spike distribution model. The LFP was either modeled as a sine wave or as an experimentally collected LPF signal filtered to a narrow band. Spikes were distributed by either a binomial distribution (correlated or not) or a Gaussian distribution around a preferred phase. Combinations of these categories gave four simulation platforms to assess these methods. In addition, we tested our bias correction on a set of experimentally collected spikes and LFP signals. Estimation methods were judged based on bias, variance, and MSE criteria.

### Sine-wave LFP with binary spike distribution model

When the PS-SFC and FB-SFC methods were applied to simulated data consisting of a sine wave with spikes placed by a binary distribution, the corrected coherence showed a reduction in bias but an increase in variance (Fig. 1). As expected, as the number of spikes increased, the mean uncorrected coherence estimate approached the asymptotic coherence and the variance of the estimate decreased (Fig. 1*A*). Since the LFP is represented as an ideal sine wave, there is no difference between the asymptotic coherence and the simulated coherence. When the correction was applied, the mean estimated coherence was consistent with the asymptotic coherence for all numbers of spikes, indicating a removal of bias (Fig. 1*B*). The bias removal was confirmed in Fig. 1*C* (*inset*), but as a consequence of the correction, the variance of the estimate was increased. Since the SFC from one spike is trivially equal to 1, the bias will be greater for estimates of lower asymptotic coherence values. This effect is evident from the difference in bias between the two coherence strengths. As can be seen most clearly in the case of <10 spikes, the bias of the uncorrected estimator for γ = 0.25 was lower than the bias for γ = 0.01. In effect, the bias correction trades bias for variance. The overall direction of this trade, however, can be summarized by the mean squared error (MSE) measurement. Importantly, the MSE was either reduced or unaffected in all cases except the case of both very few spikes and high coherence. For the majority of cases in this simulation, applying the correction reduced the bias of the coherence estimator and lowered the error, regardless of method. Because the difference between power and peak height is constant for a perfect sine wave, there was no difference between PS-SFC and FB-SFC estimates and therefore the FB-SFC results were not shown for this simulation. In the following text, we compared these conclusions to a more realistic spike distribution model.

### Sine-wave LFP with Gaussian spike distribution model

When a Gaussian distribution was used to place spikes with respect to a sine wave (Fig. 2), the SFC approached its asymptotic value with the same trend as that of the binary spike distribution for both the PS-SFC and FB-SFC methods. This is the trend predicted in the derivation of the bias correction. Again, the bias was decreased and the variance was increased when the correction was applied (Fig. 2*C*). The results are qualitatively similar to those of the case of binary spike distribution, with the exception that the variance in both uncorrected and corrected estimates was slightly reduced. This is likely explained by the fact that a Gaussian distribution allows for finer gradients of coherence when there are fewer numbers of spikes. Changing the spike distribution also had no effect on the agreement between the asymptotic and simulated coherence. Finally, the MSE was reduced except in the case of strong coherence and few spikes, similar to the binary distribution model (Fig. 2*D*). Similarly to the previous sine-wave model, the methods produced identical estimates and so only the PS-SFC results are shown.

### Experimentally collected LFP models

In a more realistic simulation, replacing the sine wave with a real LFP signal, the MSE was reduced for all conditions simulated for the PS-SFC method (Fig. 3). Data shown are for a Gaussian spike distribution model, although results were quantitatively similar for the binary distribution (data not shown). The first noteworthy feature is that the asymptotic coherence is lower than the simulated coherence. Due to this decrease, the mean values of coherence estimates for all numbers of spikes were reduced (Fig. 3*A*). In addition, for the same reason that there is greater bias for lower asymptotic coherence, the lower asymptotic coherence in this simulation increased the bias of the uncorrected estimator (Fig. 3*C*). The variance of uncorrected and corrected estimates was reduced compared with the sine-wave models. Again, correcting the estimate reduced the MSE, but to an even greater degree than that for sine-wave models. This was due to two factors: *1*) an increase in the bias of the uncorrected estimator made a greater change after applying the correction and *2*) the lower variance in the uncorrected estimate reduced the variance in the corrected estimate. For this more realistic simulation, the MSE was either reduced or unaffected under all simulated conditions after the correction was applied.

When the FB-SFC was applied to the real LFP simulation (Fig. 4), the asymptotic coherence was slightly less than the simulated coherence. This was similar to the PS-SFC result, except less pronounced. Apart from this aspect, there were no noticeable differences from when the FB-SFC was applied to the sine-wave simulations. An additional case of very high coherence (γ = 0.5) was added, to demonstrate that the bias is removed after correction for a wider range of coherence strengths. Because the bias is not as great for this highest level of coherence, correcting the bias increases MSE for <20 spikes. It should be noted that when applying either method to this model (Gaussian distributed spikes around the peaks of a well-filtered, experimentally derived LFP), the asymptotic coherence did not exactly equal the simulated coherence. This is due to the way spike field coherence was modeled in the simulations. Coherence simulated by placing spikes in a random distribution around a phase of a filtered signal is more accurately reflected by the FB-SFC method than by the PS-SFC. This is not caused by a deficiency of the PS-SFC method, but because its output is influenced by factors other than the parameters used to simulate coherence. We expand on this subject in the discussion.

### Experimentally collected neural data

Finally, applying the bias correction to real data confirms the results seen in the simulations. Of the 54 neurons recorded from three rats, 9 neurons had enough spikes that the asymptotic coherence could be estimated with an uncertainty of <1%. SFC was then estimated again from a random subset of spikes from the same neurons. For each number of spikes, the squared error of the estimated coherence was found with respect to the asymptotic coherence. The mean squared error was then found by repeating this process over 1,000 trials. Each estimate was also bias corrected and the MSE of the bias corrected estimates were compared with the noncorrected estimates (Fig. 5). When the PS-SFC method was applied to the selected neural data, correcting the bias reduced the MSE for all numbers of spikes. However, when the FB-SFC method was used, MSE was reduced or unchanged for all numbers of spikes except very small numbers (*n* <3). The reason the bias correction is not as effective in reducing MSE for the FB-SFC method is that it has a higher variance for the estimation parameters chosen. As the number of spikes in the estimate increases, the bias correction does not have as large an effect on each neuron and the MSE is not reduced as much. The reduction to the MSE for >20 spikes appears small compared with the reduction for <10 spikes (Fig. 5*A*). However, normalizing these values to the asymptotic coherence gives a sense of the relative size of this error. The percentage change reveals that the bias correction can still reduce the MSE by >5% of the asymptotic coherence value, even when 50 spikes are used (Fig. 5*B*, *inset*).

## DISCUSSION

### Bias correction

Through the course of applying the SFC estimation methods to simulated and experimental data, it was demonstrated that these methods overestimate the SFC and the degree of this bias is a function of the number of spikes used and the asymptotic SFC. Of course, given enough data, many of the issues involved in coherence estimation are no longer of practical concern, including SFC bias. However, we cannot assume that this will always—or even often—be the case. It is when the limits of the estimator are pushed that the quality of estimation needs to be examined. There certainly exist experimental paradigms where it is difficult to collect large numbers of spikes; one example is if SFC is to be estimated over short time intervals. This could occur if one is studying coherence during transient events (i.e., stimuli, ripples, or epileptic seizures) or if coherence is suspected of changing rapidly with time. Compounding the problem are neurons that fire with a very low rate. Pyramidal cells with a firing rate of <1 Hz are not uncommon in the hippocampus. In cases such as these, the bias has an appreciable effect on the estimate.

The proposed bias correction has been shown to have two major effects on the quality of the coherence estimate. In general, it reduces bias while adding variance to the estimate. This can result in either an increase or a decrease in the MSE, depending on the asymptotic coherence and the number of spikes in the STA. For the more realistic simulated models and the experimental data, the results from the conditions tested showed that the MSE is usually either reduced or unaffected, except in the case of high coherence (γ >0.25) and low numbers of spikes (*n* <4) for the FB-SFC method. We and others (Foffani et al. 2007; Fries et al. 1997) have found coherence levels are typically lower than this in neuronal data, even for well-synchronized cases. Therefore it seems unlikely that the correction would increase the error in estimation.

The bias correction had the greatest effect when the number of spikes was <100 and, beyond this, its effect was less pronounced. However, this does not mean the bias correction should not be used in cases of >100 spikes because when asymptotic coherence is low, results may still be biased. Consider a comparison between two conditions, both of zero actual coherence. In either condition, the uncorrected estimate will have an expected value of 1/*n*. If in one condition, 200 spikes are available in the observation window, the expected estimated coherence will be 0.005. For some this might be an acceptable error. However, if the other condition has only 100 spikes available in the observation window, the expected value of the uncorrected estimate will be 0.01. The result is that one estimate is twice as great as the other, even though the coherence is the same in both cases. Since the bias correction does not increase the MSE for cases of large numbers of spikes, it can be beneficial even for cases of >100 spikes.

Of course, other methods of estimating unbiased SFC-related measures have been proposed. Since the degree of bias depends on the number of spikes, the simplest approach is to fix the number of spikes used for each estimate. The obvious shortcoming of this technique is that the number of spikes in each estimate must be limited to the estimate with the fewest spikes. This means data must be disregarded from the analysis and that some estimates will have more variance than they would have had with all available spikes. Bootstrapping techniques (Vinck et al. 2010a) can address both of these problems, but there are still issues with dependence on spike distribution, dependence on parameter choice, and computational feasibility. One recent measure (Vinck et al. 2010b), called the pairwise phase consistency (PPC), addresses the bias problem by effectively fixing the number of spikes at two and averaging over the quantities calculated for all possible pairs of spikes in the sample. The expected value of the PPC was shown to be equivalent to the asymptotic limit (*n* → ∞) of a widely used and biased measure, the phase-locking value. Although this is a well-formulated solution to the bias problem of the phase-locking value, it is unknown whether it can be applied to other SFC-related measures. It would be interesting to see whether a technique similar to the PPC can be developed for unbiased estimation of PS-SFC and FB-SFC and how it would compare with the bias correction proposed herein.

### Asymptotic and simulated coherence of PS-SFC and FB-SFC

For some simulations, the FB-SFC and PS-FSC methods had asymptotic coherences values different from the simulated coherence. On the one hand, using a sine wave as a model of LFP, the asymptotic coherence for both estimation methods was equal to the simulated coherence. On the other hand, when using LFP and simulated spikes, the asymptotic coherences were different from each other and different from the simulated coherence. Finally, when applied to real data, there were differences in the asymptotic coherence estimates between the two methods as well. These discrepancies are largely due to subtle differences in the quantities these two methods measure and to how those quantities change when the data are no longer ideal. The FB-SFC measures the peak amplitude of an STA constructed from a filtered signal and the PS-SFC measures the power spectrum of the unfiltered STA. Outside of an ideal sine wave, these two quantities, although similar, are not necessarily equivalent. Therefore there will be a predisposed tendency for these two methods to produce different estimates. Additionally, it was observed that the FB-SFC asymptotic coherence was much closer to the simulated coherence than the PS-SFC. This is most likely because in these simulated data, the spikes were placed with respect to the phase of a filtered LFP signal, which is the same signal the STA is constructed from in the FB-SFC method. Essentially, the FB-SFC has the advantage of being able to more directly reflect the simulation parameters that control coherence strength. Had the coherence been simulated in a different way, this probably would not be the case.

### Factors that influence SFC estimation

As with any estimation method, it is important to be aware of which controllable factors influence the end result in addition to the parameter being estimated. It has been demonstrated that the number of spikes used to estimate the coherence causes a bias in this estimate, but even if a very large number of spikes is used, other parameters in the estimation process can still cause different kinds of bias. For the PS-SFC method, the size and shape of the STA window together probably constitute the largest factor that can affect the estimate. With an ideal sine wave, if spikes occur at the peak of the oscillation, they are just as correlated with a peak that happened many cycles ago or many cycles in the future. This will cause an STA created from sine waves to have the same average amplitude at any time lag from the center (spike). Since in real systems spikes are coherent with oscillations that are not perfectly periodic, there is an attenuation of STA power with increasing time lag from the center. This is to be expected because there is no reason a spike should be correlated with a cycle of an oscillation that happened far away in time. This general phenomenon, which occurs in many biological systems, is a consequence of the mixing condition (Brillinger 1981; Rosenblatt 1956). Therefore the closest estimate of the coherence is at the center of the STA. The practical implication is that the STA window in which the PS-SFC is calculated will influence the estimated coherence. Compounding this issue is the fact that the degree of power attenuation is frequency dependent. Coherent high-frequency oscillations attenuate within much shorter time windows of the STA center than do low-frequency oscillations (Fig. 6*A*). It follows that for each frequency there is an optimal time window with which to estimate the STA power in that frequency. The inverse of this statement is that the power spectrum for a given STA window will automatically enhance a certain frequency range while depressing surrounding frequencies (Fig. 6*B*). Furthermore, this problem is not solved by the multitaper technique with discrete prolate spheroidal sequences since additional tapers do not emphasize the exact center of the STA. Therefore the STA window size should be chosen to optimally estimate a specific frequency band when using the PS-SFC method. Comparisons between different frequency bands should require different STA windows.

Instead of using the entire STA window to estimate coherence, the FB-SFC method estimates it from the peak amplitude at the center of the STA. Since this estimate comes from a single point instead of a range, this method may produce more uncertainty in the estimate. This accounts for the increased variance from the FB-SFC method compared with that from PS-SFC. However, this feature also allows the filter bank method to estimate the coherence fairly between all frequency bands since it is not dependent on the STA window (Fig. 6*D*). The major factor that influences the estimate of coherence in this case is the fractional bandwidth of the set of filters. This is defined as the bandwidth of a filter divided by the central frequency. The FB-SFC method relies on the peak amplitude of the filtered signal being proportional to the power. This is true of a pure sine wave, but this proportionality breaks down as the bandwidth of the signal increases. Therefore larger fractional bandwidths will increase the overall value of coherence across the spectrum, while decreasing the frequency resolution. Although there is a small difference in the ratio of the spectral peaks at 6 and at 150 Hz as the fractional bandwidth increases, we attribute this to a loss of frequency resolution. In fact the peak at 6 Hz is noticeably wider and less well resolved at the highest fractional bandwidth (0.67) than that at the other fractional bandwidths. This is the fractional bandwidth that results when each band is double the bandwidth of the last band (i.e., 2–4 Hz, 4–8 Hz, 8–16 Hz, etc.) and is a common approximation of the conventional rhythms of the brain. However, when examining the broad spectrum, it is usually better to have a higher-frequency resolution since coherence can sometimes occur in narrow bands. Therefore a fractional bandwidth of <0.33 is a better choice.

### Appropriate application of methods

Our results suggest that the estimation methods examined here are influenced by factors other than the actual coherence. The major factor that is common to both methods is the number of spikes, which causes a positive bias in the SFC estimate. In addition, each method has controllable parameters that can affect the estimate. With this in mind, it is important to know the experimental circumstances under which these methods can be appropriately applied.

Probably the most common experimental paradigm is to compare a single neuron's SFC in a specific frequency band between two conditions. In this case, the minimum comparison needed is a single SFC estimate for each condition. This type of comparison could potentially be made by using a specialized statistical test, such as an adapted version of the nonparametric test for LFP–LFP coherence (Maris et al. 2007). However, it may be desirable to estimate a population-average SFC instead of single-neuron SFC and even necessary if the same neurons cannot be recorded in both conditions. Pooling trials (spike-centered LFP segments) across neurons as a possible application of the nonparametric test would be inappropriate due to differing preferred phases and firing rates of the neurons. A more logical approach is to consider the coherence values estimated from individual neurons as observations and the comparison is between a sample of observations in one condition and a sample in the other condition. The problem is that when these observations are biased estimates, no statistical test can separate the effect of the number of spikes from the coherence. The bias must be removed before this type of comparison can be made and this is solved by the bias correction.

Sometimes it is unknown which oscillation frequencies a neuron will be coherent with a priori and an examination of SFC across a broad spectrum is required. Even if statistical comparisons between frequency bands are not experimentally necessary, it is useful to know at what frequencies coherence is actually occurring and this can be accomplished only by eliminating the frequency-dependent bias in the estimate. The PS-SFC method, if used incorrectly, suffers from this type of bias as discussed in the previous section. Therefore a fair comparison of frequency bands requires recalculating the estimate for multiple STA window sizes. Alternatively, the FB-SFC method can be used for this type of comparison. Since this method does not depend on the STA window size, coherence estimates can be made fairly across the entire spectrum.

In conclusion, it is important to remember that in real conditions the SFC can only be estimated and there is always uncertainty in this estimate. A very large number of spikes will usually give an estimate with little bias, but this type of data set is not always available, as discussed earlier. In addition, it is difficult to tell just how many spikes are needed when coherence is low. The proposed bias correction addresses these problems by severing the estimate's dependence on the number of spikes.

## GRANTS

This work was supported by a grant from the Coulter Foundation.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## ACKNOWLEDGMENTS

We thank G. Foffani for helpful discussions and the Drexel University Major Research Initiative in Neuroengineering.

- Copyright © 2010 The American Physiological Society

## REFERENCES

- Bendat and Piersol, 1971.↵
- Brillinger, 1981.↵
- Courtemanche et al., 2002.↵
- Donoghue et al., 1998.↵
- Eggermont and Smith, 1995.↵
- Foffani et al., 2007.↵
- Fries et al., 2001.↵
- Fries et al., 1997.↵
- Fries et al., 2008.↵
- Goldberg et al., 2004.↵
- Gray and Singer, 1989.↵
- Harris et al., 2002.↵
- Lee et al., 2005.↵
- Maris et al., 2007.↵
- Mendell and Henneman, 1968.↵
- Mendell and Henneman, 1971.↵
- Murthy and Fetz, 1996.↵
- O'Keefe and Recce, 1993.↵
- Quyen et al., 2008.↵
- Rosenblatt, 1956.↵
- Vinck et al., 2010a.↵
- Vinck et al., 2010b.↵