## Abstract

Single neurons carry out important sensory and motor functions related to the larger networks in which they are embedded. Understanding the relationships between single-neuron spiking and network activity is therefore of great importance and the latter can be readily estimated from low-frequency brain signals known as local field potentials (LFPs). In this work we examine a number of issues related to the estimation of spike and LFP signals. We show that spike trains and individual spikes contain power at the frequencies that are typically thought to be exclusively related to LFPs, such that simple frequency-domain filtering cannot be effectively used to separate the two signals. Ground-truth simulations indicate that the commonly used method of estimating the LFP signal by low-pass filtering the raw voltage signal leads to artifactual correlations between spikes and LFPs and that these correlations exert a powerful influence on popular metrics of spike–LFP synchronization. Similar artifactual results were seen in data obtained from electrophysiological recordings in macaque visual cortex, when low-pass filtering was used to estimate LFP signals. In contrast LFP tuning curves in response to sensory stimuli do not appear to be affected by spike contamination, either in simulations or in real data. To address the issue of spike contamination, we devised a novel Bayesian spike removal algorithm and confirmed its effectiveness in simulations and by applying it to the electrophysiological data. The algorithm, based on a rigorous mathematical framework, outperforms other methods of spike removal on most metrics of spike–LFP correlations. Following application of this spike removal algorithm, many of our electrophysiological recordings continued to exhibit spike–LFP correlations, confirming previous reports that such relationships are a genuine aspect of neuronal activity. Overall, these results show that careful preprocessing is necessary to remove spikes from LFP signals, but that when effective spike removal is used, spike–LFP correlations can potentially yield novel insights about brain function.

## INTRODUCTION

Much of our mechanistic understanding of brain function comes from extracellular recordings, which provide a direct measure of electrical activity near the tip of the recording electrode. The resulting voltage signal is composed of the spikes emitted by one or more neurons and the local field potentials (LFPs), which represent the total synaptic current in the neuronal circuit from which the recordings are obtained. Typically, the LFP is extracted by low-pass filtering the wideband voltage signal, whereas spikes are identified by high-pass filtering, thresholding, and subsequent sorting. The possibility of isolating both signals from the same electrode is of great interest because it allows the experimenter to relate the responses of individual neurons to those of the larger-scale circuits in which they are embedded.

Popular methods for studying spike–LFP relationships in the time domain include the spike-triggered average (STA) of the LFP, the spike-field coherence (SFC), spike–LFP phase-locking histograms, and the prediction of spike timing from LFP features. These types of analysis have proven useful for inferring the role of feedforward and feedback circuitry in functions as diverse as perception, attention, memory, and motor control (Andersen et al. 2004; Chalk et al. 2010; Destexhe et al. 1999; Fries et al. 2008; Jacobs et al. 2007; Pesaran et al. 2002; Saleh et al. 2010). Similar methods applied to the space domain permit comparisons of single-neuron activity to columnar systems in the cortex (Nauhaus et al. 2008; Xing et al. 2009). More generally, comparisons of spiking and LFP activity facilitate comparison with other signals, such as blood oxygenation level dependent (BOLD) and electroencephalographic (EEG) activity, that can be obtained noninvasively (Goense and Logothetis 2008).

Of course, before attempting to analyze the relationship between two types of signals, one must ensure that they can be estimated independently. LFP estimation often relies on the fact that the frequency content of LFPs and individual spike waveforms differ substantially, such that LFPs can be isolated by simply low-pass filtering the voltage signal. However, to the extent that the LFPs and spikes contain overlapping frequencies, this approach will fail to completely separate the two signals. The resulting *spectral contamination* is of particular concern when estimating the causal relationships between LFPs and spikes because contamination of one signal by the other will result in spurious correlations. Spectral contamination is most problematic when spikes and LFPs are estimated from the same electrode recordings, but it can also affect recordings obtained from separate electrodes, particularly when the corresponding spike density functions are highly correlated. Although various algorithms for efficient removal of spikes from LFPs have been proposed (David et al. 2010; Galindo-Leon and Liu 2010), a systematic appraisal of the consequences of this contamination in common analytical settings has not yet been attempted.

In this work we use simulations to estimate the extent of the spectral contamination of LFPs by individual spikes and spike trains, under conditions similar to those of typical extracellular recordings. We show that such contamination introduces artifactual spike–LFP relationships into many of the popular time-domain metrics mentioned earlier. In particular, strong spike–LFP coherence and spike-triggered LFPs can be obtained from random combinations of spikes and LFPs, if the former are not removed properly. We confirm these simulation results with real recordings from the primate visual system.

To address the issues associated with spike contamination we develop a novel Bayesian method that estimates veridical LFPs by removing action potential signatures from raw voltage signals. We show that application of this method to both simulations and real data has a substantial impact on measures of the relationships between spikes and LFPs. We also test two other spike removal methods and show that the Bayesian algorithm performs better in most cases. Overall, these results show that the way in which electrode signals are processed can significantly affect the conclusions drawn about the nature of spikes, LFPs, and interactions between these two signals.

## METHODS

### Electrophysiological recordings

##### EXPERIMENTAL SETUP.

Two rhesus macaque monkeys took part in the experiments. Both underwent a sterile surgical procedure to implant a headpost and recording cylinder and, after recovery, monkeys were seated comfortably in a primate chair (Crist Instruments) and trained to fixate a small spot on a computer monitor in return for a liquid reward. Eye position was monitored at 200 Hz with an infrared camera (SR Research) and required to be within 3° of the fixation point for the reward to be dispensed. Recordings were obtained from well-isolated single units in areas MT/V5 and V6, both of which were identified based on anatomical magnetic resonance imaging scans and the physiological properties of the neurons (Fattori et al. 2009; Maunsell and Newsome 1987). All aspects of the experiments were approved by the Animal Care Committee of the Montreal Neurological Institute and were conducted in compliance with regulations established by the Canadian Council of Animal Care.

##### PROCEDURE AND VISUAL STIMULI.

On each trial, the animal acquired fixation, after which the stimulus appeared and remained stationary for 200 ms. The stimulus then moved at a constant direction and speed for 500 ms. Stimuli were displayed at 85 Hz at a resolution of 1,920 × 1,200 pixels and the viewing area subtended 70 × 42° of visual angle at a distance of 42 cm. Stimuli consisted of sinusoidal gratings of optimal spatiotemporal frequency or random-dot patches displayed on a gray background (luminance of 70.3 cd/m^{2}). Stimulus size was optimized for each cell. Motion direction was sampled in 30° steps and each stimulus was repeated five times in blockwise random order.

### Signal processing

Electrophysiological signals were recorded with a standard data acquisition system (Plexon Multichannel Acquisition Processor [MAP] System). Specifications of this system typically include a dedicated low-pass filter for LFP signals (four-pole high cut at 170 Hz) followed by digitization at 1 kHz. As we show in the following text, this type of filtering invariably introduces spurious correlations into measures of spike–LFP coherence. We therefore obtained custom hardware modifications to the MAP system that included wideband analog filters (two-pole high cut at 2.5 kHz) and a higher digitization rate (10 kHz). These modifications allowed us to detect spectral contamination of the LFPs by spikes; subsequent spike sorting and LFP analysis were performed through off-line digital filtering.

To remove line noise at 60 Hz, we assumed that the underlying signal had equal power in a small range of frequencies around 60 Hz and filtered the signal so that the power at 60 Hz was similar to that of surrounding frequencies. We transformed the wideband signal into the Fourier domain and minimized the squared error between the absolute value of the Fourier coefficients between 55 and 65 Hz and the function
*f* is frequency, *f*_{0} is the peak frequency (nominally 60 Hz), *a* is the height of the 60 Hz peak, *b* controls its width, *c* adjusts the shape of the peak, and *d* is the baseline frequency content. This functional form provides a parsimonious parameterization of the frequency content around 60 Hz. We then corrected the 60 Hz artifact by filtering the signal in the Fourier domain with the function's inverse

### Simulations

To assess the likely contribution of spectral contamination to various measures of spike–LFP coherence, we analyzed a signal that, by construction, contained no temporal relationship between spikes and LFPs. This *composite* signal was obtained by taking the frequency content of the wideband signal obtained from our electrophysiological recordings, after removing the spikes using the spike removal algorithm presented in the following text, and shuffling the phases (Fig. 1*A*). To this signal we added back waveforms taken at random from the original electrophysiological signal, preserving the timing of each spike while scaling the waveforms to achieve a desired signal-to-noise ratio (SNR). To obtain reliable estimates of the underlying spike waveforms, we first selected typical waveforms that were within 3SDs of the mean waveform at all times. We then projected this ensemble of waveforms onto their first five principal components, obtaining a bank of denoised waveforms. To achieve a desired mean firing rate, Poisson distributed spurious spikes were added to the signal.

The LFP components of this signal thus had the same amplitude spectra as those of a real neural signal, but lacked any consistent relationship between spike times and LFP phases. Consequently, any apparent relationship between the two signals would necessarily be artifactual.

### Data analysis

In an initial preprocessing step, we corrected the phase distortions introduced by analog filters by filtering the wideband signal backward in time with a digital Bessel filter that closely modeled the phase response of the analog filters used in the preamplifier (Nelson et al. 2008). We then filtered the wideband signal with digital filters similar to those used conventionally to estimate LFPs. The filters used were a fourth-order 170-Hz low-pass causal Butterworth filter (digital version of the filter used on the Plexon PreAmp separation for the LFP board) and a fourth-order two-pass Butterworth (same filter with zero phase delay). Filtering was implemented with Matlab software (The MathWorks).

For some analyses, we divided the LFP frequency spectrum into the following bands, similar to those used in previous studies (Khawaja et al. 2009; Ray et al. 2008): θ (4–8 Hz), α (8–12 Hz), β (16–24 Hz), low γ (25–55 Hz), and high γ (65–140 Hz). Since the θ, α, and β bands did not show significant spike contamination in any of our metrics, in both simulations and real data, for reasons of compactness we combined all three bands into a single “low-frequency” band (4–24 Hz).

Spike detection and sorting based on the high-pass filtered wideband signal were performed off-line using established methods (Quiroga et al. 2004). These results were compared for consistency with results from Plexon's Off-line Sorter software, which uses spike waveforms sampled at 40 kHz. Spike timing information was considered accurate if agreement between the results of the two spike sorting methods was above a criterion level of 95%; otherwise, the sorting parameters were altered until the agreement exceeded the criterion. In three cases where this level could not be reached, the data were rejected from further analysis.

SNRs were computed for each isolated unit as the peak-to-trough voltage of the mean spike waveform, divided by the square root of the mean power of the complete raw signal. We obtained the STA of the LFP by averaging local field potentials centered around the time of spikes (Dayan and Abbot 2001). We subsequently filtered the LFP with complex Morlet wavelets and used the magnitude of the coefficients as an estimate of the instantaneous power of the LFP at given frequencies (Goupillaud et al. 1984). We averaged the magnitude of the coefficients around the time of spikes, normalized by the average magnitude of the coefficients of the whole LFP signal, to obtain a normalized wavelet STA, which can reveal transient increases or decreases in the magnitude of oscillations at given frequencies around the time of spikes.

To test whether spikes preferred certain phases of the LFP, the instantaneous phase and amplitude of the LFP at 35 different frequencies (log spaced from 4 to 150 Hz) were first calculated with the Hilbert transform of different filtered versions of an LFP recording. A single phase value was then assigned at each frequency for every spike in a recording (Jacobs et al. 2007). By convention, a phase value of 180° corresponds to the positive (hyperpolarizing) peak of the cycle. Finally, histograms of the distributions of spike phase for different LFP frequencies across all spikes were compiled; 50 phase bins, ranging from 0 to 360° were used. The uniformity of the distribution of spike phase values for our frequency ranges was assessed for each case, using the Rayleigh statistical test, corrected for the number of frequency ranges (Rutishauser et al. 2010). A cell was considered phase-locked at a specific frequency range if the null hypothesis of uniformity of the phase distribution could be rejected at a corrected *P* < 0.01.

Spike-field coherence (SFC) was calculated as the power spectrum of the LFP STA divided pointwise by the sum of the power spectra of all LFP segments used to compute the STA (Fries et al. 2001). This measure can reveal phase-locking at high frequencies that is otherwise obscured in the STA of the LFP by low-frequency noise. A value of 1 at a given frequency of the SFC means that all spikes appear at the same phase for this frequency, whereas a value of 0 indicates no consistent relationship between the phase of the LFP at a given frequency and spikes.

Tuning curves were quantified for single units (spike count) and LFPs (normalized LFP band power), in the three frequency bands described earlier. We used Morlet wavelets to measure the instantaneous power at different frequencies of the LFP across time, as in the wavelet STA analysis. LFP tuning curves were derived by normalizing the mean magnitude of the Morlet coefficients over a time period that spanned 100–500 ms after the onset of stimulus motion against that in a baseline period consisting of 150 ms of spontaneous activity before the stimulus presentation. Average tuning vectors for each band were also computed for both single units and LFPs.

### Spike removal algorithm

The purpose of the algorithm is to estimate the local field potential based on a measured wideband voltage trace of length *n*. We assume that this wideband signal **y** is the superposition of a low-frequency local field potential **w**, high-frequency spike components **η*** ^{k}*, an offset μ, and white noise

**ε**(Fig. 1

*B*).

*m*is the number of sorted neurons emitting spikes. The high-frequency component of the

*k*th neuron

**η**

*is created by the convolution of the neuron's spike train*

^{k}**s**

*, assumed known (see Fig. 1*

^{k}*B*), and the neuron's spike waveform

**ϕ**

^{k}**ϕ**

*are not set to the mean waveforms determined by spike sorting or by spike-triggered averaging the wideband signal, as in previously proposed methods of eliminating spike remnants (Pesaran et al. 2002). Rather, the waveforms*

^{k}**ϕ**

*are considered as free parameters and are estimated optimally as a part of the algorithm. The offset μ is introduced to compensate for the fact the mean of*

^{k}**η**

*is likely to be negative.*

^{k}The generative model (*Eq. 3*) contains more unknowns (**w**, **ϕ*** ^{k}*, and μ) than knowns. Thus assumptions must be introduced to make the estimation of the local field potential

**w**well-posed. We first assume that the local field potential is smooth or, equivalently, that most of its power lies in the low frequencies

*N*(

**a**,

**Σ**) represents a multivariate Gaussian with mean

**a**and covariance

**Σ**. Capital gamma (

**Γ**) is a matrix that embodies our assumption of smoothness. Multiplication of a vector

**x**by this matrix (

**Γx**) produces a low-pass filtered version of the target

**x**. The matrix was chosen to be circulant, a property necessary to make estimation of the local field potential tractable; definitions and properties of circulant matrices are discussed at length in the Supplemental appendix.

^{1}Gamma (γ) controls the strength of the prior. Our second assumption is that

**ε**is generated by a white noise process,

*p*(

**ε**) =

*N*(0, σ

^{2}

**I**). Finally, we assume that spike waveforms

**ϕ**

*lie in a subspace*

^{k}**B**

**B**so that spike waveforms were described by one parameter per time sample in the interval from −1 to +2 ms relative to the peak depolarization. This duration of 3 ms was chosen based on the traditional length of snippets used in spike sorting, generally ranging from 1.5 to 3 ms.

To estimate the local field potential **w**, the spike waveforms **ϕ*** ^{k}*, and the offset μ, we used Bayesian inference to obtain maximum a posteriori (MAP) model parameters. By Bayes' theorem, the log-posterior of the model is

Here *k* is a constant factor. By taking the log of this expression and setting partial derivatives with respect to the parameters to 0, we obtain the MAP estimates of the parameters
Here **A**^{+} = (**A**′**A**)^{−1}**A**′ denotes the Penrose–Moore pseudoinverse. These equations can be understood as follows. First, the optimal offset of the model *μ̄* is simply the mean of the wideband signal with the LFP and spike contributions removed. Second, the optimal local field potential **w̄** is related to the spike-free, mean-corrected wideband signal **z** = **y** − ∑_{k}**s*** ^{k}**(

**Bφ̄**

*) −*

^{k}*μ̄*by matrix multiplication with (γ

^{2}

**Γ**+ σ

^{2}

**I**)

^{−1}γ

^{2}

**Γ**. Multiplying a vector by this matrix low-pass filters the vector and thus the model tells us that the optimal LFP is obtained by low-pass filtering the wideband signal after removing spike waveforms around the time of every spike.

The key equation is the second in the series, which governs the spike waveforms. By the properties of the pseudoinverse, (**s*** ^{k}**

**B**)

^{+}(

**s**

***

^{k}**B**)

**a**=

**a**. Now to recover a waveform

**a**from a train of waveforms (

**s**

***

^{k}**B**)

**a**, one simply needs to take the spike-triggered average of this signal normalized by the autocorrelation of the spike train. Thus the optimal spike waveform for the

*k*th neuron is given by the decorrelated STA of the wideband signal minus the LFP, the spike contributions from other neurons, and the offset. However, since the optimal LFP is mostly determined by low-pass filtering the wideband signal, the overlap between spikes from different neurons is often negligible and spike trains are uncorrelated at short time lags due to the absolute refractory period, it follows that to a first approximation the optimal waveform for a given neuron is equal to the STA of the high-pass filtered wideband signal.

In practice, the optimal spike waveform determined by the model can vary substantially from the STA of the high-pass filtered wideband signal because of the various corrections it applies. A fast and effective model-based spike removal algorithm can be derived by decoupling the *Eq. 8* system and applying the conjugate gradient descent to solve the resulting equations, as shown in the appendix. An implementation of this algorithm in Matlab is available from our website (http://apps.mni.mcgill.ca/research/cpack/lfpcode.zip). For the model to be fully specified, the strength of the prior relative to the noise γ^{2}/σ^{2} must be estimated because this determines what the model considers as signal and what it discounts as noise. In the appendix, we show how marginal likelihood (or evidence) optimization is used to estimate prior and noise levels. We also detail in the appendix how to choose an appropriate form for the prior matrix **Γ** and how to optimize performance for long recordings.

### Other spike removal algorithms

Two methods to remove spike waveforms from wideband signals were also tested in both simulated and real data. The first method, known as *mean spike removal*, simply subtracts the mean spike waveform from the wideband signal at the corresponding spike times. Here we used a 3-ms window (1 ms before the spike and 2 ms after) to recover the mean spike waveform. This was chosen for comparison with our Bayesian algorithm, which used a window of the same size. A second method uses linear *interpolation* to replace the wideband signal around the time of each spike, using the same time window as that in the other methods (3 ms).

## RESULTS

### Spectral composition of spike waveforms

Spikes are by definition fast events that are dominated by temporal frequencies far higher than those that make up the LFP signals. Thus it is natural to assume that low-pass filtering the raw voltage signal will remove the spikes and that what remains will be a pure estimate of the LFPs. This assumption can be examined by computing the power spectrum of a typical spike waveform and Fig. 2*A* shows one such waveform recorded extracellularly in our laboratory. As is typical in such recordings, the duration of the action potential is on the order of 1 ms and the corresponding waveform contains both a peak and a trough. Consequently, the power spectrum (Fig. 2*B*, *top*) is dominated by high frequencies. However, closer inspection of the power spectrum shows that there is a small amount of power at low frequencies in this fast waveform (as shown in Fig. 2*B*, *bottom*), particularly in some of the frequencies (100–200 Hz) that are often examined in LFP studies (Gold et al. 2006).

The *bottom panel* of Fig. 2*A* shows the result of applying a low-pass filter to the spike waveform shown in the *top panel*. Although this filtering reduces the amplitude of the waveform, some residual signal remains and this residual retains the basic shape of the original waveform. Moreover, the residual signal is stretched in time over an interval determined by the cutoff frequency of the low-pass filter (here on the order of 10 ms).

Beyond the low-frequency content of individual spike waveforms, the temporal structure of spike trains can contain additional power at low frequencies (Bair et al. 1994). Consequently, if spikes are not properly removed, a succession of spikes can contain more power at low frequencies than an individual waveform would indicate. These considerations suggest that the commonly used technique of low-pass filtering the raw voltage signal is not entirely sufficient to isolate the LFP signals. To estimate the magnitude of this problem, we developed a set of ground-truth simulations that allowed us to quantify the effects of spike contamination in a context where the veridical LFP signal is known.

### Spectral contamination in simulated data

The basis of the simulations was a phase-randomized wideband signal with the same frequency content as that of a real wideband signal (see methods and Fig. 1*A* for details). By combining this signal with spike waveforms taken from electrophysiological recordings (see following text), we obtained a *composite* signal similar to that typically recorded with electrophysiological methods but in which the LFPs and spike times were completely dissociated. This allowed us to perform ground-truth simulations in which any temporal relationship between spikes and LFPs was necessarily an artifact of the method of analysis. The duration of the phase-randomized signal in all our simulations was 3 min.

Figure 3 shows a short segment of the signals used in our ground-truth simulations. The *top panel* shows a phase-randomized wideband signal before (red line) and after the addition of spikes with a mean amplitude of 2 dB (blue line, offset downward to facilitate reading). The *middle panel* shows these same signals after low-pass filtering. Here the difference between the signals is much less visible, which suggests that in some analysis scenarios low-pass filtering may be sufficient to remove spike traces. The *bottom panel*, however, reveals that the deviations between the two signals (green line, scaled by a factor 3 to facilitate reading) are centered around the time of spikes (gray lines) and are highly stereotyped. These small spike remnants can bias analyses of LFPs, as we show in the next sections.

### Spike-triggered average (STA) of the LFP

As shown in Fig. 2, a typical spike waveform contains power at low frequencies, which can potentially complicate the isolation of the LFP component of the raw voltage signal. In particular, although the magnitude of the low-frequency component is relatively small (Fig. 3, *bottom*), its importance could be amplified by analytical techniques aimed at finding the components of the LFP that are correlated with spikes. One common example of such an analysis is the STA of the LFP.

The STA of the LFP is obtained by averaging the LFP signals recorded around the time of each spike. In principle this method is useful for inferring causal relationships between the two signals and it is commonly used to study local neuronal synchronization related to the discharge of a specific neuron and to compute the average spike-related modulation in specific LFP bands. The shape of the STA reveals the nature of such correlations, whereas a flat STA suggests that the two signals are unrelated.

We first examined the STA of a composite signal in which the SNR of the spike relative to the LFP was set to 2 dB and the mean spike rate was 9 spikes/s. Prior to calculation of the STA the raw composite signal was low-pass filtered to remove frequencies ≥170 Hz. However, inspection of the resulting STA in Fig. 4*A* (*top*) reveals clear residual structure, indicating that low-pass filtering is insufficient to separate spikes from LFPs. The artifact has a shape similar to that of the mean waveform, but has a longer duration, as expected from Fig. 2*A* (*bottom*). The exact shape of the artifact is determined by the properties of the low-pass filter.

For comparison we also computed the STA associated with a signal that was filtered with a single-pass Butterworth filter (Fig. 4*A*, *middle*), which is commonly implemented in analog circuits that are part of LFP recording hardware. In addition to the spurious STA, this filter also introduces frequency-dependent phase shifts, as reported previously (Nelson et al. 2008). Thus the commonly used approach of low-pass filtering the raw voltage signal appears to alter both the magnitude and timing of spike–LFP correlations.

The effects of this contamination are decomposed into different frequency bands in Fig. 4*B*. Here the STAs affected by spike contamination are shown in blue, whereas STAs obtained after the application of the spike removal algorithm introduced in the next section are shown for reference in red. The residual influence of spikes is present in all frequency bands and is particularly strong at low- and high-gamma frequencies.

Although the STA can detect linear relationships between spikes and LFPs, it is insensitive to higher-order correlations between these signals. The wavelet STA, however, can reveal transient changes in the relative instantaneous power in different frequency bands around the time of spikes, the signature of a more complex relationship between LFPs and spikes (Burns et al. 2010; also see methods). Figure 4*C* (*top*) shows that the wavelet STA is also sensitive to contamination by spikes. Here contamination is restricted to LFP frequencies ≥100 Hz at 0 ms time lag (deep red color). Note that the stripe pattern visible at low frequencies is due to chance fluctuations in the power at these frequencies and is not significant.

### Bayesian spike removal algorithm

To address the contamination issues illustrated in Fig. 4, we developed an algorithm that removes spikes from raw voltage traces. The algorithm assumes that the wideband voltage signal is generated by the sum of a low-frequency LFP component, spikes of limited duration occurring at known times, an offset, and white noise (see Fig. 1*B*, *bottom diagram*, and methods). These assumptions are embedded in a probabilistic model and the most likely LFP is then determined through Bayesian inference. The optimal parameters of the model have straightforward interpretations (see methods and the appendix): spike waveforms are essentially STAs of the high-pass filtered wideband signal, whereas the LFP is given by the low-pass filtered wideband signal after the subtraction of the contribution of spikes.

Figure 4*A* shows that the magnitude of the STA of the LFP when the composite signal is preprocessed through different methods. As mentioned earlier, a spurious STA is obtained following standard low-pass filtering (*top two panels*), whereas our spike removal algorithm greatly reduces this artifact (*bottom*). Figure 4*B* shows the STAs at different frequencies, with the LFP derived through standard low-pass filtering indicated in blue. The red lines show the STAs at different frequencies following application of our spike removal algorithm. The magnitude of the spurious correlation is now much smaller, particularly at high frequencies. Because the spike removal method removes the mean spike only at spike times, and thus will not remove all spike remnants when spike shapes for a given neuron are not completely stereotyped, some residual second-order correlations between spikes and LFPs could remain. Nevertheless, the method is effective at attenuating spurious correlations in the wavelet STA, as shown in Fig. 4*C*. Additional details on the derivation and implementation of the spike removal algorithm are given in the appendix, which also covers more challenging conditions, such as long recordings or nonstationary spike shapes that result from electrode drift.

### Phase synchronization

Another quantity commonly used to assess spike–LFP relationships is the synchronization between spikes and LFPs at different frequencies. This is typically manifested as a preference for spikes to occur at specific phases of the LFP oscillation cycle and such preferences may reveal synchronization of the local presynaptic activity. This synchronization can be an alternative way of encoding behaviorally relevant signals, particularly those related to visual attention (Fries et al. 2008), motor planning (Denker et al. 2007), motor output (Baker et al. 1999), and memory consolidation (Paz et al. 2008; Siapas et al. 2005). The frequency at which spikes are locked can also provide hints about the spatial and dynamic properties of the neural circuitry driving the neuron (Liu and Newsome 2006).

One method that is commonly used to assess the synchronization between spikes and LFPs is the spike-field coherence (SFC), obtained from the power spectrum of the STA. Figure 5*A* shows the SFC as a function of frequency for simulations in which the composite signal was preprocessed with low-pass filtering (blue) and our spike removal algorithm (red). When spikes are removed by low-pass filtering alone, prominent phase synchronization is seen at the higher frequencies and this artifactual relationship disappears following spike removal.

Figure 5*B* shows the phase-locking histogram, which measures the distribution of the instantaneous phase of the LFP at different frequencies at the time of spikes. Here it is clear that this measure is particularly sensitive to spectral contamination of the LFP. Note that the slant visible in the phase–frequency relationship is due to frequency-dependent phase shifts induced by the causal low-pass filter.

### LFP tuning curves

A final quantity of interest is the tuning curve, which captures the magnitude of the spike or LFP response as a function of an experimentally manipulated sensory or motor variable. Comparison of tuning curves obtained from the two signals provides additional insight into local cortical processing, the architecture of such neuronal circuits (Liu and Newsome 2006), the spatial extent of the LFP (Katzner et al. 2009), and links of LFP features to synaptic inputs and neuronal outputs (Khawaja et al. 2009). We therefore examined the tuning of a composite signal in which untuned LFPs were combined with the waveforms of well-tuned spikes. For the latter both the waveforms and the tuning curves were taken from a single-unit recording of a neuron found in the middle temporal area (MT) of a macaque monkey. The tuning curve for the spiking responses of this neuron for the motion direction of a visual stimulus is shown in Fig. 6*A*.

Figure 6 shows the LFP tuning curve for the high gamma band, with low-pass filtering of the raw signal and after our spike removal algorithm was applied (Fig. 6*B*). Because the LFP signal is untuned by construction, successful spike removal should yield responses that respond roughly equally to all motion directions, whereas spike contamination will manifest itself as tuning for the downward-rightward motion shown in Fig. 6*A*. The results show that both LFP tuning curves are untuned for visual motion direction, suggesting that the spikes do not introduce spurious tuning into the LFPs, even at higher frequencies (high gamma band).

### Summary of spike contamination in simulated data

Our simulation results reveal that spike contamination can introduce strong artifactual correlations between spikes and LFPs. In further stimulations we found that both spike amplitude and the firing rate of a neuron affect the magnitude of spike contamination in LFPs, with larger firing rates and spike amplitudes leading to greater contamination. Supplemental Figs. S1 and S2 and Table 1 provide summaries of the magnitude of spike contamination for different spike amplitudes and firing rates. As shown in Supplemental Fig. S1, the normalized amplitude of the STA of the LFP rises linearly as spike amplitude increases. Similarly, higher values of SNR give higher artifactual values of the SFC and raise the probability of a neuron being considered phase-locked to specific frequency ranges (Table 1). A weaker effect is observed for high neuronal firing rates (Table 1).

### Other spike removal approaches

As an alternative to simple low-pass filtering, one might attempt to estimate LFPs by removing individual spikes from the wideband signal. This has traditionally been accomplished using one of two methods. The first method involves subtracting each neuron's mean action potential waveform from the raw signal at the corresponding spike times. A second approach involves removing each action potential and interpolating the wideband signal across the resulting gap. Figure 7 compares the performance of these approaches with the Bayesian algorithm, for three of the metrics discussed earlier (LFP STAs, wavelet STAs, and phase-locking histograms). Here the simulated data are the same as those used in the previous figures.

The *top panel* of Fig. 7*A* shows that the mean spike removal method did not completely remove the residual spike correlation from the STA [a difference of 0.48 dB or 16% more residual compared with the Bayesian spike removal (green line)]. By contrast, the interpolation method (Fig. 7*A*, *middle panel*) seems to perform overly aggressive spike removal, in the process introducing an artificial positive deflection of the LFP STA around the moment of the spike (a difference of 1.68 dB or 44% more residual compared with the Bayesian spike removal). For the wavelet STA (Fig. 7*B*), the mean spike removal method (*top panel*) succeeds in attenuating the spurious correlations, as does the interpolation method (*middle panel*). Indeed the latter method slightly outperforms our Bayesian algorithm (*bottom panel*) for this particular measure (SD of the wavelet coefficients for the interpolation method is 0.0013, for the mean spike removal and the Bayesian algorithm 0.0015 and for the low-pass filter 0.0036). However, as shown in Fig. 7*C* (*middle panel*), interpolation can introduce significant spurious correlations into the phase locking histogram. For this measure mean spike removal (Fig. 7*C*, *top panel*) also fails to remove the spurious phase locking correlations, while the Bayesian algorithm shows little residual structure (*bottom panel*). These results are confirmed by the *P* values of Rayleigh statistical tests (0.008 for the mean spike removal, 0.004 for the interpolation method and 0.74 for the Bayesian algorithm, with the hypothesis of phase-locking for a value of *P* < 0.01).

Thus we find that removal of individual spikes is generally preferable to low-pass filtering the wideband signal. Of the methods for removing individual spikes, removal of the mean spike waveform most often leaves residual, artifactual spike–LFP correlations. Interpolation methods perform better by aggressively removing all signals related to each spike, but in the process they can introduce other types of spurious correlations. In most cases our Bayesian algorithm outperformed both previous approaches. Furthermore, as we show in the discussion and the appendix, the Bayesian algorithm generalizes the advantages of these previous approaches, while casting the problem in a rigorous mathematical framework within which key assumptions can be made explicit.

### Real recordings

##### SINGLE NEURONS.

The previous sections showed that low-pass filtering cannot completely remove spike waveforms from composite signals constructed to have no correlation between spikes and LFPs. The resulting spectral contamination of the LFP can introduce significant bias into such measures as the STA, the wavelet STA, the SFC, and phase-locking histograms, especially for higher-frequency bands. By comparison, we found that the tuning of the LFP to stimulus features as measured by tuning curves is not significantly affected by spikes, even when spike responses are strongly tuned. We introduced a new spike removal algorithm that performs better than low-pass filtering.

To examine the consequences of these findings for signals obtained in vivo, we recorded single-unit activity and LFPs from two alert macaque monkeys. A total of 26 isolated neurons in areas V6 (19) and MT (7) neurons were analyzed. Both areas are specialized for the processing of visual motion (Fattori et al. 2009; Maunsell and Newsome 1987). In each case, we determined the selectivity of the single-unit activity for visual motion by drifting a sinusoidal grating or a random-dot stimulus in 12 different directions spaced evenly around the circle. Cells that failed to respond to the stimulus or that lacked direction selectivity (defined in methods) were excluded from the analysis.

As shown for a typical cell in Figs. 8 and 9, the trends in the real data qualitatively match those in our simulations. In several cases, including the one illustrated in Fig. 8, the LFP was correlated with the timing of spikes and this was apparent in both STA analysis and phase synchronization. These correlations were not removed by the spike removal algorithm. Thus our results are consistent with previous findings (Gregoriou et al. 2009; Whittingstall and Logothetis 2009) that indicate a relationship between spike timing and LFP phase. However, both the range of frequencies over which this relationship is found and the magnitude of the correlations are affected by the method of preprocessing the raw signals. For example, inspection of the spike-triggered coherence before spike removal (Fig. 9*A*, blue line) would lead one to conclude that high frequencies of the LFP are robustly correlated with spikes, whereas application of the spike removal algorithm shows that in fact only low frequencies of the LFP are modulated with spikes (Fig. 9*A*, red line).

Figure 10*B* shows LFP tuning curves before and after removal of the spikes of a robustly tuned neuron (Fig. 10*A*). As expected from our simulations, LFP tuning curves were essentially unaffected by spike removal. The next section quantifies the extent of this contamination for our population of cells.

##### Population results.

We assessed the level of spectral contamination for different coherence metrics for all our neuronal recordings and these results are shown in Fig. 11. Figure 11*A* shows that the normalized amplitude of the LFP STA (SNR) for various frequency bands, measured around the time of spikes (10 ms before to 10 ms after peak depolarization), drops significantly after the removal of spikes [average of 5 dB (68%) drop for high gamma and 1 dB (20%) decrease for low gamma]. The wideband signal, which contains both high and low frequencies, is similarly affected [average of 7 dB (80%) decrease]. In contrast, the STA at low frequencies is only weakly affected by the spike removal algorithm. A similar pattern was observed for the SFC, as shown in Fig. 11*B*. There is on average a 66% decrease of SFC values for the high gamma and a 30% decrease for the low gamma frequencies.

The effect of spike removal on phase-locking histograms follows a similar pattern across frequencies, as shown in Table 2. The phase-locking metric is dramatically affected by spike remnants, with the proportion of cells classified as phase-locked to high gamma frequencies dropping by half after spike removal.

In contrast to these results, Fig. 11*C* shows the correlation between the LFP and spike tuning curves, both before and after spike removal. For both high and low frequencies, the difference in the correlation coefficient of the two tuning curves is very small (average reduction of 0.06 and 0.03, respectively), indicating that spike contamination has little effect on LFP tuning curves. This conclusion is further supported by the fact that the points are clustered along the unity line, despite wide variation in the extent to which the tuning of the two signals is correlated. These results are consistent with our simulations (Supplemental Figs. S1 and S2).

### Other approaches of spike removal on real data

For comparison with previous approaches, we also estimated the effects of mean spike removal and interpolation methods on our real data and the results are shown for an example recording in Fig. 12 (same example as that in Figs. 8 and 9). As in the simulated data, mean spike removal leads to larger LFP STAs and more structure in the phase-locking histogram (Fig. 12, *A* and *C*). In this specific example, linear interpolation performs very similarly to the Bayesian spike removal algorithm. However, across the population, both mean spike removal and interpolation yield larger estimates of spike–LFP relationships than our Bayesian algorithm (summary in Table 2). Although ground-truth is not available in the real data, our simulations from the preceding sections suggest that the greater synchronization found with other algorithms is likely to be artifactual.

## DISCUSSION

In recent years there has been an increase in the sophistication of analytical methods for interpreting neuronal activity. Although these advances are important for understanding brain circuitry, they are also easily influenced by particular assumptions about the way in which the raw signals should be processed. In this study, we have explored one issue related to such signal processing by quantifying the contamination of the LFP from action potential waveforms in both simulated and real data.

Our results show that estimates of the relationship between LFPs and spikes can be considerably affected by spike contamination. This results simply from the spectral characteristics of spike trains and individual spike waveforms: both have power at the frequencies used to study LFPs (Fig. 2*B*). Consequently, spikes and LFPs cannot be isolated from a raw voltage signal with simple frequency-domain filtering. To estimate the extent of this problem we performed simulations in which spike waveforms were added to phase-randomized LFP signals, such that the spectral profile of the LFP was maintained but there was no actual LFP–spike relationship. Nevertheless, when LFPs were estimated by low-pass filtering, our analysis showed artifactual spike–LFP STAs (Fig. 4) and spurious phase-locking of spikes to LFPs (Fig. 5). This effect was particularly strong at high frequencies. In contrast, tuning curves in response to visual stimuli were not strongly affected by spike contamination (Fig. 6), even at high frequencies (≤140 Hz), unless the spike-to-LFP power ratio was unnaturally high (>5 dB). Similar results were obtained using real recordings with direction-tuned neurons from the primate visual system (Figs. 8⇑⇑–11). To address these issues we developed a novel Bayesian spike removal algorithm and demonstrated that it is capable of removing artifactual spike–LFP correlations.

### Comparison with previous studies

Although the existence and potential importance of spectral contamination of LFPs by spikes has been recognized for some time (Pesaran et al. 2002), the issue has been addressed only rarely in recent studies. One common approach in previous work is to subtract the mean spike waveform from the raw signal at each spike time (Pesaran et al. 2002; Zanos et al. 2006).

The generative model underlying our spike removal algorithm was constructed to generalize this spike subtraction method. In fact, we show in the appendix that our spike removal method reduces to the spike subtraction method when spike waveforms do not overlap, no assumptions are made on the smoothness of the LFP, and an offset term is not included in the model. The full spike removal algorithm thus works similarly to the spike subtraction method, but corrects for spike train autocorrelations and reduces low-frequency components in spike waveforms attributable to LFPs. Each of these adjustments corrects for artifacts introduced by spike subtraction, thus enabling more effective spike removal, as we show in detail in Figs. 7 and 12 and the appendix.

Another proposed method of spike removal is to replace time samples around each spike with samples interpolated from the LFP before and after the spike (Galindo-Leon and Liu 2010; Jacobs et al. 2007; Okun et al. 2010). This approach, although again being superior to low-pass filtering, introduces the difficulty of selecting the method of interpolation, including the function used to interpolate and the length of its support on either side of spikes, without a firm theoretical foundation to justify these choices. Linear interpolation, which is commonly used on wideband signals for this purpose, not only does not completely remove spurious correlations but also introduces new ones, as shown in detail in Fig. 7. The difficulty of selecting a proper method could be alleviated by casting the interpolation as a Bayesian inference problem similar to the one considered here. Within this framework, the functional form of the interpolation would be adjusted by the statistics of the local field potential. A more fundamental issue is that when recording in areas with high transient firing rates, interpolation must sometimes be done over intervals several times longer than a single spike waveform. This is likely to introduce additional artifacts into the estimated LFP.

A recent method (David et al. 2010) uses a more sophisticated linear filtering approach to remove all spike–LFP correlations from the raw signal. In the appendix, we show that our proposed spike removal algorithm reduces to this method in the limit of very long spikes. When very long spikes are assumed, however, the method essentially removes all linear relationships between spikes and LFPs, including legitimate correlations, and thus cannot provide an estimate of the degree to which spikes and LFPs are actually correlated. In contrast, our method focuses on removing the purely artifactual parts of the LFP signal that are introduced by spike waveforms of short duration (≤3 ms), while preserving legitimate correlations that occur at longer timescales. Although the method of David et al. (2010) can be modified to accommodate the assumption of finite spikes, we show in the appendix that this leads to residual artifacts around the time of spikes due to the way that the method handles edges.

### Practical considerations for effective spike removal

Our spike removal algorithm assumes that the wideband signal derived from extracellular recordings is the sum of a low-frequency LFP signal, stereotyped spikes of finite length occurring at known times, noise, and an offset (Fig. 1*B*). For the spike removal process to be most effective, it is important that careful preprocessing be done so that the signal closely matches the assumptions of the generative model.

The form of the prior assumed on the LFP signal should be adjusted so that it closely matches the empirical autocorrelation properties of the observed signal; a practical procedure to choose the form of the prior is detailed in the appendix. Line noise should be removed before spike removal because this effect is not dealt with by the generative model (see methods). Furthermore, analog filters in the recording hardware can cause frequency-dependent phase shifts in the wideband signal (Nelson et al. 2008), which may artificially lengthen the duration of the observed spike waveforms. This issue can be avoided by compensating for phase distortions by digital filtering (Nelson et al. 2008).

Our spike removal algorithm will remove artifactual linear correlations between spikes and LFPs regardless of whether the spikes in the wideband signal are completely stereotyped. In the latter case, however, residual high-frequency artifacts of random phase will be visible around the time of spikes. These artifacts may be picked up by second-order methods such as the wavelet STA and Volterra methods (Zanos et al. 2006, 2008). It is therefore crucial to perform careful spike sorting before removing spikes. It is also important that spike timings be accurate because spike timings that are inaccurate by as little as one or two samples will cause observed waveforms to shift appreciably (Sahani 1999). Spike timings that are obtained by thresholding a high-passed signal can be inaccurate; more accurate timings can be obtained by aligning spikes to their peak depolarization (Quiroga et al. 2004; Sahani 1999). Additionally, if spike waveforms shift over time, for example because of electrode drift, the algorithm should be applied on short segments of the data to obtain the best results; we show in the appendix how to perform spike removal in chunks.

Many of these prescriptions are equally applicable to other methods of spike removal and, of course, they do not guarantee complete removal of spike artifacts from the wideband signal. Nonetheless, we found empirically that the spike removal algorithm diminished artifactual correlations to a great extent. When an interesting finding relies crucially on the quality of spike removal, it may be worthwhile to back up claims with ground-truth simulations, similar to those of Fig. 1*A*.

### Hardware considerations

Several manufacturers offer products designed to record local field potentials (often designated *LFP boards*). These boards typically low-pass filter the wideband signal from a headstage, with a cutoff usually ranging from 150 to 300 Hz, and sample the resulting signal at a typical rate of 1 kHz. Our results indicate that the signals derived from these boards are at all frequencies appropriate for analyses that are insensitive to the relative timing of spikes and LFPs, for example to derive tuning curves (Figs. 6 and 10). However, our results also imply that analyses that are sensitive to temporal relationships between spikes and LFP, including the STA, spike-triggered coherence, phase-locking histograms, and predictions of spikes from LFPs, can be corrupted by the presence of spike remnants in LFPs derived from these boards. This contamination is especially strong at gamma frequencies >60 Hz.

It may be possible to work around this contamination, either by limiting analyses to frequencies <60 Hz, removing all linear correlation between spikes and LFPs using the methods of David et al. (2010), or by recording spikes and LFPs on different electrodes. However, a safer approach is to acquire the complete wideband signal during electrophysiological recording and to perform careful postprocessing using methods similar to those outlined here. For data that have already been recorded with LFP boards, it may be possible to recover the spike-free LFP by extending our algorithm to include a model of the processing performed by the boards, but this has not yet been attempted.

### Conclusions

Numerous studies examining the relationship between spikes and LFPs have reported convincing and significant results. Important examples include studies using spikes and LFPs recorded from different electrodes (Gregoriou et al. 2009; Whittingstall and Logothetis 2009), spikes and cortical EEGs (Whittingstall and Logothetis 2009), or spikes, LFPs, and BOLD signals (Goense and Logothetis 2008) and results like those shown in Figs. 8 and 9 are consistent with some of these studies. However, our results indicate that commonly used methods of preprocessing may lead to overestimation of the association between spikes and LFPs. This is particularly the case when low-pass filtering is used as the only preprocessing step of the LFP, as has often been done in past work (e.g., Chalk et al. 2010; Litaudon et al. 2008; Rasch et al. 2008). The existence and magnitude of spurious correlations in such cases will depend on the nature of the data and the types of analyses applied. Particularly when higher frequencies of the LFP are analyzed with sophisticated nonlinear methods, the effect of even modest contaminations on the results can be significantly magnified (Zanos et al. 2006). On the other hand, successfully removing the spike waveforms from the raw recordings can reveal significant spike–LFP relationships that might otherwise have been masked (as in Fig. 9*B*).

We have introduced a spike removal algorithm that is designed to remove specific spike waveform(s) from the raw voltage signal, without affecting real spike–LFP correlations that are of great interest to the study of neural circuitry. In simulations, where ground-truth is available, our algorithm removed virtually all spike artifacts and the results of all subsequent analyses on the “spike-free” LFP did not show significant correlations between spike timing and any frequency band of the LFP. Similarly, when applied to real data the algorithm often eliminated artifactual spike–LFP relationships that would have emerged from LFPs obtained by low-pass filtering the raw voltage signal (Fig. 9). However, in other cases (Figs. 8 and 9) strong spike–LFP correlations remained after application of the spike removal algorithm and such synchronization will be an important topic for further study. Indeed, with accurate spike removal one can in principle study spike–LFP relationships at higher LFP frequencies (>150 Hz), despite their overlap with the spectral content of the spike waveforms.

## GRANTS

This work was supported by a Canadian Institutes of Health Research Grant MOP-79352 to C. C. Pack and a Montreal Neurological Institute Centre for Commercialization and Research fellowship to T. P. Zanos.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## ACKNOWLEDGMENTS

We thank Dr. Stavros Zanos for indicating the issue of spike contamination in previous local field potential studies and for providing analysis ideas, R. Cipriani for assistance in data collection, J. Coursol and C. Hunt for assistance with animal care and training, and N. Jabakhanji for help with computer programming.

## Footnotes

↵1 The online version of this article contains supplemental data.

- Copyright © 2011 The American Physiological Society