Abstract
Extracellular physiological recordings are typically separated into two frequency bands: local field potentials (LFPs) (a circuit property) and spiking multiunit activity (MUA). Recently, there has been increased interest in LFPs because of their correlation with functional magnetic resonance imaging blood oxygenation leveldependent measurements and the possibility of studying local processing and neuronal synchrony. To further understand the biophysical origin of LFPs, we asked whether it is possible to estimate their time course based on the spiking activity from the same electrode or nearby electrodes. We used “signal estimation theory” to show that a linear filter operation on the activity of one or a few neurons can explain a significant fraction of the LFP time course in the macaque monkey primary visual cortex. The linear filter used to estimate the LFPs had a stereotypical shape characterized by a sharp downstroke at negative time lags and a slower positive upstroke for positive time lags. The filter was similar across different neocortical regions and behavioral conditions, including spontaneous activity and visual stimulation. The estimations had a spatial resolution of ∼1 mm and a temporal resolution of ∼200 ms. By considering a causal filter, we observed a temporal asymmetry such that the positive time lags in the filter contributed more to the LFP estimation than the negative time lags. Additionally, we showed that spikes occurring within ∼10 ms of spikes from nearby neurons yielded better estimation accuracies than nonsynchronous spikes. In summary, our results suggest that at least some circuitlevel local properties of the field potentials can be predicted from the activity of one or a few neurons.
Introduction
The brain is usually studied at multiple scales from molecules to systems. How signals at one scale relate to those at other scales is poorly understood. Of particular importance toward understanding how perception and action are orchestrated by neuronal signals is to characterize how the activity of neuronal circuits consisting of tens to thousands of neurons arises from their component units and their interactions. A major challenge toward achieving this goal is the experimental difficulty of recording simultaneously from multiple nearby neurons and the theoretical efforts required to characterize biophysically realistic neural networks.
The extracellular voltage recorded extracellularly through microwires is typically separated into two frequency bands (Logothetis, 2002): the multiunit spiking activity (MUA) and the local field potentials (LFPs) (see Fig. 1A). MUA represents a weighted sum of the action potentials of neurons within a radius of ∼200 μm around the electrode tip (Holt and Koch, 1999; Gold et al., 2006). Lowpass filtering the extracellular signal (with a corner frequency of ∼100–300 Hz) yields the LFP. The biophysical nature of the MUA signal is much better understood than the origins of the LFP. The LFPs are thought to reflect the activity of large numbers of neurons in a sphere of one to several millimeters around the recording electrode (Mitzdorf, 1985; Juergens et al., 1999) (but see Katzner et al., 2009). Current source density analyses and simultaneous recordings of spikes and LFPs have suggested that LFPs are more strongly correlated with EPSPs, afterpotentials, and dendritic spikes than with the output action potentials of the surrounding neurons (Haberly and Shepherd, 1973; Mitzdorf, 1985; Logothetis, 2002). This has led to the notion that LFPs represent the input to and local processing within a given brain area.
Many approaches have been followed to compare spikes and LFPs (O'Keefe and Recce, 1993; Fries et al., 2001; Laurent, 2002; Pesaran et al., 2002; Bédard et al., 2004; Kreiman et al., 2006; Belitski et al., 2008; Montemurro et al., 2008; Rasch et al., 2008; Nauhaus et al., 2009, among others). Here we asked whether it is possible to estimate the time course of LFPs from spiking activity. At first glance, the large differences in terms of the spatial scales and biophysical origin between the two signals might suggest that this would be a challenging task. We use methods from “signal estimation theory” (Poor, 1994) to show that a linear filter operation on the activity of one or a few neurons can explain a significant fraction (but not all) of the LFP time course. The estimations have a spatial scale of ∼1 mm and a temporal scale of ∼200 ms. The linear filter estimation extrapolates across cortical regions and behavioral conditions. We also show that spikes that temporally coincide within ∼10 ms with spikes from nearby neurons lead to better estimations. Together, these results suggest that the activity of individual neurons can be related to at least some of the circuitlevel properties of local field potentials.
Materials and Methods
Electrophysiological recordings and data processing.
Electrophysiological data recorded from seven anesthetized monkeys (Macaca mulatta) are included in the present study. The surgical methods and recording setup have been described previously (Rasch et al., 2008). Briefly, simultaneous recordings of neuronal activity were made from primary visual cortex (V1) using 8–16 electrodes configured in 4 × 4 or 8 × 2 matrices in a grid of 1–2 mm. We study data recorded during spontaneous activity (labeled “spont” throughout the text) and data recorded while the monkey was shown a commercial movie (labeled “stim” throughout the text). In the spontaneous condition, the V1 data studied here (7 monkeys, 109 electrodes) were recorded during periods when the input screen was blank for ∼4 min. Multiple segments of spontaneous activity were recorded within each session. From each session, we considered five segments of spontaneous activity. Unless described otherwise in text, all electrodes recorded in each session were included in the analyses. This resulted in 545 (5 × 109) recorded time series of 4 min length, which we refer to as “trials” throughout the manuscript. In three sessions (two anesthetized monkeys), four electrodes were simultaneously placed in the lateral geniculate nucleus (LGN). In the visual stimulation condition (6 monkeys, 84 electrodes, 420 “trials”), a 4–6 min movie segment was shown to the monkey. The movie frames, synchronized to the monitor refresh rate (60 Hz), encompassed 7–12° of the visual field (Rasch et al., 2008). The analyses for the spont and stim conditions are based on 4 min segments in which we discarded the initial and final 30 s to avoid potential nonstationarities. Examples of the power spectral density for the LFP recordings and the spike trains are shown in supplemental Figure S1 (available at www.jneurosci.org as supplemental material). The LFP power spectral densities show a power law decay similar to those reported in other studies (Henrie and Shapley, 2005; Nauhaus et al., 2009, among many others), and the log–log plots can be fitted with a line of slope close to 2 (Milstein et al., 2009) (supplemental Fig. S1, available at www.jneurosci.org as supplemental material). For comparison purposes, in Results (see Spontaneous activity versus visual stimulation) and Figure 4A, we show results obtained based on recordings from awake monkey inferior temporal (ITC) cortex during passive viewing conditions from the study reported by Kreiman et al. (2006).
The data preprocessing steps have been described in detail previously (Rasch et al., 2008). Briefly, electrode signals were decimated to 7 kHz. The 7 kHz signal was lowpass filtered with a cutoff frequency of 220 Hz and resampled at 500 Hz to obtain the LFPs. Spike times were detected by applying a threshold to the highpassfiltered 7 kHz signal (fourthorder Butterworth, cutoff frequency of 500 Hz). Except when noted otherwise, the detection threshold was applied at 5 SDs of the noise component of the MUA signal. In supplemental Figure S9 (available at www.jneurosci.org as supplemental material), we examined the dependence of the estimation accuracies on the spike detection threshold. In supplemental Figure S9 (available at www.jneurosci.org as supplemental material), we also performed spike sorting using the algorithm described by Quiroga et al. (2004) to compare the performance of MUA against singleunit activity (SUA).
Linear estimation of timevarying signals from spike trains.
We asked whether we could infer the time course of the local field potentials from the spike trains recorded from one cell or a small number of cells near the LFP electrode. One possible approach to this problem is to apply methods of signal estimation theory (Poor, 1994). Such methods have been used in neurobiology by several groups (Bialek et al., 1991; Gabbiani and Koch, 1996; Rieke et al., 1997; Kreiman et al., 2000) to study the transmission of information by peripheral sensory neurons. Here we analyzed the data using the following algorithm (Bialek et al., 1991; Gabbiani and Koch, 1996) (see the scheme in Fig. 1A). Let be the spike train obtained after subtracting the mean firing rate x_{0} (where t_{i} are the spike occurrence times). Except in supplemental Figure S9B (available at www.jneurosci.org as supplemental material), x(t) contains spikes from multiunit activity. A linear estimate, L_{est}(t), of the local field potential, L(t), given the spike train was obtained by convolving x(t) with a filter, h(t), where T is the duration of the LFP recording, and τ is the integration variable spanning the convolution period. The Wiener–Kolmogorov (W–K) filter, h(t), is chosen in such a way as to minimize the mean squared error, ε^{2}, between the LFP and its estimate: where the integration is over the LFP recording duration T. An explicit formula for this filter is given by Poor (1994): where f_{c} is the cutoff frequency of the LFP, P_{Lx} is the Fourier transform of the crosscorrelation between the LFP and the spike train, and P_{xx} is the Fourier transform of the spike train autocorrelation.
We considered three types of W–K filters: (1) “trialspecific” filters, (2) “electrodespecific” filters, and (3) “monkeyspecific” filters (Fig. 1). For the trialspecific filter, each trial was divided into two segments of equal length. The first segment was used to estimate the W–K filter according to Equation 4. The second segment was used to estimate the LFP according to Equation 2. The performance of the filter in this test is labeled “estimation accuracy” throughout the text (see below for the computation of the estimation accuracy). We compared this estimation accuracy against the “reconstruction accuracy” obtained by using the same segment to compute the W–K filter and to reconstruct the LFP time course. The reconstruction accuracies are reported in the figures as black squares and likely contain some amount of overfitting. All the conclusions in this study are based on the estimation accuracies in which the data used to compute the filter are separate from the data used to estimate the LFP. For the electrodespecific filter, we used half of the trials recorded from each electrode (oddnumbered trials) to compute the W–K filter and used the remaining half of the trials (evennumbered trials) to estimate the LFPs (Fig. 1B). Unless stated otherwise, throughout the text, we report the estimation accuracies obtained by the electrodespecific filters. As described in Results, we observed that the W–K filters were very similar across electrodes and recording conditions. Therefore, we also evaluated the performance of a monkeyspecific filter in which the data from half of the electrodes were used to compute the W–K filter, and the data from the remaining half of the electrodes were used to estimate the LFP (Fig. 1C).
The W–K filter in Equation 4 will not be causal in general [i.e., h(t) ≠ 0 for t > 0]. Imposing the causality constraint on the filter h requires solving the causal Wiener–Hopf equation (Poor, 1994). Because h has a finite support in the time domain, causality can also be implemented by introducing a delay in the filter (Bialek et al., 1991). In the present study, we used this second method to impose causality (see Results, Causal filters) (see Fig. 8).
Data analysis.
All the data were analyzed using MATLAB (MathWorks). The spike occurrence times were resampled at 500 Hz together with the LFPs. Estimates of the LFP and spike train power spectra were obtained with a fast Fourier transform algorithm and Bartlett windowing using nfft = 2048. As expected, the estimation accuracies, and in particular the reconstruction accuracies, depend on nfft (supplemental Fig. S2A, available at www.jneurosci.org as supplemental material). Small values of nfft do not capture enough information of the LFP time course, and very large values lead to overfitting. Throughout the text, we use nfft = 2048.
The same analysis was performed to obtain estimates of the crosscorrelation between spike trains and LFPs. The cutoff frequency of the LFP was estimated by fitting the squared gain of the lowpass (Butterworth) filter transfer function to the power spectrum of the LFP, L(t). The optimal W–K filter, h(t), was obtained by deconvolving the crosscorrelation of the spike train and the LFP with the power spectrum of the spike train according to Equation 4. The LFP estimations [L_{est}(t)] were obtained by computing the convolution of the W–K filter and the spike train, x(t) (see Eq. 2), in the frequency domain with the use of a fast Fourier transform.
Performance measures.
Because the filter h is derived by minimizing the mean squared error between the LFP and the estimated LFP (see Eq. 3), a natural measure for the quality of the estimations is the mean squared error, ε^{2}. Another commonly used performance measure is the Pearson's correlation coefficient (denoted by r throughout the manuscript). Given that, in this case, both L(t) and L_{est}(t) are normalized so that they have 0 mean and SD 1, a linear relationship exists between the mean squared error and the correlation coefficient: ε^{2} = 2(1 − r) (for the derivation of this expression, see supplemental data and Fig. S1F, available at www.jneurosci.org as supplemental material). Throughout the text, we report the r values and refer to them as the estimation accuracy.
To assess the statistical significance of the estimations, we compared the results against a null hypothesis in which no correlation existed between the temporal structure of the LFPs and the spike trains. Under the null hypothesis, we generated random spike trains with the same mean firing rate as the experimental spike trains but with spike times governed by a Poisson process [x_{rand}(t)]. We then followed the same procedure in Equation 4 to compute the optimum W–K filter to estimate the LFP time course from the random spike train, h_{rand}: with x_{0}^{rand} = x_{0}. If the temporal structure of the spike trains conveys no information about the LFP time course, we would expect the correlation between the estimated LFP and the actual LFP obtained from the experimental spike trains to be close to the correlations obtained from the random spike trains. We repeated this procedure 50 times; the average estimation accuracies obtained under the null hypothesis are shown as black triangles throughout the figures. To assess the statistical significance of the LFP estimations, we performed twotailed t tests, comparing the estimation accuracies obtained under the null hypothesis against those obtained with the actual spike trains.
Results
Linear LFP estimations based on single trials
We considered 545 spontaneous activity “trials” recorded from V1 (109 electrodes, 7 monkeys), of ∼4 min length (we discarded the first and last 30 s to avoid any potential border effects). To estimate the LFP from the spike train, we computed a linear filter that minimizes the difference between the LFP and its estimate (Eqs. 2–4) (Fig. 1A). To avoid overfitting the data, the estimation filter h was constructed (as described in Materials and Methods, Linear estimation of timevarying signals from spike trains and Data analysis) with the first half of the LFP and spike train of each trial. The estimation, L_{est}(t), was computed by convolving h with the second half of each spike train. Therefore, there was no overlap between the data used to compute the filter and the data used to test the filter by estimating the LFP time course and measuring the correlation between L(t) and L_{est}(t) (see Materials and Methods, Performance measures).
Three 1 s segments of typical LFP traces and their estimations are shown in Figure 2. These three examples illustrate the range of estimation accuracies from good (Fig. 2A, where the estimation captures, to a large degree, the overall structure of the LFP time course) to poor (Fig. 2C, where the estimate is only weakly correlated with the actual LFP). We show another example electrode in Figure S3A (available at www.jneurosci.org as supplemental material), and we show more data segments for the electrode illustrated in Figure 2A in Figure S3B (available at www.jneurosci.org as supplemental material). We quantified the accuracy of the LFP estimation by computing the Pearson's correlation coefficient between the estimation and the actual LFP [referred to as the estimation accuracy (r) throughout the manuscript; see above, Performance measures]. In the examples shown in Figure 2, r has a value of 0.61 (A), 0.29 (B), and 0.01 (C). On the right in Figure 2, we show the corresponding W–K filters for time lags between −800 and +800 ms [the actual number of points in the W–K filter was nfft + 1 (2049), but there were only small deviations from 0 outside the ±800 ms range]. As illustrated by the examples in Figure 2, wider W–K filters tended to yield higher estimation accuracies. Typical features of the W–K filter include a sharp downstroke for tens of milliseconds for negative time lags followed by a slower upstroke lasting >100 ms for positive time lags (see also the example W–K filters in Figs. 3C, 4B, 8 and supplemental Figs. S3, S7, S9, available at www.jneurosci.org as supplemental material).
Figure 3A shows the estimation accuracy for each one of the 109 electrodes recorded from V1 during spontaneous activity after averaging across all trials (five trials). We found that r varies significantly across trials and electrodes. Whereas the estimation accuracy for some trials is higher than 0.6, it is almost 0 for others (see the overall distribution in supplemental Fig. S2B, available at www.jneurosci.org as supplemental material). Overall, most of the electrodes (95 of the 109 electrodes) yielded a statistically significant estimation accuracy compared with the null hypothesis obtained by creating a Poisson spike train with the same mean firing rates (p < 0.01, twotailed t test) (see above, Performance measures). The average estimation accuracy across V1 was 0.36 ± 0.15 (mean ± SD, n = 109 electrodes) (Fig. 3A). On average, the estimations obtained under the null hypothesis led to a correlation of r̄_{rand} = 0.001 ± 0.035 (mean ± SD, n = 109 electrodes) (Fig. 3A, triangles and bottom dotted line). Despite the variability across trials and electrodes, we found that the estimations were highly significant [twosample Kolmogorov–Smirnov test comparing the distribution in supplemental Figure S2B (available at www.jneurosci.org as supplemental material) against the corresponding distribution obtained under the null hypothesis, p < 10^{−10}].
To gain additional insight into what contributes to the variability in the estimation accuracies across trials and across electrodes, we plotted the estimation accuracy as a function of basic properties of the spike trains and LFPs (supplemental Fig. S4, available at www.jneurosci.org as supplemental material). The estimation accuracy increased with higher firing rates (correlation coefficient of 0.49) (Fig. S4A,D, available at www.jneurosci.org as supplemental material), with the total power of the LFP signal (correlation coefficient of 0.44) (supplemental Fig. S4B,E, available at www.jneurosci.org as supplemental material) and also with the coefficient of variation (CV) of the interspike interval distribution (correlation coefficient of 0.59) (supplemental Fig. S4C,F, available at www.jneurosci.org as supplemental material). It should be noted that these relatively large CV values are based on MUA and not SUA. Based on supplemental Figure S4 (available at www.jneurosci.org as supplemental material), subsequent analyses in the manuscript focused on those electrodes that had a firing rate of at least 5 spikes/s and a CV of at least 1 (88 of the 109 electrodes).
Given the intrinsic noise in the spike trains and LFPs, we were interested in estimating an upper bound for the estimation accuracy r. For this purpose, we asked how well this linear algorithm performs when estimating the LFP time course by estimating the LFP measured in one electrode, LFP_{1}, from the LFP recorded from a nearby electrode, LFP_{2} (supplemental Fig. S6A, available at www.jneurosci.org as supplemental material). Because the recordings were performed simultaneously, two electrodes that are nearby are likely to be measuring a similar signal from slightly different places. As expected, the LFP–LFP estimation decreased with the distance between the two electrodes (supplemental Fig. S6C, available at www.jneurosci.org as supplemental material). For a distance of 1 mm (the smallest distance available in our dataset), the mean LFP–LFP estimation was 0.55 ± 0.05 (mean ± SD, n = 109 electrodes; compare with 0.36 for the spike–LFP estimation accuracy with a distance of 0 mm).
Toward a general function to map spikes into LFPs
Is the map between spikes and LFPs specific to each trial/electrode or is there a general filter that can extrapolate across different trials or even different recording electrodes? We asked whether we could find a general filter that, when applied to a spike train in V1 recorded from a given electrode, would give us an estimate of the corresponding LFP recorded from the same electrode. We addressed this question in three steps, with increasing degrees of extrapolation. First, we constructed an electrodespecific filter to predict the LFPs of all the trials recorded from a given electrode after computing the W–K filter using different trials recorded from the same electrode (Fig. 1B, one filter per electrode as opposed to one filter per trial as used in the previous section). Second we constructed a monkeyspecific filter to predict the LFPs of all the trials and electrodes of each monkey after computing the W–K filter using different electrodes recorded from the same monkey and same brain region (Fig. 1C, one W–K filter per monkey). Finally, we assessed whether a W–K filter computed from all the electrodes in a given monkey could be used to estimate the LFPs recorded in a different monkey. Note that, in all cases, the LFP estimate for a given electrode was computed from the spike trains recorded from the same electrode (using Eq. 2). The extrapolation refers to the way in which the W–K filter is computed (using Eq. 4 and data from different trials, different electrodes, or different monkeys).
We constructed the electrodespecific filter for a given electrode by minimizing the sum of the mean squared errors in the estimations of individual trials recorded from that electrode (Fig. 1B,C). We find the filter (h_{mean}) that minimizes where N is half of the total number of trials for the electrode, L_{j} is the LFP corresponding to trial j, and x_{j} is the spike train of trial j with the mean subtracted. The minimization of the error leads to an expression for the Fourier transform of the filter (for a derivation of this expression, see supplemental data, available at www.jneurosci.org as supplemental material), where _{j}P_{Lx} and _{j}P_{xx} are the cross power spectrum and power spectrum, respectively, of the LFP and spike train in trial j. In a similar manner, we constructed a single monkeyspecific filter for each monkey by using Equations 6 and 7 but summing over all the trials and electrodes used to construct the filter (half of the electrodes for each monkey; see below).
In both cases (electrodespecific W–K filter and monkeyspecific W–K filter), we estimated the LFP in each trial by convolving the spike train in the same trial with h_{mean} (instead of using the filter h computed for each trial as in the previous section). When building an electrodespecific W–K filter, we used different trials to compute the filter and to compute the estimation: we picked half of the trials for a given electrode to build the filter according to Equation 7 and used the remaining trials from the same electrode to predict the LFP time course (Fig. 1B). When building a monkeyspecific W–K filter, we used different electrodes to compute the filter and to compute the estimation (Fig. 1C): we randomly picked half of the electrodes to build the filter according to Equation 7 and used the remaining electrodes to predict the LFP time course.
Using electrodespecific filters, the mean estimation accuracy was r̄ = 0.31 ± 0.14 (Fig. 3B, dashed line). Using monkeyspecific filters, the mean estimation accuracy was r = 0.30 ± 0.16 (Fig. 3C, dashed line). As expected, the mean estimation accuracy in this case is lower than the values reported in Figure 3A, yet it is remarkable that a general linear filter can achieve a rather accurate estimation of the LFP time course. The distribution of LFP estimations were highly significant compared with those obtained using Poisson spike trains with the same mean rate [p < 10^{−6} (electrodespecific filters) and p < 10^{−4} (monkeyspecific filters), twosample Kolmogorov–Smirnov test]. The shape of each of the monkeyspecific W–K filters is shown in Figure 3C.
The properties of the Wiener–Kolmogorov filter in the context of encoding timevarying signals by spike trains have been studied using idealized neuron models (Gabbiani, 1996; Gabbiani and Koch, 1996). Using an integrateandfire model, assuming a linear filter, neglecting the refractory period and assuming an exponentially distributed spike threshold, it has been shown that the spiketriggered average (STA) of the signal yields an estimation filter that converges onto the W–K optimal filter in the limit of low firing rates (Gabbiani and Koch, 1996). This result also holds for experimental recordings (Wessel et al., 1996). We therefore quantified how well we could estimate the LFP time course using Equation 2 from a filter built by computing the spiketriggered average of the LFP instead of the optimal Wiener–Kolmogorov filter. Supplemental Figure S7, A and B (available at www.jneurosci.org as supplemental material), compares the shape of the STA with the shape of the average electrodespecific optimal W–K filter. The STA shows a sharp negativity for negative time lags and an upswing for positive time lags. This stereotypical structure of the STA is common across V1 electrodes and resembles the shape of the optimal filter h_{mean}. The STA yielded a good estimate of the LFP time course (supplemental Fig. S7C, available at www.jneurosci.org as supplemental material). As expected, this estimate was slightly worse than the one obtained using the optimal W–K filter. Thus, to a coarse approximation and in the limit of lowfiring rates, we can think of the W–K filter as the average LFP surrounding a spike.
Spontaneous activity versus visual stimulation
The results presented thus far correspond to recordings during spontaneous activity. We asked whether the linear estimation of LFPs from spike trains would extend to conditions in which V1 neurons are activated by visual input (as opposed to spontaneous firing). To examine the influence of visual stimulation on the spike–LFP relationship, we repeated the analyses in recordings from 84 electrodes in V1 while the anesthetized monkeys were shown commercial movies (see Materials and Methods). Consistent with the previous findings, we could also linearly estimate the LFP time course during visual stimulation (Fig. 4A). The shape of the W–K filter for the visual stimulation condition was similar to the corresponding shape during spontaneous activity. The mean estimation accuracy during the stimulation condition was 0.32 ± 0.13 using electrodespecific filters. Furthermore, there was a strong correlation between the estimation accuracies during spontaneous activity and those during stimulusdriven activity (correlation coefficient of 0.64) (supplemental Fig. S5, available at www.jneurosci.org as supplemental material).
Given the similarity across conditions (spontaneous activity vs visual stimulation) for a given monkey, we asked whether the LFP estimations could extrapolate across monkeys. Using the same approach described by Equations 6 and 7, we computed the W–K filter using all the recordings from a given monkey (Eqs. 6–7) and used this filter to estimate the LFPs recorded in a different monkey (Eq. 2). Figure 4B shows that the W–K filter from one monkey can be used to estimate the LFPs recorded from the same area (V1) in a different monkey (in all cases, the LFP estimate for a given electrode is based on the spike trains from the same electrode). This observation is attributable to the stereotypical shape of the W–K filters across monkeys (Fig. 4B, bottom).
Is the relationship between spikes and LFPs specific to a given brain area or is there a universal mapping between these two scales of neural analysis? We addressed this question by considering physiological recordings from the macaque monkey ITC described by Kreiman et al. (2006). During the ITC recordings, awake monkeys passively viewed a stream of grayscale objects. Using the same approach as in Figure 3, the mean correlation between the estimated LFP and the actual LFP for 125 ITC electrodes was r̄ = 0.27 ± 0.17. Given the similarity between the V1 W–K filters and the ITC W–K filters, we asked whether we could use the V1 W–K filters to estimate the ITC LFPs based on the ITC spike trains. Using the general filter computed with the V1 data to estimate the ITC LFP from the ITC spike trains, we obtained a mean estimation accuracy of r̄ = 0.25 ± 0.12. In other words, using a general filter computed with data from V1, we obtained an estimation accuracy of 0.30 in the remaining half of the V1 electrodes and an estimation accuracy of 0.25 in the ITC data. This shows a remarkable degree of extrapolation and suggests a general relationship between spikes and LFPs across different neocortical areas.
In contrast to the extrapolation across conditions (spontaneous activity and stimulusdriven activity in V1) (Fig. 4A), monkeys (Fig. 4B), and neocortical areas (V1 and ITC) (Fig. 4A), the W–K filter showed a poor performance in the LGN recordings. In three recording sessions in two monkeys, eight additional electrodes were simultaneously positioned in the LGN. Using 95 trials recorded from the LGN, we estimated the LGN LFPs from the LGN spiking activity as described in Materials and Methods (supplemental Fig. S8, available at www.jneurosci.org as supplemental material) shows the estimation accuracy over all LGN trials using electrodespecific filters (compare with Fig. 3B). The algorithm performed substantially worse when trying to estimate LGN LFPs than when trying to estimate V1 LFPs: the average estimation accuracy for the LGN was r = 0.03. We observed that the STA of the LFP was essentially flat and noisy in the LGN, which further indicates the lack of a clear relationship between spikes and LFPs in the LGN.
Frequency analysis
The results presented above quantify the estimation of the entire LFP time course. We asked whether there were differences in the estimation accuracy across different LFP frequency bands. For this purpose, we filtered the LFP into different frequency bands between 0.1 and 100 Hz and computed a separate W–K filter to estimate the different frequency components (Fig. 5). The estimation accuracies were statistically significant compared with the null hypothesis for all frequency bands (p < 0.01, based on a twosample Kolmogorov–Smirnov test). The lower frequency band (i.e., 0.1–20 Hz) yielded the best estimation accuracy (p < 10^{−5} for the spontaneous condition and p < 10^{−5} for the stimulation condition, based on a twosample Kolmogorov–Smirnov test comparison against the null hypothesis).
These observations suggest that the spiking activity contains more information about low frequencies of the LFP than about high frequencies. This result is consistent with the results reported previously (Rasch et al., 2008). Rasch et al. showed that spikes seem to be locked at the onset of low frequency oscillations in the LFP and that the spikeLFP coherence level is higher for low frequencies than for high frequencies. Indeed, the spike–LFP coherence (Pesaran et al., 2002) can be used as a measure of the error in the estimations. High coherence in a given frequency band translates into small LFP estimation errors (for additional details about the relationship between W–K filtering and spike–LFP coherence, see supplemental data, available at www.jneurosci.org as supplemental material).
Robustness of the estimations to spike time jitter
Intuitively, the possibility of estimating the LFP time course from the spike train depends on the temporal structure of the spike train. This is further emphasized by the poor LFP estimations when using the null model consisting of a spike train with the same mean rate but random spike times. To further investigate the robustness of the LFP estimates to spike time distortions, we recomputed the W–K filter after adding temporal jitter to the spike trains. We created synthetic spike trains from the experimental ones by randomly jittering the spike times. The LFP was then estimated from these synthetic spike trains. For each spike train, the spike times were moved from their actual occurrence times by a random amount taken from a 0 mean Gaussian distribution with an SD σ_{jitter}.
The robustness of the estimations to time jittering was evaluated by plotting the mean estimation accuracy as a function of σ_{jitter} and computing the amount of temporal distortion required to cause a 50% drop in the estimation accuracy (_{50}σ_{jitter}) (Fig. 6). The value of _{50}σ_{jitter} was obtained by fitting the following equation: r = r_{0} − >, where r_{0} is the mean estimation accuracy when σ_{jitter} = 0, and _{50}σ_{jitter} and n are free parameters. Figure 6A shows the mean estimation accuracy (r̄) as a function of σ_{jitter} for an example electrode in V1. To compare the decrease in r with time jitter across electrodes, we normalized the estimation accuracy by r_{0} (Fig. 6B). The distribution of _{50}σ_{jitter} for all V1 electrodes is shown in Figure 6C; the mean value was 172 ± 89 ms.
 Download figure
 Open in new tab
 Download powerpoint
The robustness of the LFP estimate to time jitter was strongly correlated with the width of the optimal W–K filter. We found a linear relationship between the width of the optimal filter and _{50}σ_{jitter} (slope of 0.35, correlation coefficient of 0.64). The relatively long values of _{50}σ_{jitter} are consistent with the higher values of r̄ at low frequencies described in the previous section (Fig. 5).
Dependence on distance between spike trains and LFP electrode
Spikes constitute a measure of activity within a small vicinity of the electrode, on the order of 200 μm (Holt and Koch, 1999; Logothetis, 2002). In contrast, LFPs measure neural activity within a larger area (Mitzdorf, 1985; Kruse and Eckhorn, 1996; Juergens et al., 1999; Logothetis, 2002; Kreiman et al., 2006; Belitski et al., 2008; Liu et al., 2008; Nauhaus et al., 2009) (but see Katzner et al., 2009). Because recordings were performed simultaneously with multiple electrodes, we were able to study the spatial resolution of the relationship between spikes and LFPs by considering spiking activity and LFPs recorded from separate electrodes. For this purpose, we considered two electrodes separated by a distance D, and we used the spike train in one electrode to estimate the LFP in the other electrode according to Equations 2–4. In Figure 7, the estimation accuracy is plotted as a function of the distance between the spike electrode tip and the LFP electrode tip (D = 0 corresponds to the case in which LFPs and spikes were recorded from the same electrode as discussed in the previous sections). Figure 7A shows the results of this analysis for an example recording session where electrodes spanned a distance of up to 4 mm. We fitted the function r = r_{0} − , where D_{50} is a free parameter (Fig. 7A, dotted line). We show the distribution of D_{50} values in Figure 7B. To compare the decrease in r with distance across electrodes, we normalized the estimation accuracy by r_{0} (Fig. 7C). The estimation accuracy decreases by ∼50% when the two electrodes are 1 mm apart (this is the minimum electrode distance in these data). Beyond 1 mm, we found a weak dependence of the estimation accuracy with distance. This significant drop with distance was also observed when we considered only the 0.1–20 Hz frequency band of the LFP.
 Download figure
 Open in new tab
 Download powerpoint
Causal filters
In the analyses presented thus far, at any given time point t, both the spikes before t and the spikes after t contribute to the estimate of the LFP (see Eq. 2 and the shape of the W–K filters in Figs. 2, 3). We conjectured that there may be a temporal asymmetry such that the LFP estimate at time t using those spikes before t may yield a different estimation accuracy than those occurring after t. To evaluate this possibility, we constructed two different causal filters and compared their performance to the estimation accuracies from the noncausal (NC) filter used in the previous sections.
In the first case (C^{+}), we constructed a filter that was set to 0 for negative time lags [i.e., h_{C+}(α) = 0 for α < 0]. In the second case (C^{−}), we constructed a filter that was set to 0 for positive time lags [i.e., h_{C}_{−}(α) = 0 for α > 0]. In both cases, after imposing the causality constraint, the estimates were computed as described above in Materials and Methods. Figure 8 shows r for the three different filters, NC, C^{+}, and C^{−}, for 88 V1 electrodes during spontaneous activity. We found that the predictions computed with NC and C^{+} were statistically indistinguishable (p > 0.1, twosample Kolmogorov–Smirnov test). In contrast, the LFP estimations using the C^{−} filter were significantly worse than those based on either NC or C^{+} (Fig. 8) (p < 0.01). The estimations based on C^{−} were still significantly better than chance (p < 0.01, twosample Kolmogorov–Smirnov test).
Estimations based on multiple temporally clustered spike trains
It has long been assumed that local field potentials reflect the synchronous activity of ensembles of neurons (Mitzdorf, 1985; Fries et al., 2001; Pesaran et al., 2002; Berens et al., 2008). Based on this assumption, we hypothesized that spikes clustered in the time dimension (i.e., cooccurring within short time windows) across multiple electrodes could show an enhanced contribution to the LFP time course than isolated spikes (i.e., those spikes that do not cooccur with other spikes within short time windows). To test this hypothesis, we considered an electrode recording LFPs and all simultaneously recorded spike trains from nearby electrodes within a sphere of radius R. For a given time window τ, we separated those spike times that coincided within τ ms with spikes in any other electrode within the sphere (“timeclustered” spikes) and those spikes that occurred more than τ ms away from spikes in any other electrode within this sphere (“timeisolated” spikes). We used the same procedure in Equations 2–4 to estimate the W–K filter and the LFP based on these two types of spike models. In Figure 9A–C, we show that, overall, the timeclustered spikes (red circles) yielded higher estimation accuracies than the timeisolated spikes for all values of R and most values of τ. In those cases in which the timeisolated spikes yielded better estimation accuracies, we noted that there was a wide discrepancy in the number of timeisolated and timeclustered spikes (e.g., τ = 2 ms in Fig. 9A). Given the correlation between estimation accuracy and firing rate (supplemental Fig. S4A, available at www.jneurosci.org as supplemental material), we repeated the analysis restricting to those values of R, τ and those trials in which the number of spikes for the timeisolated spikes was within 20% of the number of spikes in the timeclustered spikes (Fig. 9D–F). We observed that timeclustered spikes yielded a significantly higher estimation accuracy for τ = 8 and τ = 16 ms (we could not test this hypotheses for smaller values of τ because in those cases there were always many more spikes in the timeisolated condition) (Fig. 9A–C). Overall, spikes that are temporally clustered with nearby spikes within ∼10 ms contain more information about the LFP time course than isolated spikes.
Discussion
There has been increased interest recently in local field potentials (Fries et al., 2001; Logothetis, 2002; Pesaran et al., 2002; Mehring et al., 2003; Bédard et al., 2004; Kreiman et al., 2006; Kraskov et al., 2007; Nir et al., 2007; Belitski et al., 2008; Nauhaus et al., 2009, among others) attributable to the correlation between LFPs and functional imaging measurements (Logothetis, 2002; Kayser et al., 2004; Mukamel et al., 2005), the use of LFPs as a measure of local synchrony (Fries et al., 2001; Womelsdorf et al., 2007), and the potential use of LFPs for prosthetic applications (Pesaran et al., 2002; Mehring et al., 2003). Additionally, understanding the biophysical origins of LFPs could provide insights into the computations performed in local cortical circuits.
The biophysics underlying the generation of LFPs is not clearly understood. Although several pieces of evidence suggest that the spatial extent of the LFPs is larger than the spatial extent of the spiking signals (Mitzdorf, 1985; Kruse and Eckhorn, 1996; Juergens et al., 1999; Logothetis, 2002; Bédard et al., 2004; Kreiman et al., 2006; Logothetis et al., 2007; Katzner et al., 2009; Nauhaus et al., 2009), the exact quantitative limits remain unclear. Additionally, current source density analyses and simultaneous recordings of LFPs and spikes have suggested that the LFPs are better understood in terms of synaptic potentials, afterpotentials, and other membrane potentials that show a slower temporal resolution than spikes themselves (Mitzdorf, 1985; Kamondi et al., 1998). Despite the large differences in spatial and temporal scales, our observations show that a linear filter operation on the spiking activity of one or a few neurons can explain a significant fraction of the variations in the LFP time course.
The shape of the W–K filter was highly stereotypical (Figs. 2, 3, 4B, 8 and supplemental Figs. S3, S7, S9, available at www.jneurosci.org as supplemental material; see also the spiketriggered average LFP in Fig. S7A, available at www.jneurosci.org as supplemental material): a sharp downstroke for negative time lags was followed by a slower upstroke for positive time lags. This shape is remarkably similar to the waveform shape in extracellular action potentials (see the experimental recordings as well as the computational simulations by Gold et al., 2006, their Fig. 1). The biophysics of the signals that give rise to the extracellular action potential have been characterized in multiple studies (Koch, 1999). For example, the simulations in the study by Gold et al. (2006) show that this prototypical shape can be generated by combining well known channels, most prominently fast inactivating Na^{+} channels and a variety of slower K^{+} channels. Although it is dangerous to extrapolate from the waveform of an extracellular action potential to the W–K filter in our study, it is tempting to speculate that part of the linear estimation that we report here might be accounted for by (weighted) linear and temporally smoothed averaging over many transmembrane potentials.
For an individual neuron, there can be strong nonlinearities in the generation of action potentials from the incoming synaptic currents and the integration of subthreshold dendritic events (Koch, 1999). The linear relationship that we report here concerns the estimation of an average over large numbers of dendritic events across many neurons. Furthermore, there is a topographical organization in V1 whereby neurons in the vicinity of the electrode registering the spiking activity may share similar properties. It is possible that this topography plays an important role in the structure of the LFP signal. If this conjecture holds, then we might expect that the relationship between spikes and LFPs may be more complex in areas that do not show such a strong topography as V1.
The local spiking activity analyzed here in the form of MUA consists of action potentials derived from multiple neurons. The exact number of neurons in the MUA is not known. However, in contrast with the LFPs, the spatial extent of the MUA decays rapidly with distance (Kreiman et al., 2006), suggesting that the MUA has a spatial extent of ∼200 μm and that the number of components is unlikely to be larger than hundreds of neurons. Additionally, we performed spike sorting on the MUA (supplemental Fig. S9, available at www.jneurosci.org as supplemental material). Even with the best available methods for spike sorting, we cannot guarantee that the SUA consists only of a single neuron, but the SUA signal is unlikely to consist of more than a few neurons (Lewicki, 1998). The small difference in the estimation accuracies using different spike cutoff thresholds (supplemental Fig. S9A, available at www.jneurosci.org as supplemental material) as well as the small difference in the estimation accuracies between SUA and MUA (supplemental Fig. S9B, available at www.jneurosci.org as supplemental material) further supports the idea that part of the LFP time course can be understood from highly local spiking signals. Therefore, it seems unlikely that the current estimation accuracies based on linear filtering are a consequence of the averaging of a noisy signal over tens to hundreds of neurons in the MUA. Carrying the argument to the extreme, even a single neuron might carry significant amounts of information about the LFP, a macroscopic property of the local circuit (Katzner et al., 2009).
We chose to use a linear filter because this transformation provides the simplest possible insights into the nature of the relationship between spikes and LFPs. It is conceivable that more complex nonlinear operations can account for an even higher fraction of the LFP time course. We emphasize that there remains a significant fraction of the LFP time course that is not accounted for by our linear estimates. The values reported throughout the text correspond to the mean estimation accuracies. There are cases in which the correlation between the linear estimate and the LFP was as high as ∼0.6 (see the distribution in Fig. 3 and supplemental Fig. S2B, available at www.jneurosci.org as supplemental material). This value is close to the best estimation accuracy of the LFP time course in one electrode from the LFP signal in a nearby electrode.
How generic is the map between spikes and LFPs? Several pieces of evidence argue that at least part of the LFP time course can be explained by a rather general linear function of the local spike trains. First, the estimation accuracies were quite stationary over time: for a given electrode, a filter computed in one trial performed quite well on a separate trial. Second, the estimations extrapolated across electrodes for a given monkey and even across monkeys (Fig. 4). Third, we also showed that a W–K filter computed using spikes and LFPs recorded from primary visual cortex could estimate the LFP time course in inferior temporal cortex using the inferior temporal cortex spikes (see Results, Spontaneous activity versus visual stimulation) (we further note that these are recordings performed in different labs, with different behavioral conditions, and using different tasks).
The generality of the linear W–K filter used here seems to break down outside of neocortex. First, the linear W–K filter did not work well when we attempted to estimate the LFP time course in LGN (supplemental Fig. S8, available at www.jneurosci.org as supplemental material) or in the human hippocampus (data not shown). It seems that, outside of neocortex, the relationship between these two signals is more complicated and might require nonlinear operations, extensive recordings from more nearby units, or building more complex models that incorporate other aspects of the circuit architecture. The LFP in the LGN during spontaneous activity may be generated or governed by modulatory inputs that are less related to the spiking of cells. Another nonexclusive possibility is that the special geometry of neocortex, with six layers and a columnar architecture, is an important component of the linear relationship between spikes and LFPs. As noted above, it is possible that topography also plays a role in the relationship between spikes and LFPs. These observations constrain future development of biophysical models of the origin of local field potential signals.
When we used causal filters (Fig. 8), we noted a significant asymmetry: the LFP estimations based on a filter that was constrained to be 0 for negative time lags (C^{+}) were significantly higher than those estimations based on a filter that was constrained to be 0 for positive time lags (C^{−}). Furthermore, the estimations based on the C^{+} filter were almost as high as the ones obtained when using the noncausal filter (Fig. 8).
The possibility of linearly estimating the time course of the local field potentials from spike trains during spontaneous activity extrapolates to conditions of visual stimulation (Fig. 4A). We also observed that the linear estimation holds in another brain area (inferior temporal cortex) under a passive viewing task in awake monkeys (see Results, Spontaneous activity versus visual stimulation) (Kreiman et al., 2006). Therefore, the basic map between spikes and LFPs holds for different brain areas (V1, ITC), stimulation conditions (spontaneous activity, visual stimulation), and behavioral states (anesthesia, awake). The estimation accuracies and, in particular, the estimation accuracies for different frequency bands may show a strong dependence on the experimental conditions. Several investigators have noted that the stimulus, task, and behavioral state of the animal (e.g., attention, shortterm memory) can influence the relationship between spikes and LFPs (Fries et al., 2001; Pesaran et al., 2002; Liu and Newsome, 2005; Buschman and Miller, 2007; Womelsdorf et al., 2007; Belitski et al., 2008). In particular, there are many cases in which spikes and LFPs (or particular frequency bands of the LFPs) can be decorrelated (Kreiman et al., 2006; Belitski et al., 2008; Berens et al., 2008).
It has been assumed that LFPs may reflect synchronous activity across ensembles of local neurons (Fries et al., 2001; Laurent, 2002; Logothetis, 2002; Womelsdorf et al., 2006; Montemurro et al., 2008). Our observation that temporally clustered spikes yield better estimation accuracies than temporally isolated spikes (Fig. 9) seems to be compatible with the idea that LFPs may reflect computations that are local both in space as well as in the time domain. However, the estimation accuracies were robust to significant amounts of spike time jittering (Fig. 6). More research is necessary to further our understanding of the biophysical signals that give rise to the LFPs, how many neurons are involved, and how neuronal events across local ensembles give rise to the field potential.
As in other systems, a multiscale analysis of neuronal circuits provides rich insights that cannot be achieved by focusing on a single level only. Many efforts in neuroscience aim to understand perception and behavior using only macroscopic signals such as blood oxygenation leveldependent measurements. As emphasized previously (Logothetis, 2008), it is not always trivial to interpret these measurements without a detailed understanding of the underlying architecture and cellular function. Other efforts in neuroscience correlate perception and behavior with the activity of single neurons. Circuits of neurons may show emergent properties that are not always easy to visualize by looking at individual neurons without studying their interactions (Koch and Segev, 1989). The combination of spikes and LFPs provides an ideal resolution to understand the relationship between individual neurons and local neuronal circuits.
Footnotes

This research was sponsored by the Whitehall Foundation, the Klingenstein Fund, National Institutes of Health Grant 1R21EY01971001, and the Children's Hospital Boston Faculty Development Award. We thank Fabrizio Gabbiani and Ulf Knoblich for comments on this manuscript and helpful discussions. We thank Stefano Panzeri, Rodrigo Quiroga, and Alberto Mazzoni for help with the spike sorting (supplemental Fig. S9, available at www.jneurosci.org as supplemental material). The V1 data used throughout the manuscript were recorded in the Logothetis laboratory. The inferior temporal cortex data shown in the last two columns of Figure 4A were recorded by Chou Hung and James DiCarlo.
 Correspondence should be addressed to Gabriel Kreiman, Children's Hospital, Harvard Medical School, 3 Blackfan Circle, Boston, MA 02115. gabriel.kreiman{at}tch.harvard.edu