Local field potentials reflect multiple spatial scales in V4

Local field potentials (LFP) reflect the properties of neuronal circuits or columns recorded in a volume around a microelectrode (Buzsáki et al., 2012). The extent of this integration volume has been a subject of some debate, with estimates ranging from a few hundred microns (Katzner et al., 2009; Xing et al., 2009) to several millimeters (Kreiman et al., 2006). We estimated receptive fields (RFs) of multi-unit activity (MUA) and LFPs at an intermediate level of visual processing, in area V4 of two macaques. The spatial structure of LFP receptive fields varied greatly as a function of time lag following stimulus onset, with the retinotopy of LFPs matching that of MUAs at a restricted set of time lags. A model-based analysis of the LFPs allowed us to recover two distinct stimulus-triggered components: an MUA-like retinotopic component that originated in a small volume around the microelectrodes (~350 μm), and a second component that was shared across the entire V4 region; this second component had tuning properties unrelated to those of the MUAs. Our results suggest that the LFP reflects neural activity across multiple spatial scales, which both complicates its interpretation and offers new opportunities for investigating the large-scale structure of network processing.


INTRODUCTION
Local field potentials (LFP) are low-frequency oscillations in the extracellular electric potential detectable through microelectrode recordings. Although they reflect a variety of electrical phenomena, the total synaptic current in a volume around the microelectrode is the major contributor to LFPs (Buzsáki et al., 2012). They thus offer a complementary signal to spikes, reflecting subthreshold network activity at larger spatial and longer temporal scales than are accessible through single unit recording.
The nature of the relationship between spikes and LFPs is a controversial issue that has attracted much interest (Bauer et al., 1995;Henrie and Shapley, 2005;Mukamel et al., 2005;Kreiman et al., 2006;Liu and Newsome, 2006;Nir et al., 2007;Belitski et al., 2008;Gieselmann and Thiele, 2008;Rasch et al., 2008;Katzner et al., 2009;Khawaja et al., 2009;Xing et al., 2009;Eggermont et al., 2011;Hwang and Andersen, 2011;Jia et al., 2011;Liebe et al., 2011;Lindén et al., 2011;Tsui and Pack, 2011;Zanos et al., 2011b;Lashgari et al., 2012). In a typical experimental scenario, a well-understood property of multi-unit activity (MUA) is compared and contrasted with that of the LFP. For instance, it has been established, by comparing the orientation tuning of V1 MUA and LFP activity, that the LFP activity in V1 could result from the spiking activity in a small volume around the electrode, on the order of 250 μm (Katzner et al., 2009). Xing et al. (2009) came to a similar estimate by comparing the size of MUA and LFP receptive fields (RFs) in V1. On the other hand, Kreiman et al. (2006) found that selectivity of LFPs for objects is best explained by hypothesizing an integration radius of a few millimeters.
The difference in the estimated integration radii across experiments may reflect the selection of different components of the LFP for analysis. For instance, the power in the high-frequency gamma band tends to be correlated with spiking activity (Ray and Maunsell, 2011) while the amplitude of the signal at lower frequencies has distinct tuning (Belitski et al., 2008); distinct components may have distinct integration properties (Berens et al., 2008). In addition, the physical size of the area under study, the regularity of its organization, and the correlation structure of its input (Lindén et al., 2011) can influence the properties of the LFP, and this may explain some of the discrepancies in the literature.
Many studies of the LFP have focused on V1 in particular, which is an unusual cortical area in that it is quite large, and it is known to have extremely precise columnar organization based on orientation selectivity (Ohki et al., 2005(Ohki et al., , 2006 and retinal position (Hubel and Wiesel, 1977;Blasdel and Fitzpatrick, 1984). Thus, LFPs in V1 may be unrepresentative of visual or sensory cortex as a whole. As a step toward understanding the LFP in cortex at large, we have recorded LFPs and MUAs in cortical area V4, a region that occupies an intermediate position in the visual hierarchy.
V4 recordings were carried out in two macaque monkeys, both of whom were implanted chronically with 96-electrode Utah arrays, which allowed us to relate the properties of the receptive fields to their physical location on the cortical surface. Using a sparse-noise stimulation procedure, we found that LFPs in V4 exhibit well-defined receptive fields whose positions change smoothly as a function of position on the cortical surface. However, a detailed analysis of the temporal properties of these signals revealed striking changes in RF position and size as a function of time following stimulus onset, such that the retinotopy of MUAs matched that of LFPs only at a restricted set of time lags.
These results could be explained by a model in which the LFP reflects multiple sources of input: a local, retinotopic input and a distant, shared input that had a similar effect across all of V4.
By fitting such a model to our data, we found that the local input exhibited consistent retinotopy that approximated that of the simultaneously recorded MUAs. The shared input arrived at latencies that differed from those of the retinotopic input and that differed substantially between the two animals. These results suggest that the LFP reflects neural activity across multiple spatial scales, which both complicates its interpretation and offers new opportunities for investigating the large-scale structure of network processing.

PRELIMINARY ANALYSIS
We used a sparse noise presentation paradigm to estimate the receptive fields of LFPs and multi-units (MUA). Sparse bar stimuli were flashed on a screen while the animal was rewarded for fixating a static red target ( Figure 1A). The stimuli were located on a log-polar grid and scaled in length and width proportionally to eccentricity to account for the scaling of neuronal RFs with eccentricity (Motter, 2009).
We first removed the remnants of individual spikes on each electrode, using a Bayesian spike removal algorithm (Zanos et al., 2011b). We then determined whether the amplitude of the LFP or its power in different frequency bands were modulated by the stimulus. We filtered the signal in narrow frequency bands and applied the Hilbert transform to get an estimate of the instantaneous power of the signal as a function of time (Freeman, 2007). In this preliminary analysis, we estimated the receptive fields (RFs) of LFPs using reverse correlation (De Boer and Kuyper, 1968;Marmarelis and Marmarelis, 1978). We used these RFs to predict the signal in a separate validation dataset (see section "Methods" for details). We obtained poor predictions (mean r = 0.03 for array 1, r = 0 for array 2) for all examined frequency bands (Figures 1B,C, red lines). This is likely due to the short duration of each stimulus; generally, power modulations tend to be visible after sustained stimulation (Khawaja et al., 2009).
We repeated the analysis using the amplitude of the LFP in different frequency bands. This yielded considerably better predictions (Figures 1B,C, blue lines) across the 5 lowest frequency bands examined, encompassing the range from 0.5 to 30 Hz (mean r = 0.24 for array 1, r = 0.14 for array 2). Note that shaded error bars represent ±2 s.d.; the lack of overlap in the lowest 5 frequency bands indicates that the worst fits to the amplitude of the LFPs are nevertheless better than the best fits to the power of the LFP. Thus for the following analyses, we focused on the amplitude of the LFP filtered between 0.5 and 40 Hz. Figure 1D shows a slice of a typical LFP RF, capturing the selectivity for space and orientation at a lag of 70 ms following stimulus onset. Casual inspection reveals little modulation of the FIGURE 1 | Sample stimulus and receptive field. (A) Sample stimulus frame. The stimulus was composed of sparse bars arranged on a log-polar grid to account for the magnification of the visual field with eccentricity. (B) Quality of reverse correlation fits to the amplitude of the LFPs (blue lines) and power (red lines) in different frequency bands for the first array. Delta: 0.5-4 Hz, theta: 4-8 Hz, alpha: 8-12 Hz, low beta: 12-20 Hz, high beta: 20-30 Hz, low gamma: 30-50 Hz, high gamma: 50-80 Hz. Quality of fit was evaluated by the correlation (r) of the predicted and measured responses in a validation dataset. Shaded error bars represent ±1 s.d. Power is not modulated by the stimulus; rather, the stimulus influences the low-frequency components of the amplitude of the LFPs. (C) Same as in (B), for the second array. (D) Sample LFP RF estimate. The RF of the LFP is measured at a 70 ms time lag relative to stimulus onset. Each square represents the spatial RF for a given orientation. The shape of the spatial RF varies little with orientation. The color bar indicates peak z-values estimated through bootstrapping. Here, as in all subsequent RF illustrations, a Gaussian smoothing kernel with σ = 0.7 is applied. (E) Separable RF estimate. The RF shown in (D) is approximated as separable in space (left) and orientation (right). Little information is lost in the process, and z-values for the spatial envelope of the RF are markedly increased (see legend of color bar).

Frontiers in Computational Neuroscience
www.frontiersin.org March 2013 | Volume 7 | Article 21 | 2 spatial structure of the RF with orientation, aside from a scale factor. As this pattern of results was typical of both our LFP and MUA recordings, we assumed for the rest of the analyses that orientation tuning could be treated separately from space-time selectivity. As shown in Figure 1E, the consequent reduction in free parameters led to more reliable receptive field estimates (peak z-values shown next to color bars). LFP spatial RFs were, however, strongly modulated as a function of time following stimulus onset. Figure 2A (top) illustrates the spatial RF of an LFP measured on one array at three different time lags. The LFP was responsive to a large area of the visual field at early time lags (80-100 ms) and a smaller area at later time lags (140 ms). This pattern of changes was typical of LFPs measured on this array, as shown below in more detail. By contrast, the MUA RF measured on the same electrode (Figure 2A, bottom) showed little evidence of such a change in size.
In addition to changes in size, LFP RFs frequently appeared to shift their preferred positions as a function of time. Figure 2C illustrates the receptive field of an LFP typical of the second array. Between 70 and 90 ms time lags, the RF shifted its preference toward high eccentricities. At 120 ms, the LFP responded to stimuli at a foveal location, far from the initial RF peak (peak z-value: 9.2; |z| = 4.5 corresponds to p = 0.001, corrected for multiple comparisons). Note that the polarity of the response reversed with time; while eccentric stimuli caused a low-latency negative deflection in the LFP signal, foveal stimuli caused a positive deflection at longer time lags. Again, the MUA RF ( Figure 2C, bottom) showed no evidence of such a change.
To quantify these effects, we fit each RF time slice with a Gaussian, which captured the selectivity of the RF with four parameters: preferred eccentricity, preferred angle, radial size, and angular size. We estimated the uncertainty in these parameters through bootstrapping. Figures 2B,D show the changes in RF position (top) and size (bottom) for the example LFPs (blue) and MUAs (red). LFP RF parameter changes were highly significant across time lag (blue lines; shaded error bars correspond to 95% confidence intervals). By contrast, parameters for MUAs were comparatively stable across time lag (red lines), partly as a consequence of their shorter duration.
Hence, while the LFPs in the examples were well tuned for space, their spatial RFs were not static over time. As a result, the position and size of LFP receptive fields diverged substantially from those of corresponding MUAs at some time lags.

ARRAY ANALYSIS
Given the stability of MUA receptive fields across time lags (Figure 2), we refit the data on the assumption that their RFs were separable in time and space (see section "Methods"). We then Here, the LFP RF shows a late, foveal excitatory region (120 ms) far from the initial, peripheral preference (70 ms). plotted the estimated preferred eccentricity and polar angle of each MUA according to their location on the array. Figure 3A shows the result for a single MUA; here the RF was found at roughly 15 • eccentricity along the vertical meridian in the lower visual field (bottom plot). These values were color-coded separately and displayed at the location of the recording electrode within the array (cross-shaped outlines). For eccentricity (left), blue colors corresponded to small values (foveal locations), while red corresponded to high values (eccentric locations). The same color similarly mapped the range of polar angle (right) from 270 • (lower visual field) to 360 • (right visual field), which spanned the ensemble of RFs we recovered for this array. Repeating this process for every electrode revealed the retinotopy of MUAs across 4 × 4 mms of cortex ( Figure 3B). Here eccentricity and angle for each electrode are coded as described above, with electrodes that did not yield significant RFs represented in white. As expected, preferred eccentricity and polar angle changed smoothly as a function of position on the array, and the direction of the eccentricity gradient, illustrated by a green arrow, was roughly orthogonal to that of polar angle.

Frontiers in Computational
LFP RFs exhibited a similar retinotopic organization, as illustrated in Figure 3C for the same array. As expected from the examples shown in the previous section, LFP retinotopy changed as a function of time lag. These changes were coherent across the array: at later time lags (140 ms), a greater proportion of the array responded to foveal locations (as shown by the increased representation of dark blue colors, left plot) and angles near 270 • (again shown in dark blue, right plot). Figure 3D shows that the preferred angle and eccentricity of LFPs best matched those of MUAs at a time lag of 140 ms, as measured by the root-mean-squared (RMS) discrepancy. As suggested by the example shown in Figure 2A, the mean LFP receptive field size changed dramatically as a function of time lag ( Figure 3E). LFP RFs appeared larger at earlier lags, shrinking in size to a value close to that of the mean MUA receptive field size at longer time lags (dashed line).
Corresponding results for the second array are shown in Figure 4. In this case LFP RFs formed a retinotopic map at early time lags (Figure 4B, top), consistent with that of MUAs, with the best match occurring at 80 ms ( Figure 4C). Strikingly, however, a foveal component appeared at later time lags, overtaking the retinotopy of the LFPs completely by 170 ms (Figure 4B, bottom). Thus, the majority of LFPs measured in the second array showed biphasic receptive fields similar to that illustrated in Figure 2C. The changes in retinotopy were accompanied by modest changes in mean RF size ( Figure 4D).
These results show that the LFP retinotopy changes with time lag in a concerted fashion across the cortical surface. While at some time lags the LFP retinotopy matched that of the MUAs, at others it considerably diverged. Thus, LFP RFs reflect more than the underlying retinotopy of MUAs.
Furthermore, the relationship between LFPs and MUAs was qualitatively different between the two arrays. For the first array, the retinotopy of the LFPs best matched those of MUAs at late time lags (140 ms); in the second array, LFP RFs were aligned with MUAs at early time lags (80 ms). In addition, our second array showed an array-wide foveal component unseen in the first array.

ROBUSTNESS OF RETINOTOPY
The striking differences between MUA and LFP receptive fields within an array and in the results between arrays could conceivably be a signature of a transient electrical artifact; that is, a source of noise that contaminated recordings on a particular recording day, such as line noise, reward artifact, cross-talk between electrodes, etc. To examine this, we repeated the analyses for data recorded on another day in each array. LFP RFs estimated on other recording days exhibited the same qualitative pattern of shift in retinotopy across time lag characteristic of each array. The day 1 LFP-day 2 LFP RMS discrepancy, averaged across time lags, was 0.8 for both arrays. By contrast, at the optimal temporal lag, the within-day LFP-MUA RMS discrepancy was ∼2 (Figures 3D, 4C). The tuning of the LFP is thus consistent across days.
It is also possible that signal processing exaggerated the differences between the two signals. The results presented in Figures 3B and 4A reflect MUA obtained by full-wave rectifying a bandpass filtered (750-3500 Hz) signal; we refer to this as the rMUA (cf. Xing et al., 2009). The MUA is also commonly defined by the density of threshold crossings in band-pass filtered voltage traces (Katzner et al., 2009); we term this the tMUA. These different definitions could potentially isolate different components of the signal (Supèr and Roelfsema, 2005). We thus repeated our analyses for the tMUA (see section "Methods" for details).
We found similar retinotopies with both measures of multiunit activity, with rMUA-tMUA RMS discrepancies of 0.4 and 0.9 for arrays 1 and 2, respectively. We found fewer significantly tuned electrodes with the threshold method, however, especially in the second array (tMUA: N = 57 in array 1, N = 23 in array 2; rMUA: N = 65 in array 1, N = 69 in array 2). We found that tMUA RFs were slightly (7%), but significantly smaller than rMUA receptive fields (p < 0.05 for each array, two-sided Wilcoxon rank sum test). These results are consistent with the rMUA having similar properties to the tMUA, while integrating over a slightly larger cortical area.
Thus, the changing retinotopy across time lags and the qualitatively dissimilar properties of the MUA and LFPs are unlikely due to transient noise sources or to the choice of data preprocessing for the MUA (above) or the LFP (Figures 1B,C).

TEMPORAL MIXTURE MODEL
The previous section showed that at some time lags, LFP RFs are organized in a retinotopic fashion similar to MUAs. Yet, this close correspondence between LFP and MUA retinotopy is broken at other time points. In the second array, in particular, the appearance of a foveal component at late time points ( Figure 4B) hints at the interplay between an MUA-like retinotopic component, which changes from electrode to electrode, and a component tuned for foveal locations, shared by all electrodes. Figure 5 illustrates how the interplay between these two mechanisms might account for the data. LFP RF time slices measured on two electrodes on the same array are plotted in Figure 5A. While at early time lags (left and center), the RFs are markedly different, they have similar shapes at later time points (right). These results are explained in Figure 5B by the interplay of electrode-specific components (green lines) and a second component (black lines) shared by all electrodes on the same array.
Here the changing receptive field positions result from electrodespecific response components that are stronger at early time lags and a shared component that is stronger at later time lags.
More formally, we assumed that the RFs f e τ,r,θ,o as a function of electrode number e, time lag τ, eccentricity r, angle θ and orientation o were given by: Here p r,θ,o is a component shared by all electrodes on a given array and q e r,θ,o is specific to each electrode; they are weighted differentially depending on electrode number and time lag by factors a e τ and b e τ . The shared component p r,θ,o could take on any spatial configuration. The electrode-specific component q e r,θ,o was constrained to have a Gaussian spatial RF profile and a separable orientation tuning curve.
We fit the temporal mixture model for each array by minimizing the squared error between the model and the data (see section "Methods" for details). The resulting model yielded a highly significant improvement over a baseline model without the constant component (0.57 vs. 0.45 R 2 , p 0.001 for array 1; 0.37 vs. 0.30 R 2 , p 0.001 for array 2; F-test). Based on the model fits, we reconstructed composite LFP RFs that were then fit with Gaussians at every time point to reconstruct retinotopies. These are illustrated in Figure 6 for array 1. The simulations replicated the pattern of increased representation of low eccentricity and 270 • locations at longer time lags ( Figure 6A) and the apparent decrease in receptive field size in time ( Figure 6B).
The underlying mechanism for this switch is illustrated in Figure 7. The shared component was triggered by stimuli of any orientation across a fairly broad region of space ( Figure 7A), with peak selectivity at central locations (∼10 • eccentricity, 315 • polar angle). In absolute terms, both the shared and retinotopic components were strongest at early time lags ( Figure 7B). However, because the retinotopic component decayed more slowly, it was relatively stronger at late time points (Figure 7B, bottom). It follows that at early lags, the observed RFs were both more broadly spatially tuned and biased toward representing central locations than at later lags.
Similar results are shown for the second array in Figure 8. in receptive field size with increasing lag (Figure 8B). Figure 9A shows that this shared component was strongly tuned for a foveal portion of the visual field. As with the first array, the retinotopic component peaked at early time lags (Figure 9B, middle); unlike the other array, however, the shared component manifested itself mostly at later time lags (Figure 9B, top). Together, these results explain the observed changes in retinotopic organization in terms of a gradual switch in the importance of two distinctly tuned components. The first, retinotopic component was strongest in both arrays at early time Both the median shared gain (top) and the retinotopic gain (middle) peak at early time lags. However, the retinotopic gain decays more slowly as a function of time. Therefore, the retinotopic gain is relatively larger at late time lags (bottom). This creates a shift in the representation from broadly tuned (shared component) to more tightly tuned (retinotopic component). Shaded error bars represent 95% confidence intervals for the median.
lags, while the second, shared component differed qualitatively between the two arrays. Thus, shared components may represent idiosyncratic, large-scale biases in visual representation (Jia et al., 2011), a matter we explore in more detail in the discussion.

RETINOTOPIC COMPONENT
The extracted retinotopic components are illustrated in Figures 10A,B; they were retinotopically arranged in a manner similar to MUAs (Figures 3B and 4A). This link is shown in more detail in Figures 10C,D (top row) and sizes (bottom row) of LFP and MUA RFs measured on the same electrodes. The eccentricity and polar angle of the extracted RFs were similar to those of MUAs measured at the same location, save for a cluster of observations at the bottom right of the second array ( Figure 10B). These electrodes had a retinotopic component that was foveal and thus overlapped with the shared component; this made precise estimation of their location and size problematic. Consistent with previous literature, retinotopic LFP RFs were larger than corresponding MUA receptive fields (Figures 10C,D); mean and median sizes are documented in Table 1. We estimated the integration radius of the LFP using the method introduced in (Xing et al., 2009). This involves first estimating the integration radius in visual coordinates, then translating this into A B Figure 7. The shared RF is tightly tuned for foveal locations. The gain of the shared RF grows larger at late time lags. Therefore, the RFs switch from an early retinotopic to a purely foveal late representation. cortical coordinates σ cLFP using the estimated cortical magnification factor (m), according to the formula:

FIGURE 9 | Temporal mixture model parameters-Array 2. (A) and (B) as in
Here σ vLFP and σ vMUA correspond to the mean size of LFP and MUA RFs in visual coordinates, respectively, and σ cMUA is the integration radius of the MUA in cortical coordinates. Because the retinotopy of V4 is less regular than in V1, magnification can change depending on the position on the cortical surface. We therefore estimated m by averaging the magnitude of the gradients of eccentricity and polar angle across the array (see section "Methods" for details   integration of activity proximal to the electrode, consistent with previous results in V1 (Katzner et al., 2009;Xing et al., 2009).

ORIENTATION AND TEMPORAL TUNING
Additional information about the relationship between the MUA and the LFP may be gained by comparing the orientation and temporal tuning of the two signals. We found that the LFPs were essentially untuned for orientation, with a mean circular variance (CV; Ringach et al., 2002) of 0.95 for the retinotopic component (minimum CV: 0.88) and 0.94 for the shared component across both arrays. On the other hand, some MUAs were tuned for orientation (min CV: 0.65), with a mean CV of 0.87 for significantly tuned MUAs across both arrays. This data is consistent with the idea that the LFP integrates over a larger area than the MUA, although the poor tuning prevents further analysis of the integration radius in the manner of Katzner et al. (2009). More interesting is the temporal tuning of both signals. Figure 11A illustrates the temporal filters of significantly tuned MUAs (blue lines; 50-240 ms) as a function of their position on the first array. Temporal filters are stereotyped through the array, with a rapid rise followed by a slower decay. There is some indication of suppression at late time lags (segments below the gray line). The majority of filters have a peak latency of 70 ms (Figure 11C), with a minority having a peak around 120 ms. These results are mirrored in array 2 (Figure 11B), where the filters are also highly stereotyped, although here the decay appears faster. The peak latency is also centered around 70 ms (mean: 73 ms; Figure 11D).
These results contrast strongly with the time filters of the retinotopic component of the LFP in array 1 (Figure 11E), in which the temporal filters differed dramatically across the array. Moreover, the filters are of considerably longer duration than the corresponding MUA filters, reflecting the fact that the signal is modulated at low frequencies ( Figure 1B). While the dominant polarity is negative, some filters have roughly equal positive and negative polarity phases (bottom center) or have mostly positive polarity (right middle). The peak latency occurs late compared to the MUA, and varies widely (mean: 118 ms; Figure 11G). Similar trends are visible for array 2 (Figure 11F), with filters showing large variation in shape and polarity. Peak latency is also more broadly distributed than the corresponding MUA data, and longer, with a mean of 90 ms ( Figure 11H). Thus, while MUAs and the retinotopic component of LFPs have similar retinotopy (Figure 10), the two signals diverge strongly in terms of temporal selectivity. LFPs have more sluggish dynamics, longer duration, and longer peak latencies; they are also more variable in shape and polarity than MUAs. One must therefore be careful in interpreting the LFP as a spatially smoothed version of the MUA, since the LFP does not reflect the MUA per say, but rather subthreshold activity (Buzsáki et al., 2012). While in some instances, like retinotopy, the relationship between subthreshold activity and the MUA is sufficiently well understood to make the relationship between LFP and MUA transparent (Carandini and Ferster, 2000), in the case of temporal tuning the relationship is complex, and the LFP and MUA reveal themselves as highly distinct.

GENERAL DISCUSSION
The local field potential is a complex signal which offers a window into cortical processing at larger spatial and longer temporal scales than those associated with single units. Components of this signal have been shown to correlate with attention (Fries et al., 2001(Fries et al., , 2008Gregoriou et al., 2009), cortical inhibition (Henrie and Shapley, 2005;Atallah andScanziani, 2009), arousal (Van Swinderen et al., 2004;Andretic et al., 2005), synchronicity (Gray and Singer, 1989;Mukamel et al., 2005;Nir et al., 2007), and other network phenomena.
While the LFP has proven a highly interesting signal, its interpretation has been marred by our lack of understanding of its biophysical sources and its relationship to spikes. Action potential generation and passive propagation have been well understood for several years (Hodgkin and Huxley, 1952;Koch, 1999). By contrast, LFPs have only recently been modeled in a biophysically detailed fashion (Bedard et al., 2006;Milstein et al., 2009;Lindén et al., 2011;Buzsáki et al., 2012). Modeling studies have unequivocally concluded that the LFP is an intrinsically more complex signal than spikes, reflecting a variety of distinct electrical phenomena (Buzsáki et al., 2012). Lindén et al. (2011) show that the integration radius of the LFP depends both on cortical layer and the correlation structure of the input.
It follows that the structure of the LFP and its relationship to spikes may well vary from area to area in idiosyncratic and unpredictable ways. We thus set out to estimate the receptive fields of LFPs in area V4 of two macaques and compared their properties to those of MUAs. Our results show that, in the context of a sparse noise presentation paradigm (reverse correlation), where the LFP signal is dominated by transient as opposed to sustained activity, the LFP reflects multiple sources of inputs.
In both our subjects, one component of the LFP reflected the underlying retinotopic organization of the cortical sheet (Figure 10). This component was strongest at 80-90 ms following stimulus onset (Figures 7B and 9B), and decayed slowly to baseline at 200-250 ms. It arose from the integration of activity within a patch of cortex of ∼350 um. Retinotopic signals were mixed with components shared by all electrodes on a given array. In the case of our first array, the shared component was broadly tuned for space and peaked at early time lags (Figure 7). For the second array, the shared component was tuned for foveal locations and peaked at late time lags (Figure 9).
What is the source of the shared component? We used state-ofthe-art signal processing to eliminate potential signal distortion by analog filters and spike remnants (Nelson et al., 2008;Zanos et al., 2011b; see section "Methods" for details). While it remains possible that the shared component is artifactual, its tuning properties are inconsistent with distortion caused by faulty grounding, for example. In the second array, in particular, we see that the shared component has temporal tuning properties which are very different from MUAs or early LFPs recorded on the array. Careful inspection of the shared component measured at late time lags (Figure 4B, 170 ms) shows that there is a small but visible gradient in angular selectivity from the left to the right of the array; this gradient is not captured by the temporal mixture model ( Figure 8A).
Hence, the shared components actually change across the array, albeit more modestly than the retinotopic components. We thus hypothesize that the shared components reflect largescale biases in the input to area V4. Jia et al. (2011) found that one component of low-gamma LFPs in V1 have similar orientation tuning across 4 mm of cortex, independent of the preference of local MUAs. It has been hypothesized that this reflects large scale biases in orientation representation in striate cortex, where orientations aligned with the preferred polar angle of neurons are slightly overrepresented (Freeman et al., 2011). Such large-scale biases in representation, which have no functional role per say, could vary idiosyncratically from animal to animal. We conjecture that this could explain the sharp difference between the arrays in the tuning of the shared common component.
Another potential source of discrepancy in the tuning of the shared component lies in the sampling of cortical layers. Different layers could be targeted as a function of position due to the curvature of cortex. The changing polarity of the temporal filters in Figures 11E,F is consistent with this idea. The second animal was much smaller (5-6 kg) than the first animal (9-10 kg). Consequently, the distance between the lunate and superior temporal sulci at the level of the implant was smaller in the second animal (4-5 mm vs. 6-7 mm) and the cortical surface was more curved. Thus, it is unlikely that the sample of layers is exactly the same in both animals, although this could not be verified histologically. These factors highlight that as an epiphenomenal signal (although see Anastassiou et al., 2011), the LFP is complex and noisy, and care must be taken in its interpretation.

THE INTEGRATION RADIUS OF THE LFP
These results may bear on the continuing debate regarding the integration radius of the LFP. We and others have demonstrated that different components of the LFP have tuning properties which are consistent with either local (∼300 um) or global (several mms) integration of MUA activity (Kreiman et al., 2006;Liu and Newsome, 2006;Katzner et al., 2009;Xing et al., 2009;Jia et al., 2011). We thus conclude that it is not meaningful to speak of the integration radius of the LFP, for the simple reason that the LFP reflects multiple sources of inputs with different integration scales.
We conjecture that the discrepancies in previous studies are due to signal processing and analysis choices which enhance the relative strength of one component of the LFP over another. An important distinction is that some studies (Katzner et al., 2009;Xing et al., 2009) have examined the amplitude of the LFP, similar to what is done here, while others (Kreiman et al., 2006;Liu and Newsome, 2006;Jia et al., 2011) have examined power at a selected frequency. The amplitude of the LFP is sensitive to transients in the local field potential, while power is sensitive to sustained oscillatory activity. Lashgari et al. (2012) have found that transient and sustained LFPs in V1 have markedly different tuning properties; this translates into changing relationships with MUA activity. It will be interesting to compare the retinotopy of amplitude and power components of the LFP directly in a paradigm which triggers both stereotyped deflections and oscillatory activity; natural movies may be able to evoke both phase and power modulations (cf. Figure 6, Rasch et al., 2008).
Another potential source of variability in the properties of the LFP may result from the treatment of the dimension of time. Katzner et al. (2009), for instance, analyzed the first component of the Singular Value Decomposition of the temporal-orientation tuning curve. Xing et al. (2009) instead analyzed responses at a latency corresponding to the peak deviation of the LFP. Such choices would not permit analysis of the multi-component temporal responses of the kind we have reported here (Figure 2).
It may well be the case that in V1 the major contribution to the LFP is separable with respect to time lag and highly local. In this respect, V1 may be a special case, as its retinotopy is remarkably precise (Hubel and Wiesel, 1977;Blasdel and Fitzpatrick, 1984;Ohki et al., 2005Ohki et al., , 2006; by contrast, higher-level areas have less precise retinotopy. Tonotopy in primary auditory cortex is significantly less precise than retinotopy in V1, at least in rodents (Castro and Kandler, 2010). Interestingly, LFP spectrotremporal receptive fields are much more broadly tuned than those of MUAs in both cat (Eggermont et al., 2011) and monkey (Kajikawa and Schroeder, 2011).
Given that the integration radius of the LFP varies with the correlation structure of the input (Lindén et al., 2011), part of the disagreement may reflect genuine inter-areal differences in the LFP. While this complicates the interpretation of the LFP, it may afford an opportunity to study how input correlation structures are forwarded and modified in a hierarchy, with wide implications for our understanding of encoding and decoding neural activity (Averbeck et al., 2006).

TASK
The recording methods have been described in detail previously (Zanos et al., 2011a). Briefly, we implanted chronic microelectrode Utah arrays in area V4 of two macaques (Macaca mulatta). Area V4 was identified based on stereotactic coordinates and anatomical landmarks (Ghose and Ts'O, 1997). After recovery, the monkey was seated comfortably in a primate chair (Crist Instruments) and trained to fixate for liquid reward. Eye position was monitored at 200 Hz with an infrared camera (SR Research). 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.

SIGNAL ACQUISITION AND PROCESSING
We recorded wideband signals at 10 kHz (bandpass filtered in hardware between 0.07 and 2500 Hz) over the 96 channels of each Utah array. We monitored the power spectrum of recorded wideband signals on a daily basis to minimize line noise and other artifacts. Recordings were referenced against a ground located on the skull 2-3 cm away from the array. The same recording equipment was kept in place for both animals.
The wideband signal was band-pass filtered between 750 and 3500 Hz, rectified, band-passed between 2 and 40 Hz, and downsampled to 100 Hz to form the MUA signal (Xing et al., 2009). For comparison with previous literature, we also computed an alternative MUA based on applying a low threshold (3σ) to the wideband signal bandpassed between 750 and 3500 Hz (Katzner et al., 2009). We used a detection deadtime of 1 ms, a two-sided threshold, and binned the events at 100 Hz (10 ms time bins); this gave similar results to the rectification-based method (see section "Results," Robustness of retinotopy). Action potential artifacts were removed from the wideband signal using a Bayesian method (Zanos et al., 2011b); the despiked wideband signal was then downsampled and band-pass filtered to produce the LFP (see Preliminary Analysis and Receptive Field Estimation section for filters specific to each analysis).

STIMULUS
Each animal was trained to fixate a red spot (2 • fixation window) while sparse dark bar stimuli were flashed on a uniform gray screen. The monitor was refreshed at a rate of 75 Hz; stimuli changed every odd frame (37.5 Hz); and each bar stimulus lasted for 6 monitor frames (12.5 Hz; 80 ms). Stimuli were presented in a single continuous trial, which lasted 25 min for the data presented for the first array and 30 min for the data presented for the second array. We repeated the experiment on other recording days for each array, with similar results (data not shown). The bar stimuli were placed along a 12 × 12 polar grid (Figure 1A), such that stimuli at the periphery were longer and wider than those near the fovea; bars scaled linearly with eccentricity. The length of the bars was chosen so that no bars touched when presented simultaneously; bar width was set to 0.25 times the eccentricity. Four different orientations were used. The grid was confined to the lower right corner of the screen. It spanned 120 • of polar angle and 5-40 • of eccentricity for array 1 and 3-50 • of eccentricity for array 2. On average, 7 bar stimuli were on the screen at any given time.
We split the data into a fit dataset and a validation dataset: for each 5 s block of the data, the first 4 s were assigned to the fit dataset and the last second to the validation dataset. We estimated orientation-spatial-temporal receptive fields for both band-pass filtered LFPs and their power using standard reverse correlation on the fit dataset (De Boer and Kuyper, 1968;Marmarelis and Marmarelis, 1978; Figures 1B-E). We applied a Gaussian spatial kernel (σ = 0.8) to the estimated receptive fields and predicted the signal in the validation dataset, assuming a linear model with the estimated filter.

RECEPTIVE FIELD ESTIMATION
The results of the preliminary analysis (previous section, Figures 1B,C) showed that the stimulus modulated the amplitude of the LFP but not its power, and that the modulation was concentrated at low frequencies. For the remaining analyses, we thus band-pass filtered the despiked wideband signal in the range (0.5-40 Hz) and downsampled to 100 Hz to produce the LFP signal. Inspection of the reverse correlation filters led us to a low-dimensional parameterization for each temporal slice of the LFP receptive fields: the selectivity of the RF is given by the product of an orientation filter and a Gaussian spatial envelope. Specifically, we assumed that the contribution of the stimulus presented τ epochs ago to the internal response in the kth time bin was given by: s k−τ,r,θ,o is the stimulus presented τ epochs ago. d is a bias. G a (b, c) is a Gaussian curve evaluated at a, centered around b, with width (standard deviation) c. v o is the orientation filter. We fit the model for each time slice τ = 1 to 24 by least-squares. We initialized the parameters by fitting a Gaussian to the spatial envelope of the reverse correlation estimate of the filter through least-squares; this envelope was determined by taking the first singular vector of the SVD of the reverse correlation estimate (Ahrens et al., 2008). The model for MUAs was similar, but this time we assumed that the receptive fields did not change in shape across time slices, but were simply scaled by a time-dependent gain: Here u corresponds to the weights of a separable time filter. We fit the model through least-squares. An MUA RF was deemed significantly tuned if the fit was significant at the p < 0.0001 level according to a χ 2 test (Wood, 2006). We found this criterion too lenient for LFP RFs, presumably because the correlation structure of the LFP did not follow the assumptions of the test. Instead, an LFP RF fit was deemed significant if its R 2 value was greater than the observed R 2 values on any electrode on the same array for time slices from 10 ms to 40 ms. Importantly, this simulation-based method preserves the correlation structure of both data and input, while eliminating the relationship between the two signals (Goldfine et al., 2013), and the resulting threshold corresponds to an effective p ∼ 0.01.

TEMPORAL MIXTURE FIT
The temporal mixture model illustrated in Figure 5 was as follows. The LFP RF g e τ,r,θ,o measured on electrode e was assumed to be a noisy version of the underlying RF f e τ,r,θ,o . The underlying RF was given by a mixture of a shared component p r,θ,o and a component specific to the electrode q e r,θ,o : The shared component was unconstrained while the specific component was a Gaussian in space modulated by orientation: The temporal mixture model was fit using an iterative leastsquares algorithm to minimize the mismatch between f e τ,r,θ,o and g e τ,r,θ,o , which was estimated by reverse correlation. The shared component was initialized to the mean of all RFs at t = 90 ms for array 1 and t = 200 ms for array 2. a e τ was then set by leastsquares on the assumption that b e τ = 0. Then the following steps were alternatively repeated until convergence: 1. The residual g e τ,r,θ,o − a e τ p r,θ,o was reshaped into a matrix with {e,τ} in one dimension and {r,θ,o} along the second dimension. The first singular values of this matrix were used to determine b e τ and q e r,θ,o . q e r,θ,o was then fit to Equation (6). 2. The first singular values of g e τ,r,θ,o − b e τ q e r,θ,o were used to determine a e τ and p r,θ,o .
Once the temporal mixture model was fit, we extracted f e τ,r,θ,o for each electrode and time lag and fit the reconstructed RF as an orientation filter multiplied by a Gaussian spatial envelope. The parameters determined through this process are plotted in Figures 6 and 8.

ESTIMATION OF THE INTEGRATION RADIUS OF THE LFP
We applied the method of Xing et al. (2009) to estimate the integration radius of the LFP on the cortical surface σ cLFP : σ vLFP and σ vMUA correspond to the size of LFP and MUA RFs in visual coordinates, respectively, and σ cMUA is the integration radius of the MUA in cortical coordinates. We used σ cMUA = 100 μm (Xing et al., 2009) and estimated (σ 2 vLFP − σ 2 vMUA ) by taking the 20% trimmed mean of this quantity for electrodes where we could measure both LFP and MUA receptive fields.
The cortical magnification factor m, measured in millimeters per unit of visual space, captures the change in visual coordinates that corresponds to a unit change in position on the cortical surface. Thus, the local cortical magnification factor corresponds to the inverse of the magnitude of the gradient of the visual quantity measured (log eccentricity or polar angle in our case). Unfortunately, estimating the magnitude of the gradient of the retinotopies of the MUAs directly is infeasible due to missing measurements. The measured retinotopies were too irregular to be fit reliably with a simple surface such as a plane. Instead, we obtained a smoothed estimate of log eccentricity and polar angle using Gaussian Process Regression (Rasmussen and Williams, 2006). We used the parameters suggested by the GPML for Matlab toolbox manual (Rasmussen and Williams, 2006; Gaussian likelihood, isometric squared exponential covariance, linear + constant mean function, marginal likelihood optimization for hyperparameters, exact inference).
We then computed the magnitude of the gradients of the smoothed surfaces, took their average across the surfaces, and inverted them to obtain the cortical magnification factor for log eccentricity and polar angle for each array.

ORIENTATION AND TEMPORAL SELECTIVITY
We evaluated the orientation selectivity of the retinotopic component of the LFP and the MUA by computing the circular variance of the orientation selectivity coefficients [v o in Equations (4) and (6)] as follows (Ringach et al., 2002): For the non-retinotopic component of the LFP, the same formula was used, with v o being estimated from the first singular vector of the SVD of the spatial-orientation filter. We used the time mixture parameters b τ (Equation 5) as an estimate of the temporal selectivity of the retinotopic component of the LFP (Figures 11E,F). For the MUA, we opted to take the first singular vector of (spatial-orientation)-temporal filters estimated by reverse correlation as an estimate of the temporal filters. These are more directly comparable to b τ than u τ in Equation (4), since u τ corrects for the slight auto-correlation in the stimulus while b τ , being ultimately based on a reverse correlation estimate, does not.